#coding:utf-8import mathimport copyimport numpy as npimport matplotlib.pyplot as pltisdebug = 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()