视频链接
数据集下载地址：无需下载

学习目标：

了解什么是 EM 算法
知道极大似然估计
知道 EM 算法实现流程

讲 EM 算法主要是为了后面的 HMM 做准备。

1. 初始 EM 算法

EM 算法（Expectation-Maximization algorithm，期望最大化算法）是一种迭代算法，用于在概率模型中寻找最大似然估计或最大后验估计，特别适用于模型中存在隐变量的情况。

EM 算法是一个基础算法，是很多机器学习领域算法的基础，比如隐式马尔科夫算法（HMM）等。EM 算法是一种迭代优化策略，由于它的计算方法中每一次迭代都分两步，所以算法被称为EM 算法：

一个为期望步（E 步）
另一个为极大步（M步）

EM 算法受到缺失思想影响，最初是为了解决数据缺失情况下的参数估计问题，其算法基础和收敛有效性等问题在 Dempster、Laird 和 Rubin 三人于 1977 年所做的文章《Maximum Likelihood from Incomplete Data via the EM Algorithm》中给出了详细的阐述。其基本思想是：

首先根据己经给出的观测数据，估计出模型参数的值
然后再依据上一步估计出的参数值估计缺失数据的值，再根据估计出的缺失数据加上之前己经观测到的数据重新再对参数值进行估计
然后反复迭代，直至最后收敛，迭代结束。

EM 算法计算流程：

在这里插入图片描述

2. EM 算法介绍

学习目标：

知道什么是极大似然估计
知道 EM 算法实现流程

想清晰的了解 EM 算法，我们需要知道一个基础知识：“极大似然估计”。

2.1 极大似然估计（Maximum Likelihood Estimation，MLE）

极大似然估计（Maximum Likelihood Estimation，MLE），通俗理解来说，就是利用已知的样本结果信息，反推最具有可能（最大概率）导致这些样本结果出现的模型参数值！换句话说，极大似然估计提供了一种给定观察数据来评估模型参数的方法，即：“模型已定，参数未知”。

2.1.1 问题描述

假如我们需要调查学校的男生和女生的身高分布，我们抽取 100 个男生和 100 个女生，将他们按照性别划分为两组。然后，统计抽样得到 100 个男生的身高数据和 100 个女生的身高数据。

如果我们知道他们的身高服从正态分布，但是这个分布的均值 $\mu$ 和方差 $\sigma^2$ 是不知道，这两个参数就是我们需要估计的。

问题：我们知道样本所服从的概率分布模型和一些样本，我们需要求解该模型的参数。

在这里插入图片描述

我们已知的条件有两个：

样本服从的分布模型
随机抽取的样本

我们需要求解模型的参数。即，根据已知条件，通过极大似然估计，求出未知参数。

总的来说，极大似然估计就是用来估计模型参数的统计学方法。

2.1.2 用数学知识解决现实问题

问题数学化：

样本集： $X = \{ x_1, x_2, ..., x_n \}, n=100$
概率密度 $p(x_i|\theta)$ 表示抽到第 $i$ 个男生身高的概率
由于 100 个样本之间独立同分布，所以同时抽到这 100 个男生的概率是它们各自概率的乘积，也就是样本集 $X$ 中各个样本的联合概率，用下式表示： $L(\theta) = L(x_1, x_2, ..., x_n; \theta) = \prod_{i=1}^n p(x_i; \theta), \theta \in \Theta$

这个概率反映了在概率密度函数的参数是 $\theta$ 时，得到 $X$ 这组样本的概率。

我们需要找到一个参数 $\theta$ ，使得抽到 $X$ 这组样本的概率最大，也就是说需要其对应的似然函数 $L(\theta)$ 最大。

满足条件的 $\theta$ 叫做 $\theta$ 的最大似然估计值，记为： $\hat{\theta} = \mathrm{argmax} \ L(\theta)$

2.1.3 最大似然函数估计值的求解步骤

第一步：首先写出似然函数

$L(\theta) = L(x_1, x_2, ..., x_n; \theta) = \prod_{i=1}^n p(x_i; \theta), \theta \in \Theta$

第二步：对似然函数取对数

$l(\theta) = \ln L(\theta) = \ln{\prod_{i=1}^n p(x_i; \theta)} = \sum_{i=1}^n \ln{p(x_i; \theta)}$

