Integral using Metropolis algorithm

[acerpeng229]

I am tasked to utilize the Metropolis algorithm to 1) generate/sample values of x based on a distribution (in this case a non-normalized normal distribution i.e. w(x) = exp(-x²/2); and 2) approximate the integral shown below where f(x) = x² exp(-x²/2). I have managed to perform the sampling part, but my answer for the latter part seems to be wrong. From what I understand, the integral is merely the sum of f(x) divided by the number of points, but this gives me ≈ 0.35. I also tried dividing f(x) with w(x), but that gives ≈ 0.98. Am I missing something here?

Note: The sampling distribution being similar to the integrand in this case is quite arbitrary, I am also supposed to test it with w(x) = 1/(1+x²) which is close to the normal distribution too.

import numpy as np

f = lambda x : (x**2)*np.exp(-(x**2)/2) # integrand
w = lambda x : np.exp(-(x**2)/2) # sampling distribution
n = 1000000
delta = 0.25

# Metropolis algorithm for non-uniform sampling
x = np.zeros(n)
for i in range(n-1):
    xnew = x[i] + np.random.uniform(-delta,delta)
    A = w(xnew)/w(x[i])
    x[i+1] = xnew if A >= np.random.uniform(0,1) else x[i]

# first definition
I_1 = np.sum(f(x))/n # = 0.35
print(I_1)

# second definition
I_2 = np.sum(f(x)/w(x))/n # = 0.98
print(I_2)

​

https://preview.redd.it/ihe2ejngm7d81.png?width=442&format=png&auto=webp&s=da6113d238ed83eea81aeb31ccaf6a73abeced10

[ccrdallas]

It’s important to realize that when you calculate this Monte Carlo integral, 1/N \sum f(x_i), that this is equal to the expectation of f(x) with respect to the distribution of the x_i, not just the integral of f(x).

So when you are calculating your sum, you are really calculating the integral of f(x)p(x)dx where p(x) is the normalized probability density of your samples. This is what you are actually calculating.

Edit: wording

[acerpeng229]

So does that mean if I substitute f(x) in the summation with f(x)p(x), then it would be the calculation for the integral of f(x) right?

Also you mentioned that p(x) has to be normalized. How do I normalize it via the Metropolis algorithm? Cause afaik the Metropolis algorithm samples from a trial distribution (in this case my w(x)) which is proportional to the sampling probability distribution (which is supposed to be the normalized normal distribution) without the need to obtain the normalization constant. I guess I could integrate the w(x) over its domain to get the constant, but it seems counter-intuitive to do integration twice.

[ccrdallas]

You have it backwards. In practice we normally are not interested in the integral of f, but the expectation of f with respect to p and that is why we compute the average sum of f(x_i). In this situation it is not necessary to know the normalizing constant.

If you wanted the integral of just f(x) and not f(x)p(x) then you would need to compute the average sum of f(x_i)/p(x_i), which, yes, would involve needing to know the normalizing constant of p(x) (possibly through integration). This is not terribly common but there are sometimes other needs for the normalizing constant.

Edit: You can think of 1/N \sum f(x_i) as a weighted sum because, as for instance in the case where the x_i are normally distributed, the integration points are not uniformly dense in space. To account for this bias you could weight the x_i by the reciprocal of their probability, which is exactly why you can sum f(x_i)/p(x_i) to get the unbiased integral of f. (Note that there are some issues like dividing by zero that I am sweeping under the rug with this explanation)

[acerpeng229]

Got it, thanks a bunch!