统计学习方法 第九章笔记： EM 算法

文章目录

• 三硬币模型
• EM推导背景
• EM算法步骤
• EM算法的导出
• EM算法的收敛性

如果概率模型都是观测变量，那么给定数据就可以用极大似然估计法或者贝叶斯估计发去获得模型。但是，有时候概率模型既有观测变量，又有隐变量或者潜在变量，这样就不能用这些估计方法了。本文要介绍的EM算法可以解决这类问题。

三硬币模型

为了更好的描述EM算法，先提出一个实际的问题。问题描述如下：

假设有3枚硬币，分别记作A，B和C。 这些硬币正面向上的概率分别是 \pi、p和q 。进行如下抛硬币试验：先抛硬币A, 根据其结果选出硬币B或者硬币C，正面选硬币B，反面选硬币C；然后掷选出的硬币，抛硬币的结果，出现正面记作1，出现反面记作0；独立重复10次试验。结果如下：

1，1，0，1，0，0，1，0，1，1

假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率。

我们用\theta = (\pi,p,q)表示模型的参数，y表示某次得到的0或者1的观测结果，z表示隐变量，表示某次A的结果。一般地，用Y表示观测变量数据，用Z表示隐变量的数据。在上面的问题中，Y＝(y_1， y_2,…,y_n)^T表示观测到硬币B、C的0，1序列，Z＝(z_1,z_2,…,z_n)^T表示未观测到的硬币A的0，1序列。

考虑用极大似然估计求模型参数\hat \theta=(\pi,p,q)的, 即：
\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

\hat \theta = \arg\max\limits_{\theta} L(\theta) \tag 2

直接求(2)的参数θ会比较困难，因为上式存在隐含随机变量Z。如果Z是个已知的数，那么使用极大似然估计来估算会很容易。在这种Z不确定的情形下，EM算法就派上用场了。

EM推导背景

对于式子1，我们做如下变形：
\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}
其中，对于每个实例i，用D_i表示样本实例隐含变量z的某种分布，且D_i满足条件\sum_zD_i(z)=1,D_i(z)\geq0。由于对数函数是凹函数，则根据詹森不等式可以从(3)缩放为(4)。

具体来说，问题简化为如何求解随机变量的期望。假如我们的期望表达式为：
E(x)=\sum x*p(x)
我们可以把(3)式当中的D_j(z_{i})看上式中的概率p(x)，把\frac{P(y_j,z_i|\theta) }{D_j(z_i)}看作是D_j(z_{i})的取值x，可以得到：
\sum\limits_Z D_j(z_i)\frac{P(y_j,z_i|\theta) }{D_j(z_i)}
就是\frac{P(y_j,z_i|\theta) }{D_j(z_i)}的期望。经过詹森不等式变形之后，就得到了式子(4)，于是似然函数L(θ)可以简写成：
L(θ) \geq J(z,D) \tag 5
的形式，那么：

我们可以不断优化L(θ)的下界 J(z,D)来达到优化L(θ)的目的。

现在有了两个问题：

1，式子(5)什么时候取等？

在詹森不等式中说到，当自变量X=E(X)时，即为常数的时候，等式成立。也就是：
\frac{P(y_j,z_i|\theta) }{D_j(z_i)} = c \tag 6
2，D函数怎么来的？

前面提到\sum_zD_i(z)=1，则有以下的推导过程：
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

\sum\limits_Z P(y_j,z_i|\theta) = c \\ \tag 7

将(7)代入(6)得到：
\begin{aligned} D_j(z_i) &= \frac{P(y_j,z_i|\theta)}{\sum\limits_Z P(y_j,z_i|\theta)} \\ &= \frac{P(y_j,z_i|\theta)}{ P(y_j|\theta)} \\ &= P(z_i|y_j,\theta) \end{aligned}
所以D函数本质上是隐变量z的后验概率：
D_j(z_i) = P(z_i|y_j,\theta)
到此为止，我们对式子(1)、(5)有了更多的了解，也就是有：
\begin{aligned} L(θ)& \geq J(z,D) \\ &=\sum_{j=1}^{10}\sum_{i=1}^{2}\color{green} D_j(z_i) \color{black}log\frac{P(y_j,z_i|\theta) }{\color{green}D_j(z_i)} \\ &=\sum_{j=1}^{10}\sum_{i=1}^{2} \color{green}P(z_i|y_j,\theta) \color{black}log\frac{P(y_j,z_i|\theta) }{ \color{green}P(z_i|y_j,\theta) } \end{aligned}
将上面的式子一般化可以得到：
\begin{aligned} L(\theta) &\geq \sum_Z \color{green}P(Z|Y,\theta^{(i)})\color{black}\log \frac{P(Y,Z|\theta)}{\color{green}P(Z|Y,\theta^{(i)})\color{black}}\\ &= \sum_Z \color{green}P(Z|Y,\theta^{(i)})\color{black}\log \frac{P(Y|Z,\theta)P(Z|\theta)}{\color{green}P(Z|Y,\theta^{(i)})\color{black}} \end{aligned} \tag 8
了解了上面的基础知识，下面就可以了解实际的EM算法了。

