20201008 -

0. 引言

提示：本篇文章并没有详细的说明EM算法数学推导，虽然前面通过GMM的例子能够明白大致的思想，但是在底层数学推导部分没有非常完整说明（后续有时间可能会继续添加），如果想知道数学原理的读者，就不要浪费时间再看这篇文章了，可以参考CS229的一个笔记；如果是想得到一个大致的了解，看完GMM的推导就可以了。

在学习HMM的训练算法中，其使用了EM算法，而不理解EM算法的话，实际上连很多参数都不明白，就更不用提理解了，因此必须首先学习一下EM算法。我最初学习机器学习的时候，当时很多算法都是理解了，然后能够用python的库使用了，没有具体深入，能明白大概是什么意思，这也是大多数实战书上的流程。这次就趁这个这个机会，把这部分内容来学习一下。

在学习EM算法的过程中，实际上我参考了非常多的资料，因为书上的讲解例子都是从数学的角度出发，都太过抽象，所以本次学习过程中，都是参考了一些国外的博客上的文章。

但是即使是这样，很多文章虽然表面上看起来是一样的，但是从推导上，很多地方都不一样。这也造成了，看了很多文章，反而更糊涂的结果。本篇文章讲从表层来入手，从最开始的一些内容，然后再到实际的数学公式来介绍。

（本人非数学的专业，本文仅记录自己在参考各个文章时候的思考。）

1. 初步的感受

EM算法的流程非常像k均值算法，都是先初步随机设置一个值，然后循环进行更新，EM算法只有两个步骤，分别是E和M过程，但实际上底层涉及的东西很多。

本小节的内容主要参考自文章[1]。

1.1 真实世界的案例

假设统计了一个学校的男女身高，男女的不同分布，通过随机采样两种性别的N个人，然后得到了相应的统计值，那么只需要进行估计相应的参数即可，这部分就是统计学中一些参数估计的方法。在采集的过程中，需要身高的数值，同时还需要性别。从一维数轴上来看，就是不同颜色（不同性别）的点，假设他们都符合高斯分布，然后对分布的参数进行估计。

假设我们并不知道男女的性别，仅仅只有身高的信息。那么应该如何进行估计，如何进行分析，知道两个分布参数，同时估计某个人是什么性别呢？在学习统计学的时候，一般进行参数估计都是针对单个高斯变量，通过最大似然估计就能完成。这里，性别分布也不知道，高斯分布的参数也不知道。这本质上就是一个GMM（高斯混合模型）的问题。

把问题简化一下，假设我们知道某些信息，我们该怎么解决这个参数估计或者概率分布的问题。 1）假设知道男女的标签，那么只需要通过最大似然估计的方式，即可。 $$\begin{gathered} {\hat{\mu }}_{MLE}=\frac{1}{N}\sum _i^N{x}_i=\bar{x} \\ {\hat{\theta }}_{MLE}=\sqrt{\frac{1}{N}\sum _i^N(x_i-\hat{x}{)}^2} \end{gathered}$$ 2）假设知道了分布的参数，想求解每个观察值是属于什么类别，只需要给每个观察值赋予一个混合的概率值（表示其属于某个概率分布的情况），然后同时求解最大似然方程即可。 $$\begin{gathered} P(x_i=male)=\frac{{\alpha }_1\mathcal{N}(x_i|{\mu }_1,{\Sigma }_1)}{Z} \\ P(x_i=female)=\frac{{\alpha }_2\mathcal{N}(x_i|{\mu }_2,{\Sigma }_2)}{Z} \\ Z=\sum _{j=1}^2{\alpha }_j\mathcal{N}(x_i|{\mu }_j,{\Sigma }_j) \end{gathered}$$ （上述公式摘抄自文章[2]，后续会对其进行具体分析） 从上述的内容中可以看出，假设我们知道一部分信息，我们可以推导出另外一部分信息。这就是他们很多在介绍EM算法的时候，提到的鸡和蛋的问题。在无监督学习中，通过GMM进行聚类的时候，并不知道这些信息。理解EM算法最直观的例子就是GMM的训练过程了。 问题就是，我们仅仅知道测量值的时候，我们如何进行GMM的参数估计还有相应的标签赋值呢？

