使用scipy在Python中计算多元正态分布的累积分布函数

Question

为了计算多元正态分布的累积分布函数，我按照这个例子（针对一元情况）进行操作，但无法解释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]]

如果联合累积分布函数定义为：

P (X1 ≤ x1, . . . ,Xn ≤ xn)

那么期望的输出应该是一个介于0和1之间的实数。

- statBeginner

---

我认为你不能在多元情况下使用 scipy.stats.norm。 - cel

scipy.stats库中有multivariate_normal函数（http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multivariate_normal.html），但是它没有`cdf`方法。 - Warren Weckesser

Answer 1

在大量搜索后，我认为Noah H. Silbert的这篇博客文章描述了一个可以用于计算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

- statBeginner

---

能否将一个点数组传递给 mvn.mvnun 函数？我看了代码，好像只能通过循环来实现？ - ZK Zhao

@cqcn1991，我正在寻找多元CDF，通过文件传递一个数组进行绘制。你能找到解决方案了吗？请看这里：http://stackoverflow.com/questions/37057938/bivariate-cdf-ccdf-distribution-python - Sitz Blogz

mvn.mvnun 的问题在于它不是确定性的。至少，这段代码每次都会给出不同的结果：https://pastebin.com/L0WSTRui - David Dale

这是那篇博客文章的链接（https://www.statisticalmodelcitizen.com/2018/11/19/multivariate-normal-cdf-values-in-python/），我重新发布了它（我的旧博客数据库一段时间前就已经损坏了，最近才恢复了这些文章）。虽然Genz开发的算法不是确定性的，但该代码产生的概率仅在第9位小数处有所不同。对我来说，计算多元正态积分的快速准确算法的好处远远超过它不是确定性的成本。 - Noah Motion

请注意，2D情况被视为特殊情况。如果这是与R使用的实现相同（我猜是），它是确定性的，并且可能非常接近精确：https://github.com/scipy/scipy/blob/8a64c938ddf1ae4c02a08d2c5e38daeb8d061d38/scipy/stats/mvndst.f#L988-L1031 - Benjamin Christoffersen

Answer 2

从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

文档可以在这里找到。

- adamconkey

Answer 3

如果您不关心性能（即仅偶尔执行），则可以使用multivariate_normal创建多元正态分布pdf，然后通过integrate.nquad计算cdf。

- ZK Zhao

---

请您详细说明如何使用这个？并且，这能用于寻找依赖于多元正态分布的函数的期望值吗？ - vicky113