EM算法原理详解与高斯混合模型

借助于machine learning cs229和文章【1】中的内容把EM算法的过程顺一遍，加深一下印象。
关于EM公式的推导，一般会有两个证明，一个是利用Jesen不等式，另一个是将其分解成KL距离和L函数，本质是类似的。

下面介绍Jensen EM的整个推导过程。

1. Jensen不等式

   回顾优化理论中的一些概念。设f是定义域为实数的函数，如果对于所有的实数x，$f^{''}(x)\ge 0$，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的（$H\ge 0$），那么f是凸函数。如果$f^{''}(x)\gt 0$或者$H\gt 0$，那么称f是严格凸函数。

   Jensen不等式表述如下：

   如果f是凸函数，X是随机变量，那么
   

   $$E[f(x)]\ge f(E[x])$$
   

   特别地，如果f是严格凸函数，那么

   $$E[f(x)]\gt f(E[x])$$ 当且仅当$p(X=E(X))=1$，也就是说X是常量。

   这里我们将$f(E[X])$简写为$f(EX)$。

   如果用图表示会很清晰：

   图中，实线$f$是凸函数，$X$是随机变量，有$0.5$的概率是$a$，有$0.5$的概率是$b$。（就像掷硬币一样）。$X$的期望值就是$a$和$b$的中值了，图中可以看到

   $$E[f(x)]\ge f(E[x])$$ 成立。

   当$f$是（严格）凹函数当且仅当$-f$是（严格）凸函数。

   Jensen不等式应用于凹函数时，不等号方向反向，也就是$E[f(x)]\le f(E[x])$。

---

1. EM算法

   给定的训练样本是$\{x^{(1)},...,x^{(m)}\}$，样例间独立，我们想找到每个样例隐含的类别z，能使得p(x,z)最大。p(x,z)的最大似然估计如下：

$$\begin{align} \mathcal{l}(\theta) &=\sum_{i=1}^{m}log\ p(x^{(i)};\theta)\\ &=\sum_{i=1}^{m}log\sum_z p(x^{(i)},z^{(i)};\theta) \end{align}$$
第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求 $\theta$一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。
EM是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化 $\mathcal{l}(\theta)$，我们可以不断地建立 $\mathcal{l}(\theta)$的下界（E步），然后优化下界（M步）。这句话比较抽象，看下面的。
对于每一个样例i，让 $Q_i$表示该样例隐含变量z的某种分布， $Q_i$满足的条件是 $\sum_z Q_i(z)=1,Q_i(z)\ge 1$。（如果z是连续性的，那么clip_image032[2]是概率密度函数，需要将求和符号换做积分符号）。比如要将班上学生聚类，假设隐藏变量z是身高，那么就是连续的高斯分布。如果按照隐藏变量是男女，那么就是伯努利分布了。

可以由前面阐述的内容得到下面的公式：