EM算法步骤

输入: 观测变量数据Y，隐变量数据Z，联合分布P(Y,Z|\theta)，条件分布P(Z|Y,\theta)

输出: 模型参数\theta

1，选择参数的初值\theta^{(0)}，开始迭代

2，E步：记\theta^{(i)}为第 i 次迭代参数\theta的估计值，在第i+1次迭代的E步，计算
\begin{aligned} Q(\theta, \theta^{(i)}) =& \sum_Z\color{green}P(Z|Y, \theta^{(i)})\color{red}\log P(Y,Z|\theta) \end{aligned}
3，M步：求使Q(\theta, \theta^{(i)})最大化的\theta，确定第i+1次迭代的参数估计值
\theta^{(i+1)}=\arg\max_\theta Q(\theta, \theta^{(i)})
4，重复2，3两步，直到收敛。

上面的步骤中提到了Q函数 ，它是指完全数据(Y,Z)的对数似然函数 \log P(Y, Z|\theta)关于给定观测数据Y的当前参数\theta^{(i)}下，对观测数据 Z的条件概率分布 P(Z|Y,\theta^{(i)})的期望表达式。

也就是说对于期望函数E(x)=\sum x*p(x)而言，\log P(Y, Z|\theta) 是x，P(Z|Y,\theta^{(i)}) （也就是上一节提到的D函数）是p(x)。

要注意的是，Q(\theta, \theta^{(i)})中的θ表示是要极大化的参数。EM算法，或者说期望最大化算法的两个步骤实际在：

1. 给定\theta^{(i)}后，求期望的表达式，也就是Q函数的具体形式（先优化Z）
2. 求得到的Q函数(期望表达式)取到极大值的新参数\theta^{(i+1)} （再优化θ）

EM算法的导出

为什么EM算法能近似实现对观测数据的极大似然估计？Q函数又是怎么来的？

我们知道：
L(\theta)=\log P(Y|\theta)=\log \sum_Z P(Y,Z|\theta)=\log(\sum_Z P(Y|Z,\theta)P(Z|\theta))
那么有：
\begin{aligned} L(\theta)-L(\theta^{(i)})&=\log \left(\sum_Z\color{green}P(Y|Z,\theta^{(i)})\color{black}\frac{P(Y|Z,\theta)P(Z|\theta)}{\color{green}P(Y|Z,\theta^{(i)})}\color{black}\right)-\log P(Y|\theta^{(i)})\\ &\ge\sum_Z \color{green}P(Z|Y,\theta^{(i)})\color{black}\log \frac{P(Y|Z,\theta)P(Z|\theta)}{\color{green}P(Z|Y,\theta^{(i)})\color{black}}-\log P(Y|\theta^{(i)})\\ &=\sum_Z P(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}-\color{red}\sum_ZP(Z|Y,\theta^{(i)})\color{black}\log P(Y|\theta^{(i)})\\ &=\sum_ZP(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})} \end{aligned}
上面式子中的\color{green}P(Y|Z,\theta^{(i)})就是我们之前讨论的D函数，是为了构造期望， 进而应用詹森不等式，具体是怎么来的也讨论过了，这里不再赘述。

令：
B(\theta,\theta^{(i)}) = L(\theta^{(i)})+ \sum_{Z} P(Z|Y,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^{(i)}) P(Y|\theta^{(i)})}
则有：
L(\theta)\ge B(\theta,\theta^{(i)})
也就是B(\theta,\theta^{i}) 是L(\theta) 的一个下界，又由于：
P(Y,Z|\theta)=P(Y|Z,\theta)P(Z,\theta)=P(Z|Y,\theta) P(Y|\theta)
取\theta = \theta^{(i)} 时等号成立，因为：
\begin{aligned} B(\theta^{(i)},\theta^{(i)}) &= L(\theta^{(i)})+ \sum_{Z} P(Z|Y,\theta^{(i)})log\frac{P(Y|Z,\theta^{(i)})P(Z,\theta^{(i)})}{P(Z|Y,\theta^{(i)}) P(Y|\theta^{(i)})} \\ &= L(\theta^{(i)}) + \sum_{Z} P(Z|Y,\theta^{(i)})log1 \\ &=L(\theta^{(i)}) \end{aligned}
所以，任何可以使B(\theta,\theta^{(i)}) 增大的\theta 都可以使L(\theta) 增大。所以为了使L(\theta) 尽可能增大，我们选择\theta^{(i+1)} 使得B(\theta,\theta^{i}) 尽可能大：
\begin{aligned} \theta^{(i+1)} &= \arg \max_{\theta} B(\theta,\theta^{(i)})\\ &= L(\theta^{(i)})+ \sum_{Z} P(Z|Y,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^{(i)}) P(Y|\theta^{(i)})} \\ &=\arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^{(i)})log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}) \\ &=\arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^{(i)})log(P(Y|Z,\theta)P(Z,\theta))-\arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^{(i)})log(P(Z|Y,\theta^i)P(Y|\theta^{(i)}))\\ &=arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^{(i)})log(P(Y|Z,\theta)P(Z,\theta))\\ & =arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^{(i)})log(P(Y,Z|\theta))\\ &=arg \max_{\theta}Q(\theta,\theta^{(i)}) \end{aligned}
上面式子解释了Q函数的由来，并且完成了EM的一次迭代过程，得到了\theta^{(i+1)}。上面的推导过程也说明了：EM算法是通过不断求解使得下界B(\theta,\theta^{i}) 极大化的参数θ，来近似求解对数似然函数L(θ)极大值的算法。

