The Metropolis-Hastings Sampler is the most common Markov-Chain-Monte-Carlo (MCMC) algorithm used to sample from arbitrary probability density functions (PDF). Suppose you want to simulate samples from a random variable which can be described by an arbitrary PDF, i.e., any function which integrates to 1 over a given interval. This algorithm will do just that, as illustrated by the Plot done with Matplotlib. Notice how the samples follow the theoretical PDF.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
#!/usr/bin/env python # Author: Flávio Codeço Coelho # License: GPL from math import * from RandomArray import * from matplotlib.pylab import * def sdnorm(z): """ Standard normal pdf (Probability Density Function) """ return exp(-z*z/2.)/sqrt(2*pi) n = 10000 alpha = 1 x = 0. vec =  vec.append(x) innov = uniform(-alpha,alpha,n) #random inovation, uniform proposal distribution for i in xrange(1,n): can = x + innov[i] #candidate aprob = min([1.,sdnorm(can)/sdnorm(x)]) #acceptance probability u = uniform(0,1) if u < aprob: x = can vec.append(x) #plotting the results: #theoretical curve x = arange(-3,3,.1) y = sdnorm(x) subplot(211) title('Metropolis-Hastings') plot(vec) subplot(212) hist(vec, bins=30,normed=1) plot(x,y,'ro') ylabel('Frequency') xlabel('x') legend(('PDF','Samples')) show()
A good exercise is to plug your own PDF, and see how the algorithm is capable of sampling from any PDF. Another interesting thing is to try and see how fast the algorithm converges to the PDF given.
Thank you for this recipe, I found it very helpful while learning about Metropolis-Hastings.
One thing though, this is technically just a Metropolis algorithm because the acceptance ratio does not incorporate the proposal distribution. In this case, it not a problem because you have a symmetric proposal distribution (so it cancels out), but might confuse readers.
RandomArray from Numeric has been deprecated for a while now. Interestingly, the code just runs when the import RandomArray line is commented out. I don't know which line was supposed to use this package.
^ you can use numpy libraries for the proposal distribution.
import numpy as np
in line 20
innov=np.random.uniform(-alpha, alpha, n) and similarly line 24
great post. and thanks for the comment Tim, important clarification there!