MCMC.py

import numpy as np
import numpy.random as rd
import matplotlib.pyplot as mp

plist=[]
threshold=0.8
LSave=0
def gstep(v):
    stp=rd.normal(0,.01,len(v))
    return(v+stp)

def L(pi,L1,L2,R1,R2,u):
        FL=np.array([[pi*L1[i]*L2[j] for i in range(4)] for j in range(4)])
        FR=np.array([[(1-pi)*R1[i]*R2[j] for i in range(4)] for j in range(4)])
        F=FL+FR
        lobs=sum((u*np.log(F)).reshape(16,1))[0]
        return(lobs)

u=np.array([[4,2,2,2],[2,4,2,2],[2,2,4,2],[2,2,2,4]])
L1=rd.random_sample(4)
L2=rd.random_sample(4)
R1=rd.random_sample(4)
R2=rd.random_sample(4)
L1=L1/sum(L1)
L2=L1/sum(L2)
R1=R1/sum(R1)
R2=R2/sum(R2)
pi=rd.random_sample(1)[0]
N=sum(u.reshape(16,1))[0]

lsave=L(pi,L1,L2,R1,R2,u)

for xx in range(10000):
    L1new=gstep(L1)
    L1new=L1new/sum(L1new)
    L2new=gstep(L2)
    L2new=L2new/sum(L2new)
    R1new=gstep(R1)
    R1new=R1new/sum(R1new)
    R2new=gstep(R2)
    R2new=R2new/sum(R2new)
    pnew=pi+rd.normal(0,.01)
    while pnew<0 or pnew>1:
        pnew=pi+rd.normal(0,.01)
    
    loglike=L(pnew,L1new,L2new,R1new,R2new,u)
    if loglike>lsave or (loglike<lsave and rd.uniform(0,1)>threshold):
        R1=R1new
        L1=L1new
        R2=R2new
        L2=L2new
        pi=pnew
        lsave=loglike
        plist.append(pi)
        LSave+=L1
    else:
        plist.append(pi)
        LSave+=L1

print LSave/10000        
mp.hist(plist)        
mp.show()