下图给出EM算法的直观解释。图中上方曲线为L(\theta),下方曲线为B(\theta,\theta^{(i)})。B(\theta,\theta^{(i)})为对数似然函数L(\theta)的下界。

通过图片我们总结ＥＭ算法的迭代过程如下：

(1) 初始状态：在\theta = \theta ^{(i)} 时，L(\theta^{(i)})=B(\theta^{(i)},\theta^{(i)})

(2) EM算法的更新过程：寻找\theta^{(i+1))} 使得B(\theta,\theta^{(i)})极大化。因为L(\theta)\ge B(\theta,\theta^{(i)})，所以得B(\theta,\theta^{(i)})极大化也使得L(\theta) 增加。

(3) 根据第二步找到的\theta^{(i+1))} 重新计算B(\theta,\theta^{(i)})（或者说重新计算Q函数），进行下一次迭代。

上述过程暴露一个问题：对数似然函数不断增大，但EM算法不能保证找到全局最优解。

EM算法的收敛性

我们已经证明了EM算法的合理性，接下来我们来研究一下算法是否收敛。

要证EM算法的收敛性即证:
P(Y|\theta^{(i+1)}) \ge P(Y|\theta^{(i)})
同样由：
P(Y,Z|\theta)=P(Y|Z,\theta)P(Z,\theta)=P(Z|Y,\theta) P(Y|\theta)
可以得到：
P(Y|\theta) = \frac{P(Z,Y|\theta)}{P(Z|Y,\theta)}
两边同取对数：
logP(Y|\theta) = logP(Z,Y|\theta) - logP(Z|Y,\theta)
左右两边同时对 P(Z|Y,\theta^{(i)})求期望，则有：
\begin{aligned} logP(Y|\theta) &= \sum_ZP(Z|Y,\theta^{(i)})logP(Y|\theta) \\ &= \sum_ZP(Z|Y,\theta^{(i)})logP(Z,Y|\theta) - \sum_ZP(Z|Y,\theta^{(i)})logP(Z|Y,\theta) \\ \end{aligned}
其中，\sum_ZP(Z|Y,\theta^{(i)})=1,logP(Y|\theta) 是与Z无关的变量。

根据Q函数的定义：
Q(\theta,\theta^{(i)}) = \sum_Z logP(Y,Z|\theta)P(Z|Y,\theta^{(i)})
以及这里我们令:
H(\theta,\theta^{(i)}) = \sum_ZP(Z|Y,\theta^{(i)})logP(Z|Y,\theta)
所以：
\begin{aligned} logP(Y|\theta) &= \sum_ZP(Z|Y,\theta^{(i)})logP(Z,Y|\theta) - \sum_ZP(Z|Y,\theta^{(i)})logP(Z|Y,\theta) \\ &= Q(\theta,\theta^{(i)}) - H(\theta,\theta^{(i)}) \end{aligned}
进而有：
\begin{aligned} & \quad logP(Y,\theta^{(i+1)})-logP(Y,\theta^{(i)}) \\ &= Q(\theta^{(i+1)},\theta^{(i)}-H(\theta^{(i+1)},\theta^{(i)}) - (Q(\theta^{(i)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)}))\\ &= Q(\theta^{(i+1)},\theta^{(i)})- Q(\theta^{(i)},\theta^{(i)}) -( H(\theta^{(i+1)},\theta^{(i)}) - H(\theta^{(i)},\theta^{(i)})) \end{aligned}
因为\theta^{i+1} 使得 Q(\theta,\theta^{(i)}) 增加，所以
Q(\theta^{(i+1)},\theta^{(i)})- Q(\theta^{(i)},\theta^{(i)}) \ge 0
所以，我们的目标是证明H(\theta^{(i+1)},\theta^{(i)}) - H(\theta^{(i)},\theta^{(i)}) \le 0 ，证明如下：
\begin{aligned} & \quad H(\theta^{(i+1)},\theta^{(i)}) - H(\theta^{(i)},\theta^{(i)}) \\ &=\sum_{Z}\bigg(log(\frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}) \bigg)P(Z|Y,\theta^{(i)}) \\ & \le log\bigg(\sum_{Z}\frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}P(Z|Y,\theta^{(i)}) \bigg)\\ &=logP(Z|Y,\theta^{(i+1)})=log1=0 \end{aligned}
不等号由詹森不等式得到。所以：
logP(Y,\theta^{(i+1)})-logP(Y,\theta^{(i)}) \ge 0 \\ P(Y,\theta^{(i+1)})-P(Y,\theta^{(i)}) \ge 0 \\ P(Y,\theta^{(i+1)})\ge P(Y,\theta^{(i)})
证明了EM算法的收敛性。