EM算法(讲解+程序).doc
EM算法(讲解+程序)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 Mstep计算重新估计参数，包括均值、方差、混合系数并且估计此参数下的期望值。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% 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（Na_real(1）)' for i=1：round（Na_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(Na_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=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(2pi）(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 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_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
