Implementation of the Metropolis-Hasting algorithm for solving gaussian integrals

Question

I am currently having issue with the implementation of the Metropolis-Hastings algorithm.

I am trying to use the algorithm to calculate integrals of the form

In using this algorithm, we can obtain a long chain of configurations ( in this case, each configuration is just a single numbers) such that in the tail-end of the chain the probability of having a particular configuration follows (or rather tends to) a gaussian distribution.

My code seems to be messing up with obtaining the said gaussian distributions. There is a strange dependence on the transition probablity (the probablity of picking a new candidate configuration depending on the previous configuration in the chain). However, if this transition probability is symmetric, there should be no dependence on this function at all (it only affects speed at which phase space [space of potential configurations] is explored and how quickly the chain converges to the desired distribution)!

In my case I am using a normal distribution transition function (which satisfies the need to be symmetric), with width d. For each d I use I do indeed get a gaussian distribution however the standard deviation, sigma, depends on my choice of d. The resulting gaussian should have a sigma of roughly 0.701 but I find that the value I actually get depends on the parameter d, when it shouldn't.

I am not sure where the error in this code is, any help would be greatly appreciated!

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm


'''
We want to get an exponential decay integral approx using importance sampling.
We will try to integrate x^2exp(-x^2) over the real line.
Metropolis-hasting alg will generate configuartions (in this case, single numbers) such that 
the probablity of a given configuration x^a ~ p(x^a) for p(x) propto exp(-x^2).

Once configs  = {x^a} generated, the apporximation, Q_N, of the integral, I, will be given by 
Q_N = 1/N sum_(configs) x^2
lim (N-> inf) Q_N -> I
'''

'''
Implementing metropolis-hasting algorithm
'''

#Setting up the initial config for our chain, generating first 2 to generate numpy array
x_0 = np.random.uniform(-20,-10,2)

#Defining function that generates the next N steps in the chain, given a starting config x
#Works by iteratively taking the last element in the chain, generating a new candidate configuration from it and accepting/rejecting according to the algorithm
#Success and failures implemented to see roughly the success rate of each step
def next_steps(x,N):
    i = 0
    Success = 0
    Failures = 0
    Data = np.array(x)
    d = 1.5                       #Spread of (normal) transition function
    while i < N:
        r = np.random.uniform(0,1)
        delta = np.random.normal(0,d)
        x_old = Data[-1]
        x_new = x_old + delta
        hasting_ratio = np.exp(-(x_new**2-x_old**2) )
        if hasting_ratio > r:
            i = i+1
            Data = np.append(Data,x_new)
            Success = Success +1
        else:
            Failures = Failures + 1
    print(Success)
    print(Failures)
    return Data

#Number of steps in the chain
N_iteration   = 50000

#Generating the data
Data = next_steps(x_0,N_iteration)

#Plotting data to see convergence of chain to gaussian distribution
plt.plot(Data)
plt.show()

#Obtaining tail end data and obtaining the standard deviation of resulting gaussian distribution
Data = Data[-40000:]
(mu, sigma) = norm.fit(Data)
print(sigma)


#Plotting a histogram to visually see if guassian
plt.hist(Data, bins = 300)
plt.show()

Answer 1

You need to save x even when it doesn't change. Otherwise the center values are under-counted, and more so as d increases, which increases the variance.

import numpy as np
from scipy.stats import norm


"""
We want to get an exponential decay integral approx using importance sampling.
We will try to integrate x^2exp(-x^2) over the real line.
Metropolis-hasting alg will generate configuartions (in this case, single numbers) such that 
the probablity of a given configuration x^a ~ p(x^a) for p(x) propto exp(-x^2).

Once configs  = {x^a} generated, the apporximation, Q_N, of the integral, I, will be given by 
Q_N = 1/N sum_(configs) x^2
lim (N-> inf) Q_N -> I
"""

"""
Implementing metropolis-hasting algorithm
"""

# Setting up the initial config for our chain
x_0 = np.random.uniform(-20, -10)

# Defining function that generates the next N steps in the chain, given a starting config x
# Works by iteratively taking the last element in the chain, generating a new candidate configuration from it and accepting/rejecting according to the algorithm
# Success and failures implemented to see roughly the success rate of each step
def next_steps(x, N):
    Success = 0
    Failures = 0
    Data = np.empty((N,))
    d = 1.5  # Spread of (normal) transition function
    for i in range(N):
        r = np.random.uniform(0, 1)
        delta = np.random.normal(0, d)
        x_new = x + delta
        hasting_ratio = np.exp(-(x_new ** 2 - x ** 2))
        if hasting_ratio > r:
            x = x_new
            Success = Success + 1
        else:
            Failures = Failures + 1
        Data[i] = x
    print(Success)
    print(Failures)
    return Data


# Number of steps in the chain
N_iteration = 50000

# Generating the data
Data = next_steps(x_0, N_iteration)

# Obtaining tail end data and obtaining the standard deviation of resulting gaussian distribution
Data = Data[-40000:]
(mu, sigma) = norm.fit(Data)
print(sigma)