2
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ Can you edit your post to clarify what you’d like to know? What’s your question? $\endgroup$
    – Sycorax
    Mar 15, 2021 at 13:33
  • $\begingroup$ 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). $\endgroup$
    – ArnoV
    Mar 15, 2021 at 13:41
  • 1
    $\begingroup$ 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. $\endgroup$
    – Alex2006
    Mar 15, 2021 at 14:10
  • $\begingroup$ What is your numerical estimate of the normalizing constant, just out of curiosity? $\endgroup$ Mar 15, 2021 at 14:57
  • 4
    $\begingroup$ Shouldn't the volume be $4k^2$? $\endgroup$
    – g g
    Mar 15, 2021 at 15:06

1 Answer 1

2
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ I indeed made a mistake in my volume calculation ! $\endgroup$
    – ArnoV
    Mar 15, 2021 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.