09_期望极大法EM1_统计学习方法

EM内容较多，方便阅读，分成2个部分
下接：09_期望极大法EM2_统计学习方法

EM算法是一种迭代算法，1977年由Dempster等人总结提出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（exception）；M步，求极大（maximization）。所以这一算法称为期望极大算法，简称EM算法。

一、EM算法的引入

概率模型有时既含有观测变量，又含有隐变量或潜在变量。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数；当模型含有隐变量时，就不能简单地使用这些估计方法，EM算法就是含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计法。下文以讨论极大似然估计为例，极大后验概率估计与其类似。

在朴素贝叶斯算法里面已经讲解过极大似然估计和贝叶斯估计，今天再对极大似然估计做一个详细的介绍，因为实在是太常见了。

1、极大似然估计

总体样本服从$P(X|\theta )$分布，现有$N$个样本，估计模型的参数$\theta$。那么得到的参数$\theta$应该是使现有的$N$个样本最容易出现的参数，这就是极大似然估计的思想（模型已定，参数未知）。每个样本之间是相互独立的，那么可以得到这$N$个样本的联合概率即极大似然函数，用下式表示：
$$L(\theta )=L(x_1,\cdots ,x_n|\theta )=\prod _{i=1}^{N}p(x_i|\theta ),\theta \in \Theta \text{\hspace{1em}(1)}$$
对数似然函数：
$$l(\theta )=lnL(\theta )=\sum _{i=1}^{N}lnp(x_i|\theta ),\theta \in \Theta \text{\hspace{1em}(2)}$$
求极大似然函数估计值的一般步骤：

• （1）写出似然函数；
• （2）对似然函数取对数，并整理；
• （3）求导数，令导数为0，得到似然方程；
• （4）解似然方程，得到的参数即为所求；

简单实例说明，实例来自博客。假设现在有两枚硬币1和2，,随机抛掷后正面朝上概率分别为$P_1，P_2$。为了估计这两个概率，做实验，每次取一枚硬币，连掷5下，记录下结果，如下：

可以很容易地估计出$P_1$和$P_2$，如下：

$$\begin{array}{rl}& P1=（3+1+2）/15=0.4\\ & P2=（2+3）/10=0.5\end{array}$$
上面的计算很简单很理所当然，但上面是用了极大似然估计的。可以见证明，极大似然估计证明先验概率

2、EM入场

还是上面的问题，现在我们抹去每轮投掷时使用的硬币标记，如下：

问题还是要估计两枚硬币随机抛掷后正面朝上概率$P_1，P_2$。

显然，此时我们多了一个隐变量$Z$，可以把它认为是一个5维的向量$（z_1,z_2,z_3,z_4,z_5)$，代表每次投掷时所使用的硬币，比如$z_1$，就代表第一轮投掷时使用的硬币是1还是2。但是，这个变量$Z$不知道，就无法去估计$P_1$和$P_2$，所以，我们必须先估计出$Z$，然后才能进一步估计$P_1$和$P_2$。
但要估计$Z$，我们又得知道$P_1$和$P_2$，这样我们才能用最大似然概率法则去估计$Z$，这不是鸡生蛋和蛋生鸡的问题吗？如何破？

答案就是先随机初始化一个$P_1$和$P_2$，用它来估计$Z$，然后基于$Z$，按照最大似然概率法则去估计新的$P_1$和$P_2$。这就是最开始说的E步$(P_1,P_2⟶Z)$和M步$(Z⟶P_1,P_2)$，如此反复下去，直到收敛到一个不再改变的参数$\theta$。

二、EM算法推导

1、Jensen不等式

在支持向量机里面介绍过凸函数。今天讲下Jensen不等式。

看到很多博客里面有不同的评论说凹凸函数用混了，这边为了不引起歧义在使用的时候直接带上条件，那么肯定就不会弄混了。

如果$f(x)$是凸函数($f^{\prime \prime }(x)\ge 0$)，那么有：
$$\frac{f(x_1)+f(x_2)+\cdots +f(x_n)}{n}\ge f(\frac{x_1+x_2+\cdots +x_n}{n})$$
当且仅当$x_1=x_2=\cdots =x_n$时上式取等号。

加权形式：
$$\sum _{i=1}^N{\alpha }_{i}f(x_i)\ge f(\sum _{i=1}^N{\alpha }_i{x}_i),\sum _{i=1}^N\alpha =1\text{\hspace{1em}(3)}$$
Jensen不等式证明

推得若$f(x)$是凸函数，X是随机变量，有
$$E[f(X)]\ge f(E[X])$$
当且仅当X是常量时，上式取等号。

如果$f(x)$是凹函数($f^{\prime \prime }(x)\le 0$)，那么同理有：
$$E[f(X)]\le f(E[X])\text{\hspace{1em}(4)}$$

