Why does integrating scipy.multivariate_normal give incorrect probability?

Question

I'm trying to integrate an independent bivariate normal distribution over a square region. The numerical integration does not match a Monte Carlo simulation. What's going wrong here?

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%}')

Output:

Numerical integration gives probability 22.1%

Monte Carlo gives probability 31.9%

Answer 1

The diagonals of the covariance matrix are the variances, while the scale argument of np.random.normal is the standard deviation. One way to fix your calculation is to change the cov argument of the pdf function to

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