EM算法讲解+程序.docx
《EM算法讲解+程序.docx》由会员分享，可在线阅读，更多相关《EM算法讲解+程序.docx（13页珍藏版）》请在冰点文库上搜索。
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，重新进行下面步骤，直到最后收敛才结束。
四、算法的流程图
是否收敛
否
五、实验结果
a_best=
0.80220.1978
mu_best=
2.71483.9307
4.98823.0102
cov_best=
（:
:
1）=
5.4082-0.0693
-0.06930.2184
（:
:
2）=
0.0858-0.0177
-0.01770.0769
f=
-1.6323
数据X的分布
每次迭代期望值
利用EM估计的参量值与真实值比较（红色：
真实值青绿色：
估计值）
六、参考文献
1.M.Jordan.PatternRecognitionAndMachineLearning
2.XiaoHan.EMAlgorithm
七、附录
closeall;clear;clc;
%参考书籍Pattern.Recognition.and.Machine.Learning.pdf
%http:
//www.pr-
%lwm@pr-
%2009/10/15
%%
M=2;%numberofGaussian
N=200;%totalnumberofdatasamples
th=0.000001;%convergentthreshold
K=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];
%generatethedata
x=[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
%end
figure
（1）,plot（x（1,:
）,x（2,:
）,'.'）
%这里生成的数据全部符合标准
%%%%%%%%%%%%%%%%%%参数初始化
a=[1/3,2/3];
mu=[12;21];%均值初始化完毕
cov（:
:
1）=[10;
01];
cov（:
:
2）=[10;
01];%协方差初始化
%%EMAlgorothm
%loop
count=0;
figure
（2）,holdon
while1
a_old=a;
mu_old=mu;
cov_old=cov;
rznk_p=zeros（M,N）;
forcm=1:
M
mu_cm=mu（:
cm）;
cov_cm=cov（:
:
cm）;
forcn=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）;
%Estep
%开始求rznk
rznk=zeros（M,N）;%r（Z
pikn=zeros（1,M）;%r（Z
pikn_sum=0;
forcn=1:
N
forcm=1:
M
pikn（1,cm）=a（cm）*rznk_p（cm,cn）;
%pikn_sum=pikn_sum+pikn（1,cm）;
end
forcm=1:
M
rznk（cm,cn）=pikn（1,cm）/sum（pikn）;
end
end
%求rank结束
%Mstep
nk=zeros（1,M）;
forcm=1:
M
forcn=1:
N
nk（1,cm）=nk（1,cm）+rznk（cm,cn）;
end
end
a=nk/N;
rznk_sum_mu=zeros（M,1）;
%求均值MU
forcm=1:
M
rznk_sum_mu=0;%开始的时候就是错在这里，这里要置零。
forcn=1:
N
rznk_sum_mu=rznk_sum_mu+rznk（cm,cn）*x（:
cn）;
end
mu（:
cm）=rznk_sum_mu/nk（cm）;
end
%求协方差COV
forcm=1:
M
rznk_sum_cov=zeros（K,M）;
forcn=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;
iftbreak;
end
end%while1
holdoff
f=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）,
holdon
plot（x（1,:
）,x（2,:
）,'.'）;
plot（mu_real（1,:
）,mu_real（2,:
）,'*r'）;
plot（mu_best（1,:
）,mu_best（2,:
）,'+c'）;
holdoff
升级会员 