统计学习方法 第九章笔记： 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,\dots ,y_n{)}^T$表示观测到硬币B、C的0，1序列，$Z＝(z_1,z_2,\dots ,z_n{)}^T$表示未观测到的硬币A的0，1序列。

考虑用极大似然估计求模型参数$\hat{\theta }=(\pi ,p,q)$的, 即：
$$\begin{array}{rl}L(\theta )& =\log P(Y|\theta )\\ & =\log \sum _{Z}P(Y,Z|\theta )\\ & =\sum _{j=1}^{10}log\sum _{i=1}^{2}P(y_j,z_i|\theta )\end{array}\text{\hspace{1em}(1)}$$

$$\hat{\theta }=\arg \underset{\theta }{\max }L(\theta )\text{\hspace{1em}(2)}$$

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

EM推导背景

对于式子1，我们做如下变形：
$$\begin{array}{rl}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)}(3)\\ & \ge \sum _{j=1}^{10}\sum _{i=1}^2{D}_j(z_i)log\frac{P(y_j,z_i|\theta )}{D_j(z_i)}(4)\end{array}$$
其中，对于每个实例i，用$D_i$表示样本实例隐含变量z的某种分布，且$D_i$满足条件${\sum }_z{D}_i(z)=1,D_i(z)\ge 0$。由于对数函数是凹函数，则根据詹森不等式可以从(3)缩放为(4)。

具体来说，问题简化为如何求解随机变量的期望。假如我们的期望表达式为：
$$E(x)=\sum x\ast 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 _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(\theta )\ge J(z,D)\text{\hspace{1em}(5)}$$
的形式，那么：

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

现在有了两个问题：

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

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

前面提到${\sum }_z{D}_i(z)=1$，则有以下的推导过程：
$$\begin{gathered} P(y_j,z_i|\theta )=D_j(z_i|\theta )c \\ \sum _{Z}P(y_j,z_i|\theta )=\sum _Z{D}_j(z_i|\theta )c \end{gathered}$$

$$\begin{gathered} \sum _{Z}P(y_j,z_i|\theta )=c \\ \text{\hspace{1em}(7)} \end{gathered}$$

将(7)代入(6)得到：
$$\begin{array}{rl}{D}_j(z_i)& =\frac{P(y_j,z_i|\theta )}{\sum _{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{array}$$
所以D函数本质上是隐变量z的后验概率：
$$D_j(z_i)=P(z_i|y_j,\theta )$$
到此为止，我们对式子(1)、(5)有了更多的了解，也就是有：
$$\begin{array}{rl}L(\theta )& \ge J(z,D)\\ & =\sum _{j=1}^{10}\sum _{i=1}^2{D}_j(z_i)log\frac{P(y_j,z_i|\theta )}{D_j(z_i)}\\ & =\sum _{j=1}^{10}\sum _{i=1}^{2}P(z_i|y_j,\theta )log\frac{P(y_j,z_i|\theta )}{P(z_i|y_j,\theta )}\end{array}$$
将上面的式子一般化可以得到：
$$\begin{array}{rl}L(\theta )& \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^{(i)})}\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}\end{array}\text{\hspace{1em}(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{array}{rl}Q(\theta ,{\theta }^{(i)})=& \sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y,Z|\theta )\end{array}$$
3，M步：求使$Q(\theta ,{\theta }^{(i)})$最大化的$\theta$，确定第$i+1$次迭代的参数估计值
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }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\ast 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{array}{rl}L(\theta )-L({\theta }^{(i)})& =\log \left(\sum _{Z}P(Y|Z,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(i)})}\right)-\log P(Y|{\theta }^{(i)})\\ & \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\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)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})\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)})P(Y|{\theta }^{(i)})}\end{array}$$
上面式子中的$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{array}{rl}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{array}$$
所以，任何可以使$B(\theta ,{\theta }^{(i)})$ 增大的$\theta$ 都可以使$L(\theta )$ 增大。所以为了使$L(\theta )$ 尽可能增大，我们选择${\theta }^{(i+1)}$ 使得$B(\theta ,{\theta }^i)$ 尽可能大：
$$\begin{array}{rl}{\theta }^{(i+1)}& =\arg \underset{\theta }{\max }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 \underset{\theta }{\max }\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 \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|Z,\theta )P(Z,\theta ))-\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Z|Y,{\theta }^i)P(Y|{\theta }^{(i)}))\\ & =arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y|Z,\theta )P(Z,\theta ))\\ & =arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})log(P(Y,Z|\theta ))\\ & =arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\end{array}$$
上面式子解释了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{array}{rl}logP(Y|\theta )& =\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|\theta )\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z,Y|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )\end{array}$$
其中，${\sum }_{Z}P(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 _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )$$
所以：
$$\begin{array}{rl}logP(Y|\theta )& =\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z,Y|\theta )-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )\\ & =Q(\theta ,{\theta }^{(i)})-H(\theta ,{\theta }^{(i)})\end{array}$$
进而有：
$$\begin{array}{rl}& 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{array}$$
因为${\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{array}{rl}& H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})\\ & =\sum _Z(log(\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}))P(Z|Y,{\theta }^{(i)})\\ & \le log(\sum _Z\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}P(Z|Y,{\theta }^{(i)}))\\ & =logP(Z|Y,{\theta }^{(i+1)})=log1=0\end{array}$$
不等号由詹森不等式得到。所以：
$$\begin{gathered} 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)}) \end{gathered}$$
证明了EM算法的收敛性。