考虑到$E(X)={\sum }_{X}p(x)\times x,x\in X$，则有$E(f(X))={\sum }_{X}p(x)\times f(x),x\in X$，$f(X)$是$X$的函数。

2、EM推导过程

（1）统计学习方法EM推导

第一部分内容，看看文字描述，能理解EM算法大概是怎么回事，是在解决一个什么问题。但是具体要使用，就需要数学公式推导，证明这个想法确实可行。这一部分内容看课本的时候很难理解。然后找到一篇讲解和带部分公式推导还可以的博客。但是我决定最终用来总结推导的内容还是以参考课本为主。对其中理解起来较困难的地方展开和进一步解释。

一般地，用$Z$表示隐随机变量的数据，$Y$表示观测随机变量的数据又称为不完全数据，$Y$和$Z$连在一起称为完全数据。假设给定观测数据$Y$，其概率分布是$P(Y|\theta )$，其中$\theta$是需要估计的模型参数，那么不完全数据$Y$的似然函数是$P(Y|\theta )$，对数似然$l(\theta )=lnP(Y|\theta )$；假设$Y$和$Z$的联合概率分布是$P(Y,Z|\theta )$，那么完全数据的对数似然函数是$lnP(Y,Z|\theta )$。

面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据）$Y$关于参数$\theta$的对数似然函数，即极大化
$$\begin{array}{rl}l(\theta )& =lnP(Y|\theta )=ln\sum _{Z}P(Y,Z|\theta )\\ & =ln\left(\sum _{Z}P(Y|Z,\theta )P(Z|\theta )\right)\end{array}\text{\hspace{1em}(5)}$$
注意这一极大化的主要困难是式（5）中有未观测数据并有包含和的对数。事实上，EM算法是通过迭代逐步近似极大化$l(\theta )$的。假设在第$i$次迭代后$\theta$的估值是${\theta }^{(i)}$。我们希望新估计值$\theta$能使$l(\theta )$增加，即$l(\theta )\ge l({\theta }^{(i)})$，并逐步达到极大值。
$$\begin{array}{rl}l(\theta )-l({\theta }^{(i)})& =ln\left(\sum _{Z}P(Y|Z,\theta )P(Z|\theta )\right)-lnP(Y|{\theta }^{(i)})\\ & =ln\left(\sum _{Z}P(Y|Z,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(i)})}\right)-lnP(Y|{\theta }^{(i)})\end{array}$$

$(lnx{)}^{\prime \prime }=-\frac{1}{x^2}\le 0$，为凹函数，${\sum }_{Z}P(Y|Z,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(i)})}$是$\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(i)})}$的期望，利用Jensen不等式可得凹函数有$f(E[X])\ge E[f(X)]$，则得到上式下界：
$$\begin{array}{rlr}l(\theta )-l({\theta }^{(i)})=& ln\left(\sum _{Z}P(Y|Z,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(i)})}\right)-lnP(Y|{\theta }^{(i)})& (5.1)\\ & \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})ln\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-lnP(Y|{\theta }^{(i)})& \\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})ln\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}& \end{array}$$

为什么这么凑Jensen不等式？从后面Q函数的意义可以知道为什么。

令
$$B(\theta ,{\theta }^{(i)})\widehat{=}l({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})ln\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\text{\hspace{1em}(6)}$$
则
$$l(\theta )\ge B(\theta ,{\theta }^{(i)})\text{\hspace{1em}(7)}$$
即函数$B(\theta ,{\theta }^{(i)})$是$l(\theta )$的一个下界，而且式（6）可知
$$l({\theta }^{(i)})=B({\theta }^{(i)},{\theta }^{(i)})\text{\hspace{1em}(8)}$$
因此，任何可以使$B(\theta ,{\theta }^{(i)})$增大的$\theta$，也可以使$l(\theta )$增大。为了使$l(\theta )$有尽可能的增大，选择${\theta }^{(i+1)}$使$B(\theta ,{\theta }^{(i)})$达到极大，即
$${\theta }^{(i+1)}=arg\underset{\theta }{\max }B(\theta ,{\theta }^{(i)})\text{\hspace{1em}(9)}$$
现在求${\theta }^{(i+1)}$的表达式。省去对$\theta$的极大化而言是常数的项，由式（6）、式（9），有
$$\begin{array}{rl}{\theta }^{(i+1)}& =arg\underset{\theta }{\max }\left(l({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})ln\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\right)\\ & =arg\underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})lnP(Y|Z,\theta )P(Z|\theta )\right)\\ & =arg\underset{\theta }{\max }\left(\sum _{Z}P(Z|Y,{\theta }^{(i)})lnP(Y,Z|\theta )\right)\\ & =arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\end{array}\text{\hspace{1em}(10)}$$
式（10）等价于EM算法的一次迭代，即求$Q$函数及其极大化。EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

