大家好，又见面了，我是你们的朋友全栈君。

1、总述

期望最大（Expectation Maximization）算法是一种从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。EM算法是机器学习十大算法之一，或许确实是因它在实际中的效果很好吧。

演示

2、定义

EM算法，全称Expectation Maximization Algorithm，译作最大期望化算法或期望最大算法，它是一种迭代算法，用于含有隐变量（hidden variable）的概率参数模型的最大似然估计或极大后验概率估计。

3、感性例子：

例子简介：

假设现在有两枚硬币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$，（硬币1实验15次正面6次，硬币2实验10次正面5次）如下：

$$\begin{array}{rl}{P}_1& =(3+1+2)/15=0.4\\ P_2& =(2+3)/10=0.5\end{array}$$
到这里，一切似乎很美好，下面我们加大难度。

加入隐变量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初级版。

EM初级版

我们不妨这样，先随便给$P_1$和$P_2$赋一个值，比如：
$$\begin{gathered} P_1=0.2 \\ {P}_2=0.7 \end{gathered}$$
然后，我们看看第一轮抛掷最可能是哪个硬币。
如果是硬币1，得出正正反正反的概率为（注意这里是有先后顺序的）
$$0.2\ast 0.2\ast 0.2\ast 0.8\ast 0.8=0.00512$$

如果是硬币2，得出正正反正反的概率为
$$0.7\ast 0.7\ast 0.7\ast 0.3\ast 0.3=0.03087$$
然后依次求出其他4轮中的相应概率。做成表格如下：

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

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

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

$$\begin{array}{rl}{P}_1& =(2+1+2)/15=0.33\\ P_2& =(3+3)/10=0.6\end{array}$$
设想我们是全知的神，知道每轮抛掷时的硬币就是如前面例子中的那样，$P_1$和$P_2$的最大似然估计就是0.4和0.5（下文中将这两个值称为$P_1$和$P_2$的真实值）。那么对比下我们初始化的$P_1$和$P_2$和新估计出的$P_1$和$P_2$：

初始化的P1	估计出的P1	真实的P1	初始化的P2	估计出的P2	真实的P2
0.2	0.33	0.4	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$的最大似然估计。

这里有两个问题：
1、新估计出的$P_1$和$P_2$一定会更接近真实的$P_1$和$P_2$？
答案是：没错，一定会更接近真实的$P_1$和$P_2$，数学可以证明，但这超出了本文的主题，请参阅其他书籍或文章。

2、迭代一定会收敛到真实的$P_1$和$P_2$吗？
答案是：不一定，取决于$P_1$和$P_2$的初始化值，上面我们之所以能收敛到$P_1$和$P_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	0.00512	0.03087
2	0.02048	0.01323
3	0.08192	0.00567
4	0.00512	0.03087
5	0.02048	0.01323

利用上面这个表，我们可以算出每轮抛掷中使用硬币1或者使用硬币2的概率。比如第1轮，使用硬币1的概率是：
$0.00512/(0.00512+0.03087)=0.14$
使用硬币2的概率是$1-0.14=0.86$
依次可以算出其他4轮的概率，如下：

轮数	$z_i$=硬币1	$z_i$=硬币2
1	0.14	0.86
2	0.61	0.39
3	0.94	0.06
4	0.14	0.86
5	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的概率下，去计算出现现在这种抛掷情况的可能性，就是说计算一下期望看看有多少次等价的硬币1抛出来为正，多少次抛出来为负，以第一轮为例，硬币1在第一轮中被使用的该概率是0.14，现在出现了三次正面两次反面。

第1轮的3正2反相当于
$0.14\ast 3=0.42$正
$0.14\ast 2=0.28$反
依次算出其他四轮，列表如下：

轮数	正面	反面
1	0.42	0.28
2	1.22	1.83
3	0.94	3.76
4	0.42	0.28
5	1.22	1.83
总计	4.22	7.98

$P_1=4.22/(4.22+7.98)=0.35$

可以看到，改变了$z$值的估计方法后，新估计出的$P_1$要更加接近0.4。原因就是我们使用了所有抛掷的数据，而不是之前只使用了部分的数据。

这步中，我们根据E步中求出的$z$的概率分布，依据最大似然概率法则去估计$P_1$和$P_2$，被称作M步。

例子总结

以上，我们用一个实际的小例子，来实际演示了EM算法背后的idea，共性存于个性之中，通过这个例子，我们可以对EM算法究竟在干什么有一个深刻感性的认识，掌握EM算法的思想精髓。

