从似然函数到 EM 算法 (附代码实现)

什么是 EM 算法

最大期望算法（Expectation-maximization algorithm，又译为期望最大化算法），是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐性变量。

最大期望算法经过两个步骤交替进行计算，

第一步是计算期望（E），利用对隐藏变量的现有估计值，计算其最大似然估计值； 第二步是最大化（M），最大化在 E 步上求得的最大似然值来计算参数的值。M 步上找到的参数估计值被用于下一个 E 步计算中，这个过程不断交替进行。

极大似然估计用一句话概括就是：知道结果，反推条件 θ。

1.1 似然函数

在数理统计学中，似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性。“似然性”与 “或然性” 或“概率”意思相近，都是指某种事件发生的可能性。而极大似然就相当于最大可能的意思。

比如你一位同学和一位猎人一起外出打猎，一只野兔从前方窜过。只听一声枪响，野兔应声到下，如果要你推测，这一发命中的子弹是谁打的？你就会想，只发一枪便打中，由于猎人命中的概率一般大于你那位同学命中的概率，从而推断出这一枪应该是猎人射中的。

这个例子所作的推断就体现了最大似然法的基本思想。

多数情况下我们是根据已知条件来推算结果，而最大似然估计是已经知道了结果，然后寻求使该结果出现的可能性最大的条件，以此作为估计值。

1.3 极大似然函数的求解步骤

假定我们要从 10 万个人当中抽取 100 个人来做身高统计，那么抽到这 100 个人的概率就是 (概率连乘)：

$L(\theta$ %3DL(x1%2C…%2Cx_n%7C%5Ctheta)%3D%5Cprod%7Bi%3D1%7D%5E%7Bn%7Dp(x_i%7C%5Ctheta)%2C%5Ctheta%5Cin%5Cominus>)

现在要求的就是这个 $\theta$ 值，也就是使得 $L(\theta$ >) 的概率最大化，那么这时的参数 $\theta$ 就是所求。

为了便于分析，我们可以定义对数似然函数，将其变成连加的形式：

$H(\theta$ %3DlnL(%5Ctheta)%3Dln%5Cprod%7Bi%3D1%7D%5E%7Bn%7Dp(x_i%7C%5Ctheta)%3D%5Csum%7Bi%3D1%7D%5E%7Bn%7Dlnp(x_i%7C%5Ctheta)>)

对于求一个函数的极值，通过我们在本科所学的微积分知识，最直接的设想是求导，然后让导数为 0，那么解这个方程得到的 θ 就是了（当然，前提是函数 L(θ) 连续可微）。但，如果 θ 是包含多个参数的向量那怎么处理呢？当然是求 L(θ) 对所有参数的偏导数，也就是梯度了，从而 n 个未知的参数，就有 n 个方程，方程组的解就是似然函数的极值点了，最终得到这 n 个参数的值。

求极大似然函数估计值的一般步骤：

写出似然函数；
对似然函数取对数，并整理；
求导数，令导数为 0，得到似然方程；
解似然方程，得到的参数即为所求；

1.4 EM 算法

两枚硬币 A 和 B，假定随机抛掷后正面朝上概率分别为 PA，PB。为了估计这两个硬币朝上的概率，咱们轮流抛硬币 A 和 B，每一轮都连续抛 5 次，总共 5 轮：

硬币	结果	统计
A	正正反正反	3 正 - 2 反
B	反反正正反	2 正 - 3 反
A	正反反反反	1 正 - 4 反
B	正反反正正	3 正 - 2 反
A	反正正反反	2 正 - 3 反

硬币 A 被抛了 15 次，在第一轮、第三轮、第五轮分别出现了 3 次正、1 次正、2 次正，所以很容易估计出 PA，类似的，PB 也很容易计算出来 (真实值)，如下：

PA = （3+1+2）/ 15 = 0.4 PB= （2+3）/10 = 0.5

问题来了，如果我们不知道抛的硬币是 A 还是 B 呢（即硬币种类是隐变量），然后再轮流抛五轮，得到如下结果：

硬币	结果	统计
Unknown	正正反正反	3 正 - 2 反
Unknown	反反正正反	2 正 - 3 反
Unknown	正反反反反	1 正 - 4 反
Unknown	正反反正正	3 正 - 2 反
Unknown	反正正反反	2 正 - 3 反

OK，问题变得有意思了。现在我们的目标没变，还是估计 PA 和 PB，需要怎么做呢？