（2）Andrew NG关于EM算法的推导

这篇博客里面有介绍Andrew NG关于EM算法的推导，和李航的统计学习方法不一样，可以相互参考着看，异曲同工。不用特意知道为什么那么去凑Jensen不等式。
$$\begin{array}{rl}l(\theta )& =lnP(Y|\theta )=ln\sum _{Z}P(Y,Z|\theta )\\ & =\sum _{j=1}^{N}ln\sum _{i}P(y_j,z_i|\theta )\\ & =\sum _{j=1}^{N}ln\sum _{i}Q\frac{P(y_j,z_i|\theta )}{Q}\\ & \ge \sum _{j=1}^N\sum _{i}Qln\frac{P(y_j,z_i|\theta )}{Q}\end{array}\text{\hspace{1em}(10.1)}$$
上面$\theta$极大化求解过程中含有隐变量参数以及在对数里面有求和操作，求解会非常复杂，所以上式的变形过程中使用了Jensen不等式，$lnx$满足凹函数的条件，那么只需要${\sum }_{i}Q=1$上式变形就是成立。另外要使对数似然函数$l(\theta )$最大，那么需要在上式的变形中取等号，根据Jensen不等式有$P(y_j,z_i|\theta )/Q=c,c$是常数的时候等号成立。则有

$$\begin{array}{rl}\sum _{i}P(y_j,z_i|\theta )/c=1& ⟹\sum _{i}P(y_j,z_i|\theta )=c\\ & ⟹Q=\frac{P(y_j,z_i|\theta )}{{\sum }_{i}P(y_j,z_i|\theta )}=\frac{P(y_j,z_i|\theta )}{P(y_j|\theta )}=P(z_i|y_j,\theta )=Q_j(z_i)\end{array}\text{\hspace{1em}(10.2)}$$
其中$Q_j(z_i)$为$z_i$的后验概率，求解$Q_j(z_i)$过程称为E步，即固定参数$\theta$，求出使得Jensen变形中取等号的$Q_j(z_i)$。

所以有
$$l(\theta )\ge \sum _{j=1}^N\sum _{i}P(z_i|y_j,\theta )ln\frac{P(y_j,z_i|\theta )}{P(z_i|y_j,\theta )}\text{\hspace{1em}(10.3)}$$
M步：算出$Q_j(z_i)$也就是有了隐变量的后验概率后，就可以通过极大似然估计进行求解最优的参数$\theta$作为${\theta }^{(i+1)}$进行下一轮迭代，直到迭代停止条件。

对比式（10）和式（10.3）可以发现统计学习方法里面的$Q$函数实际上是$l(\theta )$的一个变形，和Andrew NG推导里面的$l(\theta )$形式几乎是一样的，只是因为统计学习方法在M步求解极大化时把常数项参数省去了。两个推导中的$Q$函数意义不一样，只是都用了同一个名字。

按上面方式讲解EM算法一篇非常不错的博客（EM算法）The EM Algorithm。

3、统计学习方法EM算法流程

输入：观测变量数据$Y$，隐变量数据$Z$，联合分布$P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$；

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

（1）选择参数的初值${\theta }^{(0)}$，开始迭代；初值可以任意选择，但需注意EM算法对初值是敏感的。
（2）E步：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步，计算
$$Q(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})lnP(Y,Z|\theta )=E_Z[lnP(Y,Z|\theta )|Y,{\theta }^{(i)}]\text{\hspace{1em}(11)}$$

• 上式的函数$Q(\theta ,{\theta }^{(i)})$是EM算法的核心，称为$Q$函数。
• 需要计算$P(Z|Y,{\theta }^{(i)})$，即在给定观测数据$Y$和当前的参数估计${\theta }^{(i)}$下隐变量数据$Z$的条件概率分布，作为隐变量的现估计值。
• $Q$函数：完全数据的对数似然函数$lnP(Y,Z|\theta )$关于在给定观测数据$Y$和当前参数${\theta }^{(i)}$下对未观测数据$Z$的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望称为$Q$函数。

（3）M步：求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$次迭代的参数的估计值${\theta }^{(i+1)}$
$${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(12)}$$
（4）重复第（2）步和第（3）步，直到收敛。给出停止迭代的条件，一般是对较小的正数${ϵ}_1,{ϵ}_2$，若满足
$$||{\theta }^{(i+1)}-{\theta }^{(i)}||<{ϵ}_{1}or||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||<{ϵ}_2\text{\hspace{1em}(13)}$$

4、EM算法直观解释

