Solved – Calculating integral by using Monte Carlo method

monte carlonumerical integration

I am asked to calculate the value of the following integral by using Monte Carlo method. $$I=\int_{\mathbb{R^{10}}}(2\pi)^{-10/2}\exp\left(-\sum_{i=1}^{10}\frac{x_i^2}{2}+\sin\left(\sum_{i=1}^{10}x_i\right)\right)dx_1\dots dx_{10}$$

That seems to be a ten-dimensional integral with an exponential function as outer function and trigonometric function as inner function. My question is how can I calculate a value of this integral because it seems to be an indefinite integral and I thought that Monte Carlo method can be used only for calculating values for definite integrals. I already tried to calculate the integral function by using Maxima, but Maxima refused to calculate it.

Best Answer

Jan Kukacka's comment can be used to construct a Monte Carlo method to estimate this integrate. First note that, since the bounds on the integral are $\mathbb{R}^{10}$, traditional uniform sampling on this space is not possible. However, instead of a uniform sample from this space, we can obtain a sample from a different distribution and adjust the integrand according to the distribution chosen. That is if

$$\int_{\mathbb{R}^{10}}g(x)\,dx $$

is the integral needed, then find a distribution defined on $\mathbb{R}^{10}$ with pdf $f(x)$, so that we have,

\begin{align*} \int_{\mathbb{R}^{10}}g(x)\,dx & = \int_{\mathbb{R}^{10}}g(x) \dfrac{f(x)}{f(x)}\,dx\\ & = \int_{\mathbb{R}^{10}} \dfrac{g(x)}{f(x)} f(x)\,dx \\ & = E_f\left[ \dfrac{g(x)}{f(x)} \right] \,, \end{align*} where the expectation is with respected to the distribution described by $f$. The goal now is to find an appropriate distribution such that most of the mass of the distribution aligns with $g$. As Jan Kukacka's comment pointed out, the highest part of the integral is concentrated near zero, so we can pick a 0 centered distribution. The easiest distribution to work with then is the $N_{10}(0, I_{10})$ distribution.

  1. Obtain $N$ samples from $N_{10}(0, I_{10})$
  2. Calculate $g(x)/f(x)$ for each sample
  3. Calculate sample mean of all $g(x)/f(x)$ values. That is the Monte Carlo estimate.

Here is the R code that does it.

set.seed(10)

eval <- function(x){
  integrand <- (2*pi)^(-10/2) * exp(-sum(x^2)/2 + sin(sum(x))  )/ prod(dnorm(x))
  return(integrand)
}

N <- 100000
est <- numeric(length = N)
for(i in 1:N)
{
  samp <- rnorm(10)
  est[i] <- eval(samp)
}
mean(est) #1.263585
Related Question