显然，此时我们多了一个硬币种类的隐变量，设为 z，可以把它认为是一个 5 维的向量（z1,z2,z3,z4,z5)，代表每次投掷时所使用的硬币，比如 z1，就代表第一轮投掷时使用的硬币是 A 还是 B。

但是，这个变量 z 不知道，就无法去估计 PA 和 PB，所以，我们必须先估计出 z，然后才能进一步估计 PA 和 PB。
可要估计 z，我们又得知道 PA 和 PB，这样我们才能用极大似然概率法则去估计 z，这不是鸡生蛋和蛋生鸡的问题吗，如何破？

答案就是先随机初始化一个 PA 和 PB，用它来估计 z，然后基于 z，还是按照最大似然概率法则去估计新的 PA 和 PB，然后依次循环，如果新估计出来的 PA 和 PB 和我们真实值差别很大，直到 PA 和 PB 收敛到真实值为止。

我们不妨这样，先随便给 PA 和 PB 赋一个值，比如：硬币 A 正面朝上的概率 PA = 0.2 硬币 B 正面朝上的概率 PB = 0.7

然后，我们看看第一轮抛掷最可能是哪个硬币。如果是硬币 A，得出 3 正 2 反的概率为 0.2_0.2_0.2_0.8_0.8 = 0.00512 如果是硬币 B，得出 3 正 2 反的概率为 0.7_0.7_0.7_0.3_0.3=0.03087 然后依次求出其他 4 轮中的相应概率。做成表格如下：

轮数	若是硬币 A	若是硬币 B
1	0.00512，即 0.2 0.2 0.2 0.8 0.8，3 正 - 2 反	0.03087，3 正 - 2 反
2	0.02048，即 0.2 0.2 0.8 0.8 0.8，2 正 - 3 反	0.01323，2 正 - 3 反
3	0.08192，即 0.2 0.8 0.8 0.8 0.8，1 正 - 4 反	0.00567，1 正 - 4 反
4	0.00512，即 0.2 0.2 0.2 0.8 0.8，3 正 - 2 反	0.03087，3 正 - 2 反
5	0.02048，即 0.2 0.2 0.8 0.8 0.8，2 正 - 3 反	0.01323，2 正 - 3 反

按照最大似然法则：第 1 轮中最有可能的是硬币 B 第 2 轮中最有可能的是硬币 A 第 3 轮中最有可能的是硬币 A 第 4 轮中最有可能的是硬币 B 第 5 轮中最有可能的是硬币 A

我们就把概率更大，即更可能是 A 的，即第 2 轮、第 3 轮、第 5 轮出现正的次数 2、1、2 相加，除以 A 被抛的总次数 15（A 抛了三轮，每轮 5 次），作为 z 的估计值，B 的计算方法类似。然后我们便可以按照最大似然概率法则来估计新的 PA 和 PB。

PA = （2+1+2）/15 = 0.33 PB =（3+3）/10 = 0.6

就这样，不断迭代不断接近真实值，这就是 EM 算法的奇妙之处。

可以期待，我们继续按照上面的思路，用估计出的 PA 和 PB 再来估计 z，再用 z 来估计新的 PA 和 PB，反复迭代下去，就可以最终得到 PA = 0.4，PB=0.5，此时无论怎样迭代，PA 和 PB 的值都会保持 0.4 和 0.5 不变，于是乎，我们就找到了 PA 和 PB 的最大似然估计。

总结一下计算步骤：

随机初始化分布参数 θ
E 步，求 Q 函数，对于每一个 i，计算根据上一次迭代的模型参数来计算出隐性变量的后验概率（其实就是隐性变量的期望），来作为隐藏变量的现估计值：

$Q_i(z^{(i$ %7D)%3Dp(z%5E%7B(i)%7D%7Cx%5E%7B(i)%7D%3B%5Ctheta)>)
M 步，求使 Q 函数获得极大时的参数取值）将似然函数最大化以获得新的参数值

$\theta=argmax\sum_{i}\sum_{z^{(i$ %7D%7DQ_i(z%5E%7B(i)%7D)log%5Cfrac%7Bp(x%5E%7B(i)%7D%2Cz%5E%7B(i)%7D%3B%5Ctheta)%7D%7BQ_i(z%5E%7B(i)%7D)%7D>)
然后循环重复 2、3 步直到收敛。