为什么积分scipy.multivariate_normal给出不正确的概率？

问：

我正在尝试在平方区域上积分独立的二元正态分布。数值积分与蒙特卡罗模拟不匹配。这是怎么回事？

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%