#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()

EM 算法 實例