A working implementation of random walk Metropolis-Hastings, complete with a nice plotting benchmark.

metropolis-hastings.py

"""
Here is a simple implementation of Metropolis-Hastings on a 1D probability function that is not normalized.

Read more about Metropolis-Hastings, e.g., here: https://gregorygundersen.com/blog/2019/11/02/metropolis-hastings/
"""

from matplotlib import pyplot as plt
import numpy as np
from scipy import integrate
    
def test_func(x, ell):
    
    return np.exp(-np.abs(x)/ell)*np.cos(x/ell)**2
  
sigma = 0.5

n_samples = 50000
sequence = np.zeros(n_samples)

theta_t = np.random.random()
idx = 0

# The central loop of random walk Metropolis-Hastings
while idx < n_samples:
    theta_star = np.random.normal(loc=theta_t, scale=sigma) # use theta_t as the seed for the gaussian jumping distribution
                                                            # that's what makes it random walk Metropolis-Hastings
    alpha = min(exp_func(theta_star, ell)/exp_func(theta_t, ell), 1.)
    u = np.random.random()
    sequence[idx] = theta_t
    idx += 1
    if u < alpha:
        theta_t = theta_star

        
fig, ax = plt.subplots()

theta_vals = np.linspace(-10.*ell, 10.*ell, 500)

ax.hist(sequence, bins=400, density=True)
ax.plot(theta_vals,exp_func(theta_vals, ell)/integrate.romberg(exp_func, -10*ell, 10*ell, args=[ell]))

fig.tight_layout()
fig.show()