最近在上机器翻译的课程。作为机器翻译领域的统计学模型的开山模型的核心原理，姑且记录一下吧。这个算法本科学过，也考过，但是最近听课发现，突然又不会了，好多细节不清晰了。
EM算法也叫做期望最大化算法，Expectation-Maximum Algorithm。这个算法主要用于解决含有隐变量的概率估计问题。我们先不进行复杂的数学推导，先给出一个例子。
例子：翻硬币
假设我们有三枚硬币，分别命名为零号、一号和二号，三枚硬币翻到正面的概率分别为$\lambda$、$p_1$和$p_2$，规则如下：
- 当零号硬币翻到正面，就翻三次一号硬币。
- 当零号硬币翻到反面，就翻三次二号硬币。
- 最后会出现三次结果，但是不告诉你这三次结果是来自一号还是二号硬币，你需要根据这三次结果预测三个参数。
假设，做了3次实验，最后的结果是（正正正）、（正反正）、（正正反）。
现在看来，这好像没法预测参数啊，虽然我知道了硬币的结果，但是我不知道来自哪个硬币，这怎么去预测呢？我们先把这个问题建模下来，我们现在有了三个结果，可以记作$s_i,i\in\{1,2,3\}$，其中，每个 $s_i$ 代表一次实验中观察到的三次翻硬币结果。为了之后写起来方便，我们不直接记录“正反正”这种字符串，而是记录其中正面出现的次数。比如：
$$ s_1 = \text{正正正}, \quad h_1 = 3 $$
$$ s_2 = \text{正反正}, \quad h_2 = 2 $$
$$ s_3 = \text{正正反}, \quad h_3 = 2 $$
这里的 $h_i$ 表示第 $i$ 次实验中正面出现的次数。因为每次实验都会翻三次硬币，所以反面出现的次数就是 $3-h_i$。但是问题的关键在于：我们不知道每次实验到底是由一号硬币产生的，还是由二号硬币产生的。于是，我们引入一个隐变量$z_i \in {1,2}$，当 $z_i=1$ 的时候，表示第 $i$ 次实验来自一号硬币；当 $z_i=2$ 的时候，表示第 $i$ 次实验来自二号硬币。这里的 $z_i$ 就是所谓的“隐变量”。它在真实的数据生成过程中确实存在，但是我们在观测数据中看不到它。这里的思想其实很简单，既然有我们不知道的变量，就先假设我们知道它，然后我们看看有没有好办法去处理。如果我们知道 $z_i$，那问题就非常简单了。比如我们知道：
- 第一次实验来自一号硬币；
- 第二次实验来自二号硬币；
- 第三次实验来自一号硬币。
那么我们就可以直接数数：
- 零号硬币正面出现了几次，反面出现了几次；
- 一号硬币一共翻了多少次，其中正面多少次；
- 二号硬币一共翻了多少次，其中正面多少次。
然后用频率估计概率即可。比如，一号硬币一共被翻了 $6$ 次，其中正面出现了 $5$ 次，那么自然可以估计：
$$ p_1 = \frac{5}{6} $$
问题是，我们偏偏不知道 $z_i$。所以 EM 算法要解决的问题，本质上就是：
如果隐变量不知道，那我能不能先“猜”一下隐变量，再利用这个猜测去更新参数；然后再用更新后的参数重新猜隐变量，如此反复？
这句话其实就是 EM 算法的核心思想。
如果直接写似然函数
我们先看看如果直接写出观测数据的似然函数，会发生什么。对于某一次实验 $s_i$，它可能来自一号硬币，也可能来自二号硬币。如果来自一号硬币，那么概率是：$\lambda p_1^{h_i}(1-p_1)^{3-h_i}$
，这里的 $\lambda$ 表示零号硬币翻到正面的概率，也就是选择一号硬币的概率。如果来自二号硬币，那么概率是：$(1-\lambda)p_2^{h_i}(1-p_2)^{3-h_i}$，所以，对于一次观测 $s_i$，它出现的总概率应该是这两种情况加起来，也就是$\lambda p_1^{h_i}(1-p_1)^{3-h_i}+(1-\lambda)p_2^{h_i}(1-p_2)^{3-h_i}$。三次实验相互独立，所以整体似然函数是：
$$ \prod_{i=1}^3 \left[ \lambda p_1^{h_i}(1-p_1)^{3-h_i} + (1-\lambda)p_2^{h_i}(1-p_2)^{3-h_i} \right] $$
为了方便优化，我们一般取对数，得到对数似然函数：
$$ \sum_{i=1}^3 \log \left[ \lambda p_1^{h_i}(1-p_1)^{3-h_i} + (1-\lambda)p_2^{h_i}(1-p_2)^{3-h_i} \right] $$
到这里，麻烦就出现了。如果括号里面没有加法，而只是单纯的乘积，那么取对数之后就能拆开，问题会变得很简单。但是现在对数里面有一个加法$\log(a+b)$这个东西没法直接拆成$\log a + \log b$，也就是说，隐变量的存在让似然函数变得不好直接最大化。EM 算法就是为了解决这种问题出现的。
如果隐变量已知
为了看清楚 EM 算法为什么有用，我们先假装自己知道隐变量。
假设 $z_i=1$ 表示第 $i$ 次实验来自一号硬币，$z_i=2$ 表示来自二号硬币。那么完整数据应该是：$(s_i,z_i)$，这时候我们不仅知道观测结果，还知道它是由哪枚硬币生成的。完整数据的似然函数就会简单很多。如果 $z_i=1$，那么第 $i$ 次实验的概率是$\lambda p_1^{h_i}(1-p_1)^{3-h_i}$，如果 $z_i=2$，那么第 $i$ 次实验的概率是：$(1-\lambda)p_2^{h_i}(1-p_2)^{3-h_i}$，这时候我们可以写出完整数据的对数似然函数：
$$ \sum_{i=1} \left[ \log \lambda + h_i \log p_1 + (3-h_i)\log(1-p_1) \right] + \sum_{i=2} \left[ \log(1-\lambda) + h_i \log p_2 + (3-h_i)\log(1-p_2) \right] $$
这个式子就很舒服了，它不过是我们之前写的「对数似然函数」的变形，前后两部分分别表示来自于一号的概率和来自于二号的概率。如果 $z_i$ 已知，那么每一项属于哪枚硬币都是明确的。最后对参数求最大值，其实就是做频率估计，三个参数也就分别是$\frac{\text{选择一号硬币的次数}}{\text{总实验次数}}$、$\frac{\text{一号硬币中正面出现的总次数}}{\text{一号硬币被翻的总次数}}$以及$\frac{\text{二号硬币中正面出现的总次数}}{\text{二号硬币被翻的总次数}}$。但是现实是，$z_i$ 不知道。那怎么办？EM 的想法是：既然不知道，那就不要粗暴地给每个 $z_i$ 一个确定答案，而是给它一个“概率上的猜测”。也就是说，对于每次实验，我们不说它一定来自一号硬币或者一定来自二号硬币，而是说：
这次实验有多大概率来自一号硬币？有多大概率来自二号硬币？
这就是 E 步。
E 步：估计隐变量的后验概率
假设当前参数是：$\lambda^{old}, p_1^{old}, p_2^{old}$，那么对于第 $i$ 次实验，我们可以计算它来自一号硬币的概率——$P(z_i=1 \mid s_i)$。
根据贝叶斯公式，有：
$$ q_i=\frac{ \lambda^{old}(p_1^{old})^{h_i}(1-p_1^{old})^{3-h_i} }{ \lambda^{old}(p_1^{old})^{h_i}(1-p_1^{old})^{3-h_i} + (1-\lambda^{old})(p_2^{old})^{h_i}(1-p_2^{old})^{3-h_i} } $$
如果这个数值等于0.9，就表示“这次实验有 90% 的责任可以算到一号硬币头上”。在 EM 算法的相关文章里，经常会看到一个词叫 responsibility，中文一般翻译成“责任”或者“责任度”。某个观测数据到底应该由哪个隐藏成分负责？不是硬分配，而是按概率分配责任。我个人更偏向于将这个东西叫做「软计数」，因为它实际上做的事，和计数差不多。
M 步：用责任度重新估计参数
有了软计数$q_i$以后，我们就可以重新估计参数。首先，$\lambda$ 表示选择一号硬币的概率。现在每次实验不是明确属于一号硬币，而是以 $q_i$ 的比例属于一号硬币。所以新的 $\lambda$ 应该是$\frac{1}{N}\sum_{i=1}^N q_i$，其中 $N$ 是实验次数。在这个例子里，$N=3$。接下来估计 $p_1$。一号硬币的正面次数，不再是简单地把某几次实验加起来，而是要按 $q_i$ 加权：
$$ \frac{ \sum_{i=1}^N q_i h_i }{ \sum_{i=1}^N q_i \cdot 3 } $$
分子表示“一号硬币负责的正面次数”，分母表示“一号硬币负责的总翻硬币次数”。同理，二号硬币的参数更新为：
$$ \frac{ \sum_{i=1}^N (1-q_i)h_i }{ \sum_{i=1}^N (1-q_i)\cdot 3 } $$
这就是 M 步。
M 步的名字叫 Maximization，也就是最大化。它做的事情是：在当前对隐变量的估计下，找到一组更好的参数。
所以 EM 算法的两个步骤可以非常粗略地理解为：
- E 步：用当前参数估计隐变量；
- M 步：用估计出来的隐变量更新参数。
然后不断重复这两个步骤。
带入一个具体数值
为了更直观一点，我们随便初始化一组参数：
$$ \lambda=0.5,\quad p_1=0.8,\quad p_2=0.4 $$
这里的意思是：
- 零号硬币有一半概率选择一号硬币；
- 一号硬币比较容易出正面；
- 二号硬币出正面的概率较低。
现在看三次实验：
$$ h_1=3,\quad h_2=2,\quad h_3=2 $$
对于第一次实验，也就是正正正，它来自一号硬币的概率为：
$$ \frac{0.256}{0.256+0.032} \approx 0.889 $$
这很符合直觉。因为正正正更像是由正面概率较高的一号硬币产生的。
对于第二次和第三次实验，都是两正一反，所以：
$$ \frac{0.064}{0.064+0.048} \approx 0.571 $$
于是 E 步得到：
$$ q_1 \approx 0.889,\quad q_2 \approx 0.571,\quad q_3 \approx 0.571 $$
接下来做 M 步。
先更新 $\lambda$：
$$ \frac{0.889+0.571+0.571}{3} \approx 0.677 $$
再更新 $p_1$：
$$ \frac{ 0.889\times 3 + 0.571\times 2 + 0.571\times 2 }{ 3(0.889+0.571+0.571) } \approx 0.813 $$
再更新 $p_2$：
$$ \frac{ (1-0.889)\times 3 + (1-0.571)\times 2 + (1-0.571)\times 2 }{ 3((1-0.889)+(1-0.571)+(1-0.571)) } \approx 0.734 $$
所以经过一轮 EM 之后，参数从：
$$ \lambda=0.5,\quad p_1=0.8,\quad p_2=0.4 $$
变成了：
$$ \lambda\approx 0.677,\quad p_1\approx 0.813,\quad p_2\approx 0.734 $$
这个结果可能看起来有点奇怪：为什么 $p_2$ 从 $0.4$ 一下子变成了 $0.734$？原因也很简单：我们的观测数据里正面太多了。三次实验分别是 $3,2,2$ 个正面，总共 $9$ 次翻硬币里有 $7$ 次正面。所以即使一开始我们假设二号硬币不太容易出正面，数据也会把它往“更容易出正面”的方向拉。
这也提醒我们一件事：EM 是根据当前数据和当前参数做迭代的。
总结
EM算法要解决的问题是：
我想估计概率模型的参数，但是数据里有一些变量看不见。
于是它采用了一个迭代策略：
- 先随机初始化参数；
- 用当前参数估计隐变量的后验分布；
- 用隐变量的后验分布重新估计参数；
- 重复 2 和 3，直到收敛。