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

  • 相关阅读:
    魔数,常见魔数
    正则表达式 —— 括号与特殊字符
    串行总线 —— I2C、UART、SPI
    OpenCV-Python sift/surf特征匹配与显示
    OpenCV-Python 边缘检测
    python中zip()函数基本用法
    OpenCv-Python 图像滤波
    获取WMI硬件清单
    PowerShell查询sql server
    别名的应用(New-Alias)
  • 原文地址:https://www.cnblogs.com/lytwajue/p/7290070.html
走看看 - 开发者的网上家园