使用matplotlib在散点图中创建一个置信椭圆

Question

使用matplotlib在散点图中创建一个置信椭圆

25

如何使用matplotlib在散点图中创建置信椭圆？

以下代码可以成功创建散点图，但是有没有人熟悉如何在散点图上添加置信椭圆？

import numpy as np
import matplotlib.pyplot as plt
x = [5,7,11,15,16,17,18]
y = [8, 5, 8, 9, 17, 18, 25]

plt.scatter(x,y)
plt.show()

以下是来自SAS的置信椭圆参考。

http://support.sas.com/documentation/cdl/en/grstatproc/62603/HTML/default/viewer.htm#a003160800.htm

在SAS中的代码是这样的：

proc sgscatter data=sashelp.iris(where=(species="Versicolor"));
  title "Versicolor Length and Width";
  compare y=(sepalwidth petalwidth)
          x=(sepallength petallength)
          / reg ellipse=(type=mean) spacing=4;
run;

- 2964502

1

可能是多维置信区间的重复问题。 - Saullo G. P. Castro

2

@Saullo Castro，你看过SAS中的代码吗？你认为在SAS中实现的方法和你提供的链接中的方法是一样的吗？ - 2964502

2

@tester3 - 在你提供的例子中，显示的置信椭圆是针对均值而不是从同一总体中抽取的另一个样本。(这就是type=mean所指定的内容。) 我的答案，@SaulloCastro链接到的，显示了整个总体的置信椭圆(换句话说，另一个来自该总体的样本应该落在其中的区域，与SAS中的type=predicted相同)。Jamie的答案也使用了这种方法。 - Joe Kington

4个回答

29

尝试使用已接受的答案后，我发现在计算theta时它没有正确选择象限，因为它依赖于 np.arccos:

oops

查看'possible duplicate'和Joe Kington在Github上的解决方案后，我简化了他的代码如下：

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

def eigsorted(cov):
    vals, vecs = np.linalg.eigh(cov)
    order = vals.argsort()[::-1]
    return vals[order], vecs[:,order]

x = [5,7,11,15,16,17,18]
y = [25, 18, 17, 9, 8, 5, 8]

nstd = 2
ax = plt.subplot(111)

cov = np.cov(x, y)
vals, vecs = eigsorted(cov)
theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
w, h = 2 * nstd * np.sqrt(vals)
ell = Ellipse(xy=(np.mean(x), np.mean(y)),
              width=w, height=h,
              angle=theta, color='black')
ell.set_facecolor('none')
ax.add_artist(ell)
plt.scatter(x, y)
plt.show()

负斜率

- Ben

2

有人知道如何将这个推广到三维（或可能是n维）吗？ - HansSnah

@Ben 为什么 w, h = 2 * nstd * np.sqrt(vals)？那应该是 4 * np.sqrt(vals)。 - tauran

为什么当我将plt.scatter(x, y)更改为plt.scatter(np.mean(x), np.mean(y))时，它不起作用。我想要围绕平均值的椭圆形。 - Farshid Rayhan

我尝试了这个例子，但是得到了一个空白的图形。有人知道为什么吗？ - vicemagui

@tauran：椭圆的半径为sqrt(eigenvalue1)和sqrt(eigenvalue2)，而vals包含了这两个特征值：[eigenvalue1, eigenvalue2]。在绘制椭圆时，需要指定其宽度而非半径：因此为2*sqrt(vals)。现在我们有许多椭圆，但如果您想要显示置信区间，使得95%的数据点位于该椭圆区域内，那么我们需要显示2倍的标准差。因此乘以nstd（=2）。 https://en.wikipedia.org/wiki/Normal_distribution#Standard_deviation_and_coverage & https://cookierobotics.com/007/ - Jürgen Brauer

2

一旦您拥有协方差矩阵的特征分解，就无需显式计算角度：旋转部分已经免费为您编码了该信息。

cov = np.cov(x, y)
val, rot = np.linalg.eig(cov)
val = np.sqrt(val)
center = np.mean([x, y], axis=1)[:, None]

t = np.linspace(0, 2.0 * np.pi, 1000)
xy = np.stack((np.cos(t), np.sin(t)), axis=-1)

plt.scatter(x, y)
plt.plot(*(rot @ (val * xy).T + center))

你可以在平移前应用缩放来扩大你的椭圆形：

plt.plot(*(2 * rot @ (val * xy).T + center))

- Mad Physicist

1

除了被接受的答案之外：我认为正确的角度应该是：

angle=np.rad2deg(np.arctan2(*v[:,np.argmax(abs(lambda_))][::-1])))

相应的宽度（较大的特征值）和高度应为：

width=lambda_[np.argmax(abs(lambda_))]*j*2, height=lambda_[1-np.argmax(abs(lambda_))]*j*2

由于“特征值不一定有序”，因此我们需要找到对应于最大特征值的特征向量。根据规格https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html，v [:，i]是对应于eigenvalue lambda_ [i]的特征向量; 我们应该通过np.argmax(abs(lambda_))找到特征向量的正确列。

- Kerem

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jaime · Accepted Answer

以下代码绘制一个标准差、两个标准差和三个标准差大小的椭圆：

x = [5,7,11,15,16,17,18]
y = [8, 5, 8, 9, 17, 18, 25]
cov = np.cov(x, y)
lambda_, v = np.linalg.eig(cov)
lambda_ = np.sqrt(lambda_)
from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt
ax = plt.subplot(111, aspect='equal')
for j in xrange(1, 4):
    ell = Ellipse(xy=(np.mean(x), np.mean(y)),
                  width=lambda_[0]*j*2, height=lambda_[1]*j*2,
                  angle=np.rad2deg(np.arccos(v[0, 0])))
    ell.set_facecolor('none')
    ax.add_artist(ell)
plt.scatter(x, y)
plt.show()

在这里输入图片描述