本文首发于CSDN，https://blog.csdn.net/weixin_43661031/article/details/91358990， 请勿转载。

最近在做文本挖掘的时候遇到了EM算法，虽然读书的时候简单地接触过，但当时并没有深入地去了解，导致现在只记得算法的名字。既然EM算法被列为数据挖掘的十大算法之一，正好借这个机会，重新学习一下这个经典的算法。学习的过程中，我发现网上的资料大多讲解地不够细致，很多地方解释得并不明了。因此我决定抛开别人的想法，仅从数学推导本身出发，尽力理解每一个公式的含义，并将其对应到实际的实验过程当中。这篇博客记录了我对与EM算法的思考与理解，也是我人生中的第一篇博客，希望能够对于想要学习EM算法的同学有所帮助。

极大似然估计与隐变量

前面谈到我在做文本挖掘的时候遇到了EM算法，EM算法用于估计模型中的参数。提到参数估计，最常见的方法莫过于极大似然估计——在所有的候选参数中，我们选择的参数应该让样本出现的概率最大。相信看到这篇笔记的同学一定对极大似然估计非常熟悉，而EM算法可以看作是极大似然估计的一个扩充，这里就让我们用极大似然估计来解决一个简单的例子，来开始正式的讨论。

抛硬币1.

有A，B，C三枚硬币，我们想要估计A，B，C三枚硬币抛出正面的概率\(\pi_A\), \(\theta_B\), \(\theta_C\)。我们按如下流程进行实验100次：

     随机抛一次硬币A
     若硬币A抛出正面，则抛10次硬币B
     若硬币A抛出反面，则抛10次硬币C
     

记录100次实验的结果如下：
| |硬币A|硬币B/C抛出正面的次数|
|--|----|------|
|1|1 | 5|
|2|0 | 8|
|3|1 | 4|
|...|... | ... |
|100|0 | 7|
我们将上面的实验结果表述如下：
\(z^{(i)}\)表示第i次实验中，硬币A的结果，1代表正面，0代表反面；\(x^{(i)}\)表示第i次实验中，硬币B或硬币C抛出正面的个数，则参数\(\pi_A, \theta_B, \theta_C\)的极大似然估计分别为：
\[\begin{aligned} \hat\pi_A &= \frac{\sum_{i=1}^{100}\mathbb I_{z^{(i)}=1}(z^{(i)})}{100}\\[2ex] \hat\theta_B &= \frac{\sum_{i=1}^{100}\mathbb I_{z^{(i)}=1}(z^{(i)})x^{(i)}}{10\sum_{i=1}^{100}\mathbb I_{z^{(i)}=1}(z^{(i)})}\\[3ex] \hat\theta_C&=\frac{\sum_{i=1}^{100}\mathbb I_{z^{(i)}=0}(z^{(i)})x^{(i)}}{10\sum_{i=1}^{100}\mathbb I_{z^{(i)}=0}(z^{(i)})}\\ \end{aligned}\]
即硬币A，B，C各自抛出正面的次数占总次数的比例，其中\(\mathbb I(\cdot)\)为指示函数。

抛硬币2.

实验流程与1相同，但是我们不慎遗失了硬币A的记录结果，导致我们只知道随后十次抛出了多少次正面，多少次反面，却不知道实验结果来自于硬币B还是硬币C。在这种情况下，我们是否还能估计出\(\pi_A\), \(\theta_B\), \(\theta_C\)的值呢？
| |硬币A|硬币B/C抛出正面的次数|
|--|--------|-----------|
|1|？ | 5|
|2|？ | 8|
|3|？ | 4|
|...|... | ... |
|100|？ | 7|

这时候利用极大似然估计似乎行不通了， 因为这种情况下，我们不但缺失了硬币A产生的观测值，同时也不知道哪些观测值属于硬币B，哪些观测值属于硬币C。

