EM算法(讲解+程序).doc

EM算法实验报告一、 算法简单介绍EM 算法是Dempster，Laind，Rubin于1977年提出的求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行 MLE估计，是一种非常简单实用的学习算法。这种方法可以广泛地应用于处理缺损数据、截尾数据以及带有噪声等所谓的不完全数据，可以具体来说，我们可以利用EM算法来填充样本中的缺失数据、发现隐藏变量的值、估计HMM中的参数、估计有限混合分布中的参数以及可以进行无监督聚类等等。本文主要是着重介绍EM算法在混合密度分布中的应用，如何利用EM算法解决混合密度中参数的估计。二、 算法涉及的理论我们假设X是观测的数据，并且是由某些高斯分布所生成的， X是包含的信息不完整（不清楚每个数据属于哪个高斯分布）。，此时，我们用k维二元随机变量Z（隐藏变量）来表示每一个高斯分布，将Z引入后，最终得到：， ，然而Z的后验概率满足（利用条件概率计算）：但是，Znk为隐藏变量，实际问题中我们是不知道的，所以就用Znk的期望值去估计它（利用全概率计算）。 然而我们最终是计算max：最后，我们可以得到（利用最大似然估计可以计算）：三、 算法的具体描述3.1 参数初始化对需要估计的参数进行初始赋值，包括均值、方差、混合系数以及。3.2 E-Step计算利用上面公式计算后验概率，即期望。3.3 M-step计算重新估计参数，包括均值、方差、混合系数并且估计此参数下的期望值。3.4 收敛性判断将新的与旧的值进行比较，并与设置的阈值进行对比，判断迭代是否结束，若不符合条件，则返回到3.2，重新进行下面步骤，直到最后收敛才结束。四、 算法的流程图开始参数初始化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数据X的分布每次迭代期望值利用EM估计的参量值与真实值比较（红色：真实值 青绿色：估计值）六、 参考文献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-ml.cn% lwmpr-ml.cn% 2009/10/15% M=2; % number of GaussianN=200; % total number of data samplesth=0.000001; % convergent thresholdK=2; % demention of output signal% 待生成数据的参数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);% for i=1:round(N*a_real(1)% while (x(1,i)0)&(x(2,i)0)&(x(1,i)10)&(x(2,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% endfigure(1),plot(x(1,:),x(2,:),.)%这里生成的数据全部符合标准% % 参数初始化a=1/3,2/3;mu=1 2;2 1;%均值初始化完毕cov(:,:,1)=1 0; 0 1;cov(:,:,2)=1 0; 0 1;%协方差初始化% 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); 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_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_sum_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 % 求协方差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_old(:);norm(cov_old(:)-cov(:)/norm(cov_old(:); temp_f=sum(log(sum(pikn); plot(count,temp_f,r+) count=count+1; if tth 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_best);disp(f=);disp(f);figure(3),hold onplot(x(1,:),x(2,:),.);plot(mu_real(1,:),mu_real(2,:),*r);plot(mu_best(1,:),mu_best(2,:),+c);hold off