第三步：对上式求导，并令导数为 0，得到似然方程

第四步：求解似然方程，得到的参数值即为我们要求的

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

Q：为什么要对似然函数取对数？
A：【文间跳转】极大似然函数取对数的原因

2.2 EM 算法实例描述

我们目前有 100 个男生和 100 个女生的身高，但是我们不知道这 200 个数据中哪个是男生的身高，哪个是女生的身高，即抽取得到的每个样本都不知道是从哪个分布中抽取的。

这个时候，对于每个样本，就有两个未知量需要估计：

这个身高数据是来自于男生数据集合还是来自于女生数据集合？
男生、女生身高数据集的正态分布的参数（ $\mu$ 和 $\sigma^2$ ）分别是多少？

具体问题如下图所示：

在这里插入图片描述

对于具体的身高问题使用 EM 算法求解步骤如下：

在这里插入图片描述

步骤一 · 初始化参数：先初始化男生身高的正态分布的参数，如均值 $\mu=1.65$ ，方差 $\sigma^2 =0.15$ ；

步骤二 · 计算分布：计算每一个人更可能属于男生分布或者女生分布；

步骤三 · 重新估计参数：通过分为男生的 $n$ 个人来重新估计男生身高分布的参数（最大似然估计），女生分布也按照相同的方式估计出来，更新分布；

步骤四 · 迭代：这时候两个分布的概率也变了，然后重复步骤一至三，直到参数不发生变化为止。

2.3 EM 算法流程

输入：

$n$ 个样本观察数据 $x = (x_1, x_2, ..., x_n)$ ，未观察到的隐含数据 $z = (z_1, z_2, ..., z_n)$
联合分布 $\theta)$
条件分布 $p(z|x;\theta)$
最大迭代次数 $J$

算法步骤：

步骤一 · 初始化参数：随机初始化模型参数 $\theta$ 的初值 $\theta_0$ 。
步骤二 · 计算分布： $j = 1, 2, ..., J$ 开始 EM 算法迭代：
- E 步：计算联合分布的条件概率期望 $Q_i(z_i) = p(z_i |x_i; \theta_j)$ $l(\theta, \theta_j) = \sum_{i=1}^n\sum_{z_i}Q_i(z_i)\log{\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}}$
- M 步：极大化 $l(\theta, \theta_j)$ ，得到 $\theta_{j+1}$ $\theta_{j+1} = \mathrm{argmax}(\theta, \theta_j)$
- 迭代：如果 $\theta_{j+1}$ 已经已经收敛，则算法结束。否则继续进行 E 步和 M 步进行迭代。

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

3. EM 算法实例

学习目标：

通过实例了解 EM 算法实现的流程

3.1 一个超级简单的案例

假设现在有两枚硬币 1 和 2，随机抛掷后正面朝上概率分别为 $P_1$ ， $P_2$ 。为了估计这两个概率，做实验，每次取一枚硬币，连掷 5 下，记录下结果，如下表所示：

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

可以很容易地估计出 $P_1$ 和 $P_2$ ，如下：

$\begin{aligned} & P_1 = \frac{3+1+2}{15} = 1.4\\ & P_2 = \frac{2+3}{10} = 0.5 \end{aligned}$

到这里，一切似乎很美好，下面我们加大难度。

3.2 加入隐变量 $z$ 后的求解

还是上面的问题，现在我们抹去每轮投掷时使用的硬币标记，如下：

硬币	结果	统计
Unknown（未知）	正正反正反	3 正 2 反
Unknown（未知）	反反正正反	2 正 3 反
Unknown（未知）	正反反反反	1 正 4 反
Unknown（未知）	正反反正正	3 正 2 反
Unknown（未知）	反正正反反	2 正 3 反

好了，现在我们的目标没变，还是估计 $P_1$ 和 $P_2$ ，要怎么做呢？

显然，此时我们多了一个隐变量 $z$ ，可以把它认为是一个 5 维的向量 $z_1,z_2,z_3,z_4,z_5)$ ，代表每次投掷时所使用的硬币，比如 $z_1$ ，就代表第一轮投掷时使用的硬币是 1 还是 2。但是，这个变量 $z$ 不知道，就无法去估计 $P_1$ 和 $P_2$ 。所以，我们必须先估计出 $z$ ，然后才能进一步估计 $P_1$ 和 $P_2$ 。

