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

l7wslrjt  于 2023-08-05  发布在  其他
关注(0)|答案(4)|浏览(138)

如何使用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;

jc3wubiy

jc3wubiy1#

下面的代码绘制一个、两个和三个标准偏差大小的椭圆:

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()

字符串


的数据

z3yyvxxp

z3yyvxxp2#

在尝试了公认的答案之后,我发现它在计算theta时没有正确选择象限,因为它依赖于np.arccos


的数据
看看'possible duplicate'Joe Kington's solution on 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()

字符串


zaq34kh6

zaq34kh63#

一旦你有了协方差矩阵的特征分解,就不需要显式地计算Angular :旋转部分已经为您免费编码了该信息:

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))


jyztefdp

jyztefdp4#

  • 除了接受的答案:* 我认为正确的Angular 应该是:
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.htmlv[:,i]的“特征值不一定是有序的”是对应于eigenvalue lambda_[i]的特征向量;我们应该通过np.argmax(abs(lambda_))找到特征向量的正确列。

相关问题