有些同学可能会提出，虽然我们无法得到三个硬币各自产生的样本，但是我们依然可以得到每个观测值出现的概率。比如在第一次实验中， 我们抛出了5次正面5次反面，我们可以做如下思考：

  假设这5次正面由硬币B得到，那么概率应该为 \(C_{10}^5\theta_B^5(1 - \theta_B)^5\)，而这次观测值来自于硬币B，也就是硬币A抛出正面的概率为\(\pi_A\)

  假设这5次正面由硬币C得到，那么概率应该为 \(C_{10}^5\theta_C^5(1 - \theta_C)^5\)，而这次观测值来自于硬币C，也就是硬币A抛出反面的概率为\(1-\pi_A\)

  综合起来，利用条件概率公式，这个观测值出现的概率就是\(\pi_AC_{10}^5\theta_B^5(1 - \theta_B)^5+ (1- \pi_A) C_{10}^5\theta_C^5(1 - \theta_C)^5\)

因此我们可以将样本整体的概率和似然函数利用\(\pi_A\), \(\theta_B\), \(\theta_C\)表示出来，通过对似然函数求导，令其关于\(\pi_A, \theta_B, \theta_C\)的偏导数等于0，我们可以求出三个参数的值。

这个思路听上去十分合理，我们可以顺着这个思路进行数学推导，看看可以得到什么样的结果。首先我们计算样本的概率:

\[\begin{aligned} &P(x^{(1)}, x^{(2)}, ..., x^{(100)};\pi_A, \theta_B,\theta_C) \\ &=\prod_{i=1}^{100}P(x^{(i)};\pi_A, \theta_B,\theta_C) \\ &=\prod_{i=1}^{100}[P(x^{(i)}, z^{(i)}=1;\pi_A, \theta_B,\theta_C) + P(x^{(i)}, z^{(i)}=0;\pi_A, \theta_B,\theta_C)]\\ &=\prod_{i=1}^{100}[P(x^{(i)}|z^{(i)}=1;\pi_A, \theta_B,\theta_C)P(z^{(i)}=1;\pi_A, \theta_B,\theta_C) + P(x^{(i)}|z^{(i)}=0;\pi_A, \theta_B,\theta_C)P(z^{(i)}=0;\pi_A, \theta_B,\theta_C)]\\ & =\prod_{i=1}^{100}[P(x^{(i)}|z^{(i)}=1;\theta_B)P(z^{(i)}=1;\pi_A) +P(x^{(i)}|z^{(i)}=0;\theta_C)P(z^{(i)}=0;\pi_A)] \end{aligned}\]

对应的似然函数为
\[\begin{aligned} &L(\pi_A, \theta_B,\theta_C|x^{(1)}, x^{(2)}, ..., x^{(100)})\\ &=\sum_{i=1}^{100}log[P(x^{(i)}|z^{(i)}=1;\theta_B)P(z^{(i)}=1;\pi_A) +P(x^{(i)}|z^{(i)}=0;\theta_C)P(z^{(i)}=0;\pi_A)] \end{aligned}\]