4、Jensen不等式

在完善EM算法之前，首先来了解下Jensen不等式，因为在EM算法的推导过程中会用到。

Jensen不等式在优化理论中大量用到，首先来回顾下凸函数和凹函数的定义。假设$f$是定义域为实数的函数，如果对于所有的$x$，$f(x)$的二阶导数大于等于0，那么$f$是凸函数。当$x$是向量时，如果hessian矩阵$H$是半正定（即$H>=0$），那么$f$是凸函数。

如果，$f(x)$的二阶导数小于0或者$H>0$，那么f就是凹函数。

Jensen不等式描述如下：

• 如果$f$是凸函数，$X$是随机变量，则$E[f(X)]>=f(E[X])$，当$f$是严格凸函数时，则$E[f(X)]>f(E[X])$；
• 如果$f$是凹函数，$X$是随机变量，则$f(E[X])<=E[f(X)]$，当$f$是（严格）凹函数当且仅当$-f$是（严格）凸函数。

通过下面这张图，可以加深印象：

上图中，函数$f$是凸函数，$X$是随机变量，有0.5的概率是$a$，有0.5的概率是$b$。$X$的期望值$E(x)$就是$a$和$b$的中值$\frac{a+b}{2}$了，图中可以看到$E[f(X)]>=f(E[X])$成立。

5、EM思想

EM算法推导过程中，会使用到极大似然估计法估计参数，所以，首先给出一个求最大似然函数估计值的一般步骤：

极大似然估计的介绍可以参考之前的一篇blog

（1）写出似然函数；
（2）对似然函数取对数，并整理；
（3）求导数，令导数为0，得到似然方程；
（4）解似然方程，得到的参数即为所求；
关于极大似然估计的实例，这里就不再提及了，下面介绍EM算法。

给定$m$个训练样本 { x 1 , x 2 , … , x m } \{x_1,x_2,\ldots,x_m\} {
x1​,x2​,…,xm​}，假设样本间相互独立，我们想要拟合模型$p(x,z)$到数据的参数，$z_j^{(i)}$为在$x_i$的情况下不同$z$的可能 { z 1 ( i ) , z 2 ( i ) , … , z j ( i ) } \{z_1^{(i)},z_2^{(i)},\ldots,z_j^{(i)}\} {
z1(i)​,z2(i)​,…,zj(i)​}。根据分布，我们可以得到如下这个似然函数：

$\begin{array}{rl}\ell (\theta )& =\sum _{i=1}^m\log p(x_i;\theta )\\ & =\sum _{i=1}^m\log \sum _{j}p(x_i,z_j^{(i)};\theta )\end{array}$

注：之后的公式中$z_j^{(i)}$不在带上标${}^{(i)}$的原因为：之后所有的$x_i$发生都将所有的$z_j$进行求和，因为任意$x_i$的情况下$z_j$总数是不变的，即使有变化，某种$z_j$在$x_i$情况下不可能发生那即认为$p(x_i,z_j^{(i)};\theta )=0$即可。

公式解析：1、联合概率分布不同样本之间相互独立，所以可以写成概率连乘，即为发生如下样本的概率 2、在例子中所叙述的进阶版思想合理利用所有z的可能不同的z不同的可能性相加最后的期望才是最后的似然可能性。

第一步是对极大似然函数取对数第二步是对每个样本实例的每个可能的类别z求联合分布概率之和。然而，直接求这个参数$\theta$会比较困难，因为上式存在一个隐含随机变量（latent random variable）$z$。如果$z$是个已知的数，那么使用极大似然估计来估算会很容易。在这种$z$不确定的情形下，EM算法就派上用场了。

EM算法是常用的估计参数隐变量的利器。对于上述情况，由于存在隐含变量，不能直接最大化$l(\theta )$，所以只能不断地建立$l$的下界（E-step），再优化下界（M-step），依次迭代，直至算法收敛到局部最优解。这就是EM算法的核心思想，简单的归纳一下：

EM算法通过引入隐含变量,使用MLE（极大似然估计）进行迭代求解参数。通常引入隐含变量后会有两个参数，EM算法首先会固定其中的第一个参数，然后使用MLE计算第二个变量值；接着通过固定第二个变量，再使用MLE估测第一个变量值，依次迭代，直至收敛到局部最优解。

E-Step和M-Step。

• E-Step：通过observed data和现有模型估计参数估计值 missing data；
• M-Step：假设missing data已知的情况下，最大化似然函数。

