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。

这里的问题是什么？

- user192859

Answer 1

有点离题，但你应该使用以下例程（它非常快）：

最初的回答：

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并进行减法运算。最初的回答。

- Tzoiker

---

不要使用scipy.stats.mvn，它没有文档，不适合公开使用，并且将在Scipy 2.0中被移除：https://github.com/scipy/scipy/blob/main/scipy/stats/mvn.py - undefined

Answer 2

scipy的nquad仅对有界矩形域进行数值积分。你的积分收敛的事实是由于PDF的exp(-r^2)类型权重（其显式形式请参见这里）。因此，您需要在二维空间中使用Hermite quadrature。关于这个主题已经存在一些 文章，并且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

- Nico Schlömer

---

"scipy的nquad只对有界矩形域进行数值积分" 我之前不知道，谢谢提醒！ - user192859