其中\(x^{(i)}\)关于\(z^{(i)}\)的条件分布为
\[\left\{ \begin{aligned} &P(x^{(i)}|z^{(i)}=1;\theta_B) = C_{10}^{x^{(i)}}\theta_B^{x^{(i)}}(1 - \theta_B)^{10-x^{(i)}}\\ &P(x^{(i)}|z^{(i)}=0;\theta_C) = C_{10}^{x^{(i)}}\theta_C^{x^{(i)}}(1 - \theta_C)^{10-x^{(i)}}\\ \end{aligned} \right.\]

\(z^{(i)}\)的分布为
\[\left\{ \begin{aligned} &P(z^{(i)}=1;\pi_A) = \pi_A \\ &P(z^{(i)}=0;\pi_A) = 1- \pi_A \end{aligned} \right.\]

因此我们可以得到
\[\begin{aligned} &L(\pi_A, \theta_B,\theta_C|x^{(1)}, x^{(2)}, ..., x^{(100)})\\[2ex] &= \sum_{i=1}^{100}log[\pi_AC_{10}^{x^{(i)}}\theta_B^{x^{(i)}}(1 - \theta_B)^{10-x^{(i)}}+ (1- \pi_A) C_{10}^{x^{(i)}}\theta_C^{x^{(i)}}(1 - \theta_C)^{10-x^{(i)}}] \end{aligned}\]

至此，我们成功地得到了似然函数。然而观察可以发现，这个函数是由100项对数函数相加组成，每个对数函数内部包含一个求和，想通过求导并解出导数的零点几乎是不可能的。当然我们可以通过梯度下降来极小化这个函数，借助深度学习库的自动微分系统在实现上也非常容易。但是这种做法过于简单粗暴，有没有办法来优雅地解决这个问题呢？在继续讨论之前，我们先将这类问题进行一般化表述：

我们观测到随机变量\(X\)产生的m个相互独立的样本\(\left\{x^{(1)}, x^{(2)}, ..., x^{(m)}\right\}\), \(X\)的分布由联合分布\(P(X, Z;\theta)\)决定，\(Z\)是缺失数据或无法在实验中被直接观测到，称为隐变量，我们想要从样本中估计出模型参数\(\theta\)的值。在接下来的讨论中，我们假定\(Z\)的取值是离散的，于是可以得到似然函数如下:
\[\begin{aligned} L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)}) &= \sum_{i=1}^mlogP(x^{(i)};\theta)\\ &= \sum_{i=1}^mlog(\sum_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta)) \end{aligned}\tag 1\]

接下来，我们就探讨一下，如何利用EM算法解决这个问题。

EM算法的数学原理

这一部分的数学推导，主要参考了吴恩达CS229n的笔记，并且根据个人的思考和理解，尽力对公式的每一步进行详细的解释。我们先简单地介绍一下琴生不等式。

琴生不等式

琴生不等式有多种形式，下面给出其离散形式的表述和概率论中的表述:
1.若\(\phi(x)\)为严格凹函数，\(x_1, x_2, ...,x_n\)为定义域内的n个点，\(a_1, a_2, ... a_n\)是n个正实数，且满足\(\sum_{i=1}^na_i =1\), 则下述不等式成立:
\[\phi(\sum_{i=1}^na_ix_i) >= \sum_{i=1}^na_i\phi(x_i)\]

当且仅当\(x_i=x_2=...=x_n\)时，不等式取等号。

2.若\(\phi(x)\)为严格凹函数，\(X\)为实值随机变量，且期望存在，则下述不等式成立:
\[\phi(E(X)) >= E(\phi(X))\]

当且仅当\(X \equiv E(X)\)，即\(X\)为常数时，不等式取等号。

注： 这里将函数上方为凹集的函数称为凹函数， 例如\(log\)函数就是凹函数。
相信大家对琴生不等式都十分熟悉，因此这里就不做过多的说明。接下来，我们将琴生不等式应用到我们的问题中。

EM算法的数学基础

回到我们之前的问题上， 我们想要极大化下面这个函数:
\[\begin{aligned} L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)}) &= \sum_{i=1}^mlog(P(x^{(i)};\theta)) \\ &= \sum_{i=1}^mlog(\sum_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta))\tag 1 \end{aligned}\]

但是我们无法对这个函数直接求导，因此我们借助琴生不等式，对这个函数进行变换。为了让过程看上去简洁，下面只对求和中的第\(i\)项进行计算。

令\(Q^i(x)\)满足 \(\sum_{z^{(i)}}Q^i(z^{(i)}) = 1\)，且\(Q^i(z^{(i)})>=0, \forall z^{(i)}\)，则根据琴生不等式，可以得到：
\[\begin{aligned} log(\sum_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta)) &= log(\sum_{z^{(i)}}Q^i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})})\\ &\ge \sum_{z^{(i)}}Q^i(z^{(i)})log(\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})})\tag 2 \end{aligned}\]

