EMアルゴリズムの例

#coding:utf-8
import math
import copy
import numpy as np
import matplotlib.pyplot as plt

isdebug = True

#  k       ，    k=2。
#  2           Sigma，     Mu1,Mu2。
# 1000   

#      ，  6,40,20,2 
#       6，
#     20，   40
#    1000  
def ini_data(Sigma,Mu1,Mu2,k,N): 
  #         
  global X 

  #      
  global Mu 
  #           
  global Expectations 

  #1*N   ，  N   
  X = np.zeros((1,N)) 
  #         ，      
  #      ，       
  Mu = np.random.random(2) #0-1  
  print Mu
  #  1000*2   ，           
  Expectations = np.zeros((N,k)) 

  #  N      
  for i in xrange(0,N): 
    #   0.5  1   ，  0.5  2   
    if np.random.random(1) > 0.5: 
      #  40      ，      N(40,Sigma)    
      X[0,i] = np.random.normal()*Sigma + Mu1 #
    else:
      #  40      ，      N(20,Sigma)    
      X[0,i] = np.random.normal()*Sigma + Mu2 

  if isdebug:
    print "***********"
    print u"      X："
    print X


#E                 
#  ：  Sigma，  k，   N
def e_step(Sigma,k,N):
  #        
  global Expectations
  #    
  global Mu
  #  
  global X

  #       ，           
  for i in xrange(0,N):
    #  ，     
    Denom = 0
    #      ，         
    for j in xrange(0,k):
      #    
      Denom += math.exp((-1/(2*(float(Sigma**2))))*(float(X[0,i]-Mu[j]))**2)

    #      ，        
    for j in xrange(0,k):
      #  
      Numer = math.exp((-1/(2*(float(Sigma**2))))*(float(X[0,i]-Mu[j]))**2)
      #            
      Expectations[i,j] = Numer/Denom

  if isdebug:
    print "***********"
    print u"    E（Z）："
    print len(Expectations)
    #     
    print Expectations.size
    #    
    print Expectations.shape
    #         
    print Expectations


#M       
def m_step(k,N):
  #        P(k|xi)
  global Expectations
  #  
  global X
  #       
  #    
  for j in xrange(0,k):
    Numer = 0
    Denom = 0
    #     ，      
    #            
    for i in xrange(0,N):
      #       P(k|xi)xi
      Numer += Expectations[i,j]*X[0,i]
      #        Nk，Nk  P(k|xi)  
      Denom +=Expectations[i,j]
    #          uk
    Mu[j] = Numer / Denom


#    iter_num ，     Epsilon    
#    1000 ，     0.0001  
#  ：      Sigma，    Mu1，    Mu2
#   k，   N，    iter_num，     Epsilon
def run(Sigma,Mu1,Mu2,k,N,iter_num,Epsilon):
  #      
  ini_data(Sigma,Mu1,Mu2,k,N)
  print u"  <u1,u2>:", Mu

  #  1000 
  for i in range(iter_num):
    #        
    Old_Mu = copy.deepcopy(Mu)
    #E 
    e_step(Sigma,k,N)
    #M 
    m_step(k,N)

    #               
    print i,Mu

    #    
    if sum(abs(Mu-Old_Mu)) < Epsilon:
      break

if __name__ == '__main__':

  #sigma,mu1,mu2,   ，    ，    ，        
   run(6,40,20,2,1000,1000,0.0001)
   plt.hist(X[0,:],100) #      
   plt.show()