[机器学习] EM算法

  在前面的朴素贝叶斯分类器推导中，我们是基于训练样本所有属性值都是已知的这个假设，但是实际应用中往往存在训练样本的某些属性值未知，此时就需要用到EM算法来进行参数估计。EM算法的全称为最大期望算法(Expectation-Maximization algorithm)，它是基于最大似然估计理论的一种优化算法，通常用来对存在隐变量的概率模型进行参数估计。

最大似然估计

  假设我们在校园里随机找了100个男生和100个女生，来调查学生的身高分布。首先我们统计抽样得到的这100个男生的身高，假设他们的身高服从高斯分布 $N(\mu ,{\sigma }^2)$，但是分布参数 $\mu ,\sigma$ 未知，需要从抽取的样本中进行估计，为了方便将参数记为 $\theta =[\mu ,\sigma {]}^T$。
  现在样本集为 X = { x 1 , x 2 , … , x N } X=\{x_1,x_2, \dots, x_N\} X={x1​,x2​,…,xN​}，这里 $x_i$ 表示第 $i$ 个样本的身高，$N=100$ 表示样本数。用 $p(x_i|\theta )$表示从分布为 $p(x|\theta )$ 的总体中抽到样本 $i$ 的概率，由于每个样本是独立地从总体中抽取的，那么同时抽到这100个男生的概率(联合概率)可以表示为似然函数 $L(\theta )$：
$$L(\theta )=L\left(x_1,x_2,\dots ,x_N;\theta \right)=\prod _{i=1}^{N}p\left(x_i;\theta \right)\theta \in \Theta \text{\hspace{1em}(1)}$$   似然函数反映的是在参数为 $\theta$ 的时候抽到这100个样本的概率。但是现在这100个样本已知而参数 $\theta$ 未知，所以，我们的任务就是根据这100个样本来估计一个 $\theta$，使得在当前 $\theta$ 取值下从总体中抽到这100个样本的概率最大。求得的这个 $\hat{\theta }$ 叫做 $\theta$ 的最大似然估计量：
$$\hat{\theta }=\arg \max L(\theta )\text{\hspace{1em}(2)}$$ 为了方便求解，通常情况下将式(1)转化为对数似然：
$$LL(\theta )=\ln (\theta )=\ln \prod _{i=1}^{N}p\left(x_i;\theta \right)=\sum _{i=1}^N\ln p\left(x_i;\theta \right)\text{\hspace{1em}(3)}$$最后对于上式求导，导数为0的点解得的 $\theta$ 就是估计出的参数。需要注意的是，上面这个例子中 $\theta =[\mu ,\sigma {]}^T$ 包含了两个参数，因此对于多个参数的情况，需要对每个参数求偏导，偏导为0的方程构成的方程组的解就是最终参数的估计值。

总结一下求最大似然估计值的一般步骤就是：

1. 写出似然函数；
2. 对似然函数取对数得到对数似然函数；
3. 对似然函数求导(或偏导)，得到导数(或偏导数)为0的似然方程(方程组)；
4. 求解似然方程或方程组，得到参数估计值。

EM算法

  上面最大似然函数例子中我们知道了怎么从男生的100个样本中估计出男生的身高分布参数，用同样的方法可以得到女生身高分布的参数。现在来考虑另一个问题——如果我们只得到了随机抽取的100个男生和100个女生的身高值，但是并没有记录每个身高值到底是女生还是男生的。这种情况下，两种高斯分布的样本混到一块儿了，这种情况下估计男生女生的身高分布就需要用到EM算法了。
  将上面这个例子一般化，可以描述为：对于一个包含 $m$ 个独立样本的样本集 X = { x ( 1 ) , x ( 2 ) , … , x ( m ) } X=\{x^{(1)},x^{(2)}, \dots, x^{(m)}\} X={x(1),x(2),…,x(m)}，每个样本 $i$ 有一个属性 $z_i$ 是未知的(隐变量)，现在需要对概率模型的参数 $\theta$ 做最大似然估计。在详细介绍EM算法前我们先来了解一下两个需要用到的知识。

1. Jansen不等式

  设 $f(x)$ 为定义域为实数的函数，若对于所有 $x$ 都有 $f^{\prime \prime }(x)\ge 0$，则 $f(x)$ 为凸函数，若 $f^{\prime \prime }(x)>0$，$f(x)$ 为严格凸函数；反之，若对于所有 $x$ 都有 $f^{\prime \prime }(x)\le 0$，则 $f(x)$ 为凹函数，当 $f^{\prime \prime }(x)$ 恒小于0时，$f(x)$ 为严格凹函数。
Jansen不等式表述如下：

 如果 $f(x)$ 为凸函数，$X$ 为随机变量，则 $f(E[X])\le E[f(X)]$；特别地，当 $f(x)$ 为严格凸函数时， $f(E[X])=E[f(X)]$ 当且仅当 $p(x=E[X])=1$，即 $X$ 为常量。同样，如果 $f(x)$ 为凹函数，则 $f(E[X])\ge E[f(X)]$。

下图可以直观地理解一下Jansen不等式：

2. Lazy Statistician规则

若 $Y$ 是随机变量 $X$ 的函数：$Y=g(X)$，且 $g(\cdot )$ 连续，则：

1. $X$ 是离散型随机变量，分布律为 $p(X=x_k)=p_k，k=1,2,\dots$，若 ${\sum }_{k=1}^{\infty }\mathrm{g}\left(x_k\right)p_k$ 绝对收敛，那么
   $$E[Y]=E[g(X)]=\sum _{k=1}^{\infty }\mathrm{g}\left(x_k\right)p_k\text{\hspace{1em}(4)}$$