当且仅当\(\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\)为常数时，上述不等式取等号。也就是说，对于任意\(z^{(i)}\)，\(\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\)是一个与\(z^{(i)}\)无关的量。设对于任意\(z^{(i)}, \frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})} \equiv k\)，我们可以得到：
\[\begin{aligned} &\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})} =k \\[1.9ex] \Rightarrow \quad&P(x^{(i)},z^{(i)};\theta)=kQ^i(z^{(i)}) \\[2ex] \Rightarrow \quad&\sum_{z^{(i)}}P(x^{(i)},z^{(i)};\theta)=k\sum_{z^{(i)}}Q^i(z^{(i)})\\[1ex] \Rightarrow \quad&P(x^{(i)};\theta)=k \qquad\\ \Rightarrow \quad&Q^i(z^{(i)})= \frac{P(x^{(i)},z^{(i)};\theta)}{P(x^{(i)};\theta)}=P(z^{(i)}|x^{(i)};\theta)\tag 3 \end{aligned}\]

因此当\(Q^i(z^{(i)}) =P(z^{(i)}|x^{(i)};\theta)\)时，不等式\((2)\)取等号，容易验证此时\(\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})} =P(x^{(i)};\theta)\), 与\(z^{(i)}\)无关。将\((1), (2), (3)\)综合一下，我们可以得到以下结论:
\[\begin{aligned} L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)}) &= \sum_{i=1}^mlog(\sum_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta)) \\ &= \sum_{i=1}^mlog(\sum_{z^{(i)}}\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta)}P(z^{(i)}|x^{(i)};\theta)) \\ &= \sum_{i=1}^m\sum_{z^{(i)}}P(z^{(i)}|x^{(i)};\theta)log\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta)} \\ &\ge\sum_{i=1}^m\sum_{z^{(i)}}Q^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\tag 4 \end{aligned}\]

到这里为止，我们已经拥有了推导出EM算法的全部数学基础，基于\((4)\)我们可以构建出E步和M步。上面的数学推导虽然看上去略为复杂，但实际上只用到了三个知识点：
  1.琴生不等式:\[log\sum_{z^{(i)}}Q^i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\ge\sum_{z^{(i)}}Q^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\]

  2.条件概率: \[P(z^{(i)}|x^{(i)};\theta)=\frac{P(x^{(i)},z^{(i)};\theta)}{P(x^{(i)};\theta)}\]

  3.联合分布求和等于边缘分布: \[\sum_{z^{(i)}}P(x^{(i)},z^{(i)};\theta) = P(x^{(i)};\theta)\]

对上面的数学推导有疑问的同学，可以结合上面这三点，再将整个推导过程耐心地看一遍。

\(Q\)函数

大部分关于EM算法的资料，只是在数学形式上引入了\(Q\)函数，即\(Q^i(z^{(i)})\)，以满足琴生不等式的使用条件，却没有过多地解释\(Q\)函数本身。这导致了很多人完全看懂了算法的推导，却还是不理解这些数学公式究竟在做什么，甚至不明白EM算法为什么叫做EM算法。所以在给出E步和M步之前，我想先谈一谈\(Q\)函数。

我们回顾一下\(Q\)函数所满足的条件（暂时不考虑琴生不等式取等号的限制），
\[\begin{cases} \sum_{z^{(i)}}Q^i(z^{(i)})=1 \\[2ex] Q^i(z^{(i)}) \ge 0 \end{cases}\]

\(Q^i(x)\)在\(z^{(i)}\)所有可能的取值处有定义。可以看出，\(Q^i(x)\)是\(z^{(i)}\)的样本空间上任意的一个概率分布。因此，我们可以对不等式\((4)\)进行改写。首先我们可以将含有\(Q\)的求和写成期望的形式:
\[\begin{aligned} log\sum_{z^{(i)}}Q^i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})} &= log\left(\mathbb E_{Q^i}\left[ \frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\right]\right)\\[2ex] \sum_{z^{(i)}}Q^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}&=\mathbb E_{Q^i} \left[log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\right]\tag5 \end{aligned}\]

