EM算法講解程序

1、EM算法實(shí)驗(yàn)報(bào)告一、 算法簡(jiǎn)單介紹EM 算法是Dempster，Laind，Rubin于1977年提出的求參數(shù)極大似然估計(jì)的一種方法，它可以從非完整數(shù)據(jù)集中對(duì)參數(shù)進(jìn)行 MLE估計(jì)，是一種非常簡(jiǎn)單實(shí)用的學(xué)習(xí)算法。這種方法可以廣泛地應(yīng)用于處理缺損數(shù)據(jù)、截尾數(shù)據(jù)以及帶有噪聲等所謂的不完全數(shù)據(jù)，可以具體來(lái)說(shuō)，我們可以利用EM算法來(lái)填充樣本中的缺失數(shù)據(jù)、發(fā)現(xiàn)隱藏變量的值、估計(jì)HMM中的參數(shù)、估計(jì)有限混合分布中的參數(shù)以及可以進(jìn)行無(wú)監(jiān)督聚類(lèi)等等。本文主要是著重介紹EM算法在混合密度分布中的應(yīng)用，如何利用EM算法解決混合密度中參數(shù)的估計(jì)。二、 算法涉及的理論我們假設(shè)X是觀測(cè)的數(shù)據(jù)，并且是由某些高斯分布所生成的

2、， X是包含的信息不完整（不清楚每個(gè)數(shù)據(jù)屬于哪個(gè)高斯分布）。，此時(shí)，我們用k維二元隨機(jī)變量Z（隱藏變量）來(lái)表示每一個(gè)高斯分布，將Z引入后，最終得到：， ，然而Z的后驗(yàn)概率滿足（利用條件概率計(jì)算）：但是，Znk為隱藏變量，實(shí)際問(wèn)題中我們是不知道的，所以就用Znk的期望值去估計(jì)它（利用全概率計(jì)算）。 然而我們最終是計(jì)算max：最后，我們可以得到（利用最大似然估計(jì)可以計(jì)算）：三、 算法的具體描述3.1 參數(shù)初始化對(duì)需要估計(jì)的參數(shù)進(jìn)行初始賦值，包括均值、方差、混合系數(shù)以及。3.2 E-Step計(jì)算利用上面公式計(jì)算后驗(yàn)概率，即期望。3.3 M-step計(jì)算重新估計(jì)參數(shù)，包括均值、方差、混合系數(shù)并且估計(jì)此

3、參數(shù)下的期望值。3.4 收斂性判斷將新的與舊的值進(jìn)行比較，并與設(shè)置的閾值進(jìn)行對(duì)比，判斷迭代是否結(jié)束，若不符合條件，則返回到3.2，重新進(jìn)行下面步驟，直到最后收斂才結(jié)束。四、 算法的流程圖開(kāi)始參數(shù)初始化E-StepM-step是否收斂否是結(jié)束五、 實(shí)驗(yàn)結(jié)果a_best= 0.8022 0.1978mu_best= 2.7148 3.9307 4.9882 3.0102cov_best= (:,:,1) = 5.4082 -0.0693 -0.0693 0.2184(:,:,2) = 0.0858 -0.0177 -0.0177 0.0769f= -1.6323數(shù)據(jù)X的分布每次迭代期望值利用EM估

4、計(jì)的參量值與真實(shí)值比較（紅色：真實(shí)值 青綠色：估計(jì)值）六、 參考文獻(xiàn)1. M. Jordan. Pattern Recognition And Machine Learning2. Xiao Han. EM Algorithm七、 附錄close all;clear;clc;% 參考書(shū)籍% 2009/10/15% M=2; % number of GaussianN=200; % total number of data samplesth=0.000001; % convergent thresholdK=2; % demention of output signal% 待生成數(shù)據(jù)的參數(shù)a_

5、real =4/5;1/5;mu_real=3 4; 5 3;cov_real(:,:,1)=5 0; 0 0.2;cov_real(:,:,2)=0.1 0; 0 0.1;% generate the datax= mvnrnd( mu_real(:,1) , cov_real(:,:,1) , round(N*a_real(1) )' , mvnrnd(mu_real(:,2),cov_real(:,:,2),N-round(N*a_real(1)'% for i=1:round(N*a_real(1)% while (x(1,i)>0)&&(x(2,

6、i)>0)&&(x(1,i)<10)&&(x(2,i)<10)% x(:,i)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)'% end% end% % for i=round(N*a_real(1)+1:N% while (x(1,i)>0)&&(x(2,i)>0)&&(x(1,i)<10)&&(x(2,i)<10)% x(:,i)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)'% end%

7、endfigure(1),plot(x(1,:),x(2,:),'.')%這里生成的數(shù)據(jù)全部符合標(biāo)準(zhǔn)% % 參數(shù)初始化a=1/3,2/3;mu=1 2;2 1;%均值初始化完畢cov(:,:,1)=1 0; 0 1;cov(:,:,2)=1 0; 0 1;%協(xié)方差初始化% EM Algorothm% loopcount=0;figure(2),hold onwhile 1 a_old = a; mu_old = mu; cov_old= cov; rznk_p=zeros(M,N); for cm=1:M mu_cm=mu(:,cm); cov_cm=cov(:,:,cm);

8、for cn=1:N p_cm=exp(-0.5*(x(:,cn)-mu_cm)'/cov_cm*(x(:,cn)-mu_cm); rznk_p(cm,cn)=p_cm; end rznk_p(cm,:)=rznk_p(cm,:)/sqrt(det(cov_cm); end rznk_p=rznk_p*(2*pi)(-K/2);%E step %開(kāi)始求rznk rznk=zeros(M,N);%r(Z pikn=zeros(1,M);%r(Z pikn_sum=0; for cn=1:N for cm=1:M pikn(1,cm)=a(cm)*rznk_p(cm,cn);% pikn_

9、sum=pikn_sum+pikn(1,cm); end for cm=1:M rznk(cm,cn)=pikn(1,cm)/sum(pikn); end end %求rank結(jié)束% M step nk=zeros(1,M); for cm=1:M for cn=1:N nk(1,cm)=nk(1,cm)+rznk(cm,cn); end end a=nk/N; rznk_sum_mu=zeros(M,1); % 求均值MU for cm=1:M rznk_sum_mu=0;%開(kāi)始的時(shí)候就是錯(cuò)在這里，這里要置零。 for cn=1:N rznk_sum_mu=rznk_sum_mu+rznk(

10、cm,cn)*x(:,cn); end mu(:,cm)=rznk_sum_mu/nk(cm); end % 求協(xié)方差COV for cm=1:M rznk_sum_cov=zeros(K,M); for cn=1:N rznk_sum_cov=rznk_sum_cov+rznk(cm,cn)*(x(:,cn)-mu(:,cm)*(x(:,cn)-mu(:,cm)' end cov(:,:,cm)=rznk_sum_cov/nk(cm); end t=max(norm(a_old(:)-a(:)/norm(a_old(:);norm(mu_old(:)-mu(:)/norm(mu_ol

11、d(:);norm(cov_old(:)-cov(:)/norm(cov_old(:); temp_f=sum(log(sum(pikn); plot(count,temp_f,'r+') count=count+1; if t<th break; end end %while 1 hold offf=sum(log(sum(pikn); a_best=a;mu_best=mu;cov_best=cov;f_best=f;% 輸出結(jié)果disp('a_best=');disp(a_best);disp('mu_best=');disp(mu_best);disp('cov_best=');disp(cov_best