为了计算多元正态分布的 CDF,我遵循了 this 示例(对于单变量情况),但无法解释 scipy 生成的输出:
from scipy.stats import norm
import numpy as np
mean = np.array([1,5])
covariance = np.matrix([[1, 0.3 ],[0.3, 1]])
distribution = norm(loc=mean,scale = covariance)
print distribution.cdf(np.array([2,4]))
产生的输出是:
[[ 8.41344746e-01 4.29060333e-04]
[ 9.99570940e-01 1.58655254e-01]]
如果联合CDF定义为:
P (X1 ≤ x1, . . . ,Xn ≤ xn)
那么预期输出应该是 0 到 1 之间的实数。
经过大量搜索后,我认为 Noah H. Silbert 的 this 博客文章描述了标准库中唯一可用于计算 Python 中多元正态分布的 cdf 的现成代码。 Scipy 有一种方法可以做到这一点,但正如博客中提到的,很难找到。该方法基于 Alan Genz 的一篇论文。
从博客来看,这就是它的工作原理。
from scipy.stats import mvn
import numpy as np
low = np.array([-10, -10])
upp = np.array([.1, -.2])
mu = np.array([-.3, .17])
S = np.array([[1.2,.35],[.35,2.1]])
p,i = mvn.mvnun(low,upp,mu,S)
print p
0.2881578675080012
v1.1.0 版本的 scipy
multivariate_normal
现在内置了 cdf 函数:
from scipy.stats import multivariate_normal as mvn
import numpy as np
mean = np.array([1,5])
covariance = np.array([[1, 0.3],[0.3, 1]])
dist = mvn(mean=mean, cov=covariance)
print("CDF:", dist.cdf(np.array([2,4])))
CDF: 0.14833820905742245
可以在此处找到文档。
如果您不关心性能(即只是偶尔执行),那么您可以使用
multivariate_normal
创建多元正态 pdf,然后通过 integrate.nquad
计算 cdf
scipy.stats.multivariate_normal.cdf
适用于矩形区域,给出随机向量的分量小于给定向量的分量的概率。如果我们想要样本位于由 Mahalanobis 距离 确定的椭球内的概率,则可以使用 chi2.cdf
来完成(请参阅本文):
from scipy.stats import chi2
import numpy as np
mean = np.array([1, 2])
covariance = np.array([[1, 0.8],[0.8, 1]])
x = np.array([2, 3])
y = x - mean
r2 = y @ np.linalg.inv(covariance) @ y
print(chi2.cdf(r2, len(x))) # 0.4262465792625671