在Python / R中集成多元法线

如何解决在Python / R中集成多元法线

我试图通过积分找到某个区间[a，b]上的多元法线的协方差。该函数为“ int_a ^ b int_a ^ b f（x，y）（x-mu_X）（y-mu_Y）dx dy”，其中f（x，y）=（正常pdf）/（a至b上的正常cdf）。有关公式的更多信息，请参见tmvtnorm。

当我设置a = -inf和b = inf时，我正确地获得了Cov（X_1，X_2）= 0.8。 SciPy给出相同的结果。

from mpmath import mp

mp.dps = 10
mp.pretty = True

def integrate(mu=mp.matrix([0.5,0.5]),cov=mp.matrix([[1,0.8],[0.8,2]]),a=[-1,-mp.inf],b=[0.5,4]):

    def f2(x,y):
        v = mp.matrix([x - mu[0],y - mu[1]])
        pdf = -0.5 * v.T * cov **-1 * v
        return mp.exp(pdf[0])

    c = mp.quad(f2,[a[0],b[0]],[a[1],b[1]])

    def f(x,y - mu[1]])
        pdf = -0.5 * v.T * cov **-1 * v
        return mp.exp(pdf[0]) * (x - mu[0]) * (y - mu[1])

    return mp.quad(f,b[1]]) / c

print(integrate(a=[-mp.inf,b=[mp.inf,mp.inf])) # 0.8

现在，我更改积分范围，并使用R进行验证：

install.packages("tmvtnorm")
library(tmvtnorm)

mtmvnorm(mean=c(0.5,0.5),sigma=matrix(c(1,0.8,2),2,lower=c(-1,-Inf),upper=c(0.5,4))

$tmean
[1] -0.1220891161  0.0004821634

$tvar
          [,1]      [,2]
[1,] 0.1646990 0.1312234
[2,] 0.1312234 1.4575935

但是Python返回print(integrate(a=[-1,4])) # 0.4419680296而不是0.1312234。

我的错误在哪里？

更新

正如评论中指出的那样，一个问题是我使用了正态分布的均值。修复此错误后，数值稳定性出现问题（高斯求积？）。欢迎任何评论！

from mpmath import mp
from scipy.stats import norm
import numpy as np

mp.dps = 10
mp.pretty = True

# mean formula from
# https://en.wikipedia.org/wiki/Truncated_normal_distribution
def mean(a,b,mu,sigma):
    alpha = (a - mu)/sigma
    beta = (b - mu)/sigma
    numerator = norm.pdf(alpha) - norm.pdf(beta)
    denominator = norm.cdf(beta) - norm.cdf(alpha)

    return mu + (numerator/denominator) * sigma

def integrate(mu=np.array([0.5,cov=np.array([[1,0.1],[0.1,1]]),a=[-5,-5],b=[5,5]):
    sigmas = np.sqrt(cov)
    mu_0 = mean(a[0],b[0],mu[0],sigmas[0,0])
    mu_1 = mean(a[1],b[1],mu[1],sigmas[1,1])
    cov = mp.matrix(cov)

    def f2(x,y):
        v = mp.matrix([x - mu_0,y - mu_1])
        pdf = -0.5 * v.T * cov **-1 * v
        return mp.exp(pdf[0])

    c = mp.quad(f2,y - mu_1])
        pdf = -0.5 * v.T * cov **-1 * v
        return mp.exp(pdf[0]) * (x - mu_0) * (y - mu_1)

    return mp.quad(f,b[1]]) / c,mu_0,mu_1

print(integrate(a=[-5,-3],3])) # cov=0.09541323939,R says 0.09524736
print(integrate(a=[-5+10,-3+10],b=[5+10,3+10])) # cov=0.5297249671,R says 0.0001140013
print(integrate(a=[0,0],b=[np.infty,np.infty])) # cov=0.132181971,R says 0.02536199
print(integrate(a=[-np.infty,-np.infty],np.infty])) # cov=0.1,R says 0.1