为什么积分 scipy.multivariate_normal 会给出错误的概率？

Question

我正在尝试在方形区域上整合独立的二元正态分布。数值积分与蒙特卡罗模拟不匹配。这里出了什么问题？

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

sigmaX = 0.5
sigmaY = 0.8
region = 1.0  # Region of interest is a unit square centered at the mean (0,0)

# Numerical integration for answer:
def pdf(x,y):
    return multivariate_normal.pdf([x,y], mean=[0,0], cov=[[sigmaX, 0], [0, sigmaY]])
probability, err = integrate.nquad(pdf, [[-region/2.0, region/2.0], [-region/2.0, region/2.0]])

# Monte Carlo simulation for answer:
simulations = 1_000_000
X = np.random.normal(scale=sigmaX, size=simulations)
Y = np.random.normal(scale=sigmaY, size=simulations)
hits = sum(1 for s in range(simulations) if ((abs(X[s]) < region/2.0) and (abs(Y[s]) < region/2.0)))/simulations

print(f'Numerical integration gives probability {probability:.1%}\n'
      f'Monte Carlo gives probability {hits:.1%}')

输出：

数值积分给出的概率为 22.1%

蒙特卡洛给出的概率为 31.9%

Answer 1

协方差矩阵的对角线是方差，而

scale

的

np.random.normal

参数是*标准差。修复计算的一种方法是将

cov

函数的

pdf

参数更改为

cov=[[sigmaX**2, 0], [0, sigmaY**2]]