最近接触了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(θ)=lnP(Y|θ)=ln∑zP(z|θ)P(Y|z,θ)
EM算法是一种迭代算法,通过迭代的方式求取目标函数
L(θ)=lnP(Y|θ)
L(θ)=lnP(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,θ))−lnP(Y|θn)
使用Jensen不等式:
ln∑jλjyj≥∑jλjlogyj,λj≥0,∑jλj=1
ln∑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) ,也就是说第
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=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)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=argmaxθEZ|Y,θn[lnP(Y,z|θ)]=argmaxθ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[lnP(Y,z|θ)]=∑zP(z|Y,θn)lnP(Y,z|θ)
3. M步:求
Q
Q 函数的极大值点,来作为第
n+1
n+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) 在
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)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.
因此,
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] 《统计学习方法》,李航
对于原创博文:如需转载请注明出处http://www.cnblogs.com/Determined22/