depth first search

#!/usr/bin/python
"""
Author: Jeremy M. Stober
Program: METROPOLIS.PY
Date: Wednesday, July  9 2008
Description: An implementation of the metropolis-hastings
algorithm. Replicates the example from
http://www.inference.phy.cam.ac.uk/mackay/erice.pdf.
"""
 
import os, sys, getopt, pdb, string, random
import pylab
pylab.ion()
 
def drawp():
    """
    For comparing with direct draws from P.
    """
    return random.randint(0,20)
 
def drawq(x):
    return random.choice((x-1,x+1))
 
def p(x):
    """
    Compute the probability at sample points. Does not have to be normalized.
    """
    if x in range(21):
        return 1.0 / 21.0
    else:
        return 0.0
 
def q(x,y):
    """
    A proposal density.
    """
    if x in (y-1,y+1):
        return 0.5
    else:
        return 0.0
 
def flip(p):
    if random.random() < p:
        return True
    else:
        return False
 
def metropolis(x0, mtime):
    """
    An implementation of the Metropolis-Hastings algorithm with simple
    P and Q distributions.
    """
    xt = x0
 
    samples = [x0]
 
    for t in range(mtime):
        xn = drawq(xt)
        a = ( p(xn) / p(xt) ) * (q(xt, xn) / q(xn, xt) )
        if a >= 1:
            xt = xn # new state is accepted
        else:
            if flip(a):
                xt = xn # new state is accepted with probability a
 
        samples.append(xt)
        print xt,a
 
    return samples
 
if __name__ == "__main__":
    main()
 
s = metropolis(10,2400)
 
pylab.figure(0)
pylab.plot(s)
 
pylab.figure(1)
pylab.hist(s)
 
pylab.show()