缺失數(shù)據(jù)em算法python_pythonem算法的實(shí)現(xiàn)

1、缺失數(shù)據(jù)em算法python_pythonem算法的實(shí)現(xiàn)Iff數(shù)據(jù)集：偽造數(shù)據(jù)集(兩個(gè)高斯分布混合)數(shù)據(jù)集長(zhǎng)度：1000運(yùn)行結(jié)果:theParameterssetis:alpha0:0.3,mu0:0.7,sigmod0:-2.0,alpha1:0.5,mu1:0.5,sigmod1:1.0theParameterspredictis:alpha0:0.4,mu0:0.6,sigmod0:-1.7,alpha1:0.7,mu1:0.7,sigmod1:0.9importnumpyasnpimportrandomimportmathimporttimedefloadData(mu0,sigma0

2、,mu1,sigma1,alpha0,alpha1):Iff初始化數(shù)據(jù)集這里通過(guò)服從高斯分布的隨機(jī)函數(shù)來(lái)偽造數(shù)據(jù)集:parammu0:高斯0的均值:paramsigmaO:高斯0的方差:parammu1:高斯1的均值:paramsigma1:高斯1的方差:paramalpha0:高斯0的系數(shù):paramalpha1:高斯1的系數(shù):return:混合了兩個(gè)高斯分布的數(shù)據(jù)Iff#定義數(shù)據(jù)集長(zhǎng)度為1000length=1000#初始化第一個(gè)高斯分布，生成數(shù)據(jù)，數(shù)據(jù)長(zhǎng)度為length*alpha系數(shù)，以此來(lái)#滿足alpha的作用data0=np.random.normal(mu0,sigma0,int

3、(length*alpha0)#第二個(gè)高斯分布的數(shù)據(jù)data1=np.random.normal(mu1,sigma1,int(length*alpha1)初始化總數(shù)據(jù)集兩個(gè)高斯分布的數(shù)據(jù)混合后會(huì)放在該數(shù)據(jù)集中返回dataSet=#將第一個(gè)數(shù)據(jù)集的內(nèi)容添加進(jìn)去dataSet.extend(data0)添加第二個(gè)數(shù)據(jù)集的數(shù)據(jù)dataSet.extend(data1)#對(duì)總的數(shù)據(jù)集進(jìn)行打亂(其實(shí)不打亂也沒(méi)事，只不過(guò)打亂一下直觀上讓人感覺(jué)已經(jīng)混合了#讀者可以將下面這句話屏蔽以后看看效果是否有差別)random.shuffle(dataSet)#返回偽造好的數(shù)據(jù)集returndataSetdefcal

4、cGauss(dataSetArr,mu,sigmod):Iff根據(jù)高斯密度函數(shù)計(jì)算值依據(jù)：“9.3.1高斯混合模型”式9.25注：在公式中y是一個(gè)實(shí)數(shù)，但是在EM算法中(見(jiàn)算法9.2的E步)，需要對(duì)每個(gè)j都求一次yjk,在本實(shí)例中有1000個(gè)可觀測(cè)數(shù)據(jù)，因此需要計(jì)算1000次?？紤]到在E步時(shí)進(jìn)行1000次高斯計(jì)算，程序上比較不簡(jiǎn)潔，因此這里的y是向量，在numpy的exp中如果exp內(nèi)部值為向量，則對(duì)向量中每個(gè)值進(jìn)行exp,輸出仍是向量的形式。所以使用向量的形式1次計(jì)算即可將所有計(jì)算結(jié)果得出，程序上較為簡(jiǎn)潔:paramdataSetArr:可觀測(cè)數(shù)據(jù)集:parammu:均值:paramsig

