在 SciPy 中计算框区域上的多元正态积分

Question

我正在尝试评估多元正态分布以获得

1 < z1 < 2.178

和

2.178 < z2 < inf

的概率。分布的均值为零，所有方差均为 1，协方差为

1/sqrt(2)

。

在

R

中，我通过以下方式获得了正确的结果：

library(mnormt)
sadmvn(lower=c(1, 2.178), upper=c(2.178, Inf), mean=0, varcov=matrix(c(1, 1/sqrt(2), 1/sqrt(2), 1),2, 2))

在 SciPy 中，我正在尝试：

from scipy import stats
cov = np.array([[1, 1 / np.sqrt(2)], [1 / np.sqrt(2), 1]])
d = stats.multivariate_normal(mean=np.array([0, 0]), cov=cov)
d.cdf([2.178, np.inf]) - d.cdf([1, 2.178])

但这给出了 0.14 左右的值，而不是正确的结果 ~0.0082。我怎样才能正确地做到这一点？

Answer 1

您可以使用

lower_limit

 方法的 

multivariate_normal.cdf 参数来计算具有角

lower_limit

和

x

（第一个位置参数）的超矩形内的 CDF。

import numpy as np
from scipy import stats
cov = np.array([[1, 1 / np.sqrt(2)], [1 / np.sqrt(2), 1]])
d = stats.multivariate_normal(mean=np.array([0, 0]), cov=cov)
d.cdf([2.178, np.inf], lower_limit=[1, 2.178])
# 0.00820893693689204

看起来

lower_limit

的

x

和

multivariate_normal

分别对应于

lower

的

upper

和

sadmvn

。

您原来的方法不起作用的原因是，当在没有

cdf

的情况下调用

lower_limit

时，积分的下限被假定为

[-np.inf, -np.inf]

。因此，您的代码正在计算两个超矩形的 CDF 之间的差异，每个超矩形都有下限

[-np.inf, -np.inf]

。这个想法在一维中有效，但在 2+ 维度中，这与从角

[1, 2.178]

到角

[2.178, np.inf]

的单个超矩形内的 CDF 不同。

例如，在下图中，A点和B点定义的矩形内包围的面积与蓝色和红色矩形的面积之差并不相同。