这里\(\mathbb E_{Q^i}\)指的是在概率分布\(Q^i\)下，求随机变量\(\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\)和\(log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\)的期望。有同学会问，为什么我们平时求期望的时候只要写\(\mathbb E\)，并没有指明是在哪个概率分布下的期望。这是因为一般情况下，我们都清楚地知道随机变量\(X\)所服从的分布\(P\)，并且默认在分布\(P\)下求期望。

举个例子，我手上有一个硬币，抛了10次，问抛出正面次数的期望。这种情况下，大部分人会默认硬币是均匀的，也就是说抛出正面的次数\(X\)服从二项分布\(B(10, 0.5)\)，期望\(\mathbb EX=\mathbb E_{B(10, 0.5)}X=5\)。这时有人提出了质疑，他说我认为你这个硬币有问题，抛出正面的概率只有0.3，那么在他眼里， 期望\(\mathbb EX=\mathbb E_{B(10, 0.3)}X=3\)。

回到正题，我们利用等式\((5)\)改写不等式\((4)\)，可以得到:
\[\begin{aligned} L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)}) &= \sum_{i=1}^mlog\sum_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta)\\[2ex] &= \sum_{i=1}^m\mathbb E_{P(z^{(i)}|x^{(i)};\theta)}\left[ log\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta)} \right]\\[2ex] &\ge\sum_{i=1}^m \mathbb E_{Q^i(z^{(i)})} \left[log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\right]\tag6 \end{aligned}\]

这正是琴生不等式在概率论中的形式。我们可以将不等式倒过来理解：
  首先，假定随机变量\(z^{(i)}\)服从概率分布\(Q^i\)，\(Q^i\)是\(z^{(i)}\)的样本空间上的任意一个概率分布。这里\(Q^i\)可以是一组定值，也可以是关于参数\(\theta\)的函数。

  显然，当我们取不同的\(Q^i\)时，随机变量\(log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\)的期望也会随之改变。需要注意的是，由于\(P(x^{(i)},z^{(i)};\theta)\)与\(\theta\)相关，所以这里的期望不是一个数值，而是关于\(\theta\)的函数。

  当我们令\(Q^i\)为\(z^{(i)}\)的后验分布\(P(z^{(i)}|x^{(i)};\theta)\)时，上面的期望最大。这里有两点需要注意，1. 后验分布\(P(z^{(i)}|x^{(i)};\theta)\)也是一个关于参数\(\theta\)的函数。2. 由于期望是关于\(\theta\)的函数，所以这里的最大指的并非是最大值，而是最大的函数。

  若对于每一个\(i\)，我们都令\(Q^i\)为\(z^{(i)}\)的后验分布\(P(z^{(i)}|x^{(i)};\theta)\)，则上述期望之和等于我们要极大化的似然函数，即\[L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)}) = \sum_{i=1}^m\mathbb E_{P(z^{(i)}|x^{(i)};\theta)}\left[ log\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta)} \right]\]

通过上述分析，我们为寻找似然函数的极大值点\(\hat\theta_{MLE}\)提供了一个思路。我们不去极大化似然函数本身，而是去极大化\(\sum_{i=1}^m \mathbb E_{Q^i(z^{(i)})} \left[log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\right]\)。至于如何将这个思路实际应用，就要利用到EM算法中的E-step和M-step。

E-step和M-step

这一节中，我们先给出E-step和M-step的数学形式，随后在结合抛硬币的例子来解释这两步究竟在做什么。下面进入算法的流程，首先我们任意初始化\(\theta_0\)，按下述过程进行迭代直至收敛：

在第\(t\)次迭代中，
(E-step)对于每个\(i\)，令\[Q_t^i(z^{(i)})=P(z^{(i)}|x^{(i)};\theta_{t-1})\tag 7\]
(M-step)更新\(\theta\)的估计值，令\[\theta_t=\mathop{argmax} \limits_\theta\sum_{i=0}^m\sum_{z^{(i)}}Q_t^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q_t^i(z^{(i)})}\tag 8\]

EM算法从任意一点\(\theta_0\)出发，依次利用E-step优化\(Q^i\)，M-step优化\(\theta\)，重复上述过程从而逐渐逼近极大值点。而这个过程究竟是怎样的呢，就让我们一步步地揭开EM算法的面纱。

