MHclass.py

import numpy as np
import numpy.random as rd
import matplotlib.pylab as plt
class MH:
    def __init__(self,num_samples=1000,log_like=(lambda x: -x**2),start_pt=0,proposal=(lambda : 5*rd.randn())):
        self.num_samples=num_samples
        self.samples=np.empty(num_samples)
        self.acceptances=0
        self.accepted_steps=np.empty(num_samples)
        self.log_like=log_like
        self.start_point=start_pt
        self.proposal=proposal
        x=self.start_point
        j=0
        for i in range(self.num_samples):
            x1=x+self.proposal()
            a=self.log_like(x1)-self.log_like(x)
            if a>0 or (a<0 and np.log(rd.uniform(0.0,1.0))<a):
                self.acceptances+=1
                self.accepted_steps[j]=x1-x
                x=x1
                j+=1
            self.samples[i]=x
                    


def bimodal(x):
    return(np.log(np.exp(-x**2)+np.exp(-(x-3)**2)))
            
k,axs=plt.subplots(2,5,sharex=True)
xs=np.arange(-10,10,.01)
for j in range(5):
    Z=MH(10000,log_like=bimodal)
    axs[0,j].hist(Z.samples,bins=50,normed=True)
    axs[1,j].hist(Z.accepted_steps[:Z.acceptances])

plt.show()