多维置信区间

Question

我有许多元组(par1, par2)，即从多次重复实验中获得的二维参数空间中的点。

我正在寻找一种计算和可视化置信椭圆的方法（不确定这是否是正确术语）。下面是我在网上找到的一个示例图来说明我的意思：

来源: blogspot.ch/2011/07/classification-and-discrimination-with.html

因此，原则上需要将多元正态分布拟合到数据点的二维直方图中。有人能帮我吗？

- Raphael Roth

---

输入数据是什么？它是一个二维点的数组吗？您是否事先知道有2个聚类？ - Daniel

是的，我知道聚类的数量。但我还不知道输入数据的格式，我猜测它是一个nx2的数组，其中n是点的数量。 - Raphael Roth

在这种情况下，您应该首先对它们进行聚类，然后对每个聚类拟合高斯分布，最后绘制置信区间。请查看sklearn.cluster。 - Daniel

Answer 1

听起来你只是想要散点图的2-sigma椭圆？

如果是的话，考虑像这样做（参考这里的一些论文代码：https://github.com/joferkington/oost_paper_code/blob/master/error_ellipse.py）：

import numpy as np

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

def plot_point_cov(points, nstd=2, ax=None, **kwargs):
    """
    Plots an `nstd` sigma ellipse based on the mean and covariance of a point
    "cloud" (points, an Nx2 array).

    Parameters
    ----------
        points : An Nx2 array of the data points.
        nstd : The radius of the ellipse in numbers of standard deviations.
            Defaults to 2 standard deviations.
        ax : The axis that the ellipse will be plotted on. Defaults to the 
            current axis.
        Additional keyword arguments are pass on to the ellipse patch.

    Returns
    -------
        A matplotlib ellipse artist
    """
    pos = points.mean(axis=0)
    cov = np.cov(points, rowvar=False)
    return plot_cov_ellipse(cov, pos, nstd, ax, **kwargs)

def plot_cov_ellipse(cov, pos, nstd=2, ax=None, **kwargs):
    """
    Plots an `nstd` sigma error ellipse based on the specified covariance
    matrix (`cov`). Additional keyword arguments are passed on to the 
    ellipse patch artist.

    Parameters
    ----------
        cov : The 2x2 covariance matrix to base the ellipse on
        pos : The location of the center of the ellipse. Expects a 2-element
            sequence of [x0, y0].
        nstd : The radius of the ellipse in numbers of standard deviations.
            Defaults to 2 standard deviations.
        ax : The axis that the ellipse will be plotted on. Defaults to the 
            current axis.
        Additional keyword arguments are pass on to the ellipse patch.

    Returns
    -------
        A matplotlib ellipse artist
    """
    def eigsorted(cov):
        vals, vecs = np.linalg.eigh(cov)
        order = vals.argsort()[::-1]
        return vals[order], vecs[:,order]

    if ax is None:
        ax = plt.gca()

    vals, vecs = eigsorted(cov)
    theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))

    # Width and height are "full" widths, not radius
    width, height = 2 * nstd * np.sqrt(vals)
    ellip = Ellipse(xy=pos, width=width, height=height, angle=theta, **kwargs)

    ax.add_artist(ellip)
    return ellip

if __name__ == '__main__':
    #-- Example usage -----------------------
    # Generate some random, correlated data
    points = np.random.multivariate_normal(
            mean=(1,1), cov=[[0.4, 9],[9, 10]], size=1000
            )
    # Plot the raw points...
    x, y = points.T
    plt.plot(x, y, 'ro')

    # Plot a transparent 3 standard deviation covariance ellipse
    plot_point_cov(points, nstd=3, alpha=0.5, color='green')

    plt.show()

- Joe Kington

---

@JoeKington 我们不需要参考卡方概率分布表来确定我们的 nstd，即无论是 68%，90% 还是 95%。 - Srivatsan

@ThePredator - 如果你将其用作测试，那么是的。（换句话说，这与p置信水平下的另一个分布是否不同/相同？）如果你只是将其用作描述，则不是。你正确估计标准偏差和均值的置信度与你拥有的样本数量是完全独立的问题。 - Joe Kington

@JoeKington：最后一个问题，您能否详细说明使用 np.degrees(np.arctan2(*vecs[:,0][::-1])) 计算角度的原理？从这个网站上看到是 arctan(y)/(x)，但您使用了 arctan2。 - Srivatsan

