计算多元正态分布的概率

Question

我的问题涉及数学和计算机科学，但因为我需要一个高效的实现，所以我把它放在这里。

问题：

我有一个估计的二元正态分布，用Python的矩阵表示，但我之后还需要在Java中实现同样的计算。（这里是一些示例值）

mean = numpy.matrix([[0],[0]])
cov = numpy.matrix([[1,0],[0,1]])

当我接收到一个包含整数值的列向量（x,y）时，我想计算这个特定组合的概率。

value = numpy.matrix([[4],[3]])
probability_of_value_given_the_distribution = ???

从数学的角度来看，这个计算就是在我的正态分布的概率密度函数上，对3.5 < x < 4.5和2.5 < y < 3.5进行积分。

我想知道：

有没有办法避免实际实现这个过程，因为这涉及到矩阵上的表达式和双重积分？此外，如果我自己来实现，这会花费我很长时间，而且计算成本也很高。一个近似的解决方案对我来说完全可以接受。

我的思考：

在单变量正态分布中，可以简单地使用累积分布函数（甚至可以存储标准正态分布的值然后进行归一化），但不幸的是，似乎没有多元分布的封闭形式的累积分布函数。

另一种单变量的做法是使用二元近似的反函数（也就是把正态分布近似为二项分布），但扩展到多元分布时，我不知道如何考虑协方差。

我真的希望有人已经实现了这个，我需要尽快得到这个结果（正在完成我的论文），但我找不到任何相关的信息。

Answer 1

如果你想了解二元正态分布的概率密度函数，下面有几行代码可以帮你实现这个功能：

import numpy as np

def multivariate_pdf(vector, mean, cov):
    quadratic_form = np.dot(np.dot(vector-mean,np.linalg.inv(cov)),np.transpose(vector-mean))
    return np.exp(-.5 * quadratic_form)/ (2*np.pi * np.linalg.det(cov))

mean = np.array([0,0])
cov = np.array([[1,0],[0,1]])
vector = np.array([4,3])

pdf = multivariate_pdf(vector, mean, cov)

Answer 2

OpenTURNS 提供了一个高效的方法来计算多元正态分布的累积分布函数（CDF），你可以查看这个 代码。

import numpy as np

mean = np.array([0.0, 0.0])
cov = np.array([[1.0, 0.0],[0.0, 1.0]])

现在我们来创建一个多元正态分布，使用这些参数。

import openturns as ot
multinormal = ot.Normal(mean, ot.CovarianceMatrix(cov))

接下来，我们计算一下在这个范围内的概率：[3.5, 4.5] x |2.5, 3.5：

prob = multinormal.computeProbability(ot.Interval([3.5,2.5], [4.5,3.5])) 
print(prob)

计算得到的概率是

1.3701244220201715e-06