EM算法原理

前言

EM 算法，全称 Expectation Maximization Algorithm。期望最大算法是一种迭代算法，用于含有隐变量（Hidden Variable）的概率参数模型的最大似然估计或极大后验概率估计。它是一个基础算法，是很多机器学习领域算法的基础，比如隐式马尔科夫算法（HMM）， LDA 主题模型的变分推断等等。

一、EM算法的思想

我们经常会从被观察的样本数据中，找出样本的模型参数。 最常用的方法就是极大化模型分布的对数似然函数。

但是在一些情况下，我们得到的观察数据有未观察到的隐含数据，此时我们未知的有隐含数据和模型参数，因而无法直接用极大化对数似然函数得到模型分布的参数。怎么办呢？这就是 EM 算法可以派上用场的地方了。

EM 算法解决这个的思路是使用启发式的迭代方法，既然我们无法直接求出模型分布参数，那么我们可以先猜想隐含数据（ EM 算法的 Expection-Step ，又称 E步），接着基于观察数据和猜测的隐含数据一起来极大化对数似然，求解我们的模型参数（ EM 算法的 Maximization-Step ，又称 M步 )。由于我们之前的隐藏数据是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。不过没关系，我们基于当前得到的模型参数，继续猜测隐含数据（ EM 算法的 E步 ），然后继续极大化对数似然，求解我们的模型参数（ EM 算法的 M步 )。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。

从上面的描述可以看出， EM 算法是迭代求解最大值的算法，同时算法在每一次迭代时分为两步， E步 和 M步 。一轮轮迭代更新隐含数据和模型分布参数，直到收敛，即得到我们需要的模型参数。

E-Step 主要通过观察数据和现有模型来估计参数，然后用这个估计的参数值来计算上述对数似然函数的期望值；而 M-Step 是寻找似然函数最大化时对应的参数。由于算法会保证在每次迭代之后似然函数都会增加（每一次迭代都是 Maximization ），所以函数最终会收敛。

一个最直观了解 EM 算法思路的是 K-Means 算法，在 K-Means 聚类时，每个聚类簇的质心是隐含数据。我们会假设 K 个初始化质心，即 EM 算法的 E步 ；然后计算得到每个样本最近的质心，并把样本聚类到最近的这个质心，即 EM 算法的 M步 。重复这个 E步 和 M步 ，直到质心不再变化为止，这样就完成了 K-Means 聚类。

当然， K-Means 算法是比较简单的，实际中的问题往往没有这么简单。上面对 EM 算法的描述还很粗糙，我们需要用数学的语言精准描述。

二、EM算法的推导

对于 m 个样本观察数据 $x=(x^{(1)},x^{(2)},...,x^{(m)})$ 中，找出样本的模型参数 θ , 极大化模型分布的对数似然函数如下：
$$L(\theta )=\sum _{i=1}^{m}logP(x^{(i)};\theta )\text{⋯⋯ (一)}$$
$$=\sum _{i=1}^{m}log\sum _{z}P(x^{(i)},z;\theta )\text{⋯⋯ (二)}$$
第(一)步是对极大似然函数取对数，第(二)步是对每个样本的每个可能的类别 z 求联合概率分布之和。如果这个 z 是已知的数，那么使用极大似然法求解参数 $\theta$ 会很容易。但如果 z 是隐变量，我们就需要用 EM 算法来求解。
事实上，隐变量估计问题也可以通过梯度下降等优化算法，但事实由于求和项将随着隐变量的数目以指数级上升，会给梯度计算带来麻烦；而 EM 算法则可看作一种非梯度优化方法。

