EM算法作业.docx

EM算法简单介绍及应用EM算法是当存在数据缺失问题时，极大似然估计（MLE）的一种常用迭代算法，由于其每次迭代都分两步：E步，求期望（expectation）;M步，求极大（maximization），所以称之为EM算法。EM算法主要适用于以下常见的两种情况的参数估计:（1）观测到的数据不完善；（2）似然函数不是显然的或者函数的形式非常复杂导致难以用极大似然传统方法进行估计。该报告首先通过简单的实例对EM算法的原理及其计算方法进行说明，然后简单介绍了EM算法的的收敛性，最后就EM算法在GMM参数估计中的应用进行了详细的说明并通过程序实现迭代得到参数估计.一.实例分析设一次实验可能有四个结果，其发生的概率分别为其中，现进行了197次试验，四种结果的发生次数分别为75,18,70,34.求的MLE.以表示四种结果发生的次数，此时总体分布为多项分布，故其似然函数 由此式求解的MLE比较麻烦，可以考虑用EM算法添加数据，通过引入两个潜变量，使得求解比较容易。现假设第一种结果可以分成两部分，其发生的概率分别为和，令和分别表示落入这两部分的次数；再假设第三种结果分成两部分，其发生的概率分别为和，令和分别表示落入这两部分的次数。则在完全数据下的对数似然函数其对数似然为虽然在该题目中仅知道，不知道的值，但是当和已知时，得到下面根据EM算法分两步进行迭代：E步：在已有观测数据和第步估计值的条件下，求基于完全数据的对数似然函数的期望（即把其中与有关的部分积分掉）： M步：求关于的最大值，即找使得这样就完成了由到的一次迭代。重复上面两步，直至收敛即可得到的MLE.2.EM算法的收敛性算法简单、收敛稳定是EM算法的最大优点，下面的定理说明EM算法得到的估计序列是收敛的。定理1：设为观测数据的似然函数，为EM算法得到的参数估计序列，为对应的似然函数序列，则是单调递增的，即.可见在EM算法E步与M步的交替运算下，都提高了观察数据的似然函数的值。定理2：设为观测数据的对数似然函数，为EM算法的得到的参数估计序列，为对应的对数似然函数序列。（1）如果有上界，则收敛到某一值.（2）在函数与满足一定条件时，由EM算法得到的参数估计序列的收敛值是的稳定点。根据定理可以发现，得到的收敛性主要是针对对数函数值给出，而不是针对估计序列的收敛性；而且一般情况下用EM算法得到的估计值，只能保证收敛到似然函数的一个稳定点，并不能保证收敛到全局最大值点。3.EM算法的应用Gauss混合分布（GMM）的参数估计混合高斯模型是指随机变量X的概率密度函数为如下形式.其中即混合模型有M个分支，每个分支的权重为.设样本观测值为，则GMM的对数似然函数为其中每个分支的分布为我们需要估计参数，由于直接用极大似然估计不易求解，这里我们用EM算法引入潜变量当表示第个观测样本是GMM的第个分支产生的，则引入潜变量后的对数似然函数为）下面用EM算法进行参数估计.E步，求函数 根据Bayes公式得M步：极大化函数，更新参数具体为分别将对求偏导并令其为0,得到然后利用拉格朗日乘数法得到令M=2，观测数据为2 5 7 1 9 3 13 17，取均值的初始值1 2，方差的初始值1 1，权重的初始值：0.4 0.6，则15次迭代后，2个分支的参数值收敛到均值1.9574 9.5344,方差0.80782 4.6821,权重0.31799 0.68201,具体迭代过程如下:迭代次数 00.40000.60001 21 110.12050.87950.91805.16140.91805.161420.18320.81681.80288.31900.75665.088730.23670.76331.85238.76030.76054.969540.27240.72761.89719.08240.78024.863650.29390.70611.92729.28820.79304.786960.30570.69431.94299.40640.80014.738770.31180.68821.95059.46950.80384.711480.31490.68511.95409.50200.80584.696990.31650.68351.95589.51830.80684.6895100.31720.68281.95669.52650.80734.6857110.31760.68241.95709.53060.80764.6839120.31780.68221.95739.53260.80774.6830130.31790.68211.95749.53360.80784.6825140.31800.68201.95749.53410.80784.6823150.31800.68201.95749.53430.80784.6821160.31800.68201.95749.53440.80784.6821170.31800.68201.95759.53450.80784.6821下图为均值、方差在迭代过程中收敛的图像，可见它们经过15次迭代后，都趋于稳定%-EM算法在混合高斯模型中的应用-参数估计x=input(请输入观测数据：);N=size(x,2);miu0=input(请输入M个均值的初始值：);sigma0=input(请输入M个方差的初始值：);weight0=input(请输入M个高斯分布的权重的初始值：);M=size(miu0,2);%分支个数eps=1e-4;%限制误差miu(1,:)=miu0;sigma(1,:)=sigma0;weight(1,:)=weight0;flag=1;%循环控制变量，等于0终止循环n=1;while flag=1 for k=1:M Nk=0;Mk=0;Wk=0; for i=1:N u(i,k)=g(x(i),miu(n,:),sigma(n,:),weight(n,:); w(i,k)=f(x(i), miu(n,k),sigma(n,k),weight(n,k); z(i,k)= w(i,k)/u(i,k); Nk=Nk+z(i,k); Mk=Mk+x(i)*z(i,k); Wk=Wk+(x(i)-miu(n,k)2)*z(i,k); end miu(n+1,k)=Mk/Nk; sigma(n+1,k)=sqrt(Wk/Nk); weight(n+1,k)=Nk/N; end %计算迭代的误差确定终止条件 e_miu=norm(miu(n+1,:)-miu(n,:); e_sigma=norm(sigma(n+1,:)-sigma(n,:); e_weight=norm(weight(n+1,:)-weight(n,:); e_total=max(e_miu e_sigma e_weight); %判断终止条件是否满足 if e_totaleps flag=0; else n=n+1; endendmiu_disp=miu(n,:);sigma_disp=sigma(n,:);weight_disp=weight(n,:);%输出结果disp(M个分支的均值：,num2str(miu_disp);disp(M个分支的方差：,num2str(sigma_disp);disp(M个分支的权重：,num2str(weight_disp);%作图for j=1:M subplot(1,2,1); plot(miu(:,j),.); hold on; subplot(1,2,2); plot(sigma(:,j),*); hold on;end函数1（以f.m为文件名存在MATLAB根目录下）function z=f (x,miu,sigma,weight)%本函数用以计算给定参数的混合高斯分布中每个高斯分布加权后的函数值% miu,sigma为高斯分布的均值与方差.，weight为权重p=weight./(sqrt(2*pi).*sigma).*exp(-(x-miu).2./(2*sigma.2);end函数