4
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ How many samples? 100?? Most Monte Carlo converges with $1/\sqrt{n}$. $\endgroup$ – Chinny84 Jan 13 '16 at 7:40
  • $\begingroup$ Increasing sample size doesn't help. I tried increasing it to 2000. The result is the same, just with less variance. $\endgroup$ – eliang Jan 13 '16 at 7:54
  • $\begingroup$ I am getting this exact same problem, but you seem to have fixed it in your code. What was the fix? $\endgroup$ – Alan Wolfe May 14 at 19:30
2
$\begingroup$

I see that your code is "fixed" now, but that it's using independent samples each time which make it regular monte carlo integration, not markov chain monte carlo integration.

I hit the same problem you did so looked into it and it turns out that using the Metropolis algorithm for integration of a single term function like this is not straightforward.

Like me, you thought that of each sample, the $x$ component was $f(x)$ and that the $y$ component was p(x) so that you could get an estimate with $\frac{x}{y}$ aka $\frac{f(x)}{p(x)}$.

It turns out that is not correct. One reason is that the $y$ component is not from a normalized PDF and the normalization constant is unknown. For more information, check this answer out: https://stats.stackexchange.com/a/248697

I've come across three ways to use the Metropolis algorithm for integration:

  1. Using a "Harmonic Mean Estimator" which is also known as "The worst monte carlo method ever" and is not reliable as it has infinite variance.
  2. Use clever math tricks so that you don't need to know the normalization constant because it cancels out being on both the top and bottom of a division. The harmonic mean estimator does this I believe.
  3. Ultimately you need to know the normalization constant to be able to turn $f(x)$ into $p(x)$ by division of that constant. One way to calculate the normalization constant would be to count how many samples fell into a small interval $[a,b]$ and get that as a percentage of all samples gotten, calling that $C$. If you integrate the function $y=f(x)$ over that same $[a,b]$ interval, and call it $D$, the normalization constant can be estimated as $\frac{D}{C}$. A smaller interval helps for calculating $D$ but is worse for calculating $C$. An alternate idea may be to keep a histogram of the samples and use the interval of the histogram bucket that has the highest count of samples in it, with the assumption that higher counts are more accurate.
$\endgroup$

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.