1.2 宏观的理解EM算法

既然没有办法直接的进行解决，同样是分为两个步骤，在初始化过程中，将一些参数随机，例如直接给出正态分布的参数，或者给出每个观察值的标签。下面假设我们给出了每个观察值的标签。 1）通过观察值的标签，我们可以直接对这个参数分布进行估计，然后得到了相应的概率分布 2）前一步中得到的概率分布参数，重新对观察值进行标签赋值 3）重复以上两个步骤，直到收敛

这个步骤跟k-means算法的步骤非常像。但是实际上，虽然看起来很简单，但是需要思考的内容有很多，就比如说你怎么保证这个东西就一定收敛？或者说，这一定会成功呢？！

1.3 小节

前面的示例中，是从一个实际的问题出发；实际上EM算法解决最大似然估计的方法是引入了一个潜在变量，通过这个潜在变量将本来不好优化的问题转化为可以优化的部分。（但是前面的例子中没有直白的说明这个潜在变量） 前面宏观的理解非常简单，与k-means非常像。很多文章在利用GMM进行介绍的时候，也非常容易明白，但是一旦涉及到具体的公式，问题就变得抽象了。

2. EM算法在GMM模型中的应用

本小节将主要利用EM算法来解决GMM模型的问题，通过具体的数学公式，明白实际上操作流程。本小节的内容参考自文章[3]。 上述图片也是来自文章[3]，好像是一个通过matlab的程序生成的，后续的时候可以尝试也绘制一个这种图片。文章[3]要解决的问题，就是估计GMM的参数，在估计过程中，采用软聚类的形式，也就是说，给每个点赋予一个属于某个类的概率，而不是直接属于某个类。具体问题是，存在数据集$n$个点，其属于$R^d$空间，他们由一个GMM（Gaussian Mixture Model)模型生成，现在的任务就是估计出最可能产生这些点的GMM模型，也就是求解他们的参数。

2.1 参数设置

本小节主要介绍一下后文中主要使用的参数都有哪些 1）k个成分，GMM模型中，有k个高斯分布，一般这个数值都是实现指定的。 2）每个高斯分布$j，k\in \{1,2,...,k\}$有参数$\mu$和${\sigma }^2$去确定 $$P_x(x|{\mu }^j,{\sigma }_j^2)=\mathcal{N}(X;{\mu }^j,{\sigma }_j^2)=\frac{1}{(2\pi {\sigma }_j^2{)}^{\frac{d}{2}}}{e}^{\frac{-\Vert x-{\mu }^j{\Vert }^2}{2{\sigma }_j^2}}$$ 上面这些变量的含义，很好理解，但是这里有些不一样的是，他没有使用协方差矩阵，在有些公式里面可能是协方差矩阵，这里不影响理解。 3）每个高斯分布有一个权重，这个权重代表这个可能是这个高斯成分的概率 $$j\sim Multinonial(p_1,p_2,...,p_K)$$ 4）最后是所有的要去估计的参数集合 $$\Theta =p_1,p_2,...,p_K,{\mu }^1,{\mu }^2,...,{\mu }^K,{\sigma }_1^1,{\sigma }_2^2,...,{\sigma }_K^2$$

2.2 最大似然估计

最大似然估计是指选取能够满足所有观测值最大出现概率的参数值。 $$\text{P}(X_n|\Theta )=\prod _{i=1}^n\sum _{j=1}^K{p}_j\mathcal{N}(x_i;{\mu }^j,{\sigma }_j^2)$$ 而求解这个方程的过程会非常复杂，正如前面所说，很多内容都不知道。

2.3 EM过程求解

首先初始化变量$\Theta$。

2.3.1 E步骤（Expectation）