如果我们得到的观察数据有未观察到的隐含数据 $z=(z^{(1)},z^{(2)},...,z^{(m)})$ ，此时我们的极大化模型分布的对数似然函数如下：
$$L(\theta )=\sum _{i=1}^{m}logP(x^{(i)};\theta )=\sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)},z^{(i)};\theta )$$
上面这个式子是没有 办法直接求出 θ 的。因此需要一些特殊的技巧，我们首先对这个式子进行缩放如下：
$$\sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)},z^{(i)};\theta )=\sum _{i=1}^{m}log\sum _{z^{(i)}}{Q}_i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\text{⋯⋯①}$$
$$\ge \sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\text{⋯⋯②}$$
上面第①式引入了一个未知的新的分布 $Q_i(z^{(i)})$ ，第②式用到了Jensen不等式：
$$log\sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_{j}log(y_j)$$
其中：
$${\lambda }_j\ge 0,\sum _j{\lambda }_j=1$$
或者说由于对数函数（以自然数为底）是凹函数，所以有：
$$f(E(x))\ge E(f(x))$$
此时如果要满足Jensen不等式的等号，则有：
$$\frac{P(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}=c\text{c为常数}$$
由于 $Q_i(z^{(i)})$ 是一个分布，所以满足：
$$\sum _{z^{(i)}}{Q}_i(z^{(i)})=1\text{⋯⋯③}$$
从上面第①式和第③式，我们可以得到：
$$Q_i(z^{(i)})=\frac{P(x^{(i)},z^{(i)};\theta )}{{\sum }_{z^{(i)}}P(x^{(i)},z^{(i)};\theta )}=\frac{P(x^{(i)},z^{(i)};\theta )}{P(x^{(i)};\theta )}=P(z^{(i)}|x^{(i)};\theta )$$
如果 $Q_i(z^{(i)})=P(z^{(i)}|x^{(i)};\theta ))$ , 则第②式是我们的包含隐藏数据的对数似然的一个下界。如果我们能极大化这个下界，则也在尝试极大化我们的对数似然。即我们需要最大化下式：
$$\arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$$
去掉上式中为常数的部分，则我们需要极大化的对数似然下界为：
$$\arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})logP(x^{(i)},z^{(i)};\theta )$$
即通过上面的分析我们得到了 $L(\theta )\ge J(z,Q)$ 的形式（z为隐变量），那么我们就可以通过不断优化下界 $J(z,Q)$ 来使得对数似然函数 $L(\theta )$ 不断提高，如下图所示：

上式也就是我们的 EM 算法的 M步 ，那 E步 呢？注意到上式中 $Q_i(z^{(i)})$ 是一个分布，因此，$Q_i(z^{(i)})logP(x^{(i)},z^{(i)};\theta )$ 可以理解为 $logP(x^{(i)},z^{(i)};\theta )$ 基于条件概率分布 $Q_i(z^{(i)})$ 的期望。

EM算法迭代优化的另一种理解：坐标上升法（Coordinate ascent）

途中直线为迭代优化的路径，因为每次只优化一个变量（参数），所以可以看到它每走一步都是平行于坐标轴的。
EM 算法类似于坐标上升法，E 步 ：固定参数 $\theta$ ，优化 Q；M 步 ：固定 Q，优化参数 $\theta$ 。交替将极值推向最大。

但是要注意，迭代一定会收敛，但不一定会收敛到真实的全局最优参数值，因为可能会陷入局部最优。所以 EM 算法的结果很受初始值的影响。

三、EM算法流程

输入：观察数据 $x=(x^{(1)},x^{(2)},...,x^{(m)})$ ，联合分布 $p(x,z;\theta )$ , 条件分布 $p(z|x;\theta )$ , 最大迭代次数 $J$ 。

• ① 随机初始化模型参数 $\theta$ 的初值 ${\theta }_0$ 。

• ② for j from 1 to J开始 EM 算法迭代：
  a) E步 ：计算联合分布的条件概率期望：
  $$Q_i(z^{(i)})=P(z^{(i)}|x^{(i)};\theta )$$
  $$L(\theta ,{\theta }^j)=\sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})logP(x^{(i)},z^{(i)};\theta )$$
  b) M步 ：极大化 $L(\theta ,{\theta }_j)$ ,得到 ${\theta }^{j+1}$ ：
  $${\theta }^{j+1}=\arg \underset{\theta }{\max }L(\theta ,{\theta }^j)$$
  c) 如果 ${\theta }^{j+1}$ 已收敛，则算法结束。否则继续回到 步骤 a) 进行 E步 迭代。

输出：模型参数 $\theta$ 。

四、举例

为了让大家更直观的感受EM算法，在这里给大家举两个例子来说明。

例一

我们将引用 Nature Biotech 的 EM tutorial 文章中的例子。

假设有两枚硬币 A 和 B ，他们的随机抛掷的结果如下图所示：

我们很容易估计出两枚硬币抛出正面的概率：
$${\theta }_A=\frac{24}{30}=0.80$$
$${\theta }_B=\frac{9}{20}=0.45$$
现在我们加入隐变量，抹去每轮投掷的硬币标记：

Coin	Statistics
Coin *	5 H, 5 T
Coin *	9 H, 1 T
Coin *	8 H, 2 T
Coin *	4 H, 6 T
Coin *	7 H, 3 T