假设我们现在随机初始化了\(\theta_0\)，进入第一轮迭代：
(E-step)

\[\begin{aligned}Q_1^i(z^{(i)}) &=P(z^{(i)}|x^{(i)};\theta_0) \\[2ex] &= \frac{P(x^{(i)}|z^{(i)};\theta_0)P(z^{(i)};\theta_0)}{\sum_{z^{(i)}}P(x^{(i)}|z^{(i)};\theta_0)P(z^{(i)};\theta_0)} \end{aligned}\]

由于我们已经假定模型参数为\(\theta_0\)，所以此时\(Q_1^i(z^{(i)})\)不再是与\(\theta\)有关的函数，而是由一组常数构成的概率分布。结合抛硬币的例子来看，这一步是在我们已知模型参数\(\pi_{A0},\theta_{B0}, \theta_{C0}\)的基础上(虽然这是我们瞎猜的)，去推测每一次的观测值是由哪个硬币产生的，或者说我们对每一次观测值做一个软分类。比如我们根据初始化的参数，计算出\(Q_1^i(z^{(i)}=1)=0.2\)，\(Q_1^i(z^{(i)}=0)=0.8\)。可以解释为第\(i\)个观测值有20%的概率来自于硬币B，80%的概率来自于硬币C；或者说硬币A抛出了0.2个正面，0.8个反面。

(M-step)\[\theta_1=\mathop{argmax} \limits_\theta\sum_{i=0}^m\sum_{z^{(i)}}Q_1^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q_1^i(z^{(i)})}\]

考虑到\(Q_1^i(z^{(i)})\)是一组常数，我们可以舍弃常数项，进一步简化上面这个要极大化的函数
\[\begin{aligned} \theta_1 &=\mathop{argmax} \limits_\theta\sum_{i=0}^m\sum_{z^{(i)}}Q_1^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q_1^i(z^{(i)})}\\ &=\mathop{argmax} \limits_\theta\sum_{i=0}^m\sum_{z^{(i)}}Q_1^i(z^{(i)})logP(x^{(i)},z^{(i)};\theta) \end{aligned}\]

由于\(Q_1^i(z^{(i)})\)不再与\(\theta\)相关，因此上面的函数变成了对数函数求和的形式，这个函数通常来说是容易求导的，令导数等于0，我们可以求出新的参数\(\theta_1\)。我们仍旧以抛硬币为例进行解释，
\[\begin{aligned} 令&f(\pi_A, \theta_B,\theta_C) \\ &= \sum_{i=1}^m\sum_{z^{(i)}}Q_1^i(z^{(i)})logP(x^{(i)},z^{(i)};\pi_A, \theta_B, \theta_C)\\[2ex] &= \sum_{i=1}^{100}[Q_1^i(z^{(i)}=1)logP(x^{(i)},z^{(i)}=1;\pi_A, \theta_B) + Q_1^i(z^{(i)}=0)logP(x^{(i)},z^{(i)}=0;\pi_A, \theta_C)]\\[2ex] &=\sum_{i=1}^{100}[Q_1^i(z^{(i)}=1)log(P(x^{(i)}|z^{(i)}=1;\theta_B)P(z^{(i)}=1;\pi_A)) + Q_1^i(z^{(i)}=0)log(P(x^{(i)}|z^{(i)}=0;\theta_C)P(z^{(i)}=0;\pi_A))]\\[2ex] &=\sum_{i=1}^{100}[Q_1^i(z^{(i)}=1)log(\theta_B^{x^{(i)}}(1-\theta_B)^{10-x^{(i)}}\pi_A) + Q_1^i(z^{(i)}=0)log(\theta_C^{x^{(i)}}(1-\theta_C)^{10-x^{(i)}}(1-\pi_A))]\\[2ex] &=\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)(x^{(i)}log\theta_B+(10-x^{(i)})log(1-\theta_B))+\sum_{i=1}^{100}Q_1^i(z^{(i)}=0)(x^{(i)}log\theta_C+(10-x^{(i)})log(1-\theta_C))\\[2ex] &+\sum_{i=1}^{100}[Q_1^i(z^{(i)}=1)log\pi_A+Q_1^i(z^{(i)}=0)log(1-\pi_A)] \end{aligned}\]
令\(\frac{\partial f}{\partial \pi_A}=0, \frac{\partial f}{\partial \theta_B}=0, \frac{\partial f}{\partial \theta_C}=0\), 可以得到，
\[\begin{aligned} &\pi_{A1}=\frac{\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)}{100}\\[2ex] &\theta_{B1}=\frac{\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)x^{(i)}}{10\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)}\\[3ex] &\theta_{C1}=\frac{\sum_{i=1}^{100}Q_1^i(z^{(i)}=0)x^{(i)}}{10\sum_{i=1}^{100}Q_1^i(z^{(i)}=0)} \end{aligned}\]