$$\begin{align} \sum_{i=1}^{m}log p(x;\theta) &=\sum_{i=1}^{m}log\sum_{z^{(i)}} p(x^{(i)},z^{(i)};\theta) &(1)\\ &=\sum_{i=1}^{m}log\sum_z Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} &(2)\\ &\ge \sum_{i=1}^{m}\sum_z Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} &(3) \end{align}$$
（1）到（2）比较直接，就是分子分母同乘以一个相等的函数。（2）到（3）利用了Jensen不等式，考虑到 $log（x）$是凹函数（二阶导数小于0），而且
$$\sum_z Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$
就是 $\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$的期望（回想期望公式中的Lazy Statistician规则）
设Y是随机变量X的函数 $y=g(x)$（g是连续函数），那么
（1） X是离散型随机变量，它的分布律为 $P(X=x_k)=p_k,k=1,2,...$.若 $\sum_{k=1}^{\infty}g(x_k)p_k$绝对收敛，则有
$E(Y)=E[g(X)]=\sum_{k=1}^{\infty}g(x_k)p_k$
（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$
对应于上述问题，Y是 $\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$，X是 $z^{(i)}$， $Q_i(z^{(i)})$是 $p_k$，g是 $z^{(i)}$到 $\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$的映射。这样解释了式子（2）中的期望，再根据凹函数时的Jensen不等式：
$$f(E_{z^{(i)}\sim Q_i}[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}])\ge E_{z^{(i)}\sim Q_i}[f(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})})]$$
可以得到（3）。
这个过程可以看作是对 $l(\theta)$求了下界。对于 $Q_i$的选择，有多种可能，那种更好的？假设 $\theta$已经给定，那么 $l(\theta)$的值就决定于 $Q_i(z^{(i)})$和 $p(x^{(i)},z^{(i)})$了。我们可以通过调整这两个概率使下界不断上升，以逼近 $l(\theta)$的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率能够等价于 $l(\theta)$了。按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：
$$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c$$
c为常数，不依赖于 $z^{(i)}$。对此式子做进一步推导，我们知道 $\sum_z{Q_i(z^{(i)})}=1$，那么也就有 $\sum_z p(x^{(i)},z^{(i)};\theta)=c$，（多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是c），那么有下式：
$$\begin{align} Q_i(z^{(i)})&=\frac{ p(x^{(i)},z^{(i)};\theta)}{\sum_z p(x^{(i)},z^{(i)};\theta)}\\ &=\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i);\theta})}\\ &=p(z^{(i)};x^{(i)}\theta) \end{align}$$
至此，我们推出了在固定其他参数 $\theta$后， $Q_i(z^{(i)})$的计算公式就是后验概率，解决了 $Q_i(z^{(i)})$如何选择的问题。这一步就是E步，建立 $l(\theta)$的下界。接下来的M步，就是在给定 $Q_i(z^{(i)})$后，调整 $l(\theta)$，去极大化 $l(\theta)$的下界（在固定 $Q_i(z^{(i)})$后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：
循环重复直到收敛 {
（E步）对于每一个i，计算
$$Q_i(z^{(i)})=p(z^{(i)};x^{(i)}\theta)$$
（M步）计算
$$\theta=\arg \max_{\theta}\sum_{i=1}^{m}\sum_z Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$
那么究竟怎么确保EM收敛？假定 $\theta^{t}$和 $\theta^{t+1}$是EM第t次和t+1次迭代后的结果。如果我们证明了 $l(\theta^t)\le l(\theta^{t+1})$，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。下面来证明，选定 $\theta^{(t)}$后，我们得到E步
$Q_i(z^{(i)})=p(z^{(i)};x^{(i)}\theta)$
这一步保证了在给定 $\theta^{(t)}$时，Jensen不等式中的等式成立，也就是
$$l(\theta^{(t)})=\sum_{i=1}^{m}\sum_z Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$
然后进行M步，固定 $Q_i(z^{(i)})$，并将 $\theta^{(t)}$视作变量,由EM的操作可以推导出以下式子成立：
$$\begin{align} l(\theta^{(t+1)}) &\ge\sum_{i=1}^{m}\sum_z Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i(z^{(i)})} &(4)\\ &\ge \sum_{i=1}^{m}\sum_z Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i(z^{(i)})} &(5)\\ &=l(\theta)&(6) \end{align}$$
等式（4）成立是由于琴生不等式，等式（5）成立是由于我们在M步取得是max操作。
如果我们定义
$$J(Q,\theta)=\sum_{i=1}^{m}\sum_{z^{(i)}} Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$
从前面的推导中我们知道 $l(\theta)\ge J(Q,\theta)$，EM可以看作是J的坐标上升法，E步固定 $\theta$，优化 $Q$，M步固定$Q$优化 $\theta$。

1. 重新审视混合高斯模型

   我们已经知道了EM的精髓和推导过程，再次审视一下混合高斯模型。之前提到的混合高斯模型的参数$\Phi$和$\mu$，为了简单，这里在M步只给出$\Phi$和$\mu$的推导方法。