当所有参数都被确定了，那么就需要求解每个数据的后验概率（每个数据点属于每个高斯分布的概率），换句话说，就是已经看到了这些点，那么每个成分的权重是多少。 $$\text{P}(j|i)=\text{P}(y=j|X_n)=\frac{\text{P}(j)\text{P}(X_n|j)}{\text{P}(X_n)}=\frac{\text{P}(j)\text{P}(X_n|j)}{{\sum }_{j=1}^K\text{P}(j)\text{P}(X_n|j)}$$ 先验概率是$\text{P}(j)$，似然值是$\text{P}(X_n|j)$。在另一篇文章中《先验概率及后验概率等解释》记录了这些名词的解释。 这里的$\text{P}(j)$就是前面的每个成分的权重，从历时解释的角度来说，{\rm \textbf P}(j|i)就是要求，在看到这个样本点之后，其属于成分$j$的概率是多少，实质上就是要确定每个样本点属于每个高斯成分的概率，正是一个打标签的过程。类别k-means方法，就是随机中心样本点之后，求每个点的标签。

2.3.2 M步骤（Maximization）

前面的步骤已经得到了每个点属于每个高斯成分的后验概率，这部分将重新估计初始化的数值。首先重新估计$\text{P}(j)$。 $$\hat{\text{P}}(j)=\frac{{\hat{n}}_j}{n}=\frac{{\sum }_{i=1}^n\text{P}(j|i)}{n}$$ 然后，利用最大似然来估计其他高斯分布的参数。 $$\text{P}(X_n|\Theta )=\prod _{i=1}^n\sum _{j=1}^K{p}_j\mathcal{N}(x_i;{\mu }^j,{\sigma }_j^2)$$ 将其改写为$log$似然函数。 $$\begin{gathered} \log (\text{P}(X_n|\Theta ))=\ell (X_n;\Theta ) \\ \ell (X_n;\Theta )=\sum _{i=1}^{n}log\sum _{j=1}^K\hat{\text{P}}(j)\mathcal{N}(x_i;{\mu }^j,{\sigma }_j^2) \end{gathered}$$ 其中$p_j$已经被替换为前面估计出来的$\hat{\text{P}}(j)$。最后是各个参数的估计方式。 $$\begin{gathered} {\hat{\mu }}^j=\frac{1}{{\hat{n}}_j}\sum _{i=1}^n\text{P}(j|i)x_i \\ {\hat{\sigma }}_j^2=\frac{1}{{\hat{n}}_j}\sum _{i=1}^n\text{P}(j|i)\Vert x_i-{\hat{\mu }}^j{\Vert }^2 \end{gathered}$$ （正如前面介绍参数的时候所说，这里没有利用协方差矩阵的形式，仅仅是一个例子，能够明白其推到的过程）

2.3.3 重复上述步骤

通过重复上述两个步骤，直到最后数值收敛即可。

2.4 小节

本部分的内容，主要就是第一小节的公式实例推导，但前面也说过，EM算法的通过引入了一个潜在的变量来实现参数估计，这里的情况比较特殊，没有具体说明，在后面具体引入通用公式的时候再来具体说明。英语好的读者可以直接阅读文章[3]来获得最佳的体验。 利用GMM的例子来进行解释，能够知道大致的步骤，但是对于实际的问题来说，例如HMM训练算法的推导，还是不一样，因此需要一个更通用的公式来解决这个鸡和蛋的问题。

3. EM推导过程

注意：如果上述内容能够明白，就已经大致知道了EM算法的思想，由于本人能力有限，下面的推导也是看了其他的一些文章，虽然能明白他要干什么，但是对数学公式不甚理解，为了不误人子弟，请读者谨慎阅读。

其实在看这部分内容的时候，耗费了很久，很多文章的角度不同，深度也不同，我也很难明白。但是我觉得，学习这个东西的目标应该循序渐进，既然以后要用到这个东西，那么可以先从怎么使用来入手。并不是说，在sklearn中直接使用了GMM这种形式（EM被隐式的调用），而是我能够把自己的问题形式化为直接使用EM的形式，这是第一个目标，而第二个目标，就是能够理解底层的数学原理。先简单说一下EM要优化的目标。

3.1 优化目标