在下图中给出了EM算法的直观解释。图中上方曲线为$l(\theta )$，下方曲线为$B(\theta ,{\theta }^{(i)})$。$B(\theta ,{\theta }^{(i)})$为对数似然函数$l(\theta )$的下界。两个函数在点$\theta ={\theta }^{(i)}$处相等。由式（9）和（10），EM算法找到下一个点${\theta }^{(i+1)}$使函数$B(\theta ,{\theta }^{(i)})$极大化，也使函数$Q(\theta ,{\theta }^{(i)})$极大化。这时由于$l(\theta )\ge B(\theta ,{\theta }^{(i)})$，函数$B(\theta ,{\theta }^{(i)})$的增加，保证对数似然函数$l(\theta )$在每次迭代中也是增加的。EM算法在点${\theta }^{(i+1)}$重新计算$Q$函数值，进行下一次迭代。在这个过程中，对数似然函数$l(\theta )$不断增大。从图中可以推断出EM算法不能保证找到全局最优值。

5、EM算法在非监督学习中的应用

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

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

三、EM算法的收敛性

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。EM算法的最大优点是简单性和普适性。对EM算法的两个疑问：

• 1、EM算法得到的估计序列是否收敛？
• 2、如果收敛，是否收敛到全局最大值或者局部最大值？

EM算法定理1：设$P(Y|\theta )$为观测数据的似然函数，${\theta }^{(i)}(i=1,2,\cdots )$为EM算法得到的参数估计序列，$P(Y|{\theta }^{(i)})(i=1,2,\cdots )$为对应的似然函数序列，则$P(Y|{\theta }^{(i)})$是单调递增的，即
$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$$

EM算法定理2：设$l(\theta )=lnP(Y|\theta )$为观测数据的对数似然函数，${\theta }^{(i)}(i=1,2,\cdots )$为EM算法得到的参数估计序列，$l({\theta }^{(i)})(i=1,2,\cdots )$为对应的对数似然函数序列。

• （1）如果$P(Y|\theta )$有上界，则$l({\theta }^{(i)})=lnP(Y|{\theta }^{(i)})$收敛到某一值$l^{\ast }$；
• （2）在函数$Q(\theta ,{\theta }^{\prime })$与$l(\theta )$满足一定条件下，由EM算法得到的参数估计序列${\theta }^{(i)}$的收敛值${\theta }^{\ast }$是$l(\theta )$的稳定点。

定理1$(P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)}))$的证明：

由于
$$P(Y|\theta )=\frac{P(Y,Z|\theta )}{P(Z|Y,\theta )}$$
取对数有
$$lnP(Y|\theta )=lnP(Y,Z|\theta )-lnP(Z|Y,\theta )\text{\hspace{1em}(14)}$$
根据式（11）有$Q$函数
$$Q(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})lnP(Y,Z|\theta )$$
令
$$H(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})lnP(Z|Y,\theta )\text{\hspace{1em}(15)}$$
于是对数似然函数可以写成
$$lnP(Y|\theta )=Q(\theta ,{\theta }^{(i)})-H(\theta ,{\theta }^{(i)})\text{\hspace{1em}(16)}$$
在式（16）中分别取$\theta$为${\theta }^{(i+1)}$和${\theta }^{(i)}$并相减，有
$$lnP(Y|{\theta }^{(i+1)})-lnP(Y|{\theta }^{(i)})=[Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})]-[H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})]\text{\hspace{1em}(17)}$$
由于${\theta }^{(i+1)}$使$Q(\theta ,{\theta }^{(i)})$达到极大，所以有
$$Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})\ge 0\text{\hspace{1em}(18)}$$
对于式（17）等式右边第2项
$$\begin{array}{rl}H({\theta }^{(i+1)},& {\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})\\ & =\sum _Z\left(ln\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}\right)P(Z|Y,{\theta }^{(i)})\\ & \le ln\left(\sum _Z\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}P(Z|Y,{\theta }^{(i)})\right)\\ & =ln\sum _{Z}P(Z|Y,{\theta }^{(i+1)})=0\end{array}\text{\hspace{1em}(19)}$$
由式（18）和式（19）可以证明$lnP(Y|{\theta }^{(i+1)})-lnP(Y|{\theta }^{(i)})\ge 0$。证毕。

定理2证明略，EM算法的收敛性包含关于对数似然序列$l({\theta }^{(i)})$的收敛性和关于参数估计序列${\theta }^{(i)}$的收敛性两层意思，前者并不蕴涵后者。定理只能保证参数估计序列收敛到对数似然函数序列的稳定点，不能保证收敛到极大值点。所以在应用中，初值的选择变得很重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

四、高斯混合模型GMM

见 09_期望极大法EM2_统计学习方法

五、EM算法的推广

见 09_期望极大法EM2_统计学习方法

参考资料：

• 统计学习方法——李航
• 如何感性地理解EM算法？
• 从最大似然到EM算法浅解
• EM算法学习
• 高斯混合模型（GMM）