但要估计 $z$ ，我们又得知道 $P_1$ 和 $P_2$ ，这样我们才能用最大似然概率法则去估计 $z$ ，这不是鸡生蛋和蛋生鸡的问题吗，如何破？

答案就是先随机初始化一个 $P_1$ 和 $P_2$ ，用它来估计 $z$ ，然后基于 $z$ ，还是按照最大似然概率法则去估计新的 $P_1$ 和 $P_2$ ，如果新的 $P_1$ 和 $P_2$ 和我们初始化的 $P_1$ 和 $P_2$ 一样，请问这说明了什么？这说明我们初始化的 $P_1$ 和 $P_2$ 是一个相当靠谱的估计（已经收敛了，不需要再进行估计了）！

就是说，我们初始化的 $P_1$ 和 $P_2$ ，按照最大似然概率就可以估计出 $z$ ，然后基于 $z$ ，按照最大似然概率可以反过来估计出 $P_1$ 和 $P_2$ 。当与我们初始化的 $P_1$ 和 $P_2$ 一样时，说明 $P_1$ 和 $P_2$ 很有可能就是真实的值。这里面包含了两个交互的最大似然估计。

如果新估计出来的 $P_1$ 和 $P_2$ 和我们初始化的值差别很大，怎么办呢？也很简单，就是继续用新的 $P_1$ 和 $P_2$ 迭代，直至收敛。

以上就是 EM 算法的初级版。

3.2.1 EM 算法初级版

我们不妨这样，先随便给 $P_1$ 和 $P_2$ 赋一个值，比如：

$P_1 = 0.2\\ P_2 = 0.7$

然后，我们看看第一轮抛掷最可能是哪个硬币。

如果是硬币 1，得出 3 正 2 反的概率为 $0.2^3 \times 0.8^2 = 0.00512$
如果是硬币 2，得出 3 正 2 反的概率为 $0.7^3 \times 0.3^2 = 0.03087$

然后依次求出其他 4 轮中的相应概率。做成表格如下：

轮数	若是硬币 1	若是硬币 2	最大概率
1（3 正 2 反）	0.00512	0.03087	硬币 2
2（2 正 3 反）	0.02048	0.01323	硬币 1
3（1 正 4 反）	0.08192	0.00567	硬币 1
4（3 正 2 反）	0.00512	0.03087	硬币 2
5（2 正 3 反）	0.02048	0.01323	硬币 1

按照最大似然法则：

第 1 轮中最有可能的是硬币 2
第 2 轮中最有可能的是硬币 1
第 3 轮中最有可能的是硬币 1
第 4 轮中最有可能的是硬币 2
第 5 轮中最有可能的是硬币 1

我们就把上面的值作为 $z$ 的估计值。然后按照最大似然概率法则来估计新的 $P_1$ 和 $P_2$ 。

$\begin{aligned} & P_1 = \frac{2 + 1 + 2}{15} = 0.33\\ & P_2 = \frac{3 + 3}{10} = 0.6 \end{aligned}$

设想我们是全知的神，知道每轮抛掷时的硬币就是如本文第001部分标示的那样。那么， $P_1$ 和 $P_2$ 的最大似然估计就是 0.4 和 0.5（下文中将这两个值称为 $P_1$ 和 $P_2$ 的真实值）。那么对比下我们初始化的 $P_1$ 和 $P_2$ 和新估计出的 $P_1$ 和 $P_2$ ：

初始化的 $P_1$	估计出的 $P_1$	真实的 $P_1$
0.2	0.33	0.4
初始化的 $P_2$	估计出的 $P_2$	真实的 $P_2$
0.7	0.6	0.5

通过表中数据我们知道，估计的 $P_1$ 和 $P_2$ 相比于它们的初始值，更接近它们的真实值了！

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

但这里有两个问题：

新估计出的 $P_1$ 和 $P_2$ 一定会更接近真实的 $P_1$ 和 $P_2$ 吗？
迭代一定会收敛到真实的 $P_1$ 和 $P_2$ 吗？

对于问题一：没错，一定会更接近真实的 $P_1$ 和 $P_2$ （数学可以证明，但这超出了本文的主题，请参阅其他书籍或文章）。
对于问题二：不一定，取决于 $P_1$ 和 $P_2$ 的初始化值，上面我们之所以能收敛到 $P_1$ 和 $P_2$ ，是因为我们幸运地找到了好的初始化值。

