[Math] Monte Carlo Integration – Metropolis Algorithm using MATLAB

definite integralsMATLABmonte carlo

I want to compute the integral
$$ I = \int_{-\infty}^{\infty} x^2 e^{-\frac{1}{2} x^2} dx = \sqrt{2\pi} \approx
2.5066$$
using the Monte Carlo method by sampling $N$ points $\{x_i\}$ according to the normal distribution $e^{-\frac{1}{2}x^2}$. An approximation of the integral is
$$ I \approx \frac{1}{N} \sum_{i=1}^{N} x_i^2$$

I want to compute it with MATLAB. Here is my code:

targ = @(x) exp(-0.5*x^2);
proprnd = @(x) x+0.5 - 1.7*rand;
start = -0.1;
nsample = 10000;
for_sample()
histfit(x);
u = x.^2;
I = sum(u)/nsample

acc =

    0.5148


I =

    6.5270

Here, for_sample calls a .m file to create the sample points and here is its histogram, and acc is the acceptance rate of the points. enter image description here

However, you can see that the result is not $I$ but a rough approximation of $I^2$. I do not know what is wrong! So I appreciate any help from you.

Update
As mentioned, There were two mistakes, one of them is taking into account the normalization constant, and the other is that my sample does not have a zero mean. Besides putting the normalization constant, I tried to change my proprnd function to be x+0.5 - rand and it works. I don't know why?! the for_sample() function should refuse any point outside the target pdf. Can anyone see where the mistake in the following code? Also, the result is not that good; it varies from 2.4 to 2.66. 

x = [start];
count = 0;
for i = [1:1:nsample-1];
    xt = proprnd(x(i)) ;
    targ_old = targ(x(i));
    targ_new = targ(xt);
    r = targ_new / targ_old;
    if r > rand;
        x(i+1) = xt;
        count = count + 1;
    else
        x(i+1) = x(i);
    end
end
acc = count / nsample

Best Answer

In order to approximate

$$ I = \int_{-\infty}^{+\infty} x^2 \exp\Big( -\frac{1}{2}x^2 \Big) \, dx $$

rewrite it as an expectation: $I = \mathbb{E}\big[ f(X) \big]$, where $f$ is some real-valued function and $X$ a random variable whose probability distribution is defined by some density function. Here, take :

$$ f(t) = \sqrt{2\pi} t^2 $$

and $X \sim \mathcal{N}(0,1)$. As a consequence :

$$ \mathbb{E}\big[ f(X) \big] = \int_{\mathbb{R}} f(t) \varphi(t) \, dt $$

where $\varphi(t) = (\sqrt{2\pi})^{-1} \exp\big( -(1/2)t^2 \big)$ is the density function of the Gaussian distribution with zero mean and unit variance.

Results on Monte Carlo integration state that if $X_1, \ldots, X_N$ are independent samples from the $\mathcal{N}(0,1)$ distribution, then :

$$ \lim \limits_{N \to +\infty} \frac{1}{N} \sum_{i=1}^{N} f(X_i) = \mathbb{E}[ f(X) \big]. $$

In practice, take $N$ large ($N=10^5$, for example). To obtain the samples $X_1,\ldots,X_n$ use the randn function, as suggested in a comment. The Metropolis-Hastings algorithm is not needed here. One may want to use this algorithm when the density of the random variable $X$ is known only up to a normalizing constant (which, for some reason, is not tractable).

** Complement ** : By the way, there is no need for an histogram here. This won't give you information about the quality of you approximation. Instead, plot

$$ \hat{I}(N) = \frac{1}{N} \sum_{i=1}^{N} f(X_i) $$

as a function of the integer $N$. You may use the function quad to obtain a reliable approximation of $I$. Then, as $N$ goes to infinity, you'll observe that the values of $\hat{I}(N)$ come closer and closer to this value you obtained with the quad function.

Related Question