EM.py

import operator
import numpy as np
import numpy.random as rm
lsave=0
#theta=[pi,L11-L14,L21-L24,R11-R14,R21-R24]
A={}
for xx in range(500):
    u=np.array([[4,2,2,2],[2,4,2,2],[2,2,4,2],[2,2,2,4]])
    L1=rm.random_sample(4)
    L2=rm.random_sample(4)
    R1=rm.random_sample(4)
    R2=rm.random_sample(4)
    L1=L1/sum(L1)
    L2=L1/sum(L2)
    R1=R1/sum(R1)
    R2=R2/sum(R2)
    pi=rm.random_sample(1)[0]
    N=sum(u.reshape(16,1))[0]

    while True:
        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]
        uL=u*FL/F
        NL=sum(uL.reshape(16,1))[0]
        uR=u*FR/F
        NR=sum(uR.reshape(16,1))[0]

        pistar=sum(uL.reshape(16,1))[0]/N
        L1star=sum(uL.T)/NL
        L2star=sum(uL)/NL
        R1star=sum(uR.T)/NR
        R2star=sum(uR)/NR

        if np.abs(lobs-lsave)<.000001:
            A[lobs]=F
            break

        R1=R1star
        R2=R2star
        L2=L2star
        L1=L1star
        pi=pistar
        lsave=lobs

print 40*A[sorted(A.iterkeys(),reverse=True)[0]]