numpy 使用scipy.stats.multivariate_normal从多元正态分布中绘制样本并计算样本概率

7rtdyuoh  于 2023-11-18  发布在  其他
关注(0)|答案(1)|浏览(96)

我想做一些可能非常简单,但给我带来困难的事情。尝试从多元正态分布中抽取N样本,并计算每个随机抽取样本的概率。在这里,我尝试使用scipy,但我也愿意使用np.random.multivariate_normal。这是最简单的。

>>> import numpy as np
>>> from scipy.stats import multivariate_normal

>>> num_samples = 10
>>> num_features = 6
>>> std = np.random.rand(num_features)

# define distribution
>>> mvn = multivariate_normal(mean = np.zeros(num_features), cov = np.diag(std), allow_singular = False, seed = 42)

# draw samples
>>> sample = mvn.rvs(size = num_samples); sample

# determine probability of each drawn sample
>>> prob = mvn.pdf(x = sample)

# print samples
>>> print(sample)
[[ 0.04816243 -0.00740458 -0.00740406  0.04967142 -0.01382643  0.06476885]
...
 [-0.00977815  0.01047547  0.03084945  0.10309995  0.09312801 -0.08392175]]

# print probability all samples
[26861.56848337 17002.29353025  2182.26793265  3749.65049331
 42004.63147989  3700.70037411  5569.30332186 16103.44975393
 14760.64667235 19148.40325233]

字符串
这让我感到困惑,原因有几个:

  • 对于rvs采样函数:我没有在docs中使用关键字参数meancov,因为在mvn = multivariate_normal(mean = np.zeros(num_features), cov = np.diag(std), allow_singular = False, seed = 42)中定义一个meancov的分布,然后在rvs调用中重复该定义似乎很奇怪。
  • 对于mvn.pdf调用,概率密度显然是>1,这对于连续多元正态分布来说并非不可能,但我想将这些数字转换为特定点的近似概率。

谢谢你,谢谢

agxfikkp

agxfikkp1#

我没有在文档中使用关键字arguments mean和cov.我错过了什么吗?
不,你所做的是允许的。发行版的设计允许调用带参数的方法(正如你在文档中所读到的那样),也允许“冻结”带参数的发行版和调用不带参数的方法。这些是等价的:

mean = np.zeros(num_features)
cov = np.diag(std)

mvn = multivariate_normal(mean=mean, cov=cov, seed=42)
sample = mvn.rvs(size=num_samples)
pdf = mvn.pdf(sample)

sample2 = multivariate_normal.rvs(mean=mean, cov=cov, size=num_samples, random_state=42)
pdf2 = multivariate_normal.pdf(sample2, mean=mean, cov=cov)

np.testing.assert_equal(sample2, sample)  # passes
np.testing.assert_equal(pdf2, pdf)  # passes

字符串
我想把这些数字转换成那个特定点的近似概率。我怎么做呢?.我想计算样本值的特定范围内的概率。
您可以定义一个边长为eps的超立方体,并计算该超立方体内的累积密度(使用SciPy 1.10.0+)。

eps = 0.01
mvn.cdf(sample - eps/2, lower_limit=sample + eps/2)
# array([2.87121214e-14, 1.81736055e-14, 2.33269634e-15, 4.00857084e-15,
#        4.48976867e-14, 3.95613589e-15, 5.95304832e-15, 1.72140983e-14,
#        1.57778144e-14, 2.04685939e-14])


你可以通过将概率密度乘以超立方体的体积得到近似相同的结果:

vol = eps**num_features
pdf * vol
# array([2.87145307e-14, 1.81751442e-14, 2.33280494e-15, 4.00830854e-15,
#        4.49021911e-14, 3.95598175e-15, 5.95348449e-15, 1.72142965e-14,
#        1.57788643e-14, 2.04692967e-14])


如果你喜欢超球面区域,你可以乘以超球面的体积,而不是超立方体的体积。对于一个6维空间,超球面的直径为epsvol = np.pi**3/6 * (eps/2)**6

相关问题