只会接近，但不一定能达到！

3.2.2 EM 算法进阶版

我们思考下，上面的方法还有没有改进的余地？

我们是用最大似然概率法则估计出的 $z$ 值，然后再用 $z$ 值按照最大似然概率法则估计新的 $P_1$ 和 $P_2$ 。也就是说，我们使用了一个最可能的 $z$ 值，而不是所有可能的 $z$ 值。

如果考虑所有可能的 $z$ 值，对每一个 $z$ 值都估计出一个新的 $P_1$ 和 $P_2$ ，将每一个 $z$ 值概率大小作为权重，将所有新的 $P_1$ 和 $P_2$ 分别加权相加，这样的 $P_1$ 和 $P_2$ 应该会更好一些。

那么，所有的 $z$ 值有多少个呢？

显然，有 $2^5 = 32$ 种，那么需要我们进行 32 次估值？

答案是并不需要，我们可以用期望来简化运算。

轮数	若是硬币 1	若是硬币 2
1（3 正 2 反）	0.00512	0.03087
2（2 正 3 反）	0.02048	0.01323
3（1 正 4 反）	0.08192	0.00567
4（3 正 2 反）	0.00512	0.03087
5（2 正 3 反）	0.02048	0.01323

在 EM 算法初始版中，我们直接根据概率大小确定使用了哪枚硬币，如 0.00512 < 0.03087，我们认为使用了硬币 2。显然这样过于草率了，也不太科学。因此我们可以使用一个更加科学的方式，比如我们确定使用哪枚硬币的概率。

利用上面这个表，我们可以算出每轮抛掷中使用硬币 1 或者使用硬币 2 的概率。

比如第 1 轮，使用硬币 1 的概率是：

$\begin{aligned} & z_{使用硬币1} = \frac{0.00512}{0.00512 + 0.03087} = 0.14\\ & z_{使用硬币2} = 1 - 0.14 = 0.86 \end{aligned}$

因此我们可以算出其他 4 轮的概率，如下表所示：

轮数	$z_i=硬币 1$	$z_i=硬币 2$
1（3 正 2 反）	0.14	0.86
2（2 正 3 反）	0.61	0.39
3（1 正 4 反）	0.94	0.06
4（3 正 2 反）	0.14	0.86
5（2 正 3 反）	0.61	0.39

上表中的右两列表示期望值。看第一行，0.86 表示从期望的角度看，这轮抛掷使用硬币 2 的概率是 0.86。相比于前面的方法，我们按照最大似然概率，直接将第 1 轮估计为用的硬币 2。此时的我们更加谨慎，我们只说，有 0.14 的概率是硬币 1，有 0.86 的概率是硬币 2，不再是非此即彼。这样我们在估计 $P_1$ 或者 $P_2$ 时，就可以用上全部的数据，而不是部分的数据，显然这样会更好一些。

这一步，我们实际上是估计出了 $z$ 的概率分布，这步被称作 E 步。

结合下表：

硬币	结果	统计
Unknown（未知）	正正反正反	3 正 2 反
Unknown（未知）	反反正正反	2 正 3 反
Unknown（未知）	正反反反反	1 正 4 反
Unknown（未知）	正反反正正	3 正 2 反
Unknown（未知）	反正正反反	2 正 3 反

我们按照期望最大似然概率的法则来估计新的 $P_1$ 和 $P_2$ ：

以 $P_1$ 估计为例，第 1 轮的 3 正 2 反相当于：

$0.14 \times 3 = 0.42 \Rightarrow 正\\ 0.14 \times 2 = 0.28 \Rightarrow 反$

依次算出其他四轮，列表如下：

轮数	正面	反面
1（3 正 2 反）	$0.14 \times 3 = 0.42$	$0.14 \times 2 = 0.28$
2（2 正 3 反）	$0.61 \times 2 = 1.22$	$0.61 \times 3 = 1.83$
3（1 正 4 反）	$0.94 \times 1 = 0.94$	$0.94 \times 4 = 3.76$
4（3 正 2 反）	$0.14 \times 3 = 0.42$	$0.14 \times 2 = 0.28$
5（2 正 3 反）	$0.61 \times 2 = 1.22$	$0.61 \times 3 = 1.83$
总计	4.22	7.98