机器学习之EM算法

EM算法是一种迭代算法，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计，我们经常会从样本观察数据中，找出样本的模型参数。 最常用的方法就是极大化模型分布的对数似然函数。

但是在一些情况下，我们得到的观察数据有未观察到的隐含数据，由于我们有未知的隐含数据和模型参数，因而无法直接用极大化对数似然函数得到模型分布的参数。怎么办呢？这就是EM算法可以派上用场的地方了。

EM算法解决这个的思路是使用启发式的迭代方法，既然我们无法直接求出模型分布参数，那么我们可以先猜想隐含数据（隐变量）「EM算法的E步」，接着基于观察数据和猜测的隐含数据一起来极大化对数似然，求解我们的模型参数*「EM算法的M步然后继续极大化对数似然，求解我们的模型参数（EM算法的M步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。

一、Jensen不等式与Lazy Statistician规则

Jensen不等式的意义是：如果f是凸函数，X是随机变量，那么函数的期望大于等于期望的函数，即：

$E(f(X))\ge f(E(X))$

特别地，如果f是严格凸函数，那么$E(f(X))=f(E(X))$当且仅当X是常数时取等号。

因为当X为常数时$f(X)$为一定值，$E(f(X))=f(X)$,$f(E(X))=f(X)$

如下如所示：

在这个表达式中，若t相当于$x_1$的概率，(1-t)相当于$x_2$的概率：

$tf(x_1)+(1-t)f(x_2)\ge f(t{x}_1+(1-t)x_2),t\in \left(0,1\right)$

二、EM算法的推导

假设我们有一个样本集
$\left\{x^{(1)},x^{(2)},...,x^{(m)}\right\}$,包含m个独立的样本。但每一个样本i对应的$z^{(i)}$是未知的，也就是隐变量。故我们需要估计概率模型$p(x,z)$的参数$\theta$,但是由于里面包含隐变量z，所以很难用最大似然求解，但是我们要是想办法知道了z，那么就很容易求解了。

对于参数估计，我们本质上还是想获得一个使似然函数最大的$\theta$：

$\mathbf{\theta }=\underset{\mathbf{\theta }}{\mathrm{argmax}}{\sum }_{\mathbf{i}=1}^{\mathbf{m}}\mathbf{l}\mathbf{o}\mathbf{g}\mathbf{P}({\mathbf{x}}^{\mathbf{i}}\mid \mathbf{\theta })$

如果我们得到的观察数据有未观察到的隐含数据$z=(z^{(1)},z^{(2)},...,z^{(m)})$,那么由边缘分布可以写成联合分布之和，我们得到极大化模型分布的对数似然函数如下：

$\mathbf{\theta }=\underset{\mathbf{\theta }}{\mathrm{argmax}}{\sum }_{\mathbf{i}=1}^{\mathbf{m}}\mathbf{l}\mathbf{o}\mathbf{g}\mathbf{P}({\mathbf{x}}^{\mathbf{i}}\mid \mathbf{\theta })=\underset{\mathbf{\theta }}{\mathrm{argmax}}{\sum }_{\mathbf{i}=1}^{\mathbf{m}}\mathbf{l}\mathbf{o}\mathbf{g}{\sum }_{{\mathbf{z}}^{\mathbf{i}}}\mathbf{P}({\mathbf{x}}^{\mathbf{i}},{\mathbf{z}}^{\mathbf{i}}\mid \mathbf{\theta })$

对于每一个样例$x^i$,让$Q_i$表示该样例隐含变量z的某种分布，$Q_i$满足条件${\sum }_z{Q}_i(z)=1,Q_i(z)\ge 0$。（如果z是连续性的，那么$Q_i$是概率密度函数，需要将求和符号换做积分符号）。比如要将班上学生聚类，假设隐藏变量z是身高，那么就是连续的高斯分布。如果按照隐藏变量是男女，那么就是伯努利分布了。

下面的公式阐述了这个思想：

上面公式推导中多了一个未知的变量而已，我也可以分别对未知的θ和z分别求偏导，再令其等于0，求解出来不也一样吗？但是可以看到里面有“和的对数”，求导后形式会非常复杂（自己可以想象下$log(f1(x)+f2(x)+f3(x)+\dots )$复合函数的求导），所以很难求解得到未知参数z和θ,所以在（2）式中分子分母乘以一个相等的函数，但还是有"和的对数"，这里就用到了我们前面讲的$Jensen$不等式，将"和的对数"转变成了"对数的和"，并且等号变成了不等号，因此：

${\sum }_{i=1}^{m}logP(x^i\mid \theta )$就这么巧妙的有了一个下界，并且如果满足$Jensen$不等式的等号，则有：

