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
    integral = sum(f_vals) / num_points
    acceptance_rate = accepted_points / num_points
    return integral, acceptance_rate

Any pointers are appreciated.

  • 1
    $\begingroup$ 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 $\endgroup$ – Tyberius Jan 24 at 20:05
  • $\begingroup$ @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. $\endgroup$ – monstergroup42 Jan 24 at 20:13

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.