这三个参数的解释是显而易见的。我们在E-step中对每个观测值进行了软分类，\(\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)\)可以看成是硬币A抛出正面的次数，所以\(\frac{\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)}{100}\)是\(\pi_A\)的极大似然估计；\(10\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)\)是我们抛硬币B的次数，\(\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)x^{(i)}\)是硬币B抛出正面的次数，所以\(\frac{\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)x^{(i)}}{10\sum_{i=1}^{100}Q_1^i(z^{(i)}=1)}\)是\(\theta_B\)的极大似然估计；对于\(\theta_C\)我们有相同的解释。

我们将这个结果与抛硬币1中极大似然估计的结果相比较可以发现，之前结果中的指示函数\(\mathbb I(\cdot)\)变成了这里的\(Q(\cdot)\)，在指示函数下，某个观测值要么来自于硬币B，要么来自于硬币C，因此也称为硬分类。而在\(Q\)函数下，某个观测值可以一部分来自于硬币B，一部分来自于硬币C，因此也称作软分类。

将上述两步综合起来，EM算法可以总结如下：我们首先初始化模型的参数，我们基于这个参数对每一个隐变量进行分类，此时相当于我们观测到了隐变量。有了隐变量的观测值之后，原来含有隐变量的模型变成了不含隐变量的模型，因此我们可以直接使用极大似然估计来更新模型的参数，再基于新的参数开始新一轮的迭代，直到参数收敛。接来下我们就讨论为什么参数一定会收敛。

EM算法的收敛性

前面写了太多的公式，但是这一部分我不打算给出收敛性的数学推导。其实数学上证明EM算法的收敛性很容易，只需要证明每一轮迭代之后，参数的似然函数递增，即\(L(\theta_{t+1}|x^{(1)}, x^{(2)}, ..., x^{(m)}) \ge L(\theta_{t}|x^{(1)}, x^{(2)}, ..., x^{(m)})\)，再结合单调有界收敛定理便可证明。

下面我从最直观的角度，展示一下EM算法迭代的全过程。
首先我们看一下\(L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)})\)与\(\begin{aligned}\sum_{i=1}^m\sum_{z^{(i)}}Q_t^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q_t^i(z^{(i)})}\end{aligned}\)之间的关系。我们之前说过，当\(Q^i(z^{(i)})\)自由变化时，\(\begin{aligned}\sum_{i=1}^m\sum_{z^{(i)}}Q^i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q^i(z^{(i)})}\end{aligned}\)的集合构成了一族函数，而当我们固定\(\theta\)为任意\(\theta^\star\)，并令\(Q^i(z^{(i)})=P(z^{(i)}|x^{(i)};\theta^\star)\)时，\(\begin{aligned}\left\{\sum_{i=1}^m\sum_{z^{(i)}}P(z^{(i)}|x^{(i)};\theta^\star)log\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta^\star)}, \forall \theta^\star \in \Theta\right\}\end{aligned}\)便构成了函数族中一个特殊的子集，注意此时\(P(z^{(i)}|x^{(i)};\theta^\star)\)被\(\theta^\star\)固定，不再与\(\theta\)有关。这一族子集具有如下性质：
\[\begin{aligned} L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)}) &\ge \sum_{i=1}^m\sum_{z^{(i)}}P(z^{(i)}|x^{(i)};\theta^\star)log\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta^\star)} \\ L(\theta^\star|x^{(1)}, x^{(2)}, ..., x^{(m)}) &= \sum_{i=1}^m\sum_{z^{(i)}}P(z^{(i)}|x^{(i)};\theta^\star)log\frac{P(x^{(i)},z^{(i)};\theta^\star)}{P(z^{(i)}|x^{(i)};\theta^\star)} \end{aligned}\]

