# Numerically integrate gaussian pdf

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$$ ?

• Can you edit your post to clarify what you’d like to know? What’s your question?
– Sycorax
Mar 15, 2021 at 13:33
• I've added a question, but I think it's clear from my post that I want a method which converges to the true integral value, not a method which outputs values close to $0$ (which is what I have at present). Mar 15, 2021 at 13:41
• The uniform sampling approach of Monte Carlo integration is likely a pitfall here. I would suggest to use importance sampling rather than the more naive Monte Carlo sampling in order to use your sampled values more efficiently and choose a bivariate Gaussian importance function. Mar 15, 2021 at 14:10
• What is your numerical estimate of the normalizing constant, just out of curiosity? Mar 15, 2021 at 14:57
• Shouldn't the volume be $4k^2$?
– g g
Mar 15, 2021 at 15:06

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

• I indeed made a mistake in my volume calculation ! Mar 15, 2021 at 16:37