是否会收敛：由于算法保证了每次迭代之后，似然函数都会增加，所以函数最终会收敛(最后有推到)。

6、EM推导

下面来推导EM算法：

对于每个实例$x_i$，用$Q_i$表示样本实例隐含变量$z$的某种分布，且$Q_i$满足条件${\sum }_j{Q}_i(z_j)=1,Q_i(z_j)>=0,\forall j$，如果$Q_i$是连续性的，则$Q_i$表示概率密度函数，需要将求和符号换成积分符号。

对于上面的式子，做如下变换：

$$\begin{array}{rl}\sum _i\log p(x_i;\theta )& =\sum _i\log \sum _{j}p(x_i,z_j;\theta )(1)\\ & =\sum _i\log \sum _j{Q}_i(z_j)\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}(2)\\ & \ge \sum _i\sum _j{Q}_i(z_j)\log \frac{p(x_i,z_j;\theta )}{Q_i(z_j)}(3)\end{array}$$

上面三个式子中，式（1）是根据联合概率密度下某个变量的边缘密度函数求解的（这里把$z$当作是随机变量）。对每一个样本i的所有可能类别$z$求等式右边的联合概率密度函数和，也就得到等式左边为随机变量$x$的边缘概率密度。由于对式（1）直接求导非常困难，所以将其分子分母都乘以一个相等的函数$Q_i(z_j)$，转换为式（2）。而在式（2）变为式（3）的过程，采用的是上面提到的Jensen不等式。分析过程如下：

首先，把（1）式中的log函数体看成是一个整体，由于log(x)的二阶导数为$-\frac{1}{x^2}$，小于0，为凹函数。所以使用Jensen不等式时，应用第二条准则：$f(E[X])>=E[f(x)]$。

到这里，问题简化为如何求解随机变量的期望。还记得当年读大学的时候，概率论中的随机变量的期望计算方法么，不记得也没关系，下面这张图比较详细：

因此，结合上面的知识点，我们可以把(2)式当中的$Q_i(z_j)$看成相应的概率$p_{z_j}$，把$\frac{p(x^i,z_j;\theta )}{Q_i(z_j)}$看作是$z_j$的函数$g(z_j)$，类似地，根据期望公式$E(x)=\sum x\ast p(x)$可以得到：
$$E(z)=\sum _{j}p(z_j)\ast g(z_j)=\sum _j{Q}_i(z_j)[\frac{p(x^i,z_j;\theta )}{Q_i(z_j)}]$$
其实就是$\frac{p(x^i,z_j;\theta )}{Q_i(z_j)}$的期望，再根据凹函数对应的Jensen不等式性质：
$$f(E_{z_j\sim Q_i}[\frac{p(x^i,z_j;\theta )}{Q_i(z_j)}])\ge E_{z_j\sim Q_i}[f(\frac{p(x^i,z_j;\theta )}{Q_i(z_j)})]$$
因此便得到了公式（3）。OK，现在我们知道上面的式（2）和式（3）两个不等式可以写成：似然函数$L(\theta )>=J(z,Q)$的形式（$z$为隐含变量），那么我们可以通过不断的最大化$J$的下界，来使得$L(\theta )$不断提高，最终达到它的最大值。使用下图会比较形象：

这里来说说上图的内在含义：首先我们固定 $\theta$，调整$Q(z)$使下界$J(z,Q)$上升至与$L(\theta )$在此点$\theta$处相等（绿色曲线到蓝色曲线），然后固定$Q(z)$，调整$\theta$使下界$J(z,Q)$达到最大值（${\theta }_t$到${\theta }_{t+1}$），然后再固定$\theta$，调整Q(z)……直到收敛到似然函数$L(\theta )$的最大值处的$\theta$。

这里有两个问题：

1. 什么时候下界$J(z,Q)$与$L(\theta )$在此点$\theta$处相等？
2. 为什么一定会收敛？

首先来解释下第一个问题（1.什么时候下界$J(z,Q)$与$L(\theta )$在此点$\theta$处相等？）。在Jensen不等式中说到，当自变量$X=E(X)$时，即为常数的时候，等式成立。而在这里，为：
$$\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}=c$$
对该式做个变换，并对所有的$z$求和，得到
$$\sum _{j}p(x_i,z_j;\theta )=\sum _j{Q}_i(z_j)c$$
因为前面提到${\sum }_j{Q}_i(z_j)=1$（概率之和为1），所以可以推导出：
$$\sum _{j}p(x_i,z_j;\theta )=c$$
因此也就可以得到下面的式子：

