EM算法实验报告一、算法简单介绍EM 算法是Dempster，Laind，Rubin于1977年提出的求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行MLE估计，是一种非常简单实用的学习算法。
这种方法可以广泛地应用于处理缺损数据、截尾数据以及带有噪声等所谓的不完全数据，可以具体来说，我们可以利用EM算法来填充样本中的缺失数据、发现隐藏变量的值、估计HMM中的参数、估计有限混合分布中的参数以及可以进行无监督聚类等等。
本文主要是着重介绍EM算法在混合密度分布中的应用，如何利用EM算法解决混合密度中参数的估计。
二、算法涉及的理论我们假设X是观测的数据，并且是由某些高斯分布所生成的，X是包含的信息不完整（不清楚每个数据属于哪个高斯分布）。
，此时，我们用k维二元随机变量Z（隐藏变量）来表示每一个高斯分布，将Z 引入后，最终得到：，，页脚内容1然而Z的后验概率满足（利用条件概率计算）：但是，Z nk为隐藏变量，实际问题中我们是不知道的，所以就用Z nk的期望值去估计它（利用全概率计算）。
然而我们最终是计算max：最后，我们可以得到（利用最大似然估计可以计算）：页脚内容2三、算法的具体描述3.1 参数初始化对需要估计的参数进行初始赋值，包括均值、方差、混合系数以及。
3.2 E-Step计算利用上面公式计算后验概率，即期望。
3.3 M-step计算重新估计参数，包括均值、方差、混合系数并且估计此参数下的期望值。
3.4 收敛性判断将新的与旧的值进行比较，并与设置的阈值进行对比，判断迭代是否结束，页脚内容3若不符合条件，则返回到3.2，重新进行下面步骤，直到最后收敛才结束。
页脚内容4四、算法的流程图页脚内容5五、实验结果a_best=0.8022 0.1978mu_best=2.71483.93074.9882 3.0102cov_best=(:,:,1) =5.4082 -0.0693-0.0693 0.2184(:,:,2) =页脚内容60.0858 -0.0177-0.0177 0.0769f=-1.6323数据X的分布页脚内容7每次迭代期望值利用EM估计的参量值与真实值比较（红色：真实值青绿色：估计值）页脚内容8六、参考文献1.M. Jordan. Pattern Recognition And Machine Learning2.Xiao Han. EM Algorithm七、附录close all;clear;clc;% 参考书籍Pattern.Recognition.and.Machine.Learning.pdf页脚内容9% % lwm@% 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];页脚内容10% 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)<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% endfigure(1),plot(x(1,:),x(2,:),'.')页脚内容11%这里生成的数据全部符合标准%% %%%%%%%%%%%%%%%% 参数初始化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 1a_old = a;mu_old = mu;cov_old= cov;页脚内容12rznk_p=zeros(M,N);for cm=1:Mmu_cm=mu(:,cm);cov_cm=cov(:,:,cm);for cn=1:Np_cm=exp(-0.5*(x(:,cn)-mu_cm)'/cov_cm*(x(:,cn)-mu_cm));rznk_p(cm,cn)=p_cm;endrznk_p(cm,:)=rznk_p(cm,:)/sqrt(det(cov_cm));endrznk_p=rznk_p*(2*pi)^(-K/2);%E step%开始求rznkrznk=zeros(M,N);%r(Zpikn=zeros(1,M);%r(Zpikn_sum=0;页脚内容13for cn=1:Nfor cm=1:Mpikn(1,cm)=a(cm)*rznk_p(cm,cn);% pikn_sum=pikn_sum+pikn(1,cm);endfor cm=1:Mrznk(cm,cn)=pikn(1,cm)/sum(pikn);endend%求rank结束% M stepnk=zeros(1,M);for cm=1:Mfor cn=1:Nnk(1,cm)=nk(1,cm)+rznk(cm,cn);end页脚内容14enda=nk/N;rznk_sum_mu=zeros(M,1);% 求均值MUfor cm=1:Mrznk_sum_mu=0;%开始的时候就是错在这里，这里要置零。
for cn=1:Nrznk_sum_mu=rznk_sum_mu+rznk(cm,cn)*x(:,cn);endmu(:,cm)=rznk_sum_mu/nk(cm);end% 求协方差COVfor cm=1:Mrznk_sum_cov=zeros(K,M);for cn=1:N页脚内容15rznk_sum_cov=rznk_sum_cov+rznk(cm,cn)*(x(:,cn)-mu(:,cm))*(x(:,cn)-mu(:,cm))';endcov(:,:,cm)=rznk_sum_cov/nk(cm);endt=max([norm(a_old(:)-a(:))/norm(a_old(:));norm(mu_old(:)-mu(:))/norm(mu_old(:));nor m(cov_old(:)-cov(:))/norm(cov_old(:))]);temp_f=sum(log(sum(pikn)));plot(count,temp_f,'r+')count=count+1;if t<thbreak;end页脚内容16end %while 1hold 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);页脚内容17disp('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页脚内容18。