在文章[2]中，文章的前半部分也是利用GMM来进行说明，来看一看他的公式。 $$\begin{gathered} \Theta =\{{\alpha }_1,{\alpha }_2,..,{\alpha }_K,{\mu }_1,{\mu }_2,...,{\mu }_K,{\Sigma }_1,{\Sigma }_2,...,{\Sigma }_K\} \\ \underset{\Theta }{\max }P(X|\Theta )=\underset{\Theta }{\max }\prod _{i}P(x_i|\Theta )=\underset{{\alpha }_j,{\mu }_j{\Sigma }_j}{\max }\prod _i⟮\sum _{j=1}^K{\alpha }_j\mathcal{N}(x_i|{\mu }_j,{\Sigma }_j)⟯ \end{gathered}$$ 转化为$log$似然函数 $$\begin{gathered} \underset{\Theta }{\max }logP(X|\Theta )=\underset{\Theta }{\max }log(\prod _{i}P(x_i|\Theta ))=\underset{\Theta }{\max }\sum _{i}log(P(x_i|\Theta )) \\ =\underset{{\alpha }_i,{\mu }_j,{\Sigma }_j}{\max }\sum _{i}log(\sum _{j=1}^K{\alpha }_j\mathcal{N}(x_i|{\mu }_j,{\Sigma }_j)) \end{gathered}$$ 这个公式跟前面讲解GMM的时候的公式是一致的。虽然从前面的讲解中，我们知道了GMM的参数如何用EM算法进行估计，但是EM算法本身的描述更为形式化，或者说更通用，他可以应用于很多优化问题，本质上的想法就是要去优化一个似然函数，或者说$log-$似然函数，如 $$\underset{\Theta }{\max }log(X|\Theta )=\underset{\Theta }{\max }log(\prod _{i}P(x_i|\Theta ))=\underset{\Theta }{\max }\sum _{i}log(P(x_i|\Theta ))$$ 那么也就是说，将问题更改为，如果这个概率的分布并不是混合高斯分布的情况下，甚至于更通用的场景下，应该怎么来解决这个优化问题。这里文章[2]介绍了一些数学基础才开始将EM算法，我这里直接就开始将EM，然后穿插着再把这些数学基础带上，这里最重要的就是这个Jensen’s inequality。

3.2 隐含变量

隐含变量是一个非常重要的概念，但是经过这些文中的描述之后，我有点不明白这其中的问题了。有些文章直接从隐含变量的角度来说明，其要最大似然估计的公式中就带有着隐含变量，因为求导过程复杂，所以采用了EM算法；而有些文章则是说是为了求解这个最大似然函数，然后构造了一个隐含变量，以此帮助进行优化。感觉很多内容都不统一，这也导致我看了越来越多的文章之后，越来越感觉非常疑惑。但是有些文章也指出，**EM算法最初的提出，就是为了解决含有隐含变量的最大似然函数的求解 **，这句话我感觉也挺正确的。总之，这里要说明的问题就是，因为这些概念的不明确，导致我从这些博客从这些文章中得到的思想是非常不一样的，也就导致了不能理解。看来得多看一些比较权威的文章或者是书籍才行。

通俗一点来讲，在前面的GMM中，一个指示样本$x$属于某个类别的变量，就被称为隐含变量，用公式写出来，就是： $$log(X|\theta )=\sum _{i}log(\sum _{z_i}p(x_i,z_i|\theta ))$$ 也就是说， 通过遍历隐含变量的空间，计算整个边缘概率，就是原始的概率结果。

3.3 EM的形式化推导

首先，EM算法就是要求解出来一个下限值，然后通过迭代的方式，来持续获取到最大的下限。要使用jensen不等式的内容来辅助推导。这部分可以看文章[2]的一个动图来加深理解。 意思就是说，找到一个函数，先固定参数，然后调整函数来让他达到最大的下限；而第二步之后，固定函数，调整参数再次达到最大下限。重复这个步骤，最终得到相应的结果。

