学习统计学习方法-EM算法及其推广

  最近在论文阅读中看到了EM算法，于是翻开了我的《统计学习方法》。之前看第一版的时候看到过EM算法，当时只是大致学习了一下，并没有仔细深入，现在细节忘了许多，回顾一章。
参考博客:统计学习笔记九----EM算法
很仔细的文章:详细解读EM算法
预备知识：极大似然估计

---

---

1.EM算法简介

  EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大后验概率估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation)；M步，求极大(maximization)，故算法全称期望极大算法(expectation maximization algorithm)。
  概率模型通常情况下含有观测变量，隐变量或潜在变量(latent variable)。当模型只有观测变量，没有隐变量时，可以直接用极大似然估计法或者贝叶斯估计法估计模型参数。当模型含有隐变量的时候，就不能简单地使用这些估计法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或者极大后验概率估计法。例如常被用来学习高斯混合模型（Gaussian mixture model，简称GMM）的参数；隐式马尔科夫算法（HMM）、LDA主题模型的变分推断等等。
:博客1已经很好的阐述了书中的三枚硬币的例子。

2.EM算法细节

  由于EM是机器学习十大算法之一，所以现有的资料博客对EM算法讲解已经很细节。所以本文只是为了将我的个人理解记录下来，如有错误，希望指正。首先EM算法从钢镚儿正反故事开始说起。

算法例子

  Nature Biotech在他的论文《Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm?. Nature biotechnology, 26(8), 897.》中用了一个很简单的例子解释了EM算法。
  有两枚硬币，分别为A和B。可能分布不均匀，所以投掷的时候正反面不一定相等。投掷这两枚硬币，开始的时候咱们知道投掷的是哪枚硬币，每枚硬币投掷一轮，每轮10次，共投掷了5 epochs。如图a所示

硬币	结果	统计
A	正 正 反 正 反	3正-2反
B	反 反 正 正 反	2正-3反
A	正 反 反 反 反	1正-4反
B	正 反 反 正 正	3正-2反
A	反 正 正 反 反	2正-3反

根据上表投掷硬币的结果：
A的估计概率为(第1，3，5轮)：$P(A)=\frac{(3+1+2)}{15}=0.4$；
B的估计概率为(第2，4轮)：$P(B)=\frac{(2+3)}{10}=0.5$。
如果不知道投掷的是哪枚硬币呢？还是上个例子:

硬币	结果	统计
Unknow	正 正 反 正 反	3正-2反
Unknow	反 反 正 正 反	2正-3反
Unknow	正 反 反 反 反	1正-4反
Unknow	正 反 反 正 正	3正-2反
Unknow	反 正 正 反 反	2正-3反

如下图所示：

  其中Unknown为隐含变量$Z=(z_1,z_2,z_3,z_4,z_5)$,其中$z_i$代表第i个硬币是A或B。由于即不知道$Z$,也不知道$P(A),P(B)$。所以陷入循环，不知道怎么入手，所以需要先初始化一手，即首先假设$P(A),P(B)$的值，然后最大似然概率估计Z，之后通过Z按照最大似然法则估计$P(A)$和$P(B)$。
  初始化$P(A)$和$P(B)$。不妨设

• A硬币正面朝上的概率：$P(A)=0.2$
• B硬币正面朝上的概率：$P(B)=0.7$
  按照$P(A)$和$P(B)$，极大似然估计Z

轮数	结果	统计	硬币 A	硬币B
1	正 正 反 正 反	3正-2反	$0.2\times 0.2\times 0.8\times 0.2\times 0.8=0.00512$	$0.7\times 0.7\times 0.3\times 0.7\times 0.3=0.03087$
2	反 反 正 正 反	2正-3反	$0.8\times 0.8\times 0.2\times 0.2\times 0.8=0.02048$	$0.3\times 0.3\times 0.7\times 0.7\times 0.3=0.01323$
3	正 反 反 反 反	1正-4反	$0.2\times 0.8\times 0.8\times 0.8\times 0.8=0.08192$	$0.7\times 0.3\times 0.3\times 0.3\times 0.3=0.00567$
4	正 反 反 正 正	3正-2反	$0.2\times 0.8\times 0.8\times 0.2\times 0.2=0.00512$	$0.7\times 0.3\times 0.3\times 0.7\times 0.7=0.03087$
5	反 正 正 反 反	2正-3反	$0.8\times 0.2\times 0.2\times 0.8\times 0.8=0.02048$	$0.3\times 0.7\times 0.7\times 0.3\times 0.3=0.01323$