碰到这种情况，我们该如何估计 ${\theta }_A$ 和 ${\theta }_B$ 的值？
我们多了一个隐变量 $Z=(z_1,z_2,z_3,z_4,z_5)$ ，代表每一轮所使用的硬币，我们需要知道每一轮抛掷所使用的硬币这样才能去估计 ${\theta }_A$ 和 ${\theta }_B$ 的值，但是估计隐变量 $Z$ 我们又需要知道 ${\theta }_A$ 和 ${\theta }_B$ 的值，才能用极大似然估计法去估计出 $Z$ 。这就陷入了一个鸡生蛋和蛋生鸡的问题。
其解决方法就是先随机初始化 ${\theta }_A$ 和 ${\theta }_B$ ，然后去估计 $Z$ ， 然后基于 $Z$ 按照最大似然概率去估计新的 ${\theta }_A$ 和 ${\theta }_B$ ，循环直至收敛。

接下来我们随机初始化 ${\theta }_A=0.6$ 和 ${\theta }_B=0.5$ ，按照上述的解决方案先计算第一轮的数据：
对于第一轮来说，如果使用的是硬币 A ，得出的 5 正 5 反 的概率为：$0.6^5\times 0.4^5$ ；如果使用的是硬币 B ，得出的 5 正 5 反 的概率为：$0.5^5\times 0.5^5$ 。我们可以算出使用的是硬币 A 和硬币 B 的概率分别为：
$$P_A=\frac{0.6^5\times 0.4^5}{0.6^5\times 0.4^5+0.5^5\times 0.5^5}=0.45$$
$$P_B=\frac{0.5^5\times 0.5^5}{0.6^5\times 0.4^5+0.5^5\times 0.5^5}=0.55$$
从期望的角度来看，对于第一轮抛掷，使用硬币 A 的概率是 0.45 ，使用硬币 B 的概率是 0.55 。同理，其他轮的概率计算结果为：

No	Coin A	Coin B
1	0.45	0.55
2	0.80	0.20
3	0.73	0.27
4	0.35	0.65
5	0.65	0.35

这一步我们实际上是估计出了 $Z$ 的概率分布，这步就是 E-Step（E步） 。

结合硬币 A 的概率和上面的投掷结果，我们利用期望可以求出硬币 A 和硬币 B 的贡献。以第二轮硬币 A 为例子，计算方式为：
$$H:0.80\times 9=7.2$$
$$T:0.80\times 1=0.8$$
同理，于是我们可以得到：

No	Coin A	Coin B
1	2.2 H, 2.2 T	2.8 H, 2.8 T
2	7.2 H, 0.8 T	1.8 H, 0.2 T
3	5.9 H, 1.5 T	2.1 H, 0.5 T
4	1.4 H, 2.1 T	2.6 H, 3.9 T
5	4.5 H, 1.9 T	2.5 H, 1.1 T
Total	21.3 H, 8.6 T	11.7 H, 8.4 T

然后再用极大似然估计来估计新的 $P_A$ 和 $P_B$ ：
$$P_A=\frac{21.3}{21.3+8.6}=0.71$$
$$P_B=\frac{11.7}{11.7+8.4}=0.58$$
这步就对应了 M-Step（M步） ，重新估计出了期望值。
如此反复迭代，我们就可以算出最终的参数值。
上述讲解对应下图：

例二

如果说 例一 需要繁琐的计算，而这对于你来说可能没那么直观，那我们接下来就举一个更直观的例子：

现在一个班里有 50 个男生和 50 个女生，且男女生分开。我们假定男生的身高服从正态分布：$\mathcal{N}({\mu }_1,{\sigma }_1^2)$ ，女生的身高则服从另一个正态分布：$\mathcal{N}({\mu }_2,{\sigma }_2^2)$ 。这时候我们可以用极大似然估计（MLE），分别通过这 50 个男生和 50 个女生的样本来估计这两个正态分布的参数。

但现在我们让情况稍微复杂一点，就是这 50 个男生和 50 个女生混在一起了。我们拥有 100 个人的身高数据，却不知道这 100 个人哪个是男生或者哪个是女生。
这个时候情况就有点尴尬，因为通常来说，我们只有知道了精确的男女身高的正态分布参数我们才能知道这 100 个人中某个人更有可能是男生还是女生。但从另一方面去考量，我们只有知道了每个人是男生还是女生才能尽可能准确地估计男女各自身高的正态分布的参数。
这个时候有人就想到我们必须从某一点开始，并用迭代的办法去解决这个问题：我们先设定男生身高和女生身高分布的几个参数（初始值），然后根据这些参数去判断每一个样本（学生）是男生还是女生，之后根据标注后的样本再反过来重新估计参数。之后再多次重复这个过程，直至收敛。

五、参考文献

① EM算法原理总结 - 刘建平
② 机器学习 - 周志华 【162-165】
③ https://www.zhihu.com/question/27976634
④ https://blog.csdn.net/zouxy09/article/details/8537620
⑤ Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm?. Nature biotechnology, 26(8), 897.
⑥ 一文详尽系列之EM算法 - Datawhale