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

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%