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

Question

我想做一些可能很简单但给我带来困难的事情。尝试从多元正态分布中抽取

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

  中不使用关键字参数

  mean

  和 cov，因为在

  mean

  中使用

  cov

  和

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

  定义分布似乎很奇怪然后在

  rvs

  调用中重复该定义。我是不是错过了什么？
• 对于

  mvn.pdf

  调用，概率密度显然是 >>>1，这对于连续多元正态分布来说并非不可能，但我想将这些数字转换为该特定点的近似概率。我该怎么做？

谢谢！

Answer 1

我不使用文档中的关键字参数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

我想将这些数字转换为该特定点的近似概率。我该怎么做？...我想计算样本值的特定 epsilon 内的概率。

您可以定义一个边长为

eps

、以每个点为中心的超立方体，并评估该超立方体内的累积密度。

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

如果您更喜欢超球面区域，您可以乘以超球体的体积而不是超立方体的体积。对于以

eps

作为超球面直径的 6 维空间，

vol = np.pi**3/6 * (eps/2)**6

。