机器学习（7）：EM（期望最大化）算法

最近自己会把自己个人博客中的文章陆陆续续的复制到CSDN上来，欢迎大家关注我的 个人博客，以及我的github。

本文主要讲解有关EM算法的相关知识，除了给出算法的基本思想，还会给出该算法的由来及公式推导，这又涉及到极大似然估计和Jensen不等式的相关内容。

本文主要是依据李航老师的《统计学习方法》和邹博老师的机器学习教学视频总结编写的。文中所用到的有关机器学习的基本概念和方法可以参考本人博客中该系列之前的文章，或者直接上网搜索相关的内容。以下文章所列出的公式以及公式的推导读者们最好是在草稿本上自己推导一遍。由于本人水平所限，文章中难免有错误和不当之处，欢迎大家多多批评指正！

一、EM算法概述

EM算法又称期望最大化算法，它是一种迭代算法，可以用来求含有隐变量的概率模型的极大似然估计。所谓隐变量就是不能直接观测的变量。比如有三枚硬币，要根据第一枚硬币的正反来决定投掷第2、3枚硬币中的哪一枚。观测变量是第2、3枚硬币的正反，隐变量是第一枚硬币的正反。

如果概率模型的变量都是观测变量，则给定数据，可以直接用极大似然估计法来求得最优的模型参数。但是如果模型含有隐变量时，极大似然估计法就不再适用了，这时可以选择用EM算法来解决。

学习EM算法必须先了解两个知识点——极大似然估计和Jensen不等式。其中极大似然估计已经在文章极大似然估计和最小二乘法讲过了，下面就先来看一下Jensen不等式的相关内容。

二、Jensen不等式

Jensen不等式有个大前提：函数$f(x)$是凸函数。以下性质都是在该前提的基础上成立的。

1. 对于二维平面

Jensen不等式的形式为：
$$f(\theta x+(1-\theta )y)\le \theta f(x)+(1-\theta )f(y)$$
这也是凸函数的一种定义。

2. 对于n维平面

$$f({\theta }_1{x}_1+{\theta }_2{x}_2+...+{\theta }_n{x}_n)\le {\theta }_{1}f(x_1)+{\theta }_{2}f(x_2)+...+{\theta }_{n}f(x_n)$$

其中${\theta }_1+{\theta }_2+...+{\theta }_n=1$且${\theta }_i\ge 0$

3. 一般形式

在上式中有${\theta }_1+{\theta }_2+...+{\theta }_n=1$且${\theta }_i\ge 0$，所以可以把${\theta }_i$看作是随机变量$x_i$发生的概率，那么$f({\theta }_1{x}_1+{\theta }_2{x}_2+...+{\theta }_n{x}_n)=f(Ex)$，而${\theta }_{1}f(x_1)+{\theta }_{2}f(x_2)+...+{\theta }_{n}f(x_n)=Ef(x)$，所以有：
$$f(Ex)\le Ef(x)$$
即函数的期望大于期望的函数，当且仅当$x_i=c$，c为常数时取等号。

三、EM算法

在极大似然估计和最小二乘法一文中提到过，可以用负对数似然函数作为模型的损失函数，想要让损失函数（负对数似然函数）取最小值，则对数似然函数应该取最大值。

若样本满足概率分布$p(x_i;\theta )$，则令目标函数为对数似然，并求它的最大值：
$$l(\theta )=lnL(\theta )=\prod _{i=1}^n\log p(x_i;\theta )=\sum _{i=1}^n\log p(x_i;\theta )$$

$$=\sum _{i=1}^n\log \sum _{z_i}p(x_i,z_i;\theta )$$

其中$z_i$是隐变量，联合概率密度$p(x_i,z_i;\theta )$对隐变量积分就得到了边缘概率密度$p(x_i;\theta )$。

之前说了，由于含有隐变量，用极大似然估计估计对其求导的方法并不适用，所以采用了EM算法，迭代的更新参数，重复该过程直到收敛到局部最大值。

