EM算法原理分析

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

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

是不是一目了然？简单解释一下：假设X是一个随机变量，有0.5的概率落在a点，有0.5的概率落在b点，因此X的期望E[X]$E[X]$便落在a,b 的中点处。根据f是凸函数，我们可以在图上画出f(a),f(b),f(E[X])$f(a),f(b),f(E[X])$的位置，而E[f(X)]$E[f(X)]$则落在f(a),f(b)$f(a),f(b)$的中点处。
由上图可知，因为f是凸函数，所以有E[f(X)]≥f(EX)$E[f(X)]\ge f(EX)$。同理，如果f是凹函数，则有E[f(X)]≤f(EX)$E[f(X)]\le f(EX)$。
EM算法
假设我们有m个独立样本(独立性假设){x(1),...,x(m)}$\{x^{(1)},...,x^{(m)}\}$,给定以下似然函数：

我们希望求出模型p(x,z)$p(x,z)$的参数θ$\theta$. 然而，由于存在隐藏变量z$z$,θ$\theta$的求解是很困难的，如果能够提前得到z$z$,那么最大似然估计将变得简单起来。（请记住这一点，因为后面的EM算法的E步其实就相当于给z做了一个先验假设，然后再做优化）
对于每个样本i，假设Qi$Q_i$是关于z的分布（∑zQi(z)=1,Qi(z)≥0$\sum _z{Q}_i(z)=1,Q_i(z)\ge 0$）,因此可得到下列不等式：

最后一步怎么得来的呢？其实就是用到了Jensen不等式。特别的，f(x)=log x$f(x)=logx$是一个凹函数，因为f′′(x)=−1x2<0$f^{{}^{″}}(x)=\frac{-1}{x^2}<0$.因此有E[f(x)]≤f(E(x))$E[f(x)]\le f(E(x))$，其中自变量x为p(x(i),z(i);θ)Qi(z(i))$\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$,代入得：
.
综上便可得到上文所述不等式。
那么，不等式什么时候取等号呢？其实上文的定理已经提到了，当自变量为常数时等号成立，对应到我们得不等式中，即：
.
事实上，我们知道∑zQi(z)=1$\sum _z{Q}_i(z)=1$,因此我们可以得到下面得推导：

也就是说，我们可以简单设置Qi$Q_i$为在参数θ$\theta$下给定x(i)$x^{(i)}$时，关于z(i)$z^{(i)}$的后验分布。
因此，我们可以得到EM算法的迭代过程如下：
循环以下两步直到收敛{
（E-step）对于每个样本i,

（M-step）

}
那么我们怎么知道EM算法是否收敛呢？我们假设θ(t)和θ(t+1)$\theta (t)和\theta (t+1)$为迭代过程中的参数，那么我们只要证明l(θ(t))≤l(θ(t+1))$l(\theta (t))\le l(\theta (t+1))$,那么就可以得到EM算法是在不断优化，直至收敛。顺着这个思想，我们假设Qi(z(i)):=p(z(i)|x(i);θ)$Q_i(z^{(i)}):=p(z^{(i)}|x^{(i)};\theta )$,此时

参数θ(t+1)${\theta }^{(t+1)}$通过最大化等式右边的式子获得，因此：

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

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

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

其中，ϕj${\varphi }_j$是系数，ϕj≥0,∑kj=1ϕj=1${\varphi }_j\ge 0,\sum _{j=1}^k{\varphi }_j=1$;p(y|θj)$p(y|{\theta }_j)$是高斯分布密度，θj=(μj,σ2j)=((μj,Σj)${\theta }_j=({\mu }_j,{\sigma }_j^2)=(({\mu }_j,{\mathrm{\Sigma }}_j)$,

称为第j个分模型。
一般混合模型可以由任意概率分布密度代替上式中的高斯分布密度，我们这里只介绍最常用的高斯混合模型。

E-step：计算

即w(i)j$w_j^{(i)}$是针对第i个样本，在参数为ϕ,μ,Σ$\varphi ,\mu ,\mathrm{\Sigma }$已知样本特征x(i)$x^{(i)}$的情况下，属于第j个分模型的概率。
M-step：最大化以下式子优化参数ϕ,μ,Σ$\varphi ,\mu ,\mathrm{\Sigma }$：

首先我们关于μl${\mu }_l$最大化以上式子。将L对μl${\mu }_l$求导，得到：

令导数等于零，可得到μl${\mu }_l$的更新规则如下：

至于Σ$\mathrm{\Sigma }$的更新跟μl${\mu }_l$类似，不再赘述。下面讲一下ϕ$\varphi$的更新。
通过观察式子，我们可以把无关变量去掉，得到：

另一方面，因为ϕj=p(z(i)=j;ϕ)${\varphi }_j=p(z^{(i)}=j;\varphi )$，所以有约束条件∑kj=1ϕj=1$\sum _{j=1}^k{\varphi }_j=1$.因此，我们使用拉格朗日乘子β$\beta$将有约束问题转换成无约束问题，如下：

值得注意的是，这里并没有把约束条件ϕj>0${\varphi }_j>0$加上，这是为什么呢？别急，下文会提到。
对以上式子求导，得到：

令导数等于零，可得到ϕj${\varphi }_j$的更新规则如下：

使用约束条件∑kj=1ϕj=1$\sum _{j=1}^k{\varphi }_j=1$，我们可以得到−β=∑mi=1∑kj=1w(i)j=∑mi=11=m(使用条件w(i)j=Qi(z(i)=j),从而∑kj=1w(i)j=1)$-\beta =\sum _{i=1}^m\sum _{j=1}^k{w}_j^{(i)}=\sum _{i=1}^{m}1=m(使用条件w_j^{(i)}=Q_i(z^{(i)}=j),从而\sum _{j=1}^k{w}_j^{(i)}=1)$，因此，我们可以进一步化简得到：

我们可以看到，ϕj${\varphi }_j$恒大于零，默认满足约束条件ϕj>0${\varphi }_j>0$。

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

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