MH.py

import numpy as np
import numpy.random as rd
import matplotlib.pylab as plt

#likelihood function
#proposal density
#Number of samples to draw (number of steps)
#starting point

#Returns --
#samples returned
#non-zero steps taken
#acceptance ratio = non-zero steps/samples



def like(x):
    try:
        z=-x+np.log(x)
    except:
        z=0
    return(z)
    
xs=np.arange(-1,10,.01)
ys=xs*np.exp(-xs)
Nsteps=10000
x0=0
accept=0
n=0

answers=np.empty(Nsteps)
accepted_steps=np.empty(Nsteps)
F=plt.figure()
stepsizes=[np.exp(i*np.log(2)) for i in range(-10,10)]
accepted=np.empty(len(stepsizes))


for stepsize in stepsizes:
    x=2
    j=0
    for i in range(Nsteps):
        step=rd.uniform(-stepsize,stepsize)
        x1=x+step
        a=like(x1)-like(x)
        if a>0:
            x=x1
            accept+=1.0
            accepted_steps[j]=step
            j+=1    
        else:
            if np.log(rd.uniform(0.0,1.0))<a:
                x=x1
                accept+=1.0
                accepted_steps[j]=step
                j+=1
        answers[i]=x
    accepted[n]=accept/Nsteps
    n+=1
    accept=0

#plt.subplot(2,1,1)
#plt.hist(answers,normed=True,bins=50,label="Mean is %1.3f\nVar is %1.3f\nRate is %0.2f" % (np.mean(answers),np.var(answers),accept/Nsteps))
#plt.plot(xs,ys)
#plt.xlim([-1,10])
#plt.ylim([-.5,1])
#plt.legend()
#plt.subplot(2,1,2)
#plt.hist(accepted_steps[:accept],normed=True,bins=20,label="Average Step is %2.2f" % (np.mean(accepted_steps[:accept])))
plt.plot(stepsizes,1/accepted,label="Accepted")
plt.legend()
plt.show()