也就是说，每一个\(\begin{aligned}\sum_{i=1}^m\sum_{z^{(i)}}P(z^{(i)}|x^{(i)};\theta^\star)log\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta^\star)}\end{aligned}\)都是\(L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)})\)的一个下界，且在唯一的一点\(\theta=\theta^\star\)处相切，\(L(\theta|x^{(1)}, x^{(2)}, ..., x^{(m)})\)构成了这一族曲线的包络线。下面的动图直观地展示了这一点：
EM算法正是在这一族曲线上进行反复迭代，逐渐逼近似然函数的极大值点。因此EM算法也可以看成是如下的流程：
  1.选定一个初始值\(\theta_0\)，
  2.选定在\(\theta_0\)与似然函数相切的函数\(\begin{aligned}\sum_{i=1}^m\sum_{z^{(i)}}P(z^{(i)}|x^{(i)};\theta^0)log\frac{P(x^{(i)},z^{(i)};\theta)}{P(z^{(i)}|x^{(i)};\theta^0)}\end{aligned}\)
  3.求出上述函数的极大值\(\theta_1\)，将\(\theta_0\)更新为\(\theta_1\)，重复上述过程直至收敛

从图中可以看出，E-step的是从一个函数跳跃到另一个函数，在\(x\)不变的情况下，直接对\(y\)进行优化。而M-step则是沿着某一个函数曲线进行优化。到此为止，我们已经从不同的角度深度剖析了EM算法的内涵。最后我补充两点来结束本文的讨论。

小结

本文从最基本的极大似然估计开始，讨论了模型中存在隐变量时，极大似然估计所遇到的困难，从而引入了EM算法。我们将琴生不等式应用于似然函数中，并详细解释了\(Q\)函数的意义，进而给出了E-step和M-step，和这两步在抛硬币例子中的实际应用。最后我们从函数图像的角度，直观地展示了EM算法收敛的过程。相信耐心读完本文的同学， 对EM算法一定有了一个比较透彻的理解。关于EM算法，我在这里在最后补充两点：
  1.EM算法并不能给我们提供任何额外的信息，它只是给我们提供了一条巧妙的路径去逼近似然函数的极大值。即使我们不用EM算法，而是借助计算机直接对似然函数求导，同样可以获得问题的解。对于抛硬币问题而言，EM算法收敛之后，依旧不能告诉我们某个观测值到底是来自B还是C，而是只能给出一个概率。
  2.EM算法无法保证求出的解是全局最优解，它只能给我们一个局部最优解，至于求出的是哪个局部最优解，就与函数的形状和给定的初值有关了。在样本规模不太大的情况下，可以对参数空间进行网格搜索，选取最优的那个解。至于为什么动图中可以求出全局最优，是因为图中的函数只有一个极大值点，是一个凸优化问题。对于任意给定的函数求最大值，是一个至今仍没有解决的问题。

这篇博客里除了部分数学推导参考了吴恩达的cs229n笔记， 其余的内容基本是自己思考过程的总结，故难免有疏漏和不当之处，欢迎大家留言与我讨论。

ps.断断续续写了三周的时间才完工，把自己的想法以一种容易理解的方式用文字表达出来的确是一个烧脑的过程，但是写完之后自己对于算法的理解更深了一层，所带来的成就感也值得一行行敲公式的辛苦。这是我的第一篇博客，希望能将这个习惯坚持下去吧。