toy_mcmc.py

#!/usr/bin/env python

""" A toy mcmc sampler for a loaded dice.
    Author: Valentin Haenel <esc@zetatech.org>

    Adapted from Matlab code presented in the course:
    'Probabilistic and Bayesian Modelling in Machine Learning and AI'
    given by Prof. Manfred Opper at TU-Berlin in Summer 2009

"""

import numpy as np

# the density we would like to sample from
target_density = [0.15, 0.15, 0.15, 0.15, 0.15, 0.25]
# total number of samples
n_samples = 100000

# the initial state
x = 1
# the collected states of the chain
states = []

for i in range(n_samples):
    # get the state transition, either -1 or 1
    s = np.sign(np.random.rand()-0.5)
    # alternatively instead of incrementing, or decrementing
    # we could switch to any of the other states
    # transition to the otehr states with prob 1/5
    #s = np.ceil(np.random.rand()*5)
    # propose a new state
    y = (x+s)%6
    # calculate the probability ratio
    r = target_density[int(y)]/target_density[int(x)]
    # if the ratio is greater than 1, accept new
    if r >= 1:
        x = y
    # if the ratio is less than 1, but larger
    # than a sample from the uniform distribution [0,1]
    # accept also
    # note that: the smaler r, the less likely this is to happen
    # meaning the acceptance probability is equal to the
    # probability ratio
    elif np.random.rand() < r:
        x = y
    # if not accepted, rejected, and we keep the current state
    states.append(x)

# make a histogram to show the samples from the target density
num = [0,0,0,0,0,0]
for s in states:
    num[int(s)] += 1
print num