5、mod:方差:return:整個(gè)可觀測(cè)數(shù)據(jù)集的高斯分布密度(向量形式)Iff計(jì)算過(guò)程就是依據(jù)式9.25寫的,沒(méi)有別的花樣result=(1/(math.sqrt(2*math.pi)*sigmod*2)*np.exp(-1*(dataSetArr-mu)*(dataSetArr-mu)/(2*sigmod*2)返回結(jié)果returnresultdefE_step(dataSetArr,alpha0,mu0,sigmod0,alpha1,mu1,sigmod1):IffEM算法中的E步依據(jù)當(dāng)前模型參數(shù)，計(jì)算分模型k對(duì)觀數(shù)據(jù)y的響應(yīng)度:paramdataSetArr:可觀測(cè)數(shù)據(jù)y:paramalph

6、aO:高斯模型0的系數(shù):parammuO:高斯模型0的均值:paramsigmodO:高斯模型0的方差:paramalpha1:高斯模型1的系數(shù):parammu1:高斯模型1的均值:paramsigmod1:高斯模型1的方差:return:兩個(gè)模型各自的響應(yīng)度Iff#計(jì)算yO的響應(yīng)度#先計(jì)算模型0的響應(yīng)度的分子gamma0=alpha0*calcGauss(dataSetArr,mu0,sigmod0)模型1響應(yīng)度的分子gamma1=alpha1*calcGauss(dataSetArr,mu1,sigmod1)#兩者相加為E步中的分布sum=gamma0+gamma1各自相除，得到兩個(gè)模型的

7、響應(yīng)度gamma0=gamma0/sumgamma1=gamma1/sum返回兩個(gè)模型響應(yīng)度returngamma0,gamma1defM_step(muo,mu1,gamma0,gamma1,dataSetArr):依據(jù)算法9.2計(jì)算各個(gè)值#這里沒(méi)什么花樣，對(duì)照書本公式看看這里就好了mu0_new=np.dot(gamma0,dataSetArr)/np.sum(gamma0)mu1_new=np.dot(gamma1,dataSetArr)/np.sum(gamma1)sigmod0_new=math.sqrt(np.dot(gamma0,(dataSetArr-muo)*2)/np.su

8、m(gamma0)sigmod1_new=math.sqrt(np.dot(gamma1,(dataSetArr-mu1)*2)/np.sum(gamma1)alpha0_new=np.sum(gamma0)/len(gamma0)alpha1_new=np.sum(gamma1)/len(gamma1)#將更新的值返回returnmu0_new,mu1_new,sigmod0_new,sigmod1_new,alpha0_new,alpha1_newdefEM_Train(dataSetList,iter=500):Iff根據(jù)EM算法進(jìn)行參數(shù)估計(jì)算法依據(jù)“9.3.2高斯混合模型參數(shù)估計(jì)的EM

9、算法”算法9.2:paramdataSetList:數(shù)據(jù)集(可觀測(cè)數(shù)據(jù)):paramiter:迭代次數(shù):return:估計(jì)的參數(shù)Iff#將可觀測(cè)數(shù)據(jù)y轉(zhuǎn)換為數(shù)組形式，主要是為了方便后續(xù)運(yùn)算dataSetArr=np.array(dataSetList)步驟1：對(duì)參數(shù)取初值，開(kāi)始迭代alpha0=0.5mu0=0sigmod0=1alpha1=0.5mu1=1sigmod1=1開(kāi)始迭代step=0while(stepiter):#每次進(jìn)入一次迭代后迭代次數(shù)加1step+=1#步驟2:E步：依據(jù)當(dāng)前模型參數(shù)，計(jì)算分模型k對(duì)觀測(cè)數(shù)據(jù)y的響應(yīng)度gamma0,gamma1=E_step(dataSetA

10、rr,alpha0,mu0,sigmod0,alpha1,mu1,sigmod1)#步驟3:M步mu0,mu1,sigmod0,sigmod1,alpha0,alpha1=M_step(mu0,mu1,gamma0,gamma1,dataSetArr)#迭代結(jié)束后將更新后的各參數(shù)返回returnalpha0,mu0,sigmod0,alpha1,mu1,sigmod1if_name_=_main_:start=time.time()#設(shè)置兩個(gè)高斯模型進(jìn)行混合，這里是初始化兩個(gè)模型各自的參數(shù)#見(jiàn)“9.3EM算法在高斯混合模型學(xué)習(xí)中的應(yīng)用”#alpha是“9.3.1高斯混合模型”定義9.2中的系數(shù)

11、a#muO是均值usigmod是方差o#在設(shè)置上兩個(gè)alpha的和必須為1,其他沒(méi)有什么具體要求，符合高斯定義就可以alphaO=0.3#系數(shù)amuO=-2#均值usigmodO=0.5#方差oalphal=0.7#系數(shù)amu1=0.5#均值usigmodi=1#方差o初始化數(shù)據(jù)集dataSetList=loadData(mu0,sigmod0,mu1,sigmod1,alpha0,alpha1)#打印設(shè)置的參數(shù)print()print(theParameterssetis:)print(alpha0:%.1f,mu0:%.1f,sigmod0:%.1f,alpha1:%.1f,mu1:%.1f,sigmod1:%.1f%(alpha0,alpha1,mu0,mu1,sigmod0,sigmod1)#開(kāi)始EM算法，進(jìn)行參數(shù)估計(jì)alpha0,mu0,sigmod0,alpha1,mu1,sigmod1=EM_Train(dataSetList)打印參數(shù)預(yù)測(cè)結(jié)果print()print(theParameterspredictis:)print(alpha0:%.1f,mu0:%.1f,si