E步很简单，按照一般EM公式得到：

$$\begin{align} w^{(i)}_j & =Q_i(z^{(i)}=j)=P(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)\\ & =\frac{\phi_j p(x^{(i)}|z^{(i)}=j;\phi,\mu,\Sigma)}{\Sigma_{k}^{j=1}\phi_j p(x^{(i)}|z^{(i)}=j;\phi,\mu,\Sigma)} \end{align}$$
简单解释就是每个样例i的隐含类别 $z^{(i)}$为j的概率可以通过后验概率计算得到。 $\phi$是每个类的概率。
在M步中，我们需要在固定 $Q_i(z^{(i)})$后最大化最大似然估计，也就是
$$\begin{align} &\sum_{i=1}^{m}\sum_{z^{(i)}} Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_i(z^{(i)})} \\ =&\sum_{i=1}^{m}\sum_{j=1}^{k} Q_i(z^{(i)}=j)log\frac{p(x^{(i)}|z^{(i)}=j;\phi,\mu,\Sigma)p(z^{(i)}=j;\phi)}{Q_i(z^{(i)})} \\ =&\sum_{i=1}^{m}\sum_{j=1}^{k} w^{(i)}_{j}log\frac {\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma^{-1}_{j}(x^{(i)}-\mu_j)).\phi_j} {w^{(i)}_j} \end{align}$$
这是将 $z^{(i)}$的k种情况展开后的样子，未知参数 $\phi_j$, $\mu_j$和 $\Sigma_j$。
固定 $\phi_j$,和 $\Sigma_j$，对 $\mu_j$求导得
$$\begin{align} &\nabla_{\mu_j}\sum_{i=1}^{m}\sum_{z^{(i)}} Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_i(z^{(i)})} \\ =&-\nabla_{\mu_j}\sum_{i=1}^{m}\sum_{j=1}^{k} w^{(i)}_{j}\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma^{-1}_{j}(x^{(i)}-\mu_j))\\ =&\frac{1}{2}\sum_{i=1}^{m}w^{(i)}_{j}\nabla_{\mu_j}(2\mu^T_j\Sigma^{-1}_{j}x^{(i)}-\mu^T_j\Sigma^{-1}\mu_j)\\ =&\sum_{i=1}^{m}w^{(i)}_{j}(\Sigma^{-1}_{j}x^{(i)}-\Sigma^{-1}\mu_j) \end{align}$$
等于0时，得到
$\mu_j=\frac{\sum_{i=1}^{m}w^{(i)}_{j}x^{(i)}}{\Sigma_{i=1}^{m}w^{(i)}_{j}}$
这就是我们之前模型中的 $\mu$的更新公式。
然后推导 $\phi$的更新公式。看之前得到的
$$\sum_{i=1}^{m}\sum_{j=1}^{k} w^{(i)}_{j}log\frac {\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma^{-1}_{j}(x^{(i)}-\mu_j)).\phi_j} {w^{(i)}_j}$$
分子和分母上与 $\phi$无关的常数都可以通过log提取出来，，实际上需要优化的公式是：
$$\sum_{i=1}^{m}\sum_{j=1}^{k} w^{(i)}_{j}log \phi_j$$
需要知道的是， $\phi_j$还需要满足一定的约束条件就是 $\Sigma_{j=1}^{k}\phi_j=1$。
这个优化问题我们很熟悉了，直接构造拉格朗日乘子。
$$L(\phi)=\sum_{i=1}^{m}\sum_{j=1}^{k} w^{(i)}_{j}log \phi_j+ \beta(\sum_{j=1}^{k}\phi_j-1)$$
还有一点就是 $\phi_j\gt 0$，但这一点会在得到的公式里自动满足。
求导得，
$\frac{\partial}{\partial \phi_j}L(\phi)=\sum_{i=1}^{m}\frac{w^{(i)}_{j}}{\phi_j}+\beta$
等于0，得到
$\phi_j=\frac{\sum_{i=1}^{m}w^{(i)}_{j}}{-\beta}$
也就是说 $\phi_j\propto \sum_{i=1}^{m}{w^{(i)}_{j}}$再次使用 $\sum_{j=1}^{k}\phi_j=1$，得到
$-\beta=\sum_{i=1}^{m}\sum_{j=1}^{k} w^{(i)}_{j}=\Sigma_{i=1}^{m}1=m$
这样就神奇地得到了 $\beta$。
那么就顺势得到M步中 $\phi_j$的更新公式：
$\phi_j=\frac{1}{m}\sum_{i=1}^{m}{w^{(i)}_{j}}$
$\Sigma$的推导也类似，不过稍微复杂一些，毕竟是矩阵。结果在之前的混合高斯模型中已经给出。
$$\Sigma=\frac{\sum_{i=1}^{m}w^{(i)}_{j}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum_{i=1}^{m}w^{(i)}_{j}}$$