2. $X$ 是连续型随机变量，概率密度函数为 $f(x)$，若 ${\int }_{-\infty }^{\infty }g(x)f(x)dx$绝对收敛，那么
   $$E[Y]=E[g(X)]={\int }_{-\infty }^{\infty }g(x)f(x)dx\text{\hspace{1em}(5)}$$

3. EM算法详细过程

  回顾完Jansen不等式和Lazy Statistician规则之后，我们来详细的介绍EM算法的具体过程。给定的训练集为 X = { x ( 1 ) , x ( 2 ) , … , x ( m ) } X=\{x^{(1)},x^{(2)}, \dots, x^{(m)}\} X={x(1),x(2),…,x(m)}，每个样本间相互独立，每个样本的隐变量 $z_i$ 未知，要估计概率模型的参数 $\theta$，我们首先建立似然函数，这里直接用对数似然表示：
$$\ell (\theta )=\sum _{i=1}^m\ln p\left(x^{(i)};\theta \right)\text{\hspace{1em}(6)}$$根据边际概率分布和联合概率分布的关系，上式可以写成：
$$\ell (\theta )=\sum _{i=1}^m\ln p\left(x^{(i)};\theta \right)=\sum _{i=1}^m\ln \sum _{z^{(i)}}p\left(x^{(i)},z^{(i)};\theta \right)\text{\hspace{1em}(7)}$$由于隐变量 $z$ 的存在，直接对上式进行最大似然估计得到 $\theta$ 的估计值式很困难的。EM算法的中心思想就是，首先建立 $\ell (\theta )$ 的下界(E步)，接下来最大化这个下界(M步)，不断地重复这个过程直到收敛。

  具体来说，首先我们假设隐变量 $z^{(i)}$ 服从某种分布 $Q_i(z^{(i)})$，满足 ${\sum }_z{Q}_i(z^{(i)})=1,Q_i(z^{(i)})\ge 0$，那么
$$\begin{array}{rl}\ell (\theta )& =\sum _i^m\ln \sum _{z^{(i)}}p\left(x^{(i)},z^{(i)};\theta \right)\\ & =\sum _i^m\ln \sum _{z^{(i)}}{Q}_i\left(z^{(i)}\right)\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}\\ & ⩾\sum _i^m\sum _{z^{(i)}}{Q}_i\left(z^{(i)}\right)\ln \frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}\end{array}\text{\hspace{1em}(8)}$$ 从第2步到第3步的不等式由Jansen不等式和Lazy Statistician规则可以得到，其中 $\ln (x)$ 为凹函数。

  式(8)可以看作是对 $\ell (\theta )$ 求了下界，此时若 $\theta$ 给定，那么 $\ell (\theta )$ 取决于 $Q_i(z^{(i)})$ 和 $p(x^{(i)},z^{(i)})$。这样，我们通过调整 $Q_i(z^{(i)})$ 和 $p(x^{(i)},z^{(i)})$ 使得 $\ell (\theta )$ 的下界不断上升，方便起见我们将下界表示为
$$J(Q,z)=\sum _i^m\sum _{z^{(i)}}{Q}_i\left(z^{(i)}\right)\ln \frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}$$ 最大化下界就是得到使式(8)中不等式变成等式的 $J(Q,z)$，这样$J(Q,z)$ 就等价于似然函数 $\ell (\theta )$。那么怎么让式(8)中的等号成立呢？根据Jansen不等式，当且仅当随机变量为常数的时候等式才成立，即：
$$\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}=C\text{\hspace{1em}(9)}$$ 这里 $C$ 表示常数，不依赖于 $z^{(i)}$，根据上式可以进一步得到
$$\begin{gathered} p\left(x^{(i)},z^{(i)};\theta \right)=C\left(Q_i\left(z^{(i)}\right)\right) \\ \sum _{z}p\left(x^{(i)},z^{(i)}\theta \right)=C\left(\sum _z{Q}_i\left(z^{(i)}\right)\right) \end{gathered}$$ 又因为 ${\sum }_z{Q}_i(z^{(i)})=1$，所以
$$\sum _{z}p\left(x^{(i)},z^{(i)};\theta \right)=C$$ 进一步推导得到
$$\begin{array}{rl}{Q}_i(z^{(i)})& =\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{{\sum }_{z}p\left(x^{(i)},z^{(i)};\theta \right)}\\ & =\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{p\left(x^{(i)};\theta \right)}\\ & =p(z^{(i)}|x^{(i)};\theta )\end{array}\text{\hspace{1em}(10)}$$ 上式就是在固定了 $\theta$ 之后如何拉升似然函数 $\ell (\theta )$ 的下界使得式(8)中不等式的等号成立，可以看到这里的 $Q_i(z^{(i)})$ 是个后验概率，这就是EM算法的E步。 接下来就是M步，在求得 $Q_i(z^{(i)})$ 后调整 $\theta$ 使得 $\ell (\theta )$ 最大化。总结一下EM算法就是以下框架：

 重复以下步骤直到算法收敛（达到最大迭代步骤或是两次迭代的 $\ell (\theta )$ 相差非常小）：
  E步：$$Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta )$$
  M步：$$\theta =\arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}{Q}_i\left(z^{(i)}\right)\log \frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}$$

可以证明在两次迭代 $t$ 和 $t+1$ 过程中 $\ell ({\theta }^{(t)})\ge \ell ({\theta }^{(t+1)})$，具体证明过程这里不讲了，有兴趣可以参考。https://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html

总结

EM算法有很多应用，最典型的就是K-means聚类和混合高斯模型。

参考资料

文章参考整理自以下三篇博客，在此致谢
从最大似然到EM算法浅解
（EM算法）The EM Algorithm
EM算法