EM算法

转自：http://www.cnblogs.com/Determined22/p/5776791.html

最近接触了pLSA模型，由于该模型中引入了主题作为隐变量，所以需要使用期望最大化（Expectation Maximization）算法求解。

      本文简述了以下内容：

      为什么需要EM算法

      EM算法的推导与流程

      EM算法的收敛性定理

      使用EM算法求解三硬币模型

为什么需要EM算法

      数理统计的基本问题就是根据样本所提供的信息，对总体的分布或者分布的数字特征作出统计推断。所谓总体，就是一个具有确定分布的随机变量，来自总体的每一个iid样本都是一个与总体有相同分布的随机变量。

      参数估计是指这样一类问题——总体所服从的分布类型已知，但某些参数未知：设  Y1,...,YN Y1,...,YN 是来自总体  Y Y 的iid样本，记  Y=(y1,...,yN) Y=(y1,...,yN) 是样本观测值，如果随机变量  Y1,...,YN Y1,...,YN 是可观测的，那么直接用极大似然估计法就可以估计参数  θ θ 。

      但是，如果里面含有不可观测的隐变量，使用MLE就没那么容易。EM算法正是服务于求解带有隐变量的参数估计问题。

EM算法的推导与流程

      下面考虑带有隐变量  Z Z (观测值为  z z )的参数估计问题。将观测数据（亦称不完全数据）记为  Y=(y1,...,yN) Y=(y1,...,yN) ，不可观测数据记为  Z=(z1,...,zN) Z=(z1,...,zN) ，  Y Y 、 Z Z 合在一起称为完全数据。那么观测数据的似然函数为

l(θ)=∏j=1NP(yj|θ)=∏j=1N∑zP(z|θ)P(yj|z,θ) l(θ)=∏j=1NP(yj|θ)=∏j=1N∑zP(z|θ)P(yj|z,θ)

      其中求和号表示对  z z 的所有可能取值求和。

      为了省事，表述成这个形式：

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(θ)=ln⁡P(Y|θ)=ln⁡∑zP(z|θ)P(Y|z,θ)

      EM算法是一种迭代算法，通过迭代的方式求取目标函数  L(θ)=lnP(Y|θ) L(θ)=ln⁡P(Y|θ) 的极大值。因此，希望每一步迭代之后的目标函数值会比上一步迭代结束之后的值大。设当前第  n n 次迭代后参数取值为  θ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,θ))−ln⁡P(Y|θn)

      使用Jensen不等式：

ln∑jλjyj≥∑jλjlogyj,λj≥0,∑jλj=1 ln⁡∑jλjyj≥∑jλjlog⁡yj,λ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)ln⁡P(z|θ)P(Y|z,θ)P(z|Y,θn)

      第二项有

−lnP(Y|θn)=−∑zP(z|Y,θn)lnP(Y|θn) −ln⁡P(Y|θn)=−∑zP(z|Y,θn)ln⁡P(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)ln⁡P(z|θ)P(Y|z,θ)P(z|Y,θn)−∑zP(z|Y,θn)ln⁡P(Y|θn)=∑z[P(z|Y,θn)ln⁡P(z|θ)P(Y|z,θ)P(z|Y,θn)−P(z|Y,θn)ln⁡P(Y|θn)]=∑zP(z|Y,θn)ln⁡P(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)ln⁡P(z|θ)P(Y|z,θ)P(Y|θn)P(z|Y,θn)

      从而有  L(θ)≥l(θ|θn) L(θ)≥l(θ|θn) ，也就是说第  n n 次迭代结束后，  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]，自己做了点注释）可以解释：

      需要注意的是，下界的曲线当然是随着迭代的进行而变化的：在第  i i 次迭代结束后，总是有不等式  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=arg⁡maxθ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=arg⁡maxθ∑zP(z|Y,θn)ln⁡[P(z|θ)P(Y|z,θ)]=arg⁡maxθ∑zP(z|Y,θn)ln⁡P(Y,z|θ)

      回顾一下随机变量的期望的表达式：

E[Z]=∑kP(Z=zk)zk E[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)zk E[Z|Y=y]=∑kP(Z=zk|Y=y)zk

      所以：

θn+1=argmaxθEZ|Y,θn[lnP(Y,z|θ)]=argmaxθQ(θ|θn) θn+1=arg⁡maxθEZ|Y,θn[ln⁡P(Y,z|θ)]=arg⁡maxθQ(θ|θn)

      上式定义了一个函数  Q(θ|θn) Q(θ|θn) ，称为  Q Q 函数。

      上式完整表明了EM算法中的一步迭代中所需要的两个步骤：E-step，求期望；M-step，求极大值。有了上面的铺垫，下面介绍EM算法的流程：

      输入：观测数据  Y Y ，不可观测数据  Z Z；

      输出：参数  θ θ ；

      步骤：1. 给出参数初始化值  θ0 θ0 ；

               2. E步：记  θn θn 为第  n n 次迭代后参数的估计值。在第  n+1 n+1 次迭代的E步，求期望（  Q Q 函数）

Q(θ|θn)=EZ|Y,θn[lnP(Y,z|θ)]=∑zP(z|Y,θn)lnP(Y,z|θ) Q(θ|θn)=EZ|Y,θn[ln⁡P(Y,z|θ)]=∑zP(z|Y,θn)ln⁡P(Y,z|θ)

               3. M步：求  Q Q 函数的极大值点，来作为第  n+1 n+1 次迭代所得到的参数估计值  θn+1 θn+1 

θn+1=argmaxθQ(θ|θn) θn+1=arg⁡maxθ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)=ln⁡P(Y|θn) 收敛到某一值  L∗ L∗ ；

                 (2) 在  Q Q 函数与  L(θ) L(θ) 满足一定条件下，由EM算法得到的参数估计序列  θn θn的收敛值  θ∗ θ∗ 是   L(θ) L(θ) 的稳定点。

      定理2中第二点的“条件”在大多数情况下都满足。只能保证收敛到稳定点，不能保证收敛到极大值点，因此EM算法受初值的影响较大。

使用EM算法求解三硬币模型

      参考资料[2]给出了三硬币模型的描述：

      假设有三枚硬币A、B、C，这些硬币正面出现的概率分别是  π π 、 p p 、 q q 。进行如下掷硬币试验： 先掷A，如果A是正面则再掷B，如果A是反面则再掷C。对于B或C的结果，如果是正面则记为1，如果是反面则记为0。进行  N N 次独立重复实验，得到结果。现在只能观测到结果，不能观测到掷硬币的过程，估计模型参数  θ=(π,p,q) θ=(π,p,q) 。

      在这个问题中，实验结果是可观测数据  Y=(y1,...,yn) Y=(y1,...,yn) ，硬币A的结果是不可观测数据  Z=(z1,...,zn) Z=(z1,...,zn) 且  z z 只有两种可能取值1和0。

      对于第  j j 次试验，

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−yj 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,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)ln⁡P(Y,z|θ)=∑j=1N{∑zP(z|yj,θn)ln⁡P(yj,z|θ)}=∑j=1N{P(z=1|yj,θn)ln⁡P(yj,z=1|θ)+P(z=0|yj,θn)ln⁡P(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.

      因此， Q Q 函数表达式为：

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，求  Q Q 函数的极大值：

      令  Q Q 函数对参数求导数，并等于零。

∂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] 《统计学习方法》，李航

[3] 梯度下降和EM算法：系出同源，一脉相承 这是最近扫到的一篇，还没有看，感觉不错，有时间看一下～