关于EM证明的补充，一般情况下可以通过Jensen不等式得到$\sum_{i=1}^{m}log p(x;\theta)$的下界来求解，但是将该公式展开似乎更直观一些。
首先是将$p(x;\theta)$拆分开来，我们注意到$p(x,z;\theta)=p(x;\theta)p(z|x;\theta)$这里有联合概率，x的生成概率，z的后验概率，取$log$然后做一下变换

$$\begin{align} ln p(x;\theta)&=ln p(x,z;\theta)- ln p(z|x;\theta）\\ &=ln p(x,z;\theta)-ln q(z)-[ ln p(z|x;\theta)-ln q(z)]\\ &=ln \frac{p(x,z;\theta)}{q(z)}-ln \frac{ln p(z|x;\theta)}{q(z)} \end{align}$$
同时乘以 $q(z)$并对 $q(z)$求和得到
$$\Sigma_z q(z)lnp(x;\theta)=\Sigma_z [q(z)(ln \frac{p(x,z;\theta)}{q(z)}-ln\frac{ln p(z|x;\theta)}{q(z)})]$$
由于z与 $p(x;\theta)$独立，并且 $\sigma_z(q)=1$,所以又
$$ln p(x;\theta)=\Sigma_z q(z)n \frac{p(x,z;\theta)}{q(z)}- \Sigma_z q(z)ln\frac{ln p(z|x;\theta)}{q(z)}$$
写作
$$ln p(x;\theta)=\mathcal{L}(q;\theta)+KL(q||p)$$
看到这里应该很明白了，在NG的推导中放缩掉的是这里的KL距离，最后得到的是这里的 $\mathcal{L}(q;\theta)$,当 $q(z)$和 $p(z;x,\theta)$相等时KL距离为0，前面的推导的等式成立。

E-step
假设当前的参数为 $\theta_{old}$,则Estep可以被描述为：固定 $\theta_{old}$找一个分布 $q(z)$,使得 $L(q,\theta_{old})$最大
由于 $ln p(x;\theta_{old}$与z无关则使 $L(q,\theta_{old})$最大姐使KL(q||p)最小（=0），也就是说 $q(z)=p(z;x,\theta_{old})$

回想高斯混合分布，用 $p(z;x,\mu,\sigma,\pi)$去求 $E[z]$，这个时候 $KL（q||p）=0$
M-step：固定 $q(z)$找一组参数 $\theta_{new}$,使得 $L(q,\theta_{new})$最大，
$ln p(x;\theta)$的增大可能来自两部分 $L(q,\theta_{new})$和 $KL(q||p)$（因为此时 $p(z;x,\theta_{new})$）而 $q(z;x,\theta_{old}),p\neq q$

这里可以看到，L的下届提高了，$ln(p(x;\theta)$的下届也提高了。

参考：
【1】
http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html#2831962
【2】http://wenku.baidu.com/view/d5e6973a87c24028915fc361.html
【3】http://cs229.stanford.edu/materials.html
【4】http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/