EM算法原理分析

EM算法主要用于含有隐藏变量的参数估计问题。
在将EM算法之前，先讲一下Jensen不等式。
定理：假设f是一个凸函数，X是随机变量，即：

$$E[f(X)]\ge f(EX)$$
此外，如果f是严格凸的，当且仅当$X=E[X]=常数（不再是随机变量）时E[f(X)]=f(EX)$.
不理解的可以看下面的图：

是不是一目了然？简单解释一下：假设X是一个随机变量，有0.5的概率落在a点，有0.5的概率落在b点，因此X的期望$E[X]$便落在a,b 的中点处。根据f是凸函数，我们可以在图上画出$f(a),f(b),f(E[X])$的位置，而$E[f(X)]$则落在$f(a),f(b)$的中点处。
由上图可知，因为f是凸函数，所以有$E[f(X)]\ge f(EX)$。同理，如果f是凹函数，则有$E[f(X)]\le f(EX)$。
EM算法
假设我们有m个独立样本(独立性假设)$\{x^{(1)},...,x^{(m)}\}$,给定以下似然函数：
$$l(\theta)=\sum_{i=1}^mlog \ p(x;\theta)=\sum_{i=1}^mlog\sum_{z}p(x,z;\theta)$$
我们希望求出模型$p(x,z)$的参数$\theta$. 然而，由于存在隐藏变量$z$,$\theta$的求解是很困难的，如果能够提前得到$z$,那么最大似然估计将变得简单起来。（请记住这一点，因为后面的EM算法的E步其实就相当于给z做了一个先验假设，然后再做优化）
对于每个样本i，假设$Q_i$是关于z的分布（$\sum_zQ_i(z)=1,Q_i(z)\ge0$）,因此可得到下列不等式：
$$\begin{align} l(\theta)&=\sum_{i=1}^mlog \ p(x^{(i)};\theta)\\ &=\sum_{i=1}^mlog\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta) \\ & = \sum_{i=1}^mlog\sum_{z^{(i)}}Q_i(z^{(i)}) {p(x^{(i)},z^{(i)};\theta)\over Q_i(z^{(i)})}\\ & \ge \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})log{p(x^{(i)},z^{(i)};\theta)\over Q_i(z^{(i)})} \end{align}$$
最后一步怎么得来的呢？其实就是用到了Jensen不等式。特别的，$f(x)=log \ x$是一个凹函数，因为$f^{''}(x)={-1\over x^2}<0$.因此有$E[f(x)]\le f(E(x))$，其中自变量x为${p(x^{(i)},z^{(i)};\theta)\over Q_i(z^{(i)})}$,代入得：
$$f(E_{z^{(i)}\sim Q_i}[{p(x^{(i)},z^{(i)};\theta)\over Q_i(z^{(i)})}])\ge E_{z^{(i)}\sim Q_i}[f({p(x^{(i)},z^{(i)};\theta)\over Q_i(z^{(i)})})]$$ .
综上便可得到上文所述不等式。
那么，不等式什么时候取等号呢？其实上文的定理已经提到了，当自变量为常数时等号成立，对应到我们得不等式中，即：
$${p(x^{(i)},z^{(i)};\theta)\over Q_i(z^{(i)})}=c$$ .
事实上，我们知道$\sum_zQ_i(z)=1$,因此我们可以得到下面得推导：
$$\begin{align} Q_i(z^{(i)})&={p(x^{(i)},z^{(i)};\theta)\over \sum_zp(x^{(i)},z;\theta)} \\ & = {p(x^{(i)},z^{(i)};\theta)\over p(x^{(i)};\theta)}\\ & = p(z^{(i)}|x^{(i)};\theta) \end{align}$$
也就是说，我们可以简单设置$Q_i$为在参数$\theta$下给定$x^{(i)}$时，关于$z^{(i)}$的后验分布。
因此，我们可以得到EM算法的迭代过程如下：
循环以下两步直到收敛{
（E-step）对于每个样本i,
$$Q_i(z^{(i)}):=p(z^{(i)}|x^{(i)};\theta).$$
（M-step）
$$\theta:=argmax_\theta \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})log{p(x^{(i)},z^{(i)};\theta)\over Q_i(z^{(i)})}.$$
}
那么我们怎么知道EM算法是否收敛呢？我们假设$\theta(t)和\theta(t+1)$为迭代过程中的参数，那么我们只要证明$l(\theta(t))\le l(\theta(t+1))$,那么就可以得到EM算法是在不断优化，直至收敛。顺着这个思想，我们假设$Q_i(z^{(i)}):=p(z^{(i)}|x^{(i)};\theta)$,此时
$$l(\theta^{(t)})=\sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})log{p(x^{(i)},z^{(i)};\theta^{(t)})\over Q_i(z^{(i)})}$$
参数$\theta^{(t+1)}$通过最大化等式右边的式子获得，因此：
$$\begin{align} \theta^{(t+1)}&\ge \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})log{p(x^{(i)},z^{(i)};\theta^{(t+1)})\over Q_i(z^{(i)})}\\ &\ge \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})log{p(x^{(i)},z^{(i)};\theta^{(t)})\over Q_i(z^{(i)})}\\ &=l(\theta^{(t)}) \end{align}$$

当$\theta$为$\theta^{(t+1)}$时，${p(x^{(i)},z^{(i)};\theta^{(t+1)})\over Q_i^{(t)}(z^{(i)})}$不一定为常数了，所以等号不一定成立，因此上述第一个式子为大于等于。
至于第二个不等式，由EM算法的M步可知，$\theta^{(t+1)}$是通过最大化上一步的函数值得到的，即：

