Numerically integrate gaussian pdf

Question

Suppose I want to numerically integrate the function $g(\mathbf{x}) = \exp\bigg(-\frac12 \mathbf{x}^\mathsf{T}\mathbf{\Lambda} \mathbf{x}\,\bigg)$ to obtain the normalization constant $$\int_{\mathcal{X}}g(\mathbf{x})\mathrm{d}\mathbf{x}= \int_{\mathcal{X}}\exp\left(- \frac{1}{2}\mathbf{x}^\mathsf{T}\mathbf{\Lambda}\mathbf{x}\right)\,\mathrm{d}\mathbf{x} = (2\pi)^{d/2}/\det(\mathbf{\Lambda})$$ with $d=2$ and $\mathcal{X} = \mathbb{R}^d = \mathbb{R}^2$. My precision matrix is

$$ \mathbf{\Lambda} = \begin{bmatrix} \lambda_{11} & \lambda_{12} \\ \lambda_{21} & \lambda_{22} \end{bmatrix} = \begin{bmatrix} 7.722 & 2.774\\\ 2.774 & 1.078 \end{bmatrix} $$

My attempt was to use naive monte-carlo integration by sampling uniformly $N$ vectors $\mathbf{x}_1, \dots \mathbf{x}_N$ on the interval $[-k, +k]$ and see the behavior of the integral approximation $$I_k = \int_{-k}^{k}g(\mathbf{x})\mathrm{d}\mathbf{x} \approx \dfrac{4k}{N}\sum_{i=1}^N g(\mathbf{x}_i)$$ where $4k$ is the volume of the square $[-k,+k]^2$.

Problem

However, even for $N=10^6$, my estimates $I_k$ are near zero as I uniformly sample on large intervals a function which takes small values. The formula for the gaussian integral gives a target value near $7.92$.

My question is: What should I change in order to obtain estimates $I_k$ which slowly converge to $7.92$ ?

Answer 1

The normalization constant is about 0.126=1/7.92. That's what you need to multiple $g(x)$ by to make it a density (i.e. integrate to 1). As stated in one of the comments, you need $4k^2$. You didn't say what $k$ you were using, but I used $k=14$ and it worked OK.

library(MASS) #ginv used to get inverse of precisin matrix
k=14
N=10000
set.seed(123)
x=matrix(2*k*runif(2*N)-k,ncol=2)
sigmainv=matrix(c(7.722,2.774,2.774,1.078),2)
sigma=ginv(sigmainv)
y=(sigmainv%*%t(x))
z=rep(0,N)
for (i in 1:N) z[i]=sum(x[i,]*y[,i])

1/(mean(exp(-0.5*z))*(4*k^2))  #estimate of normalizing constant
(2*pi)^(-1)/sqrt(det(sigma))   #actual normalizing constant
sqrt(det(sigmainv))/(2*pi)     #same thing