绘制显示粒子百分比的等高线

Question

我想要制作的是类似于这个情节的东西：

这是一个等高线图，代表两个数据集中包含的 68%、95%、99.7% 的粒子。

到目前为止，我已经尝试实现高斯 KDE 估计，并在轮廓上绘制这些粒子高斯。

这里添加文件https://www.dropbox.com/sh/86r9hf61wlzitvy/AABG2mbmmeokIiqXsZ8P76Swa?dl=0

from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
import numpy as np

# My data
x = RelDist
y = RadVel

# Peform the kernel density estimate
k = gaussian_kde(np.vstack([RelDist, RadVel]))
xi, yi = np.mgrid[x.min():x.max():x.size**0.5*1j,y.min():y.max():y.size**0.5*1j]
zi = k(np.vstack([xi.flatten(), yi.flatten()]))



fig = plt.figure()
ax = fig.gca()


CS = ax.contour(xi, yi, zi.reshape(xi.shape), colors='darkslateblue')
plt.clabel(CS, inline=1, fontsize=10)

ax.set_xlim(20, 800)
ax.set_ylim(-450, 450)
ax.set_xscale('log')

plt.show()

制作这个：

]2

其中 1) 我不知道如何在 gaussain kde 中控制 bin 数，2) 等高线标签全为零，3) 我不知道如何确定百分位数。

感谢任何帮助。

Answer 1

摘自此 example in the matplotlib 文档

您可以将数据 zi 转换为百分比比例 (0-1)，然后绘制等高线图。

您还可以在调用 plt.contour() 时手动确定计数图的级别。

下面是一个包含 2 个随机生成的正态双变量分布的示例：

delta = 0.025
x = y = np.arange(-3.0, 3.01, delta)
X, Y = np.meshgrid(x, y)
Z1 = plt.mlab.bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = plt.mlab.bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = 10* (Z1- Z2)

#transform zi to a 0-1 range
Z = Z = (Z - Z.min())/(Z.max() - Z.min())

levels =  [0.68, 0.95, 0.997] 
origin = 'lower'
CS = plt.contour(X, Y, Z, levels,
              colors=('k',),
              linewidths=(3,),
              origin=origin)

plt.clabel(CS, fmt='%2.3f', colors='b', fontsize=14)

使用您提供的数据，代码同样有效：

from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
import numpy as np

RadVel = np.loadtxt('RadVel.txt')
RelDist = np.loadtxt('RelDist.txt')
x = RelDist
y = RadVel

k = gaussian_kde(np.vstack([RelDist, RadVel]))
xi, yi = np.mgrid[x.min():x.max():x.size**0.5*1j,y.min():y.max():y.size**0.5*1j]
zi = k(np.vstack([xi.flatten(), yi.flatten()]))

#set zi to 0-1 scale
zi = (zi-zi.min())/(zi.max() - zi.min())
zi =zi.reshape(xi.shape)

#set up plot
origin = 'lower'
levels = [0,0.1,0.25,0.5,0.68, 0.95, 0.975,1]

CS = plt.contour(xi, yi, zi,levels = levels,
              colors=('k',),
              linewidths=(1,),
              origin=origin)

plt.clabel(CS, fmt='%.3f', colors='b', fontsize=8)
plt.gca()
plt.xlim(10,1000)
plt.xscale('log')
plt.ylim(-200,200)

Answer 2

@Tkanno 的回答在编程上是正确的，但并没有完全按照问题中的要求去做。

kde returns 样本根据建模分布的可能性。因此，等高线图是对样本概率的限制。 0.1 等高线图将显示样本根据建模分布出现的机会小于 10% 的限制。现在通过按照 Tkanno 的建议对 z 值进行归一化，现在绘制的是相对概率，因此在 Tkanno 的回答中，0.1 等高线图是一个极限，超过该极限，样本出现的可能性比最有可能的样本低 10 倍。

您可以通过绘制二维直方图、使用最频繁的 bin 归一化并绘制具有相同水平的等高线来绘制与 Tkanno 提出的非常相似的等高线图（但尚未平滑）。

这不能与包含90%数据的极限相提并论。我认为包含给定部分数据的等高线图的获取要复杂一些（参见 https://stats.stackexchange.com/questions/68105/contours-containing-a-given-fraction-of-x-y-points 和包图的解决方案）。显然在 R 中有一个包图的实现，也许有人 has/will 为 python 做了它。

为了说明解决问题的难度，可以想到一个有 100 个点的数据集。任何包含 95 分（不包括 5 分）的卷实际上都可以回答这个问题。可能隐含地问的是包含95个点的最小体积（因此代表最高的可能性或密度），这是一个组合优化问题。

绘制显示粒子百分比的等高线

Plotting contour lines that show percentage of particles

python

arrays

numpy

matplotlib

contour