$$\theta^{(t+1)}:=argmax_\theta \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})log{p(x^{(i)},z^{(i)};\theta^{(t)})\over Q_i(z^{(i)})}.$$
再把$\theta^{(t+1)}$迭代回去，得到的函数值肯定会大于等于上一步的函数值，因此第二个不等式成立。
综上，通过EM算法，我们总可以得到$l(\theta^{(t+1)})\ge l(\theta^{(t)})$，从而不断优化，直到收敛，收敛条件是函数值增长小于等于阈值（阈值自己设定）时，停止迭代。

二 高斯混合模型（Gaussian Misture Model, GMM）
EM算法的一个重要应用就是高斯混合模型的参数估计。
高斯混合模型（Gaussian Misture Model, GMM）是指具有如下形式的概率分布模型：

$$p(y|\theta)=\sum_{j=1}^k\phi_jp(y|\theta_j)$$
其中，$\phi_j$是系数，$\phi_j\ge 0,\sum_{j=1}^k\phi_j=1$;$p(y|\theta_j)$是高斯分布密度，$\theta_j=(\mu_j,\sigma_j^2)=((\mu_j,\Sigma_j)$,
$$p(y|\theta_j)={1\over (2\pi)^{1\over 2}\sigma_j}exp(-{(y-\mu_j)^2\over 2\sigma_j^2})$$
称为第j个分模型。
一般混合模型可以由任意概率分布密度代替上式中的高斯分布密度，我们这里只介绍最常用的高斯混合模型。

E-step：计算

$$w_j^{(i)}=Q_i(z^{(i)}=j)=P(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma).$$
即$w_j^{(i)}$是针对第i个样本，在参数为$\phi,\mu,\Sigma$已知样本特征$x^{(i)}$的情况下，属于第j个分模型的概率。
M-step：最大化以下式子优化参数$\phi,\mu,\Sigma$：

$$\begin{align} L&=\sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})log{p(x^{(i)},z^{(i)};\phi,\mu,\Sigma)\over Q_i(z^{(i)})} \\ & = \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)}=j)log{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)\over Q_i(z^{(i)})=j} \\ & = \sum_{i=1}^m \sum_{z^{(i)}} w_j^{(i)}log{{1\over (2\pi)^{1\over 2}|\Sigma_j|^{1\over 2}}exp(-{1\over 2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j))\cdot \phi_j\over w_j^{(i)}} \end{align}$$

首先我们关于$\mu_l$最大化以上式子。将L对$\mu_l$求导，得到：

$$\begin{align} {\partial L \over \partial \mu_l}&=\nabla_{\mu_l}\sum_{i=1}^m \sum_{j=1}^k w_j^{(i)}log{{1\over (2\pi)^{1\over 2}|\Sigma_j|^{1\over 2}}exp(-{1\over 2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j))\cdot \phi_j\over w_j^{(i)}} \\ & =\nabla_{\mu_l} \sum_{i=1}^m \sum_{j=1}^k w_j^{(i)} {1\over 2} (x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)\\ & = {1\over 2}\sum_{i=1}^mw_l^{(i)}\nabla_{\mu_l} 2\mu_l^T\Sigma_l^{-1}x^{(i)}-\mu_l^T\Sigma_l^{-1}\mu_l\\ & = \sum_{i=1}^mw_l^{(i)}(\Sigma_l^{-1}x^{(i)}-\Sigma_l^{-1}\mu_l) \end{align}$$

令导数等于零，可得到$\mu_l$的更新规则如下：

$$\mu_l:={\sum_{i=1}^m w_l^{(i)}x^{(i)}\over \sum_{i=1}^m w_l^{(i)}}.$$
至于$\Sigma$的更新跟$\mu_l$类似，不再赘述。下面讲一下$\phi$的更新。
通过观察式子，我们可以把无关变量去掉，得到：
$$L(\phi)=\sum_{i=1}^m \sum_{j=1}^kw_j^{(i)}log \phi_j.$$
另一方面，因为$\phi_j=p(z^{(i)}=j;\phi)$，所以有约束条件$\sum_{j=1}^k\phi_j=1$.因此，我们使用拉格朗日乘子$\beta$将有约束问题转换成无约束问题，如下：
$$L(\phi)=\sum_{i=1}^m \sum_{j=1}^kw_j^{(i)}log \phi_j+\beta(\sum_{j=1}^k\phi_j-1)$$
值得注意的是，这里并没有把约束条件$\phi_j>0$加上，这是为什么呢？别急，下文会提到。
对以上式子求导，得到：
$${\partial L(\phi) \over \partial \phi_j}=\sum_{i=1}^m{w_j^{(i)} \over \phi_j}+\beta$$
令导数等于零，可得到$\phi_j$的更新规则如下：
$$\phi_j:={\sum_{i=1}^mw_j^{(i)} \over -\beta}.$$
使用约束条件$\sum_{j=1}^k\phi_j=1$，我们可以得到$-\beta=\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}=\sum_{i=1}^m1=m(使用条件w_j^{(i)}=Q_i(z^{(i)}=j),从而\sum_{j=1}^kw_j^{(i)}=1)$，因此，我们可以进一步化简得到：
$$\phi_j:={1 \over m}\sum_{i=1}^mw_j^{(i)}.$$
我们可以看到，$\phi_j$恒大于零，默认满足约束条件$\phi_j>0$。

再简单说明一下我理解的EM算法与Kmeans算法的联系与区别：
联系：Kmeans算法可以看作EM算法的一个特例，Kmeans中的簇即为EM算法中的隐藏变量；
区别：Kmeans中每一个数据点都只属于一个簇中，属于硬分隔；
而EM算法使用后验概率的方法，相当于一个数据点分到每一个簇都有一个概率，概率和为1.

参考：吴恩达CS229 Lecture notes “The EM algorithm”
《统计学习方法》（李航著）