Welcome, guest | Sign In | My Account | Store | Cart

Metropolis-Hastings Sampler (Python recipe) by Flávio Codeço Coelho
ActiveState Code (http://code.activestate.com/recipes/414200/)

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.

Python, 43 lines
Download
 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.

Tags: algorithms

3 comments

Tim Vieira 14 years, 3 months ago  # | flag

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.

mudd1 10 years, 4 months ago  # | flag

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.

suhail dhawan 9 years, 7 months ago  # | flag

^ 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!