Self-Study: What is wrong with my Metropolis Monte Carlo integration implementation?

I am trying to integrate $$e^{-x^2}$$ from $$0$$ to $$1$$ using Metropolis sampling. While plain Monte Carlo integration gives the right result, Metropolis sampling grossly underestimates it. I have tried different probability densities for the sampling. I cannot figure out where I am making a mistake in my code. My code (in Python) is as follows.

import numpy as np

def func(x):
return np.exp(-x * x)

def weight(x):
return np.exp(-x)

def mcmc_int(func, weight, num_points=1000):
current_point = 0.0
f_vals = []
accepted_points = 0
for point in range(num_points):
delta = np.random.random()
new_point = current_point + delta
acceptance_prob = min(1, weight(new_point) / weight(current_point))
if acceptance_prob > np.random.random():
current_point = new_point
accepted_points += 1
f_vals.append(func(current_point))
integral = sum(f_vals) / num_points
acceptance_rate = accepted_points / num_points
return integral, acceptance_rate


Any pointers are appreciated.

• How many points are you using? What is the acceptance rate with this weight function? If you are trying to integrate from 0 to 1, why is there no restriction on the value of current_point or new_point? Since you only ever add delta and it can be any value from 0 upto 1, current_point will likely increase considerably outside the range 0 to 1 after a small number of iterations. @monstergroup42 – Tyberius Jan 24 at 20:05
• @Tyberius You are right, new_point needs to be restricted between 0 and 1. It works now. With 1000 points, the acceptance rate is around 80% for this weight function. However I have to append func(current_point)/weight(current_point) to f_vals instead of just func(current_point) as I wrote in the original code. – monstergroup42 Jan 24 at 20:13