# Small Markov chain Monte Carlo implementation

I've written a small markov chain monte carlo function that takes samples from a posterior distribution, based on a prior and a binomial (Bin(N, Z)) distribution.

I'd be happy to have it reviewed, especially perhaps, regarding how to properly pass functions as arguments to functions (as the function prior_dist() in my code). In this case, I'm passing the function uniform_prior_distribution() showed below, but it's quite likely I'd like to pass other functions, that accept slightly different arguments, in the future. This would require me to rewrite mcmc(), unless there's some smart way around it...

def mcmc(prior_dist, size=100000, burn=1000, thin=10, Z=3, N=10):
import random
from scipy.stats import binom
#Make Markov chain (Monte Carlo)
mc = [0] #Initialize markov chain
while len(mc) < thin*size + burn:
cand = random.gauss(mc[-1], 1) #Propose candidate
ratio = (binom.pmf(Z, N, cand)*prior_dist(cand, size)) / (binom.pmf(Z, N, mc[-1])*prior_dist(mc[-1], size))
if ratio > random.random(): #Acceptence criteria
mc.append(cand)
else:
mc.append(mc[-1])
#Take sample
sample = []
for i in range(len(mc)):
if i >= burn and (i-burn)%thin == 0:
sample.append(mc[i])
sample = sorted(sample)
#Estimate posterior probability
post = []
for p in sample:
post.append(binom.pmf(Z, N, p) * prior_dist(p, size))
return sample, post, mc

def uniform_prior_distribution(p, size):
prior = 1.0/size
return prior

• is it possible to replace binom.pmf with some more simple implementation Commented Jul 27, 2015 at 4:41

## Function Passing

To pass functions with possibly varying argument lists consider using argument unpacking.

By using a dictionary at the call site, we can use the arguments that are relevant to the chosen function whilst discarding those which aren't.

def uniform_prior_distribution(p, size, **args): #this args is neccesary to discard parameters not intended for this function
prior = 1.0/size
return prior

def foo_prior_distribution(p, q, size, **args):
return do_something_with(p, q, size)

def mcmc(prior_dist, size=100000, burn=1000, thin=10, Z=3, N=10):
import random
from scipy.stats import binom
#Make Markov chain (Monte Carlo)
mc = [0] #Initialize markov chain
while len(mc) < thin*size + burn:
cand = random.gauss(mc[-1], 1) #Propose candidate
args1 = {"p": cand, "q": 0.2, "size": size}
args2 = {"p": mc[-1], "q": 0.2, "size": size}
ratio = (binom.pmf(Z, N, cand)*prior_dist(**args1)) / (binom.pmf(Z, N, mc[-1])*prior_dist(**args2))


Depending on how you use the mcmc function you may want to pass the args1 and args2 dicts directly as parameters to mcmc.

However, this approach with only work if the common parameters (in this case p and size) can remain constant between different functions. Unfortunately I don't know enough about prior distributions to know if this is the case.

## Other points of note

Use

sample.sort()


sample = sorted(sample)


It's clearer and slightly more efficient (it saves a copy).

for i in range(len(mc)):
if i >= burn and (i-burn)%thin == 0:
sample.append(mc[i])


Use the more pythonic:

for i,s in enumerate(mc):
if i >= burn and (i-burn)%thin == 0:
sample.append(s)

def mcmc(prior_dist, size=100000, burn=1000, thin=10, Z=3, N=10):
import random
from scipy.stats import binom


Don't import inside functions. It goes against python convention and is inefficient

    #Make Markov chain (Monte Carlo)
mc = [0] #Initialize markov chain
while len(mc) < thin*size + burn:
cand = random.gauss(mc[-1], 1) #Propose candidate
ratio = (binom.pmf(Z, N, cand)*prior_dist(cand, size)) / (binom.pmf(Z, N, mc[-1])*prior_dist(mc[-1], size))
if ratio > random.random(): #Acceptence criteria
mc.append(cand)
else:
mc.append(mc[-1])


Why is this in else clause? It only make sense to use else clauses if you are using break.

    #Take sample
sample = []
for i in range(len(mc)):
if i >= burn and (i-burn)%thin == 0:
sample.append(mc[i])


Actually something like sample = mc[burn::thin] will have the same effect as this loop

    sample = sorted(sample)
#Estimate posterior probability
post = []
for p in sample:
post.append(binom.pmf(Z, N, p) * prior_dist(p, size))


I'd suggest post = [binom.pmf(Z, N, p) * prior_dist(p, size) for p in sample]

    return sample, post, mc

def uniform_prior_distribution(p, size):
prior = 1.0/size
return prior