按照最大似然法则可知：
第1轮中最有可能的是硬币B
第2轮中最有可能的是硬币A
第3轮中最有可能的是硬币A
第4轮中最有可能的是硬币B
第5轮中最有可能的是硬币A
  由此可知，第2，3，5轮更可能是A，一共被抛了15次，共2+1+2次正面，所以按照最大似然概率法则估计新的$P(A)=\frac{2+1+2}{15}=0.33$；同理，$P(B)=\frac{3+3}{10}=0.6$。对比已知正反面的实验来看

初始化P(A)	估计出的P(A)	真实P(A)	初始化P(B)	估计出的P(B)	真实P(B)
0.2	0.33	0.4	0.7	0.6	0.5

  可以看出估计的$P(A)$和$P(B)$距离真实数据更近，之后不断迭代可以得到最优解-即初始化的$P(A)$和$P(B)$和估计的$P(A)$和$P(B)$相等。
  通过以上的例子可以对EM算法有了大致的了解：
E步：估计隐变量Z的过程
M步：由估计的Z值，采用最大似然概率法则估计新的$P(A)$和$P(B)$。
用新得到的$P(A)$和$P(B)$估计新的隐变量Z。
循环迭代，直到$P(A)$和$P(B)$不再变化为止。

推导过程

  通过抛硬币，我们大致理解了EM算法。接下来进行公式计算。从总体中抽取n个样本，形成的样本集为X，样本集分布为$p(x|\theta )$，联合概率为：$$L(\theta )=L(x_1,x_2,x_n;\theta )=\prod _{i=1}^{n}p(x_i;\theta ),\theta \in \Theta$$极大似然估计为：$$\begin{array}{rl}\ell (\theta )& =\sum _{i=1}^m\log p(x;\theta )\\ & =\sum _{i=1}^m\log \sum _{z}p(x,z;\theta )\end{array}$$z为隐藏变量，不易求出，所以需要求解z。$$\begin{array}{rl}\sum _i\log p(x^{(i)};\theta )& =\sum _i\log \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )(1)\\ & =\sum _i\log \sum _{z^{(i)}}{Q}_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}(2)\\ & \ge \sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}(3)\end{array}$$根据极大似然估计，使1式最大就可，对$\theta$求导找梯度为0的点，但是求导太过于复杂，对数的和求导，很难得到z和$\theta$。
  这时就该想，如果把对数求和变成对数乘积就可以解决了，所以目标就是对数的乘积。分子分母同时乘以隐变量Z的分布函数$Q_i(z^{(i)})$，其概率之和为1.即${\sum }_z{Q}_i(z^{(i)})=1$，即可以得到上式2，但是还是有和的对数，无法求解，所以采用Jensen不等式得到3式。
  先说一下为什么要2式到3式目的：2式不好计算，所以需要找到2式的下界，用下界的最大值来拟合2式，这样就可以间接解决了2式的极大值。3式就是这个作用；
  怎么求3式呢？这里引入Jensen不等式。

Jensen不等式

设f是定义域为实数的函数
1. 如果对于所有的实数x，f(x)的二次导数$f^{\prime \prime }(x)\ge 0$，那么f是凸函数。
2. 当x是向量时，如果其hessian矩阵H是半正定的$H\ge 0$，那么f是凸函数。
3. 如果$f^{\prime \prime }(x)\ge 0$或者$H\ge 0$，那么称f是严格凸函数，如下图所示。
这里有一个误区，区别凸函数的理解，凹凸性正好相反。

Jensen不等式性质
  如果f是凸函数，X是随机变量，那么：$E[f(X)]\ge f(E[X])$，通俗的说法是函数的期望大于等于期望的函数。
  特别地，如果f是严格凸函数，当且仅当$P(X=EX)=1$，即X是常量时，上式取等号，即$E[f(X)]=f(E[X])$。
  最简单的例子，上图灰色实线f是凸函数，即 $f^{\prime \prime }(x)\ge 0$，X是随机变量，有0.5的概率是a，0.5的概率是b。X的期望就是a,b的均值，通过上图可知$E[f(X)]\ge f(E[X])$。相反，当f是严格凹函数时，不等号相反$E[f(X)]\le f(E[X])$。
  回到式2，$f(x)=\log x$是凹函数($f^{\prime }(x)=-\frac{1}{x^2}<0$)
