Integrating sine with Monte Carlo / Metropolis algorithm 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?
 A: 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:


*

*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.

*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.

*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.

