EM算法(讲解+程序)

EM算法实验报告一、算法简单介绍EM算法是Dempster，Laind，Rubin于1977年提出的求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行MLE估计，是一种非常简单实用的学习算法。这种方法可以广泛地应用于处理缺损数据、截尾数据以及带有噪声等所谓的不完全数据，可以具体来说，我们可以利用EM算法来填充样本中的缺失数据、发现隐藏变量的值、估计HMM中的参数、估计有限混合分布中的参数以及可以进行无监督聚类等等。本文主要是着重介绍EM算法在混合密度分布中的应用，如何利用EM算法解决混合密度中参数的估计。二、算法涉及的理论我们假设X是观测的数据，并且是由某些高斯分布所生成的，X是包含的信息不完整（不清楚每个数据属于哪个高斯分布）。，此时，我们用k维二元随机变量Z（隐藏变量）来表示每一个高斯分布，将Z引入后，最终得到：，，然而Z的后验概率满足（利用条件概率计算）：但是，Znk为隐藏变量，实际问题中我们是不知道的，所以就用Znk的期望值去估计它（利用全概率计算）。然而我们最终是计算max：最后，我们可以得到（利用最大似然估计可以计算）：三、算法的具体描述3.1参数初始化对需要估计的参数进行初始赋值，包括均值、方差、混合系数以及。3.2E-Step计算利用上面公式计算后验概率，即期望。3.3M-step计算重新估计参数，包括均值、方差、混合系数并且估计此参数下的期望值。3.4收敛性判断将新的与旧的值进行比较，并与设置的阈值进行对比，判断迭代是否结束，若不符合条件，则返回到3.2，重新进行下面步骤，直到最后收敛才结束。四、算法的流程图结束参数初始化是否收敛开始M-stepE-Step是否五、实验结果a_best=0.80220.1978mu_best=2.71483.93074.98823.0102cov_best=(:,:,1)=5.4082-0.0693-0.06930.2184(:,:,2)=0.0858-0.0177-0.01770.0769f=-1.6323-5051022.533.544.555.566.5数据X的分布0510152025-5.5-5-4.5-4-3.5-3-2.5-2-1.5每次迭代期望值-5051022.533.544.555.566.5利用EM估计的参量值与真实值比较（红色：真实值青绿色：估计值）六、参考文献1.M.Jordan.PatternRecognitionAndMachineLearning2.XiaoHan.EMAlgorithm七、附录closeall;clear;clc;%参考书籍[e-mail]%2009/10/15%%M=2;%numberofGaussianN=200;%totalnumberofdatasamplesth=0.000001;%convergentthresholdK=2;%dementionofoutputsignal%待生成数据的参数a_real=[4/5;1/5];mu_real=[34;53];cov_real(:,:,1)=[50;00.2];cov_real(:,:,2)=[0.10;00.1];%generatethedatax=[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)))'];%fori=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%%fori=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,:),'.')%这里生成的数据全部符合标准%%%%%%%%%%%%%%%%%%参数初始化a=[1/3,2/3];mu=[12;21];%均值初始化完毕cov(:,:,1)=[10;01];cov(:,:,2)=[10;01];%协方差初始化%%EMAlgorothm%loopcount=0;figure(2),holdonwhile1a_old=a;mu_old=mu;cov_old=cov;rznk_p=zeros(M,N);forcm=1:Mmu_cm=mu(:,cm);cov_cm=cov(:,:,cm);forcn=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);%Estep%开始求rznkrznk=zeros(M,N);%r(Zpikn=zeros(1,M);%r(Zpikn_sum=0;forcn=1:Nforcm=1:Mpikn(1,cm)=a(cm)*rznk_p(cm,cn);%pikn_sum=pikn_sum+pikn(1,cm);endforcm=1:Mrznk(cm,cn)=pikn(1,cm)/sum(pikn);endend%求rank结束%Mstepnk=zeros(1,M);forcm=1:Mforcn=1:Nnk(1,cm)=nk(1,cm)+rznk(cm,cn);endenda=nk/N;rznk_sum_mu=zeros(M,1);%求均值MUforcm=1:Mrznk_sum_mu=0;%开始的时候就是错在这里，这里要置零。forcn=1:Nrznk_sum_mu=rznk_sum_mu+rznk(cm,cn)*x(:,cn);endmu(:,cm)=rznk_sum_mu/nk(cm);end%求协方差COVforcm=1:Mrznk_sum_cov=zeros(K,M);forcn=1:Nrznk_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(:));norm(cov_old(:)-cov(:))/norm(cov_old(:))]);temp_f=sum(log(sum(pikn)));plot(count,temp_f,'r+')count=count+1;iftthbreak;endend%while1holdofff=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),holdonplot(x(1,:),x(2,:),'.');plot(mu_real(1,:),mu_real(2,:),'*r');plot(mu_best(1,:),mu_best(2,:),'+c');holdoffclearallclcx1(1)=0.2000;x2(1)=0.8000;y1(1)=2.0000;y2(1)=1.0000;z1(1)=0.0000;z2(1)=2;N=1000;forj=1:Nifrand(1)0.4a(j)=random('Normal',0,1,1,1)elsea(j)=random('Normal',2,2,1,1)endendfori=1:50W1=0;W2=0;W3=0;W4=0;forj=1:Np1=(normpdf(a(j),y1(i),z1(i)))*x1(i);p2=(normpdf(a(j),y2(i),z2(i)))*x2(i);w1(j)=p1/(p1+p2);w2(j)=p1/(p1+p2);W1=w1(j)+W1;W2=w2(j)+W2;W3=(w1(j)*a(j)+W3);W4=(w2(j)*a(j)+W4);endx1(i+1)=W1/N;x2(i+1)=W2/N;y1(i+1)=(W3/(N*x1(i+1)));y2(i+1)=(W4/(N*x2(i+1)));W5=0;W6=0;forj=1:NW5=(w1(j)*(a(j)-y1(i+1).^2)+W5);W6=(w2(j)*(a(j)-y2(i+1).^2)+W6);endz1(i+1)=sqrt(W5/(N*x1(i+1)));z2(i+1)=sqrt(W6/(N*x2(i+1)));end