如何用Scipy的multivariate_normal或integrate计算指定多元正态概率？

能否用scipy.stats.multivariate_normal简便计算？ #

不能直接实现。原因是scipy.stats.multivariate_normal.cdf仅支持计算轴对齐矩形区域的联合概率（即$P(Z_1 \leq x_1, Z_2 \leq x_2, \dots, Z_n \leq x_n)$这类形式），而我们需要的概率是在矩形区域$Z_1 < c_1, \dots, Z_n < c_n$内，额外满足$\sum_{i=1}^n Z_i > 0$的子区域概率——这个子区域由线性约束$\sum Z_i >0$切割矩形得到，并非标准的轴对齐矩形，因此无法通过multivariate_normal的现有方法直接计算。

用scipy.integrate实现的方法 #

由于各$Z_i$独立，联合概率密度函数（PDF）是单个正态PDF的乘积：
$$
f(z_1, z_2, \dots, z_n) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}} \exp\left( -\frac{(z_i - \mu_i)^2}{2} \right)
$$

我们需要在以下积分区域计算该PDF的积分：

• $z_1 \in (-\infty, c_1), z_2 \in (-\infty, c_2), \dots, z_{n-1} \in (-\infty, c_{n-1})$
• $z_n \in \left( \max\left(-\infty, -\sum_{i=1}^{n-1} z_i\right), c_n \right)$

代码实现示例 #

1. 通用n维情况（使用nquad）

import numpy as np
from scipy.integrate import nquad
from scipy.stats import norm

def joint_pdf(z, mus):
    # z为长度n的变量数组，mus为对应均值数组
    return np.prod([norm.pdf(zi, mu, 1) for zi, mu in zip(z, mus)])

def compute_probability(n, mus, cs):
    # 前n-1个变量的积分限：(-∞, ci)
    limits = [(-np.inf, ci) for ci in cs[:-1]]
    
    # 第n个变量的积分限：依赖前n-1个变量的和
    def zn_limit(*args):
        sum_prev = sum(args)
        lower = -sum_prev
        upper = cs[-1]
        return (lower, upper)
    
    limits.append(zn_limit)
    
    # 执行n维积分
    result, error = nquad(joint_pdf, limits, args=(mus,))
    return result, error

# 示例调用：n=3，均值[1, 0, -0.5]，阈值[2, 1.5, 1]
n = 3
mus = [1, 0, -0.5]
cs = [2, 1.5, 1]
prob, err = compute_probability(n, mus, cs)
print(f"计算得到的概率：{prob:.4f}，积分误差：{err:.4e}")

2. 低维特殊情况（比如n=2，用dblquad更直观）

from scipy.integrate import dblquad
from scipy.stats import norm

def compute_2d_prob(mu1, mu2, c1, c2):
    # 被积函数：z1和z2的联合PDF
    def pdf(z2, z1):
        return norm.pdf(z1, mu1, 1) * norm.pdf(z2, mu2, 1)
    
    # z1范围：(-∞, c1)；z2范围：(-z1, c2)
    result, error = dblquad(pdf, -np.inf, c1, lambda z1: -z1, lambda z1: c2)
    return result, error

# 示例调用
prob_2d, err_2d = compute_2d_prob(1, 0, 2, 1.5)
print(f"二维情况概率：{prob_2d:.4f}，积分误差：{err_2d:.4e}")

注意事项 #

• 当n较大（比如n>4）时，n维积分的计算速度会显著变慢，精度也可能下降，此时可考虑蒙特卡洛采样等近似方法替代。
• 可通过调整nquad或dblquad的epsabs、epsrel参数，平衡计算速度与精度。

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