Python绘制概率分布的百分比轮廓线

Question

给定一个具有未知函数形式的概率分布（以下为示例），我希望绘制“基于百分位”的轮廓线，即对应于积分为10％，20％，...，90％等区域的轮廓线。

## example of an "arbitrary" probability distribution ##
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt
import numpy as np

X, Y = np.mgrid[-3:3:100j, -3:3:100j]
z1 = bivariate_normal(X, Y, .5, .5, 0., 0.)
z2 = bivariate_normal(X, Y, .4, .4, .5, .5)
z3 = bivariate_normal(X, Y, .6, .2, -1.5, 0.)
z = z1+z2+z3
plt.imshow(np.reshape(z.T, (100,-1)), origin='lower', extent=[-3,3,-3,3])
plt.show()

我已经尝试了多种方法，包括使用matplotlib中的默认轮廓函数，使用scipy中的stats.gaussian_kde方法，甚至可能从分布中生成随机点样本，然后估计内核。但是没有一个方法提供解决方案。

- neither-nor

---

你的问题表述不清。有无数种方法可以将你的示例图片分割，以便每个分割面积的一半都为50%。你想要哪种分割方式？听起来你想要等高线 - 但只是那些对应于10％，20％，...，90％积分区域的等高线。是这样吗？ - Timothy Shields

@TimothyShields 感谢您的澄清。您更好地表述了我想要的内容。 - neither-nor

Answer 1

观察p（x）在p（x）≥t的轮廓内的积分，并解决所需的t值：

import matplotlib
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt
import numpy as np

X, Y = np.mgrid[-3:3:100j, -3:3:100j]
z1 = bivariate_normal(X, Y, .5, .5, 0., 0.)
z2 = bivariate_normal(X, Y, .4, .4, .5, .5)
z3 = bivariate_normal(X, Y, .6, .2, -1.5, 0.)
z = z1 + z2 + z3
z = z / z.sum()

n = 1000
t = np.linspace(0, z.max(), n)
integral = ((z >= t[:, None, None]) * z).sum(axis=(1,2))

from scipy import interpolate
f = interpolate.interp1d(integral, t)
t_contours = f(np.array([0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]))
plt.imshow(z.T, origin='lower', extent=[-3,3,-3,3], cmap="gray")
plt.contour(z.T, t_contours, extent=[-3,3,-3,3])
plt.show()

- Timothy Shields

---

谢谢您提供的完美解决方案！然而，我完全不知道如何理解这一行代码：integral = ((z >= t[:, None, None]) * z).sum(axis=(1,2))。另外，是否有办法用0.9、0.8、0.7等标签来标记等高线？ - neither-nor

@neither-nor：z 是一个代表分布概率 p(x) 的二维数组。t 是一个一维数组，包含从 0 到 z.max() 不同的阈值。mask = (z >= t[:, None, None]) 是一个形状为 t.shape + z.shape 的三维数组，其中每个 mask[i] 都是一个二维数组，由 0/1 值组成，其中 1 表示在轮廓线 p(x) >= t[i] 的内部。integral = (mask * z).sum(axis=(1,2)) 是一个一维数组，包含这些区域的积分值，其中 integral[i] 是 p(x) 在轮廓线 p(x) >= t[i] 区域上的积分值。 - Timothy Shields

Answer 2

你可以像这样做：

from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import numpy as np

X, Y = np.mgrid[-3:3:100j, -3:3:100j]

sigma = 0.5

z = bivariate_normal(X,Y,.5, .5, 0., 0.)
z1 = bivariate_normal(0, 1 * sigma, sigma, sigma, 0.0, 0.0)
z2 = bivariate_normal(0, 2 * sigma, sigma, sigma, 0.0, 0.0)
z3 = bivariate_normal(0, 3 * sigma, sigma, sigma, 0.0, 0.0)

plt.imshow(z, interpolation='bilinear', origin='lower', extent=[-3,3,-3,3])
contour = plt.contour(z,[z1,z2,z3],origin='lower',extent=[-3,3,-3,3],colors='yellow')
plt.show()

这将会得到：

- Srivatsan

---

这些等高线对应于分布的高度，而不是其中的体积量。 - neither-nor