博客1长于理论推导
https://zhuanlan.zhihu.com/p/149810914
EM算法的十重境界
https://www.zhihu.com/question/40797593/answer/275171156
博客2也长于理论推导
NLP —— 图模型(零):EM算法简述及简单示例(三硬币模型)
最近接触了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 inequality):
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/
标签: Maths/Statistics, ML, PGM
« 上一篇: 机器学习 —— 类不平衡问题与SMOTE过采样算法
» 下一篇: DL4NLP——词表示模型(一)表示学习;syntagmatic与paradigmatic两类模型;基于矩阵的LSA和GloVe