接着前面的内容继续来说，通过加入了隐含变量，或者说我本来就有隐含变量，那么似然值应该是如下过程： $$P(x_i|\theta )=\sum _{j=1}^{K}P(x_i|t_i=j,\theta )P(t_i=j|\theta )=\sum _{j=1}^{K}P(x_i,t_i=j|\theta )$$ 其中，$t_i$是一个隐含变量，其解释了$x_i$属于哪个聚类的类别。那么根据jensen不等式的内容， $$log(\sum _{j=1}^K{\beta }_j{\nu }_j)\ge \sum _{j=1}^{K}log({\beta }_j{\nu }_j)$$ 其中，各项内容分别是 $$\begin{gathered} {\beta }_j=P(t_i=j|\theta ) \\ {\nu }_j=P(x_i|t_i=j,\theta ) \end{gathered}$$ 前面已经定义了，本来要求解得最优化目标是 $$\underset{\Theta }{\max }log(X|\Theta )=\underset{\Theta }{\max }log(\prod _{i}P(x_i|\Theta ))=\underset{\Theta }{\max }\sum _{i}log(P(x_i|\Theta ))$$ 那么讲上述内容带入进去， $$\begin{gathered} \sum _{i}log(P(x_i|\theta ))=\sum _{i}log(\sum _{j=1}^{K}P(x_i|t_i=j,\theta )P(t_i=j|\theta ))\ge \\ \sum _i\sum _{j=1}^{K}log({\beta }_j{\nu }_j)=\sum _i\sum _{j=1}^{K}log(P(t_i=j|\theta )P(x_i|t_i=j,\theta )) \end{gathered}$$ 但是公式右边是一个$\theta$的函数，还要有一个分布$q$ $$\begin{gathered} log(P(x_i|\theta ))=log(\sum _{j=1}^{K}P(x_i|t_i=j,\theta )P(t_i=j|\theta )) \\ =log(\sum _{j=1}^{K}P(x_i,t_i=j|\theta )) \\ =log(\sum _{j=1}^K\frac{q(t_i=j)}{q(t_i=j)}P(x_i,t_i=j|\theta )) \\ =log(\sum _{j=1}^{K}q(t_i=j)\frac{P(x_i,t_i=j|\theta }{q(t_i=j)}) \end{gathered}$$ 再次利用jensen不等式。 $$log(\sum _{j=1}^{K}q(t_i=j)\frac{P(x_i,t_i=j|\theta )}{q(t_i=j)})\ge \sum _{j=1}^{K}log(q(t_i=j)\frac{P(x_i,t_i=j|\theta )}{q(t_i=j)})=L(\theta ,q)$$ 这部分内容定义了一个下限，保证所有的值都小于等于最大似然值。 那么下面的步骤就是最大化这个下限。

1）Expection Step

$$q^{k+1}=\arg \max L({\theta }^k,q^k)$$ 固定参数$\theta$，求解最大的下限值。而上述公式可以推导为 $$q^{k+1}(t_i)=\arg \max L({\theta }^k,q^k)=P(t_i|x_i,\theta )$$ 这个公式就是前面GMM时候求解的$P(j|x)$

2）Maximization Step

$${\theta }^{k+1}=\arg \max L({\theta }^k,q^{k+1}）$$ 固定$q$然后最大的这个下限，这里文章[2]的公式就不对了，感觉他少推到了一些东西，最后就得到了一个期望值。

3.4 小节

这部分还是没有从公式的角度，或者从数学的角度上弄明白，就感觉不知道他到底想说什么，也没有一个承上启下的过程。

再来记录一下： 因为直接对含有因变量的似然函数进行求解很难，所以采用一种另外的方法，也就是通过求解他的紧下界，也就是必须存在等于的情况；那么此时通过jensen不等式的内容，可以直接求解。而具体的实现步骤，就是通过固定某些参数，然后最大化另外的，持续这个步骤，直到收敛。

jensen不等式的内容： $$\mathrm{E}[f(x)]\ge f(\mathrm{E}X)$$ 其中$f$是一个凸函数（$f^{{}^{\prime \prime }}(x)\ge 0$），$X$是一个随机变量，也就是他的期望值在某些函数下，不等式成立。

我觉得，我大致上明白了具体的含义，通过参考了两个部分的内容：

EMEM详细推导

同时我觉得文章[1]在后面写的挺好的，就是必须有一定的基础；同时她在最后对比了这种方式与梯度下降算法的异同，贴一个图。

同时文章The Expectation-Maximization Algorithm的解释，看起来更为形式化，可以作为参考。

参考

[1]The EM Algorithm Explained [2]Gaussian Mixture Models and Expectation-Maximization (A full explanation) [3]Expectation-Maximization Algorithm Step-by-Step