Integrating sine with Monte Carlo / Metropolis algorithm

422 Views Asked by At

I'm learning Monte Carlo / Metropolis algorithm, so I made up a simple question and write some code to see if I really understand it. The question is simple: integrating sine over 0 to PI. The integral can be calculated analytically. If my code is correct, the integral should be 2.

According to Monte Carlo, we can approximate an integral with N samples:

$$ \int _{a}^{b} f(x)dx \approx \frac {1}{N} \sum _{i=1}^{N} \frac {f(X_i)}{p(X_i)} $$

where in my sine integral case:

$$ f(x) = sin(x)\\ a = 0\\ b = \pi\\ $$

I'm using uniform distribution to sample $X \in [0, \pi]$, so

$$ p(x) = \frac {1}{\pi} $$

Metropolis algorithm also requires an "acceptance probability", which is the probability that whether we should accept a transition from the current location $X$ to a new location $X'$. The acceptance probability is defined as:

$$ a(X, X') = min(1, \frac {f(X')}{f(X)}) $$

My code is at http://bl.ocks.org/eliangcs/6e8b45f88fd3767363e7. Every time you refresh your browser, it makes 100 samples and shows the integral solution. But it always seems to give me a way larger value ($\approx 2.5$) instead of the correct solution: 2. Why? I guess I made a mistake on the PDF $p(x) = 1/\pi$, which does not consider the acceptance probability $a(X, X')$. If that is the case, how should I adjust $p(x)$ then?