前言
EM算法是一种迭代算法,用于含有隐变量的概率模型参数的极大似然估计,或极大后验概率估计。
EM算法的每次迭代由两步组成:
EM算法引入
首先介绍一个使用EM算法的例子:三硬币模型
假设有3枚硬币,分别记作A,B,C,这些硬币正面出现的概率分别是π,p,q,进行如下试验:
先掷硬币A根据其结果选出硬币B或者硬币C,正面选硬币B,反面选硬币C;然后掷出选出的硬币,出现正面记为1,反面记为0;独立地重复n次试验(这里n取10),观测结果如下:
1101001101
假设只能观测到掷硬币的结果,不能观测掷硬币的过程,问如何估计三个硬币正面出现的概率,即三硬币模型的参数。
三硬币模型可写为:
p(x|θ)=∑zp(x,z|θ)=∑zp(z|θ)p(x|z,θ)=πpx(1−p)1−x+(1−π)qx(1−q)1−x
这里,随机变量
x是观测变量,表示一次试验观测的结果1或0;随机变量
z是隐变量,表示未观测到的掷硬币
A的结果;
θ=(π,p,q)是模型参数。
下面将观测数据表示为X=(x1,x2,...,xn)T,未观测数据表示为Z=(z1,z2,...,zn)T,则观测数据的似然函数为:
p(X|θ)=∑zp(Z|θ)p(X|Z,θ)
即:
p(X|θ)=∏j=1n[πpxj(1−p)1−xj+(1−π)qxj(1−q)1−xj]
考虑求模型参数
θ=(π,p,q)的极大似然估计,取对数似然:
θMLE=argmaxθlog p(X|θ)
这个问题没有解析,只有通过迭代的方法来求解。EM算法就是可以用于求解这个问题的一种迭代算法。
EM算法通过迭代求L(θ)=log p(X|θ)的极大似然估计。每次迭代包含两步:
E步:求期望;M步:求极大化。
EM算法
输入: 观测变量数据X,隐变量数据Z,联合分布p(X,Z|θ),条件分布p(Z|X,θ);
输出:模型参数θ
(1)选择参数的初值θ(0),开始迭代;
(2)E步:记θ(i)为第i次迭代参数θ的估计值,在第i+1次迭代的E步,计算:
Q(θ,θ(i))=Ez[log p(X,Z|θ)|X,θ(i)]=∑zp(X,Z|θ)p(Z|X,θ(i))
这里,p(Z|X,θ(i))表示在给定观测数据X和当前的参数估计θ(i)下隐变量数据Z的条件概率分布;
(3)M步:求使Q(θ,θ(i))极大化的θ,确定第i+1次迭代的参数估计值θ(i+1)
θ(i+1)=argmaxθQ(θ,θ(i))
(4)重复第(2)步和第(3)步,直到收敛。
注:
函数Q(θ,θ(i))是EM算法的核心,称为Q函数。
参数的初值可以任意选择,但EM算法对初值是敏感的
Q(θ,θ(i))中θ表示要极大化的参数,θ(i)表示当前估计值,每次迭代实际在求Q函数及其极大化
M步求Q函数的极大化,得到θ(i+1),完成一次迭代
第(4)步给出停止迭代的条件,一般是对较小的正数E1,E2,若满足:
||θ(i+1)−θ(i)||⩽E1,或||Q(θ(i+1),θ(i))−Q(θ(i),θ(i))||⩽E2
则停止迭代。
导出EM算法
上面我们给出了EM算法,那么为什么EM算法能近似实现对观测数据的极大似然估计呢?下面通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法。
由上面内容可以知道,极大化观测数据X关于参数θ的对数似然函数,即极大化:
L(θ)=log p(X|θ)=log ∑zp(X,Z|θ)=log(∑zp(X|Z,θ)p(Z|θ))
注意到:上式中有未观测数据并包含和的对数
事实上,EM算法通过迭代逐步近似极大化L(θ)的。
假设在第i次迭代后θ的估计值是θ(i),我们希望新估计值θ能使L(θ)增加,即L(θ)>L(θ(i)),并逐步达到极大值,为此考虑两者的差:
L(θ)−L(θ(i))=log(∑zp(X|Z,θ)p(Z|θ))−log p(X|θ(i))
利用Jensen不等式,得到其下界:
log∑jλjyi⩾∑jλjlog yj
其中λj⩾0,∑jλj=1
L(θ)−L(θ(i))=log(∑zp(X|Z,θ(i))p(X|Z,θ)p(Z|θ)p(X|Z,θ(i)))−log p(X|θ(i))⩾∑zp(Z|X,θ(i))logp(X|Z,θ)p(Z|θ)p(Z|X,θ(i))−log p(X|θ(i))=∑zp(Z|X,θ(i))logp(X|Z,θ)p(Z|θ)p(Z|X,θ(i))p(X|θ(i))
令:
B(θ,θ(i))=L(θ(i))+∑zp(Z|X,θ(i))logp(X|Z,θ)p(Z|θ)p(Z|X,θ(i))p(X|θ(i))
则:
L(θ)⩾B(θ,θ(i))
即函数
B(θ,θ(i))是
L(θ)的一个下届,因此任何使
B(θ,θ(i))增大的
θ都会使
L(θ)的下届增大,也即使
L(θ)增大,为了使
L(θ)尽可能的大,自然选择使
B(θ,θ(i))达到最大的
θ(i+1),即:
θ(i+1)=argmaxθB(θ,θ(i))
现在求θ(i+1)的表达式,省去对θ的极大化而言的常数的项,有:
θ(i+1)=argmaxθ(L(θ(i))+∑zp(Z|X,θ(i))logp(X|Z,θ)p(Z|θ)p(Z|X,θ(i))p(X|θ(i)))=argmaxθ(∑zp(Z|X,θ(i))log (p(X|Z,θ)p(Z|θ)))=argmaxθ(∑zp(Z|X,θ(i))log (p(X,Z|θ))=argmaxθQ(θ,θ(i))
EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
上图给出EM算法的直观解释。图中上方曲线为L(θ),下方曲线为B(θ,θ(i)),图中两个函数在点θ=θ(i)处相等,但此时L(θ)并没用达到极大值,这时使B(θ,θ(i))极大化,由于L(θ)⩾B(θ,θ(i)),所以L(θ)在每次迭代中也是增加的。
在这个过程中,L(θ)不断增大,从图中可以看出EM算法不能保证找到全局最优解。
使用迭代必须保证p(X|θ(i))单调递增:
log p(X|θ(i+1))⩾log p(X|θ(i))
证明:
由:
p(X)=p(X,Z)p(Z|X)
得到:
E[log p(X|θ)]=E[log p(X,Z|θ)−log p(Z|X,θ)]
求期望,连续变量即求积分,离散变量即求和,上面我们用的都是求和,这个使用求积分,是为了让大家理解上面求和符号怎么来的。
上式左端:
=∫zlog p(X|θ)p(Z|X,θ(i))dz=log p(X|θ)
上式右端:
=∫zlog p(X,Z|θ)p(Z|X,θ(i))dz−∫zlog p(Z|X,θ)p(Z|X,θ(i)dz
上式减号左边的式子与上面
θ(i+1)中的式子一样,我们将其设为
Q(θ,θ(i)),减号后面一项设为
H(θ,θ(i))
那么对数似然函数(等式)可写为:
log p(X|θ)=Q(θ,θ(i))−H(θ,θ(i))
因为求得是极大似然函数估计,那么Q(θ(i+1),θ(i))⩾Q(θ,θ(i));
假如对于所有的θ都有H(θ(i),θ(i))⩾H(θ,θ(i))
=> 那么H(θ(i),θ(i))⩾H(θ(i+1),θ(i))
这样一来,对数似然函数就会逐渐增大。
那么我们又该如何证明对于所有的θ都有H(θ(g),θ(g))⩾H(θ,θ(g))
证明:
H(θ(i),θ(i))−H(θ,θ(i))⩾0=∫zlog p(Z|X,θ(i))p(Z|X,θ(i))dz−∫zlog p(Z|X,θ)p(Z|X,θ(i))dz=∫zlog(p(Z|X,θ(i))p(Z|X,θ))p(Z|X,θ(i))dz=−∫zlog(p(Z|X,θ)p(Z|X,θ(i)))p(Z|X,θ(i))dz
利用
Jensen不等式(函数的期望大于等于期望的函数):
⩾−log∫zp(Z|X,θ)p(Z|X,θ(i))p(Z|X,θ(i))dz=0
总结
EM算法对初值比较敏感,初值的选择非常重要,通常选取几个不同的初值进行迭代,然后对得到的各个估计值加以比较,从中选择最好的。
一般条件下EM算法是收敛的,但不能保证收敛到全局最优
参考资料:
李航《统计学习方法》