如何在SciPy中正确计算指定箱型区域的多元正态分布概率？

需要计算二元正态分布的概率：随机变量满足 (1 < z_1 < 2.178) 且 (2.178 < z_2 < \infty)。该分布均值为0，所有方差均为1，协方差为 (1/\sqrt{2})。

已在R中通过以下代码得到正确结果（约0.0082）：

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中尝试的代码得到约0.14的错误结果：

from scipy import stats
import numpy as np

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

原代码对多维累积分布函数（CDF）的逻辑理解有误：stats.multivariate_normal.cdf([x1, x2]) 计算的是联合概率 (P(z_1 \leq x1, z_2 \leq x2))，直接用 d.cdf([2.178, np.inf]) - d.cdf([1, 2.178]) 无法正确表示目标区间的概率——前者等价于 (P(z_1 \leq 2.178))，后者是 (P(z_1 \leq 1, z_2 \leq 2.178))，两者的差并非目标区间的概率。

方法1：使用scipy.stats.mvn.mvnun（推荐） #

该函数专门用于计算多元正态分布在矩形区域内的概率，功能与R的sadmvn完全对应：

import numpy as np
from scipy.stats import mvn

# 定义分布参数
mean = np.array([0, 0])
cov = np.array([[1, 1 / np.sqrt(2)], [1 / np.sqrt(2), 1]])

# 定义区间边界：lower为各维度下限，upper为各维度上限
lower = np.array([1, 2.178])
upper = np.array([2.178, np.inf])

# 计算概率
prob, _ = mvn.mvnun(lower, upper, mean, cov)
print(prob)  # 输出约为0.0082

方法2：通过联合CDF组合计算 #

利用联合概率的逻辑推导，拆分目标区间的概率：
目标概率 = (P(z_1 \leq 2.178, z_2 > 2.178) - P(z_1 \leq 1, z_2 > 2.178))
其中 (P(z_j > a, z_i \leq b) = P(z_i \leq b) - P(z_i \leq b, z_j \leq a))，代码实现如下：

from scipy import stats
import numpy as np

cov = np.array([[1, 1 / np.sqrt(2)], [1 / np.sqrt(2), 1]])
d = stats.multivariate_normal(mean=np.array([0, 0]), cov=cov)

# 计算P(z1<=2.178, z2>2.178)
p1 = d.cdf([2.178, np.inf]) - d.cdf([2.178, 2.178])
# 计算P(z1<=1, z2>2.178)
p2 = d.cdf([1, np.inf]) - d.cdf([1, 2.178])
# 目标概率
prob = p1 - p2
print(prob)  # 输出约为0.0082

内容的提问来源于stack exchange，提问作者clog14