Approximating 1D integral with Metropolis - Hastings Markov Chain Monte Carlo

Question

I've been asked to approximate the integral of a one dimensional unnormalised posterior with a flat prior, using a Metropolis Hastings Markov Chain Monte Carlo, I realise that this isn't a practical endeavour, i think its more of an instructional exercise.

At the moment, I am creating a Markov chain of sampling points, $\{x_i\}$, for the distribution using the MH algorithm, and then approximating the integral of the posterior as:

\begin{equation} I \approx (\text{Prior interval}) * \frac{1}{N}\sum_i f(x_i) \end{equation}

where $f(x_i)$ is the likelihood evaluated at the point $x_i$ in the chain, and N is the length of the Markov Chain.

I have tested this algorithm in Python for numerous simple posteriors, and it converges on completely the wrong answer in every case. I've been over the code many times, and I'm almost certain there's nothing wrong with it. I also read that you need to discard the first half of the chain, and discard everything but every, say, 10th point. I've tried this and it makes no difference.

The code is shown below. It tries to integrate $\sin x$ between $0$ and $\pi$. For 20000 iterations it gives $2.45448671306 \pm 0.00502441977146$, which is obviously very wrong (answer should be 2) .The graph below shows the convergence of the algorithm

def f(x):
    if 0. <= x <= np.pi:
        return np.sin(x)
    else:
        return 0


n = 20000 #number of iterations
sigma = 0.5

x = uniform(0, np.pi) #initial x value

chain = [] #stores Markov chain 
chain.append(f(x)) #initialises chain



accept = 0 #number of accepted jumps

#generates an array of random x values from norm distribution
rands = normal(0, sigma, n) 

#Metropolis - Hastings algorithm
for i in range(1,n):
    can = x + rands[i]  #candidate for jump
    aprob = min([1., f(can)/f(x)]) #acceptance probability
    u = uniform(0,1) #rand number between 0 and 1
    if u < aprob:
        x = can
        accept += 1
    chain.append(f(x))

I = (np.pi)*sum(chain)*(1/len(chain)) #integral approximation

#uncertainty on approx.
uncertainty = (np.pi) / np.sqrt(len(chain)) * np.std(chain) 

#integral approx at each iteration
conv = np.array([(np.pi)*sum(chain[:x+1])*(1/(x+1)) for x in range(len(chain))])

print('MCMC approx: ', I, '\t','uncertainty: ', uncertainty)

Answer 1

"I've been over the code many times, and I'm almost certain there's nothing wrong with it."

There is an issue with the code, actually, in that you fail to replicate the current value of the Markov chain x when rejecting the candidate value can, while you should: thus the conditional update

   if u < aprob:
        x = can
        accept += 1
    chain.append(f(x))

should be corrected as done in the R code below:

    if (u < aprob){
        x = can
        accept = accept+1}
    chain=c(chain,f(x))}

Furthermore, there is a conceptual (and common) error in averaging the $f(x_i)$'s when the $x_i$'s are simulated from the posterior distribution, $g(x)\propto f(x)\mathbb{I}_{(0,\pi)}(x)$. Indeed, the expectation of a term in this average is$$\int_0^\pi f(x)g(x)\,\text{d}x=\int_0^\pi f(x)\times\frac{f(x)\mathbb{I}_{(0,\pi)}(x)}{\int _0^\pi f(x)\,\text{d}x}\,\text{d}x=\frac{\int_0^\pi f^2(x) \,\text{d}x}{\int _0^\pi f(x)\,\text{d}x}$$

Using the posterior sample as you do in the Monte Carlo approximation hence amounts to using the likelihood twice, i.e. integrating $f^2$ rather than $f$.

---

The Harmonic Mean Estimator

One immediate way to deduce an approximation of the integral$$\mathcal{I}=\dfrac{1}{\pi}\int_0^\pi f(x)\text{d}x=\dfrac{\cos(\pi)-\cos(0)}{\pi}=\frac{2}{\pi}$$ which is the marginal likelihood, is to use the sample $(x_i)$ from the posterior inside an harmonic mean estimate: $$\hat{\mathcal{I}}=1\big/\frac{1}{n}\sum_{i=1}^n \frac{1}{f(x_i)}$$This is a consequence of the identity$$\mathbb{E}[f(X)^{-1}]=\int_0^\pi f(x)^{-1}\frac{f(x)\mathbb{I}_{(0,\pi}(x)}{\mathcal{I}\pi}\text{d}x=\int_0^\pi\frac{\mathbb{I}_{(0,\pi}(x)}{\mathcal{I}\pi}\text{d}x=\frac{1}{\mathcal{I}}$$ The estimator returns an unbiased estimate of the inverse integral, but its major drawback is that it often has an infinite variance, hence is unreliable no matter what the number of simulations. Radford Neal has called this approach the worst Monte Carlo method ever and I concur. Only use it when you want to check that the value of the integral is what another approach tells you it should be. As in this exercise.

f <- function(x){
    if ((0<x)&(x<pi)){
        return(sin(x))}else{
        return(0)}}

n = 2000 #number of iterations
sigma = 0.5
x = runif(1,0,pi) #initial x value
chain = fx = f(x)   

#generates an array of random x values from norm distribution
rands = rnorm(n,0, sigma) 

#Metropolis - Hastings algorithm
for (i in 2:n){
    can = x + rands[i]  #candidate for jump
    fcan=f(can)
    aprob = fcan/fx #acceptance probability
    if (runif(1) < aprob){
        x = can
        fx = fcan}
    chain=c(chain,fx)}

I = pi*length(chain)/sum(1/chain) #integral harmonic approximation

Here is the corresponding result for n=1e5:

> I
[1] 2.000591

which coincides with the true value. Hence the harmonic mean did not hit too small a value of $\sin(x)$ to induce large variations in the average. But, nonetheless, $1/\sin(x)$ has an infinite integral on $(0,\pi)$ and the resulting harmonic mean estimator has an infinite variance.

An explanation for the surprisingly good performances of this "worst Monte Carlo method ever" is that sigma = 0.5 does not often return values of $\theta$ that are close to $0$ or to $\pi$. If instead one uses sigma = 0.01 in the above code, the result is quite different:

> I
[1] 3.039536

and highly variable

> I
[1] 0.5901059

Hence it cannot be trusted and should not be used.

Answer 2

The Metro-Hastings sampling procedure is over-sampling higher values due to the accept/reject criteria.

The average of the sample values returned by
sum(chain) * (1/len(chain))
is too high.

In typical method for Monte Carlo integration, random samples are drawn uniformly across the [0, pi] interval.

You can still use a Markov Chain to create "jumps". For example, changing

for i in range(1,n):
can = x + rands[i]  #candidate for jump
aprob = min([1., f(can)/f(x)]) #acceptance probability
u = uniform(0,1) #rand number between 0 and 1
if u < aprob:
    x = can
    accept += 1
chain.append(f(x))

to

for i in range(1, n):
x += rands[i]  #candidate for jump
if f(x) > 0:
    chain.append(f(x))

will get a much better approximation