最近接触了pLSA模型，由于该模型中引入了主题作为隐变量，所以需要使用期望最大化（Expectation Maximization）算法求解。
本文简述了以下内容：
为什么需要EM算法
EM算法的推导与流程
EM算法的收敛性定理
使用EM算法求解三硬币模型
为什么需要EM算法
数理统计的基本问题就是根据样本所提供的信息，对总体的分布或者分布的数字特征作出统计推断。所谓总体，就是一个具有确定分布的随机变量，来自总体的每一个iid样本都是一个与总体有相同分布的随机变量。
参数估计是指这样一类问题——总体所服从的分布类型已知，但某些参数未知：设 Y1,...,YNY1,...,YN 是来自总体 YY 的iid样本，记 Y=(y1,...,yN)Y=(y1,...,yN) 是样本观测值，如果随机变量 Y1,...,YNY1,...,YN 是可观测的，那么直接用极大似然估计法就可以估计参数 θθ 。
但是，如果里面含有不可观测的隐变量，使用MLE就没那么容易。EM算法正是服务于求解带有隐变量的参数估计问题。
EM算法的推导与流程
下面考虑带有隐变量 ZZ (观测值为 zz )的参数估计问题。将观测数据（亦称不完全数据）记为 Y=(y1,...,yN)Y=(y1,...,yN) ，不可观测数据记为 Z=(z1,...,zN)Z=(z1,...,zN) ， YY 、ZZ 合在一起称为完全数据。那么观测数据的似然函数为
l(θ)=∏j=1NP(yj|θ)=∏j=1N∑zP(z|θ)P(yj|z,θ)l(θ)=∏j=1NP(yj|θ)=∏j=1N∑zP(z|θ)P(yj|z,θ)
其中求和号表示对 zz 的所有可能取值求和。
为了省事，表述成这个形式：
l(θ)=P(Y|θ)=∑zP(z|θ)P(Y|z,θ)l(θ)=P(Y|θ)=∑zP(z|θ)P(Y|z,θ)
对数似然：
L(θ)=lnP(Y|θ)=ln∑zP(z|θ)P(Y|z,θ)L(θ)=lnP(Y|θ)=ln∑zP(z|θ)P(Y|z,θ)
EM算法是一种迭代算法，通过迭代的方式求取目标函数 L(θ)=lnP(Y|θ)L(θ)=lnP(Y|θ) 的极大值。因此，希望每一步迭代之后的目标函数值会比上一步迭代结束之后的值大。设当前第 nn 次迭代后参数取值为 θnθn ，我们的目的是使 L(θn+1)>L(θn)L(θn+1)>L(θn) 。那么考虑：
L(θ)−L(θn)=ln(∑zP(z|θ)P(Y|z,θ))−lnP(Y|θn)L(θ)−L(θn)=ln(∑zP(z|θ)P(Y|z,θ))−lnP(Y|θn)
使用Jensen不等式：
ln∑jλjyj≥∑jλjlogyj,λj≥0,∑jλj=1ln∑jλjyj≥∑jλjlogyj,λj≥0,∑jλj=1
因为 ∑zP(z|Y,θn)=1∑zP(z|Y,θn)=1 ，所以 L(θ)−L(θn)L(θ)−L(θn) 的第一项有
ln(∑zP(z|θ)P(Y|z,θ))=ln(∑zP(z|Y,θn)P(z|θ)P(Y|z,θ)P(z|Y,θn))≥∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(z|Y,θn)ln(∑zP(z|θ)P(Y|z,θ))=ln(∑zP(z|Y,θn)P(z|θ)P(Y|z,θ)P(z|Y,θn))≥∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(z|Y,θn)
第二项有
−lnP(Y|θn)=−∑zP(z|Y,θn)lnP(Y|θn)−lnP(Y|θn)=−∑zP(z|Y,θn)lnP(Y|θn)
则 L(θ)−L(θn)L(θ)−L(θn) 的下界为
L(θ)−L(θn)≥∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(z|Y,θn)−∑zP(z|Y,θn)lnP(Y|θn)=∑z[P(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(z|Y,θn)−P(z|Y,θn)lnP(Y|θn)]=∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(Y|θn)P(z|Y,θn)L(θ)−L(θn)≥∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(z|Y,θn)−∑zP(z|Y,θn)lnP(Y|θn)=∑z[P(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(z|Y,θn)−P(z|Y,θn)lnP(Y|θn)]=∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(Y|θn)P(z|Y,θn)
定义一个函数 l(θ|θn)l(θ|θn) ：
l(θ|θn)≜L(θn)+∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(Y|θn)P(z|Y,θn)l(θ|θn)≜L(θn)+∑zP(z|Y,θn)lnP(z|θ)P(Y|z,θ)P(Y|θn)P(z|Y,θn)
从而有 L(θ)≥l(θ|θn)L(θ)≥l(θ|θn) ，也就是说第 nn 次迭代结束后， L(θ)L(θ) 的一个下界为 l(θ|θn)l(θ|θn)。另外，有等式 L(θn)=l(θn|θn)L(θn)=l(θn|θn) 成立。
我们的目的是使下一次迭代后得到的目标函数值能够大于当前的值： L(θn+1)>L(θn)L(θn+1)>L(θn)，即 L(θn+1)>l(θn|θn)L(θn+1)>l(θn|θn) 。
而在当前， L(θ)L(θ) 的下界为 l(θ|θn)l(θ|θn) ，因此，任何能让 l(θ|θn)l(θ|θn) 增大的 θθ ，也可以让 L(θ)L(θ) 增大。
也就是说，能满足 l(θn+1|θn)>l(θn|θn)l(θn+1|θn)>l(θn|θn) 的 θn+1θn+1 ，一定更能满足L(θn+1)>l(θn|θn)=L(θn)L(θn+1)>l(θn|θn)=L(θn) 。
通过下图（来源：参考资料[1]，自己做了点注释）可以解释：
需要注意的是，下界的曲线当然是随着迭代的进行而变化的：在第 ii 次迭代结束后，总是有不等式 L(θ)≥l(θ|θi)L(θ)≥l(θ|θi) 和等式 L(θi)=l(θi|θi)L(θi)=l(θi|θi) 成立。
换句话说，EM算法通过优化对数似然在当前的下界，来间接优化对数似然。
ok，那么现在问题就从找满足 L(θn+1)>L(θn)L(θn+1)>L(θn) 的 θn+1θn+1 ，
转变成了找满足 l(θn+1|θn)>l(θn|θn)l(θn+1|θn)>l(θn|θn) 的 θn+1θn+1 。如何找到这样一个 θn+1θn+1 ？直接找 l(θ|θn)l(θ|θn) 的最优解呗：
θn+1=argmaxθl(θ|θn)θn+1=argmaxθl(θ|θn)
把 l(θ|θn)l(θ|θn) 中的几个与 θθ 无关的量去掉，从而有
θn+1=argmaxθ∑zP(z|Y,θn)ln[P(z|θ)P(Y|z,θ)]=argmaxθ∑zP(z|Y,θn)lnP(Y,z|θ)θn+1=argmaxθ∑zP(z|Y,θn)ln[P(z|θ)P(Y|z,θ)]=argmaxθ∑zP(z|Y,θn)lnP(Y,z|θ)
回顾一下随机变量的期望的表达式：
E[Z]=∑kP(Z=zk)zkE[Z]=∑kP(Z=zk)zk
E[g(Z)]=∑kP(Z=zk)g(zk)E[g(Z)]=∑kP(Z=zk)g(zk)
E[Z|Y=y]=∑kP(Z=zk|Y=y)zkE[Z|Y=y]=∑kP(Z=zk|Y=y)zk
所以：
θn+1=argmaxθEZ|Y,θn[lnP(Y,z|θ)]=argmaxθQ(θ|θn)θn+1=argmaxθEZ|Y,θn[lnP(Y,z|θ)]=argmaxθQ(θ|θn)
上式定义了一个函数 Q(θ|θn)Q(θ|θn) ，称为 QQ 函数。
上式完整表明了EM算法中的一步迭代中所需要的两个步骤：E-step，求期望；M-step，求极大值。有了上面的铺垫，下面介绍EM算法的流程：
输入：观测数据 YY ，不可观测数据 ZZ；
输出：参数 θθ ；
步骤：1. 给出参数初始化值 θ0θ0 ；
2. E步：记 θnθn 为第 nn 次迭代后参数的估计值。在第 n+1n+1 次迭代的E步，求期望（ QQ 函数）
Q(θ|θn)=EZ|Y,θn[lnP(Y,z|θ)]=∑zP(z|Y,θn)lnP(Y,z|θ)Q(θ|θn)=EZ|Y,θn[lnP(Y,z|θ)]=∑zP(z|Y,θn)lnP(Y,z|θ)
3. M步：求 QQ 函数的极大值点，来作为第 n+1n+1 次迭代所得到的参数估计值 θn+1θn+1
θn+1=argmaxθQ(θ|θn)θn+1=argmaxθQ(θ|θn)
4. 重复上面两步，直至达到停机条件。
EM算法的收敛性定理
定理1：观测数据的似然函数 P(Y|θ)P(Y|θ) 通过EM算法得到的序列 P(Y|θn)(n=1,2,...)P(Y|θn)(n=1,2,...) 单调递增：P(Y|θn+1)≥P(Y|θn)P(Y|θn+1)≥P(Y|θn) 。
定理2：(1) 如果 P(Y|θ)P(Y|θ) 有上界，则 L(θn)=lnP(Y|θn)L(θn)=lnP(Y|θn) 收敛到某一值 L∗L∗ ；
(2) 在 QQ 函数与 L(θ)L(θ) 满足一定条件下，由EM算法得到的参数估计序列 θnθn的收敛值 θ∗θ∗ 是 L(θ)L(θ) 的稳定点。
定理2中第二点的“条件”在大多数情况下都满足。只能保证收敛到稳定点，不能保证收敛到极大值点，因此EM算法受初值的影响较大。
使用EM算法求解三硬币模型
参考资料[2]给出了三硬币模型的描述：
假设有三枚硬币A、B、C，这些硬币正面出现的概率分别是 ππ 、pp 、qq 。进行如下掷硬币试验：
先掷A，如果A是正面则再掷B，如果A是反面则再掷C。对于B或C的结果，如果是正面则记为1，如果是反面则记为0。进行 NN 次独立重复实验，得到结果。现在只能观测到结果，不能观测到掷硬币的过程，估计模型参数 θ=(π,p,q)θ=(π,p,q) 。
在这个问题中，实验结果是可观测数据 Y=(y1,...,yn)Y=(y1,...,yn) ，硬币A的结果是不可观测数据 Z=(z1,...,zn)Z=(z1,...,zn) 且 zz 只有两种可能取值1和0。
对于第 jj 次试验，
P(yj|θ)=∑zP(yj,z|θ)=∑zP(z|θ)P(yj|z,θ)=P(z=1|θ)P(yj|z=1,θ)+P(z=0|θ)P(yj|z=0,θ)={πp+(1−π)q,π(1−p)+(1−π)(1−q),if yj=1;if yj=0.=πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yjP(yj|θ)=∑zP(yj,z|θ)=∑zP(z|θ)P(yj|z,θ)=P(z=1|θ)P(yj|z=1,θ)+P(z=0|θ)P(yj|z=0,θ)={πp+(1−π)q,if yj=1;π(1−p)+(1−π)(1−q),if yj=0.=πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj
所以有
P(Y|θ)=∏j=1NP(yj|θ)=∏j=1N(πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj)P(Y|θ)=∏j=1NP(yj|θ)=∏j=1N(πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj)
1. E-step，求期望（Q函数）：
Q(θ|θn)=∑zP(z|Y,θn)lnP(Y,z|θ)=∑j=1N{∑zP(z|yj,θn)lnP(yj,z|θ)}=∑j=1N{P(z=1|yj,θn)lnP(yj,z=1|θ)+P(z=0|yj,θn)lnP(yj,z=0|θ)}Q(θ|θn)=∑zP(z|Y,θn)lnP(Y,z|θ)=∑j=1N{∑zP(z|yj,θn)lnP(yj,z|θ)}=∑j=1N{P(z=1|yj,θn)lnP(yj,z=1|θ)+P(z=0|yj,θn)lnP(yj,z=0|θ)}
先求 P(z|yj,θn)P(z|yj,θn) ，
P(z|yj,θn)=⎧⎩⎨πnpyjn(1−pn)1−yjπnpyjn(1−pn)1−yj+(1−πn)qyjn(1−qn)1−yj=μj,n1−μj,nif z=1;if z=0.P(z|yj,θn)={πnpnyj(1−pn)1−yjπnpnyj(1−pn)1−yj+(1−πn)qnyj(1−qn)1−yj=μj,nif z=1;1−μj,nif z=0.
再求 P(yj,z|θ)=P(z|θ)P(yj|z,θ)P(yj,z|θ)=P(z|θ)P(yj|z,θ) ，
P(yj,z|θ)={πpyj(1−p)1−yj(1−π)qyj(1−q)1−yjif z=1;if z=0.P(yj,z|θ)={πpyj(1−p)1−yjif z=1;(1−π)qyj(1−q)1−yjif z=0.
因此，QQ 函数表达式为：
Q(θ|θn)=∑j=1N{μj,nln[πpyj(1−p)1−yj]+(1−μj,n)ln[(1−π)qyj(1−q)1−yj]}Q(θ|θn)=∑j=1N{μj,nln[πpyj(1−p)1−yj]+(1−μj,n)ln[(1−π)qyj(1−q)1−yj]}
2. M-step，求 QQ 函数的极大值：
令 QQ 函数对参数求导数，并等于零。
∂Q(θ|θn)∂π=∑j=1N{μj,nln[πpyj(1−p)1−yj]+(1−μj,n)ln[(1−π)qyj(1−q)1−yj]∂π}=∑j=1N{μj,npyj(1−p)1−yjπpyj(1−p)1−yj+(1−μj,n)−qyj(1−q)1−yj(1−π)qyj(1−q)1−yj}=∑j=1N{μj,n−ππ(1−π)}=(∑Nj=1μj,n)−nππ(1−π)∂Q(θ|θn)∂π=∑j=1N{μj,nln[πpyj(1−p)1−yj]+(1−μj,n)ln[(1−π)qyj(1−q)1−yj]∂π}=∑j=1N{μj,npyj(1−p)1−yjπpyj(1−p)1−yj+(1−μj,n)−qyj(1−q)1−yj(1−π)qyj(1−q)1−yj}=∑j=1N{μj,n−ππ(1−π)}=(∑j=1Nμj,n)−nππ(1−π)
∂Q(θ|θn)∂π=0∴πn+1⟹π=1n∑j=1Nμj,n=1n∑j=1Nμj,n∂Q(θ|θn)∂π=0⟹π=1n∑j=1Nμj,n∴πn+1=1n∑j=1Nμj,n
∂Q(θ|θn)∂p=∑j=1N{μj,nln[πpyj(1−p)1−yj]+(1−μj,n)ln[(1−π)qyj(1−q)1−yj]∂p}=∑j=1N{μj,nπ(yjpyj−1(1−p)1−yj+pyj(−1)(1−yj)(1−p)1−yj−1)πpyj(1−p)1−yj+0}=∑j=1N{μj,n(yj−p)p(1−p)}=(∑Nj=1μj,nyj)−(p∑Nj=1μj,n)p(1−p)∂Q(θ|θn)∂p=∑j=1N{μj,nln[πpyj(1−p)1−yj]+(1−μj,n)ln[(1−π)qyj(1−q)1−yj]∂p}=∑j=1N{μj,nπ(yjpyj−1(1−p)1−yj+pyj(−1)(1−yj)(1−p)1−yj−1)πpyj(1−p)1−yj+0}=∑j=1N{μj,n(yj−p)p(1−p)}=(∑j=1Nμj,nyj)−(p∑j=1Nμj,n)p(1−p)
∂Q(θ|θn)∂p=0∴pn+1qn+1⟹p=∑Nj=1μj,nyj∑Nj=1μj,n=∑Nj=1μj,nyj∑Nj=1μj,n=∑Nj=1(1−μj,n)yj∑Nj=1(1−μj,n)∂Q(θ|θn)∂p=0⟹p=∑j=1Nμj,nyj∑j=1Nμj,n∴pn+1=∑j=1Nμj,nyj∑j=1Nμj,nqn+1=∑j=1N(1−μj,n)yj∑j=1N(1−μj,n)
既然已经得到了三个参数的迭代式，便可给定初值，迭代求解了。
参考资料：
[1] The Expectation Maximization Algorithm: A short tutorial - Sean Borman
[2] 《统计学习方法》，李航
对于原创博文：如需转载请注明出处http://www.cnblogs.com/Determined22/
扫一扫
657
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
