#coding:gbk
import math
import copy
import numpy as np
import matplotlib.pyplot as plt
isdebug = False
# 指定k個高斯分布參數，這里指定k=2。注意2個高斯分布具有相同均方差Sigma，分別為Mu1,Mu2。
def ini_data(Sigma,Mu1,Mu2,k,N):
global X
global Mu
global Expectations
X = np.zeros((1,N))
Mu = np.random.random(2)
Expectations = np.zeros((N,k))
for i in xrange(0,N):
if np.random.random(1) > 0.5:
X[0,i] = np.random.normal()*Sigma + Mu1
else:
X[0,i] = np.random.normal()*Sigma + Mu2
if isdebug:
print "***********"
print u"初始觀測數據X："
print X
# EM算法：步驟1，計算E[zij]
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 Expectations
# EM算法：步驟2，求最大化E[zij]的參數Mu
def m_step(k,N):
global Expectations
global X
for j in xrange(0,k):
Numer = 0
Denom = 0
for i in xrange(0,N):
Numer += Expectations[i,j]*X[0,i]
Denom +=Expectations[i,j]
Mu[j] = Numer / Denom
# 算法迭代iter_num次，或達到精度Epsilon停止迭代
def run(Sigma,Mu1,Mu2,k,N,iter_num,Epsilon):
ini_data(Sigma,Mu1,Mu2,k,N)
print u"初始:", Mu
for i in range(iter_num):
Old_Mu = copy.deepcopy(Mu)
e_step(Sigma,k,N)
m_step(k,N)
print i,Mu
if sum(abs(Mu-Old_Mu)) < Epsilon:
break
if __name__ == '__main__':
run(6,40,20,2,1000,1000,0.0001)
plt.hist(X[0,:],50)
plt.show()
本代碼用于模擬k=2個正態分布的均值估計。其中ini_data(Sigma,Mu1,Mu2,k,N)函數用于生成訓練樣本，此訓練樣本時從兩個高斯分布中隨機生成的，其中高斯分布a均值Mu1=40、均方差Sigma=6，高斯分布b均值Mu2=20、均方差Sigma=6，生成的樣本分布如下圖所示。由于本問題中實現無法直接沖樣本數據中獲知兩個高斯分布參數，因此需要使用EM算法估算出具體Mu1、Mu2取值。
圖 1 ?樣本數據分布
在圖1的樣本數據下，在第11步時，迭代終止，EM估計結果為：
Mu=[ 40.55261688 ?19.34252468]
附：
