Python绘制概率分布百分位等高线

2024-04-25 14:53:02 发布

您现在位置:Python中文网/ 问答频道 /正文

给定一个未知函数形式的概率分布(如下示例),我喜欢绘制“基于百分位”的等高线,即那些与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()

enter image description here 我研究了多种方法,从在matplotlib中使用默认轮廓函数,方法包括高斯统计在scipy中,甚至可能从分布中生成随机点样本,然后估计一个核。它们似乎都没有提供解决方案。在


Tags: 方法函数importmatplotlibasnpplt形式
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()

它给出了:

enter image description here

在轮廓p(x)≥t内观察p(x)的积分,求出所需的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()

enter image description here

相关问题 更多 >