Python 中的多元正态分布积分

Question

我正在尝试在 python 中整合多元分布。为了测试它，我用二元正态分布构建了这个玩具示例。我使用 nquad() 以便稍后将其扩展到两个以上的变量。这是代码：

import numpy as np
from scipy import integrate
from scipy.stats import multivariate_normal


def integrand(x0, x1, mean, cov):
    return multivariate_normal.pdf([x0, x1], mean=mean, cov=cov)

mean = np.array([100, 100])
cov = np.array([[20, 0], [0, 20]])

res, err = integrate.nquad(integrand,
                           [[-np.inf, np.inf], [-np.inf, np.inf]],
                           args=(mean, cov))

print(res)

我得到的结果是9.559199162933625e-10。显然，这是不正确的。应该是（接近）1.

这里有什么问题？

Answer 1

scipy 的 nquad 仅在有界矩形域上进行数值积分。您的积分完全收敛的事实是由于 PDF 的 exp(-r^2) 类型权重（参见 here for its explicit form). Hence, you need Hermite quadrature in 2D. Some articles exist on this topic, and quadpy（我的一个项目）实现了这些。

您首先需要将积分转化为包含精确权重 exp(-r**2) 的形式，其中 r**2 为 x[0]**2 + x[1]**2。然后你削减这个重量并将它送入 quadpy 的 e2r2 正交：

import numpy
import quadpy


def integrand(x):
    return 1 / numpy.pi * numpy.ones(x.shape[1:])


val = quadpy.e2r2.integrate(
    integrand,
    quadpy.e2r2.RabinowitzRichter(3)
    )

print(val)

1.0000000000000004

Answer 2

有点跑题了，不过你应该改用下面的例程（速度相当快）：

from scipy.stats.mvn import mvnun
import numpy as np

mean = np.array([100, 100])
cov = np.array([[20, 0], [0, 20]])
mvnun(np.array([-np.inf, -np.inf]), np.array([np.inf, np.inf]), mean, cov)

或使用 multivariate_normal.cdf 并做减法。

Python 中的多元正态分布积分

Integration of Multivariate Normal Distribution in Python

python

numerical-integration