所以可以得到3式。 $$\begin{array}{rl}\sum _i\log p(x^{(i)};\theta )& =\sum _i\log \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )(1)\\ & =\sum _i\log \sum _{z^{(i)}}{Q}_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}(2)\\ & \ge \sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}(3)\end{array}$$
其中2式的${\sum }_{z^{(i)}}{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)})}$的期望。
原因： 设Y是随机变量X的函数(g是连续函数),$Y=g(X)$,那么
（1） X是离散型随机变量，它的分布律为$P(X=x_k)=p_k$，k=1,2,…。若绝对收敛，则有${\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$
考虑到$E(X)=\sum x\times p(x)$,f(X)是X的函数，则$E(f(X))=\sum f(x)\times p(x)$，又${\sum }_z{Q}_i(z^{(i)})=1$，所以就可以得到公式3的不等式了：$$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)})})]$$这样就将对数的加法变为对数的乘法，就可以求导分解了。但是有个问题，式2和式3是不等号，式2的最大值不是式3的最大值，目的是求取式2的最大值，似然函数$L(\theta )\ge J(z,Q)$，我们可以通过不断的最大化这个下界J，来使得L(θ)不断提高，最终达到L(θ)的最大值。可以这样：

如上图所示，我们固定$\theta$，调整$Q(z)$使下界$J(z,Q)$上升至与$L(\theta )$在此点$\theta$处相等（绿色曲线到蓝色曲线），然后固定Q(z)，调整θ使下界$J(z,Q)$达到最大值（${\theta }^t$到${\theta }^{t+1}$），然后再固定$\theta$，调整$Q(z)$……直到收敛到似然函数$L(\theta )$的最大值处的$\theta \ast$
  根据Jensen不等式，当自变量X是常数的时候，等式成立。对应2式的$X=\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}=c$,c为常数，与$z^{(i)}$无关。所以上式可以变换为：$$\sum _{z}p(x^{(i)},z^{(i)};\theta )=\sum _z{Q}_i(z^{(i)})c$$由于${\sum }_z{Q}_i(z^{(i)})=1$,所以${\sum }_{z}p(x^{(i)},z^{(i)};\theta )=c$$$\begin{array}{rl}{Q}_i(z^{(i)})& =\frac{p(x^{(i)},z^{(i)};\theta )}{c}\\ & =\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{array}$$
  在固定参数θ后，使下界拉升的Q(z)的计算公式就是条件概率，解决了Q(z)如何选择的问题。这一步就是E步，建立L(θ)的下界。
  接下来的M步，就是在给定$Q(z)$后,调整$\theta$,去极大化$L(\theta )$的下界$J(z,Q)$,毕竟在固定$Q(z)$后，下界还可以调整的更大。

EM算法总结

   期望最大EM算法是一种从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。
   当我们不知道隐变量z的具体取值时（比如是硬币A还是硬币B），也就不好判断硬币A或硬币B正面朝上的概率$\theta$，既然这样，那就首先随机初始化分布参数$\theta$，之后循环重复直到收敛
（E步，求Q函数）对于每一个$i$，计算根据上一次迭代的模型参数来计算出隐性变量的后验概率（其实就是隐性变量的期望），来作为隐藏变量的现估计值：$Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta )$
(M步，求使Q函数获得极大时的参数取值)将似然函数最大化以获得新的参数值 $$\theta =\underset{\theta }{\mathrm{argmax}}\sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$$就这样，$Q(z)$求出来代入到$\theta$，$\theta$求出来又反代回$Q(z)$，如此不断的迭代，就可以得到使似然函数$L(\theta )$最大化的参数$\theta$了。
收敛特性
   假定和是EM第 t 次和 t+1 次迭代后的结果（上标 t 表示第 t 次迭代）。如果我们证明了，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。
E步：$Q_i^{(t)}(z^{(i)})=p(z^{(i)}|x^{(i)};{\theta }^{(t)})$
根据Jensen不等式中等式成立，则$$\ell ({\theta }^{(t)})=\sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t)})}{Q_i(z^{(i)})}$$然后进行M步，固定$Q_i^{(t)}(z^{(i)})$，并将${\theta }^{(t)}$作为变量，对$\ell ({\theta }^{(t)})$求导，得到${\theta }^{(t+1)}$。$$\begin{array}{rl}\ell ({\theta }^{(t+1)})& \ge \sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t+1)})}{Q_i(z^{(i)})}\\ & \ge \sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t)})}{Q_i(z^{(i)})}\\ & =\ell ({\theta }^{(t)})\end{array}$$