如果概率模型都是观测变量,那么给定数据就可以用极大似然估计法或者贝叶斯估计发去获得模型。但是,有时候概率模型既有观测变量,又有隐变量或者潜在变量,这样就不能用这些估计方法了。本文要介绍的EM算法可以解决这类问题。
三硬币模型
为了更好的描述EM算法,先提出一个实际的问题。问题描述如下:
假设有3枚硬币,分别记作A,B和C。 这些硬币正面向上的概率分别是 π \pi π、 p p p和 q q q 。进行如下抛硬币试验:先抛硬币A, 根据其结果选出硬币B或者硬币C,正面选硬币B,反面选硬币C;然后掷选出的硬币,抛硬币的结果,出现正面记作1,出现反面记作0;独立重复10次试验。结果如下:
1,1,0,1,0,0,1,0,1,1
假设只能观测到掷硬币的结果,不能观测掷硬币的过程。问如何估计三硬币正面出现的概率。
我们用 θ = ( π , p , q ) \theta = (\pi,p,q) θ=(π,p,q)表示模型的参数,y表示某次得到的0或者1的观测结果,z表示隐变量,表示某次A的结果。一般地,用Y表示观测变量数据,用Z表示隐变量的数据。在上面的问题中, Y = ( y 1 , y 2 , … , y n ) T Y=(y_1, y_2,…,y_n)^T Y=(y1,y2,…,yn)T表示观测到硬币B、C的0,1序列, Z = ( z 1 , z 2 , … , z n ) T Z=(z_1,z_2,…,z_n)^T Z=(z1,z2,…,zn)T表示未观测到的硬币A的0,1序列。
考虑用极大似然估计求模型参数 θ ^ = ( π , p , q ) \hat \theta=(\pi,p,q) θ^=(π,p,q)的, 即:
L ( θ ) = log P ( Y ∣ θ ) = log ∑ Z P ( Y , Z ∣ θ ) = ∑ j = 1 10 l o g ∑ i = 1 2 P ( y j , z i ∣ θ ) (1) \begin{aligned} L(\theta) &= \log P(Y|\theta) \\ &= \log \sum\limits_{Z}P(Y,Z|\theta) \\&= \sum_{j=1}^{10}log\sum_{i=1}^{2} P(y_j,z_i|\theta) \end{aligned} \tag 1 L(θ)=logP(Y∣θ)=logZ∑P(Y,Z∣θ)=j=1∑10logi=1∑2P(yj,zi∣θ)(1)
θ ^ = arg max θ L ( θ ) (2) \hat \theta = \arg\max\limits_{\theta} L(\theta) \tag 2 θ^=argθmaxL(θ)(2)
直接求(2)的参数θ会比较困难,因为上式存在隐含随机变量Z。如果Z是个已知的数,那么使用极大似然估计来估算会很容易。在这种Z不确定的情形下,EM算法就派上用场了。
EM推导背景
对于式子1,我们做如下变形:
L ( θ ) = ∑ j = 1 10 l o g ∑ i = 1 2 P ( y j , z i ∣ θ ) = ∑ j = 1 10 l o g ∑ i = 1 2 D j ( z i ) P ( y j , z i ∣ θ ) D j ( z i ) ( 3 ) ≥ ∑ j = 1 10 ∑ i = 1 2 D j ( z i ) l o g P ( y j , z i ∣ θ ) D j ( z i ) ( 4 ) \begin{aligned}L(\theta) &= \sum_{j=1}^{10}log\sum_{i=1}^{2} P(y_j,z_i|\theta) \\&=\sum_{j=1}^{10}log\sum_{i=1}^{2}D_j(z_i)\frac{P(y_j,z_i|\theta) }{D_j(z_i)} \qquad (3) \\&\geq \sum_{j=1}^{10}\sum_{i=1}^{2}D_j(z_i)log\frac{P(y_j,z_i|\theta) }{D_j(z_i)} \qquad (4)\end{aligned} L(θ)=j=1∑10logi=1∑2P(yj,zi∣θ)=j=1∑10logi=1∑2Dj(zi)Dj(zi)P(yj,zi∣θ)(3)≥j=1∑10i=1∑2Dj(zi)logDj(zi)P(yj,zi∣θ)(4)
其中,对于每个实例i,用 D i D_i Di表示样本实例隐含变量z的某种分布,且 D i D_i Di满足条件 ∑ z D i ( z ) = 1 , D i ( z ) ≥ 0 \sum_zD_i(z)=1,D_i(z)\geq0 ∑zDi(z)=1,Di(z)≥0。由于对数函数是凹函数,则根据詹森不等式可以从(3)缩放为(4)。
具体来说,问题简化为如何求解随机变量的期望。假如我们的期望表达式为:
E ( x ) = ∑ x ∗ p ( x ) E(x)=\sum x*p(x) E(x)=∑x∗p(x)
我们可以把(3)式当中的 D j ( z i ) D_j(z_{i}) Dj(zi)看上式中的概率 p ( x ) p(x) p(x),把 P ( y j , z i ∣ θ ) D j ( z i ) \frac{P(y_j,z_i|\theta) }{D_j(z_i)} Dj(zi)P(yj,zi∣θ)看作是 D j ( z i ) D_j(z_{i}) Dj(zi)的取值x,可以得到:
∑ Z D j ( z i ) P ( y j , z i ∣ θ ) D j ( z i ) \sum\limits_Z D_j(z_i)\frac{P(y_j,z_i|\theta) }{D_j(z_i)} Z∑Dj(zi)Dj(zi)P(yj,zi∣θ)
就是 P ( y j , z i ∣ θ ) D j ( z i ) \frac{P(y_j,z_i|\theta) }{D_j(z_i)} Dj(zi)P(yj,zi∣θ)的期望。经过詹森不等式变形之后,就得到了式子(4),于是似然函数L(θ)可以简写成:
L ( θ ) ≥ J ( z , D ) (5) L(θ) \geq J(z,D) \tag 5 L(θ)≥J(z,D)(5)
的形式,那么:
我们可以不断优化L(θ)的下界 J ( z , D ) J(z,D) J(z,D)来达到优化L(θ)的目的。
现在有了两个问题:
1,式子(5)什么时候取等?
在詹森不等式中说到,当自变量X=E(X)时,即为常数的时候,等式成立。也就是:
P ( y j , z i ∣ θ ) D j ( z i ) = c (6) \frac{P(y_j,z_i|\theta) }{D_j(z_i)} = c \tag 6 Dj(zi)P(yj,zi∣θ)=c(6)
2,D函数怎么来的?
前面提到 ∑ z D i ( z ) = 1 \sum_zD_i(z)=1 ∑zDi(z)=1,则有以下的推导过程:
P ( y j , z i ∣ θ ) = D j ( z i ∣ θ ) c ∑ Z P ( y j , z i ∣ θ ) = ∑ Z D j ( z i ∣ θ ) c P(y_j,z_i|\theta) = D_j(z_i|\theta) c \\ \sum\limits_Z P(y_j,z_i|\theta) = \sum\limits_Z D_j(z_i|\theta) c P(yj,zi∣θ)=Dj(zi∣θ)cZ∑P(