$\frac{P(x^i,z^i\mid \theta )}{Q_i(z^i)}=c，c为常数$

下面推导下$z^i$的分布$Q_i(z^i)$:

应用到一个简单的小公理：$\frac{1}{2}=\frac{1}{2}+\frac{2}{4}=\frac{1+2}{2+4}=\frac{3}{6}$

$\frac{P(x^i,z^i\mid \theta )}{Q_i(z^i)}=\frac{{\sum }_{z}P(x^i,z^i\mid \theta )}{{\sum }_z{Q}_i(z^i)}=\frac{{\sum }_{z}P(x^i,z^i\mid \theta )}{1}=c$

得：

$c={\sum }_{z}P(x^i,z^i\mid \theta )$

故：

$Q_i(z^i)=\frac{P(x^i,z^i\mid \theta )}{c}=\frac{P(x^i,z^i\mid \theta )}{{\sum }_{z}P(x^i,z^i\mid \theta )}=\frac{P(x^i,z^i\mid \theta )}{P(x^i\mid \theta )}=P(z^i\mid x^i,\theta )$

而且：

${\sum }_z{Q}_i(z)\frac{P(x^i,z^i\mid \theta )}{Q_i(z^i)}$就是$\frac{P(x^i,z^i\mid \theta )}{Q_i(z^i)}$在分布为$Q_i(z^i)$下的期望，这就是我们的E步（求期望）。

上面的关于E步的公式推导，都是建立在$Jensen$不等式取等号的前提下完成的，所以E步求得的期望越大，那么我们要优化的似然估计$logP(x^i\mid \theta )$也就越大。

因此我们的M步就是要求出使得该期望最大的参数$\theta$,因此我们要最大化下式：

$\underset{\mathbf{\theta }}{\mathrm{argmax}}{\sum }_{\mathbf{i}=1}^{\mathbf{m}}{\sum }_{{\mathbf{z}}^{\mathbf{i}}}{\mathbf{Q}}_{\mathbf{i}}({\mathbf{z}}^{\mathbf{i}})\mathbf{l}\mathbf{o}\mathbf{g}\frac{\mathbf{P}({\mathbf{x}}^{\mathbf{i}},{\mathbf{z}}^{\mathbf{i}}\mid \mathbf{\theta })}{{\mathbf{Q}}_{\mathbf{i}}({\mathbf{z}}^{\mathbf{i}})}$

去掉常数项，注意

$Q_i(z^i)=P(z^i\mid x^i,{\theta }^j)$

是第j次迭代后得到的${\theta }^j$,所以为常数，将其去掉后得：

$\underset{\mathbf{\theta }}{\mathrm{argmax}}{\sum }_{\mathbf{i}=1}^{\mathbf{m}}{\sum }_{{\mathbf{z}}^{\mathbf{i}}}{\mathbf{Q}}_{\mathbf{i}}({\mathbf{z}}^{\mathbf{i}})\mathbf{l}\mathbf{o}\mathbf{g}\mathbf{P}({\mathbf{x}}^{\mathbf{i}},{\mathbf{z}}^{\mathbf{i}}\mid \mathbf{\theta })$

这就是M步！

三、EM算法流程

现在我们总结下EM算法的流程！

输入：观察数据$x=\left\{x^{(1)},x^{(2)},...,x^{(m)}\right\}$,联合分布为$P(x,z\mid \theta )$,条件分布为$P(z\mid x,\theta )$,最大迭代次数为$J$。

1. 随机初始化模型参数$\theta$的初值${\theta }^0$.

2. for j from 1 to J开始EM算法迭代：

   1. E步：计算联合分布的条件概率期望：

   $Q_i(z^i)=P(z^i\mid x^i,{\theta }^j)$

   $L(\theta ,{\theta }^j)={\sum }_{i=1}^m{\sum }_{z^i}{Q}_i(z^i)logP(x^i,z^i\mid \theta )$

   2. M步：极大化$L(\theta ,{\theta }^j)$,得到${\theta }^{j+1}$:

   ${\theta }^{j+1}=\underset{\theta }{\mathrm{argmax}}L(\theta ,{\theta }^j)$

   3. 如果${\theta }^j$已收敛，则算法停止。否则继续回到E步，进行迭代。

输出：模型参数$\theta$.

四、EM算法收敛性分析

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法，EM算法的最大优点是简单性和普适性，我们很自然地要问：EM算法得到的估计序列是否收敛？如果收敛，是否收敛到全局最大值或局部极大值？下面证明其收敛性！

假设$P(X\mid \theta )$为观测数据的似然函数，其中$X=\left\{x^{(1)},x^{(2)},...,x^{(m)}\right\}$,${\theta }^i(i=1,2,...)$为EM算法得到的参数估计序列，$P(X\mid {\theta }^j)(i=1,2,...)$为对应的似然函数序列，那么我们需要只需证明下式成立就可以：

