EM算法作业

EM 算法简单介绍及应用
EM 算法是当存在数据缺失问题时，极大似然估计（MLE ）的一种常用迭代算法，由于其每次迭代都分两步：E 步，求期望（expectation ）;M 步，求极大（maximization ），所以称之为EM 算法。

EM 算法主要适用于以下常见的两种情况的参数估计:（1）观测到的数据不完善；（2）似然函数不是显然的或者函数的形式非常复杂导致难以用极大似然传统方法进行估计。

该报告首先通过简单的实例对EM 算法的原理及其计算方法进行说明，然后简单介绍了EM 算法的的收敛性，最后就EM 算法在GMM 参数估计中的应用进行了详细的说明并通过程序实现迭代得到参数估计. 一.实例分析
其中)1,0(∈θ，现进行了197次试验，四种结果的发生次数分别为以4321,,,y y y y 表示四种结果发生的次数，此时总体分布为多项分布，故其似然函数 由此式求解θ的MLE 比较麻烦，可以考虑用EM 算法添加数据，通过引入两个潜变量21,z z ，
令1z 和11z y -分别表示落入这两部分的次数；再假设第三种结果分成两部分，其发生的概率分别为
4θ和4
1
，令2z 和23z y -分别表示落入这两部分的次数。

则在完全数据),(z y 下的对数似然函数 其对数似然为
虽然在该题目中仅知道y ，不知道z 的值，但是当y 和θ已知时，得到 下面根据EM 算法分两步进行迭代：
E 步：在已有观测数据y 和第i 步估计值)(i θ的条件下，求基于完全数据的对数似然函数的期望（即把其中与z 有关的部分积分掉）：
M 步：求),()
(i y Q θθ关于θ的最大值)1(+i θ，即找)1(+i θ使得
这样就完成了由)(i θ到)1(+i θ的一次迭代。

重复上面两步，直至收敛即可得到θ的MLE. 2.EM 算法的收敛性 算法简单、收敛稳定是EM 算法的最大优点，下面的定理说明EM 算法得到的估计序列是收敛的。

定理1：设)(P θx 为观测数据的似然函数，)(i θ为EM 算法得到的参数估计序列，)
()
(i x p θ为对应的似然函数序列，则)()
(i x p θ
是单调递增的，即）
（)()1()(i i x p x p θθ≥+. 可见在EM 算法E 步与M 步的交替运算下，都提高了观察数据的似然函数的值。

定理2：设)(ln θx p 为观测数据的对数似然函数，)(i θ为EM 算法的得到的参数估计序列，
)(ln )(i x p θ为对应的对数似然函数序列。

（1）如果)(p θx 有上界，则)(ln )(i x p θ收敛到某一值*L .
（2）在函数),('
θθQ 与)(ln θx p 满足一定条件时，由EM 算法得到的参数估计序列的)(i θ收敛值*θ是)(ln θx p 的稳定点。

根据定理可以发现，得到的收敛性主要是针对对数函数值给出，而不是针对估计序列的收敛性；而且一般情况下用EM 算法得到的估计值)(i θ，只能保证收敛到似然函数的一个稳定点，并不能保证收敛到全局最大值点。

3.EM 算法的应用—Gauss 混合分布（GMM ）的参数估计 混合高斯模型是指随机变量X 的概率密度函数为如下形式
1,
)|()|(1
M
1
==∑∑==M
k k
k k k k x p x p α
θαθ.
其中),,2121M M θθθαααθK K ，，（=即混合模型有M 个分支，每个分支的权重为k α. 设样本观测值为},,{21N x x x X K =，则GMM 的对数似然函数为 其中每个分支的分布为
我们需要估计参数),,2121M M θθθαααθK K ，，（=
，由于直接用极大似然估计不易求解，这里我们用EM 算法引入潜变量.2,1M}2,1{},,{Y 21N i y y y y i N K K K =∈=，且 当时k y i =表示第i 个观测样本是GMM 的第k 个分支产生的，则引入潜变量后的对数似然函数为∑∏====N
i y i y y i
N i i
i i i
x p y x p y x L 1
1
(ln(),(ln
)),(ln(）θα
θθ）
下面用EM 算法进行参数估计. E 步，求Q 函数 根据Bayes 公式得
M 步：极大化Q 函数，更新参数））1()
(|(max arg -=t t Q θθθθ
具体为分别将））
1(|(-t Q θ
θ对2
,k k σμ求偏导并令其为0,得到
然后利用拉格朗日乘数法得到
令M=2，观测数据为[2 5 7 1 9 3 13 17]，取均值的初始值[1 2]，方差的初始值[1 1]， 权重的初始值：[0.4 0.6]，则15次迭代后，2个分支的参数值收敛到
:
趋于稳定
%-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_total<eps
flag=0;
else
n=n+1;
end
end
miu_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’为文件名存在MA TLAB根目录下）
function z=f (x,miu,sigma,weight)
%本函数用以计算给定参数的混合高斯分布中每个高斯分布加权后的函数值% miu,sigma为高斯分布的均值与方差.，weight为权重
p=weight./(sqrt(2*pi).*sigma).*exp(-(x-miu).^2./(2*sigma.^2));
end
函数2（以‘g.m’为文件名存在MA TLAB根目录下）
function q=g(x,miu,sigma,weight)
%本函数用以计算给定参数的高斯分布的函数值
% miu,sigma为高斯分布的均值与方差.weight为权重
M=size(miu,2);
q=0;
for k=1:M
q=q+f(x,miu(k),sigma(k),weight(k));
end
end。