@ThePredator - arctan2 返回完整的角度（可以在任意一个四象限内）。arctan 的输出被限制在第一和第四象限之间（介于-pi / 2和pi / 2之间）。您可能会注意到，arctan 只需要一个参数。因此，它无法区分第一和第四象限中的角度以及第二和第三象限中的类似角度。这是许多其他编程语言共享的惯例，其中 C 定义了它们。 - Joe Kington

@ThePredator 正确的实现方式已由Syrtis Major在另一个回答中提供。我想这个答案应该被修改，以包括这个sigma定义不对应95％的说明。 - Zappel

Answer 2

请参考文章如何绘制协方差误差椭圆。

以下是Python代码实现：

import numpy as np
from scipy.stats import norm, chi2

def cov_ellipse(cov, q=None, nsig=None, **kwargs):
    """
    Parameters
    ----------
    cov : (2, 2) array
        Covariance matrix.
    q : float, optional
        Confidence level, should be in (0, 1)
    nsig : int, optional
        Confidence level in unit of standard deviations. 
        E.g. 1 stands for 68.3% and 2 stands for 95.4%.

    Returns
    -------
    width, height, rotation :
         The lengths of two axises and the rotation angle in degree
    for the ellipse.
    """

    if q is not None:
        q = np.asarray(q)
    elif nsig is not None:
        q = 2 * norm.cdf(nsig) - 1
    else:
        raise ValueError('One of `q` and `nsig` should be specified.')
    r2 = chi2.ppf(q, 2)

    val, vec = np.linalg.eigh(cov)
    width, height = 2 * sqrt(val[:, None] * r2)
    rotation = np.degrees(arctan2(*vec[::-1, 0]))

    return width, height, rotation

< p > Joe Kington的答案中，标准差的意义错误。通常我们使用1、2 sigma表示68%和95%的置信水平，但是他答案中的2 sigma椭圆不包含总分布的95%概率。正确的方法是使用卡方分布来估计椭圆大小，如文章所示。

- Syrtis Major

---

他回答中显示的椭圆不是2西格玛椭圆。它是一个3西格玛椭圆，并且包含了大约与3西格玛椭圆相同数量的点。 - senderle

我认为这种差异是因为这个答案描述了一个置信椭圆，而Joe的答案描述了一个N-sigma误差椭圆。这两者之间的区别在这里有所解释。 - Gabriel

由于PDF文件的链接已经失效，我在此更新我的评论。这是新链接。 - Gabriel

Answer 3

我稍微修改了上面的一个例子，用于绘制错误或置信区间轮廓。现在我认为它给出了正确的轮廓。

之前它给出了错误的轮廓，因为它将scoreatpercentile方法应用于联合数据集（蓝色 + 红色点），而应该分别应用于每个数据集。

修改后的代码如下：

import numpy
import scipy
import scipy.stats
import matplotlib.pyplot as plt

# generate two normally distributed 2d arrays
x1=numpy.random.multivariate_normal((100,420),[[120,80],[80,80]],400)
x2=numpy.random.multivariate_normal((140,340),[[90,-70],[-70,80]],400)

# fit a KDE to the data
pdf1=scipy.stats.kde.gaussian_kde(x1.T)
pdf2=scipy.stats.kde.gaussian_kde(x2.T)

# create a grid over which we can evaluate pdf
q,w=numpy.meshgrid(range(50,200,10), range(300,500,10))
r1=pdf1([q.flatten(),w.flatten()])
r2=pdf2([q.flatten(),w.flatten()])

# sample the pdf and find the value at the 95th percentile
s1=scipy.stats.scoreatpercentile(pdf1(pdf1.resample(1000)), 5)
s2=scipy.stats.scoreatpercentile(pdf2(pdf2.resample(1000)), 5)

# reshape back to 2d
r1.shape=(20,15)
r2.shape=(20,15)

# plot the contour at the 95th percentile
plt.contour(range(50,200,10), range(300,500,10), r1, [s1],colors='b')
plt.contour(range(50,200,10), range(300,500,10), r2, [s2],colors='r')

# scatter plot the two normal distributions
plt.scatter(x1[:,0],x1[:,1],alpha=0.3)
plt.scatter(x2[:,0],x2[:,1],c='r',alpha=0.3)

- Rodrigo

Answer 4

我猜你想要计算置信区间。

我对此不是很了解，但作为一个起点，我建议你检查Python的sherpa应用程序。至少在他们2011年Scipy演讲中，作者提到你可以用它确定和获取置信区间（但你可能需要有数据模型）。

请查看Sherpa演讲的视频和幻灯片。

希望这能帮到你。

- gcalmettes

---

我也看了Sherpa文档，但实际上我不知道这是什么 :) - Raphael Roth