$P(X\mid {\theta }^{i+1})\ge P(X\mid {\theta }^i)$

证明：由于

$P(X\mid \theta )=\frac{P(X,Z\mid \theta )}{P(Z\mid X,\theta )}$

取对数得

$logP(X\mid \theta )=logP(X,Z\mid \theta )-logP(Z\mid X,\theta )$

上面公式两端都取在条件概率$P(Z\mid X,{\theta }^j)$分布下的期望

${\sum }_{z}logP(X\mid \theta )P(Z\mid X,{\theta }^j)={\sum }_{z}logP(X,Z\mid \theta )P(Z\mid X,{\theta }^j)-{\sum }_{z}logP(Z\mid X,\theta )P(Z\mid X,{\theta }^j)$

整理得

$logP(X\mid \theta )={\sum }_{z}logP(X,Z\mid \theta )P(Z\mid X,{\theta }^j)-{\sum }_{z}logP(Z\mid X,\theta )P(Z\mid X,{\theta }^j)$

令

$L(\theta ,{\theta }^j)={\sum }_{z}logP(X,Z\mid \theta )P(Z\mid X,{\theta }^j)$

$H(\theta ,{\theta }^j)={\sum }_{z}logP(Z\mid X,\theta )P(Z\mid X,{\theta }^j)$

于是对数似然函数可以写成

$logP(X\mid \theta )=L(\theta ,{\theta }^j)-H(\theta ,{\theta }^j)$

在上式中将$\theta$分别取${\theta }^j$和${\theta }^{j+1}$并相减，有

$logP(X\mid {\theta }^{j+1})-logP(X\mid {\theta }^j)$

$=[L({\theta }^{j+1},{\theta }^j)-L({\theta }^j,{\theta }^j)]-[H({\theta }^{j+1},{\theta }^j]-H({\theta }^j,{\theta }^j)]$

因此，我们只需证明上式的右端非负即可，而上式的第一项是由${\theta }^{j+1}$使$L(\theta ,{\theta }^j)$达到极大，所以有

$L({\theta }^{j+1},{\theta }^j)-L({\theta }^j,{\theta }^j)\ge 0$

其第二项证明如下

$H({\theta }^{j+1},{\theta }^j）-H({\theta }^j,{\theta }^j)$

$={\sum }_{z}logP(Z\mid X,{\theta }^{j+1})P(Z\mid X,{\theta }^j)-{\sum }_{z}logP(Z\mid X,{\theta }^j)P(Z\mid X,{\theta }^j)$

$={\sum }_z\left(log\frac{P(Z\mid X,{\theta }^{j+1})}{P(Z\mid X,{\theta }^j)}\right)P(Z\mid X,{\theta }^j)$

$\le log\left({\sum }_z\frac{P(Z\mid X,{\theta }^{j+1})}{P(Z\mid X,{\theta }^j)}P(Z\mid X,{\theta }^j)\right)$

$=log\left({\sum }_{z}P(Z\mid X,{\theta }^{j+1})\right)=0$

因此

$logP(X\mid {\theta }^{j+1})-logP(X\mid {\theta }^j)\ge 0$

从上面的推导可以看出，EM算法可以保证收敛到一个稳定点，但是却不能保证收敛到全局的极大值点，因此它是局部最优的算法，当然，如果我们的优化目标$L(\theta ,{\theta }^j)$是凸的，则EM算法可以保证收敛到全局最大值，这点和梯度下降法这样的迭代算法相同。

五、EM算法总结

监督学习是由训练数据 { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x N , y N ) } \left \{ (x_1,y_1),(x_2,y_2),...,(x_N,y_N) \right \} {(x1​,y1​),(x2​,y2​),...,(xN​,yN​)}学习条件概率分布$P(Y\mid X)$或决策函数$Y=f(X)$作为模型，用于分类、回归、标注等任务，这时训练数据中的每个样本点由输入和输出对组成。

但有时训练数据只有输入没有对应的输出 { ( x 1 , ∗ ) , ( x 2 , ∗ ) , . . . , ( x N , ∗ ) } \left \{ (x_1,*),(x_2,*),...,(x_N,*) \right \} {(x1​,∗),(x2​,∗),...,(xN​,∗)}，从这样的数据学习模型称为非监督学习问题。EM算法可以用于生成模型的非监督学习。生成模型由联合概率分布$P(X,Y)$表示，可用认为非监督学习训练数据上联合概率分布产生的数据。X为观测数据，Y为未观测数据（隐变量）。

另外，在应用中，初值的选择变得非常重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。