EM算法的基本思想如上图所示，$\theta$是参数，紫色的曲线是目标函数$p(x|\theta )$，则对于一个特定的${\theta }_i$（图中绿色直线的位置）来说，找一个目标函数在${\theta }_i$处的下界函数，假设该下界函数在${\theta }_{i+1}$（图中红色直线的位置）处取得最大值。如果用${\theta }_{i+1}$更新原来的${\theta }_i$，因为是下界函数，参数更新后下界函数的值变大了，则目标函数的值肯定也变大了。不断重复该过程，就可以收敛到一个局部最大值，注意不是全局最大值。

所以从之前直接求目标函数的解析解的方式，变成了迭代的求下界函数并更新参数，从而使目标函数值更大的过程。

下面继续进行公式的推导，假设$Q_i$是$z_i$的某一个分布，则
$$l(\theta )=\sum _{i=1}^n\log \sum _{z_i}p(x_i,z_i;\theta )=\sum _{i=1}^n\log \sum _{z_i}{Q}_i(z_i)\frac{p(x_i,z_i;\theta )}{Q_i(z_i)}$$

$$\ge \sum _{i=1}^n\sum _{z_i}{Q}_i(z_i)\log \frac{p(x_i,z_i;\theta )}{Q_i(z_i)}$$

下面来解释一下$\ge$成立的原因。我们可以把$\frac{p(x_i,z_i;\theta )}{Q_i(z_i)}$整体看作是一个随机变量，不妨记作$Y$，则$\log {\sum }_{z_i}{Q}_i(z_i)\frac{p(x_i,z_i;\theta )}{Q_i(z_i)}$就可以看作是关于$Y$的期望的函数值，即$\log EY$，

而${\sum }_{z_i}{Q}_i(z_i)\log \frac{p(x_i,z_i;\theta )}{Q_i(z_i)}$可被看作是函数值的期望，即$E\log (Y)$，由Jensen不等式可知，以上不等式成立。

到此为止，我们找到了目标函数在参数$\theta$处的下界函数的表达式，为了进一步简化，不妨设目标函数和下界函数在$\theta$处的函数值是相等的。

由上述可知，当且仅当$Y$为常数时等式成立。也就是$Q_i(z_i)$正比于$p(x_i,z_i;\theta )$，又因为$Q_i(z_i)$是概率分布，所以${\sum }_z{Q}_i(z_i)=1$。下面凑出一个满足这两个约数条件的$Q_i(z_i)$：
$$Q_i(z_i)=\frac{p(x_i,z_i;\theta )}{{\sum }_{j}p(x_i,z_j;\theta )}$$

$$=\frac{p(x_i,z_i;\theta )}{p(x_i;\theta )}$$

$$=p(z_i|x_i;\theta )$$

第一个等式相当于做了一个归一化，分母是所有的$z_i$相加的和，分子是一个特定的$z_i$，这就保证了${\sum }_z{Q}_i(z_i)=1$，同时$p(x_i,z_i;\theta )$和$Q_i(z_i)$的比值是$p(x_i;\theta )$，是一个定值。

找到下界函数之后，只需要求下界函数的最大值，然后更新对应的参数$\theta$就好了。

---

EM算法可以分为E步和M步两个主要步骤，其中E步是估计参数的值，M步是使下界函数最大化。其算法流程为：

(1). 初始化参数$\theta$

(2). E步：根据上一次迭代的模型参数来计算新的下界函数：
$$Q_i(z_i)=p(z_i|x_i;\theta )$$
(3). M步：将下界函数最大化以获得新的参数值：
$$\theta =arg\underset{\theta }{\max }\sum _i\sum _{z_i}{Q}_i(z_i)\log \frac{p(x_i,z_i;\theta )}{Q_i(z_i)}$$
(4). 重复(2)、(3)两步直到收敛到局部最小值。