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

有点跑题了，但您应该改用以下例程（它非常快）：

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 做减法。

Answer 2

scipy 的 nquad 仅在有界矩形域上进行数值积分。您的积分完全收敛的事实是由于 PDF 的exp(-r^2)-type 权重（有关其显式形式，请参阅here）。因此，您需要二维的Hermite quadrature。 Some articles 在这个主题上存在，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