EM算法(講解+程序)

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

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

3、參數(shù)下的期望值。3.4 收斂性判斷將新的與舊的值進行比較，并與設置的閾值進行對比，判斷迭代是否結束，若不符合條件，則返回到3.2，重新進行下面步驟，直到最后收斂才結束。四、 算法的流程圖開始參數(shù)初始化E-StepM-step是否收斂否是結束五、 實驗結果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ù)據X的分布每次迭代期望值利用EM估

4、計的參量值與真實值比較（紅色：真實值 青綠色：估計值）六、 參考文獻1. M. Jordan. Pattern Recognition And Machine Learning2. Xiao Han. EM Algorithm七、 附錄close all;clear;clc;% 參考書籍Pattern.Recognition.and.Machine.Learning.pdf% http:/www.pr-% lwmpr-% 2009/10/15% M=2; % number of GaussianN=200; % total number of data samplesth=0.000001;

5、% convergent thresholdK=2; % demention of output signal% 待生成數(shù)據的參數(shù)a_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)&

6、#39;% for i=1:round(N*a_real(1)% 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% 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

7、)<10)% x(:,i)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)'% end% endfigure(1),plot(x(1,:),x(2,:),'.')%這里生成的數(shù)據全部符合標準% % 參數(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;

8、 rznk_p=zeros(M,N); for cm=1:M mu_cm=mu(:,cm); cov_cm=cov(:,:,cm); 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 %開始求rznk rznk=zeros(M,N);%r(Z pikn=zeros(1,M);%r(Z pikn_

9、sum=0; for cn=1:N for cm=1:M pikn(1,cm)=a(cm)*rznk_p(cm,cn);% pikn_sum=pikn_sum+pikn(1,cm); end for cm=1:M rznk(cm,cn)=pikn(1,cm)/sum(pikn); end end %求rank結束% 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_su

10、m_mu=0;%開始的時候就是錯在這里，這里要置零。 for cn=1:N rznk_sum_mu=rznk_sum_mu+rznk(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=m

11、ax(norm(a_old(:)-a(:)/norm(a_old(:);norm(mu_old(:)-mu(:)/norm(mu_old(:);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;% 輸出結果disp('a_best=');disp(a_best);disp('mu_best=');disp(mu_best);disp('cov_best=');disp(cov_b