$$\begin{array}{rl}{Q}_i(z_j)& =\frac{p(x_i,z_j;\theta )}{{\sum }_{j}p(x_i,z_j;\theta )}\\ & =\frac{p(x_i,z_j;\theta )}{p(x_i;\theta )}\\ & =p(z_j|x_i;\theta )\end{array}$$

至此，我们推出了在固定参数$\theta$后，使下界拉升的$Q(z)$的计算公式就是后验概率（条件概率），一并解决了$Q(z)$如何选择的问题。此步就是EM算法的E-step，目的是建立$L(\theta )$的下界。接下来的M-step，目的是在给定$Q(z)$后，调整$\theta$，从而极大化$L(\theta )$的下界$J$（在固定$Q(z)$后，下界还可以调整的更大）。到此，可以说是完美的展现了EM算法的E-step & M-step，完整的流程如下：

第一步，初始化分布参数$\theta$；
第二步，重复E-step 和 M-step直到收敛：

• E步骤：根据参数的初始值或上一次迭代的模型参数来计算出的隐性变量的后验概率（条件概率），其实就是隐性变量的期望值。作为隐藏变量的现有估计值：

$$Q_i(z_j)=p(z_j|x_i;\theta )$$

• M步骤：最大化似然函数从而获得新的参数值：

$$\theta =arg\underset{\theta }{\max }\sum _i\sum _j{Q}_i(z_j)\log \frac{p(x_i,z_j;\theta )}{Q_i(z_j)}$$

通过不断的迭代，然后就可以得到使似然函数$L(\theta )$最大化的参数$\theta$了。

上面多次说到直至函数收敛，那么该怎么确保EM收敛呢？(②为什么一定会收敛？)，下面进入证明阶段.

假定${\theta }^{(t)}$和${\theta }^{(t+1)}$是EM第t次和t+1次迭代后的结果。如果我们证明了$l({\theta }^{(t)})<=l({\theta }^{(t+1)})$，也就是说极大似然估计单调增加，那么最终我们就会得到极大似然估计的最大值。下面来证明，选定$\theta (t)$之后，我们得到E步：
$$Q_i^{(t)}(z_j)=p(z_j|x^i;{\theta }^{(t)})$$
这一步保证了在给定$\theta (t)$时，Jensen不等式中的等式成立，也就是
$$l({\theta }^{(t)})=\sum _i\sum _j{Q}_i^{(t)}(z_j)\log \frac{p(x_i,z_j;{\theta }^{(t)})}{Q_i(z_j)}$$
然后进行M步，固定$Q_i^{(t)}{z}_j$，并将${\theta }^{(t)}$试作变量，对上面的式子求导，得到${\theta }^{(t+1)}$，这样经过一些推导会有以下式子成立：
$$\begin{array}{rl}l({\theta }^{(t+1)})& \ge \sum _i\sum _j{Q}_i^{(t)}(z_j)\log \frac{p(x_i,z_j;{\theta }^{(t+1)})}{Q_i^{(t)}(z_j)}(4)\\ & \ge \sum _i\sum _j{Q}_i^{(t)}(z_j)\log \frac{p(x_i,z_j;{\theta }^{(t)})}{Q_i^{(t)}(z_j)}(5)\\ & =l({\theta }^{(t)})(6)\end{array}$$
在公式（4）中，得到${\theta }^{(t+1)}$，只是最大化$l({\theta }^{(t)})$，也就是$l({\theta }^{(t+1)})$的下界，并没有使等式成立，等式成立只有在固定$\theta$，并按E步得到$Q_i$时才能成立。

这样就证明了$l(\theta )$会单调增加。如果要判断收敛情况，可以这样来做：一种收敛方法是$l(\theta )$不再变化，还有一种就是变化幅度很小,即根据$l(\theta {)}^{(t+1)}-l(\theta {)}^{(t)}$的值来决定。

EM算法类似于坐标上生法（coordinate ascent）：E步：固定$\theta$，优化Q；M步：固定$Q$，优化$\theta$；交替将极值推向最大。

7、应用

• 混合高斯模型（Mixtures of Gaussians）
• K-means聚类算法

8、参考文献

如何感性地理解EM算法

EM-最大期望算法

http://cs229.stanford.edu/notes/cs229-notes8.pdf

wikipedia维基百科

JerryLead博客-（EM算法）The EM Algorithm

从最大似然到EM算法浅解

Rachel Zhang-EM算法原理