【统计学习方法】第9章 EM算法及其推广

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

1、EM算法的引入

EM算法

不完全数据：观测随机变量$Y$。
完全数据：观测随机变量$Y$和隐随机变量$Z$。

$Q$函数：完全数据的对数似然函数$\log P\left(Y,Z|\theta \right)$关于在给定观测数据$Y$和当前参数${\theta }^{\left(i\right)}$下对未观测数据$Z$的条件概率分布$P\left(Z|Y,{\theta }^{\left(i\right)}\right)$的期望
$$\begin{array}{rl}& Q\left(\theta ,{\theta }^{\left(i\right)}\right)=E_Z\left[\log P\left(Y,Z|\theta \right)|Y,{\theta }^{\left(i\right)}\right]\end{array}$$

EM算法的导出

含有隐变量$Z$的概率模型，目标是极大化观测变量$Y$关于参数$\theta$的对数似然函数，即$$\begin{array}{rl}& \max L\left(\theta \right)=\log P\left(Y|\theta \right)\\ & =\log \sum _{Z}P\left(Y,Z|\theta \right)\\ & =\log \left(\sum _{Z}P\left(Y|Z,\theta \right)P\left(Z|\theta \right)\right)\end{array}$$

对数似然函数$L\left(\theta \right)$与第$i$次迭代后的对数似然函数估计值$L\left({\theta }^{\left(i\right)}\right)$的差$$\begin{array}{rl}& L\left(\theta \right)-L\left({\theta }^{\left(i\right)}\right)＝\log \left(\sum _{Z}P\left(Y|Z,\theta \right)P\left(Z|\theta \right)\right)-\log P\left(Y|{\theta }^{\left(i\right)}\right)\\ & =\log \left(\sum _{Z}P\left(Z|Y,{\theta }^{\left(i\right)}\right)\frac{P\left(Y|Z,\theta \right)P\left(Z|\theta \right)}{P\left(Z|Y,{\theta }^{\left(i\right)}\right)}\right)-\log P\left(Y|{\theta }^{\left(i\right)}\right)\\ & \ge \sum _{Z}P\left(Z|Y,{\theta }^{\left(i\right)}\right)\log \frac{P\left(Y|Z,\theta \right)P\left(Z|\theta \right)}{P\left(Z|Y,{\theta }^{\left(i\right)}\right)}-\log P\left(Y|{\theta }^{\left(i\right)}\right)\\ & =\sum _{Z}P\left(Z|Y,{\theta }^{\left(i\right)}\right)\log \frac{P\left(Y|Z,\theta \right)P\left(Z|\theta \right)}{P\left(Z|Y,{\theta }^{\left(i\right)}\right)P\left(Y|{\theta }^{\left(i\right)}\right)}\end{array}$$

令$$\begin{array}{rl}& B\left(\theta ,{\theta }^{\left(i\right)}\right)=L\left({\theta }^{\left(i\right)}\right)+\sum _{Z}P\left(Z|Y,{\theta }^{\left(i\right)}\right)\log \frac{P\left(Y|Z,\theta \right)P\left(Z|\theta \right)}{P\left(Z|Y,{\theta }^{\left(i\right)}\right)P\left(Y|{\theta }^{\left(i\right)}\right)}\end{array}$$

则$$\begin{array}{rl}& L\left(\theta \right)\ge B\left(\theta ,{\theta }^{\left(i\right)}\right)\end{array}$$

即函$B\left(\theta ,{\theta }^{\left(i\right)}\right)$ 是$L\left(\theta \right)$ 的一个下界。选择${\theta }^{\left(i\right)}$使$B\left(\theta ,{\theta }^{\left(i\right)}\right)$达到极大，即$$\begin{array}{rl}& {\theta }^{\left(i+1\right)}＝\arg \max B\left(\theta ,{\theta }^{\left(i\right)}\right)\\ & =\arg \max \left(L\left({\theta }^{\left(i\right)}\right)+\sum _{Z}P\left(Z|Y,{\theta }^{\left(i\right)}\right)\log \frac{P\left(Y|Z,\theta \right)P\left(Z|\theta \right)}{P\left(Z|Y,{\theta }^{\left(i\right)}\right)P\left(Y|{\theta }^{\left(i\right)}\right)}\right)\\ & =\arg \max \left(\sum _{Z}P\left(Z|Y,{\theta }^{\left(i\right)}\right)\log \left(P\left(Y|Z,\theta \right)\right)P\left(Z|\theta \right)\right)\\ & =\arg \max \left(\sum _{Z}P\left(Z|Y,{\theta }^{\left(i\right)}\right)\log P\left(Y,Z|\theta \right)\right)\\ & =\arg \max Q\left(\theta ,{\theta }^{\left(i\right)}\right)\end{array}$$

EM算法：

• 输入：观测随机变量数据$Y$，隐随机变量数据$Z$，联合分布$P\left(Y,Z|\theta \right)$，条件分布$P\left(Y｜Z，\theta \right)$；
• 输出：模型参数$\theta$

1. 初值${\theta }^{\left(0\right)}$

2. $E$步：$$\begin{gathered} Q\left(\theta ,{\theta }^{\left(i\right)}\right)=E_Z\left[\log P\left(Y,Z|\theta \right)|Y,{\theta }^{\left(i\right)}\right] \\ =\sum _Z\log P\left(Y,Z|\theta \right)\cdot P\left(Z|Y,{\theta }^{\left(i\right)}\right) \end{gathered}$$

3. $M$步：$$\begin{array}{rl}& {\theta }^{\left(i+1\right)}=\arg \max Q\left(\theta ,{\theta }^{\left(i\right)}\right)\end{array}$$

4. 重复2. 3.，直到收敛。

2、EM算法的推广

F函数的极大算法

$F$函数：隐变量$Z$的概率分布为$\tilde{P}\left(Z\right)$，关于分布$\tilde{P}$与参数$\theta$的函数$$\begin{array}{rl}& \\ & F\left(\tilde{P},\theta \right)=E_{\tilde{P}}\left[\log P\left(Y,Z|\theta \right)\right]+H\left(\tilde{P}\right)\end{array}$$
其中，$H\left(\tilde{P}\right)=-E_{\tilde{P}}\left[\log \tilde{P}\left(Z\right)\right]$是分布$\tilde{P}\left(Z\right)$的熵。

对于固定的$\theta$，极大化$F$函数
$$\begin{array}{rl}& \\ & \underset{\tilde{P}}{\max }F\left(\tilde{P},\theta \right)\\ & s.t.\sum _Z{\tilde{P}}_{\theta }\left(Z\right)=1\end{array}$$

引入拉格朗日乘子$\lambda$，构造拉格朗日函数$$\begin{array}{rl}& \\ & L=E_{\tilde{P}}\left[\log P\left(Y,Z|\theta \right)\right]-E_{\tilde{P}}\left[\log \tilde{P}\left(Z\right)\right]+\lambda \left(1-\sum _Z\tilde{P}\left(Z\right)\right)\\ & =\sum _Z\log P\left(Y,Z|\theta \right)\tilde{P}\left(Z\right)-\sum _Z\log P\left(Z\right)\tilde{P}\left(Z\right)+\lambda -\lambda \sum _Z\tilde{P}\left(Z\right)\end{array}$$

将其对$\tilde{P}\left(Z\right)$求偏导，得$$\begin{array}{rl}& \\ & \frac{\partial L}{\partial \tilde{P}\left(Z\right)}=\log P\left(Y,Z|\theta \right)-1-\log P\left(Z\right)-\lambda \end{array}$$

令其等于0，得$$\begin{array}{rl}& \lambda ＝\log P\left(Y,Z|\theta \right)-1-\log P\left(Z\right)\\ & \frac{P\left(Y,Z|\theta \right)}{{\tilde{P}}_{\theta }\left(Z\right)}=e^{1+\lambda }\\ & \sum _{Z}P\left(Y,Z|\theta \right)＝e^{1+\lambda }\sum _Z{\tilde{P}}_{\theta }\left(Z\right)\end{array}$$

由于${\sum }_Z{\tilde{P}}_{\theta }\left(Z\right)=1$，得
$$\begin{array}{rl}& P\left(Y\right)=e^{1+\lambda }\end{array}$$

代回，得
$$\begin{array}{rl}& {\tilde{P}}_{\theta }\left(Z\right)=P\left(Z|Y,\theta \right)\end{array}$$

则$$\begin{array}{rl}& F\left(\tilde{P},\theta \right)=E_{\tilde{P}}\left[\log P\left(Y,Z|\theta \right)\right]+H\left(\tilde{P}\right)\\ & =\sum _Z\log P\left(Y,Z|\theta \right)\tilde{P}\left(Z\right)-\sum _Z\log P\left(Z\right)\tilde{P}\left(Z\right)\\ & =\sum _Z\tilde{P}\left(Z\right)\log \frac{P\left(Y,Z|\theta \right)}{\tilde{P}\left(Z\right)}\\ & =\sum _Z\tilde{P}\left(Z\right)\log \frac{P\left(Z|Y,\theta \right)P\left(Y|\theta \right)}{\tilde{P}\left(Z\right)}\\ & =\log P\left(Y|\theta \right)\sum _Z\tilde{P}\left(Z\right)\\ & =\log P\left(Y|\theta \right)\\ & =L\left(\theta \right)\end{array}$$

对于使$F\left(\tilde{P},\theta \right)$达到最大值的参数${\theta }^{\ast }$，有$$\begin{array}{r}L\left({\theta }^{\ast }\right)=F\left({\tilde{P}}_{{\theta }^{\ast }},{\theta }^{\ast }\right)=F\left({\tilde{P}}^{\ast },{\theta }^{\ast }\right)\end{array}$$

即，如果$F\left(\tilde{P},\theta \right)$在${\tilde{P}}^{\ast },{\theta }^{\ast }$达到局部极大值（全局最大值）,则$L\left({\theta }^{\ast }\right)$在${\tilde{P}}^{\ast },{\theta }^{\ast }$也达到局部极大值（全局最大值）。

由${\tilde{P}}_{\theta }\left(Z\right)=P\left(Z|Y,\theta \right)$，对固定的${\theta }^{\left(i\right)}$，$$\begin{array}{r}{\tilde{P}}^{\left(i+1\right)}\left(Z\right)={\tilde{P}}_{{\theta }^{\left(i\right)}}\left(Z\right)=P\left(Z|Y,{\theta }^{\left(i\right)}\right)\end{array}$$

使$F\left(\tilde{P},{\theta }^{\left(i\right)}\right)$极大化，
则$$\begin{array}{rl}& F\left({\tilde{P}}^{\left(i+1\right)},\theta \right)＝E_{{\tilde{P}}^{\left(i+1\right)}}\left[\log P\left(Y,Z|\theta \right)\right]+H\left({\tilde{P}}^{\left(i+1\right)}\right)\\ & =\sum _{Z}logP\left(Y,Z|\theta \right)P\left(Z|Y,{\theta }^{\left(i\right)}\right)+H\left({\tilde{P}}^{\left(i+1\right)}\right)\\ & =Q\left(\theta ,{\theta }^{\left(i\right)}\right)+H\left({\tilde{P}}^{\left(i+1\right)}\right)\end{array}$$

固定${\tilde{P}}^{\left(i+1\right)}$，求${\theta }^{\left(i\right)}$使$F\left({\tilde{P}}^{\left(i+1\right)},\theta \right)$极大化，得
$$\begin{array}{r}{\theta }^{\left(i+1\right)}=\arg \underset{\theta }{\max }F\left({\tilde{P}}^{\left(i+1\right)},\theta \right)=\arg \underset{\theta }{\max }Q\left(\theta ,{\theta }^{\left(i\right)}\right)\end{array}$$

即，由$EM$算法与$F$函数的极大－极大算法的到的参数估计序列${\theta }^{\left(i\right)},i=1,2,\cdots ,$是一致的。

GEM算法

$GEM$算法：

• 输入：观测数据$Y$，$F$函数；
• 输出：模型参数$\theta$

1. 初值${\theta }^{\left(0\right)}$
2. 第$i+1$次迭代，第1步：记${\theta }^{\left(i\right)}$为参数$\theta$的估计值，${\tilde{P}}^{\left(i\right)}$为函数$\tilde{P}$的估计。求${\tilde{P}}^{\left(i+1\right)}$使$\tilde{P}$极大化$F\left({\tilde{P}}^{\left(i+1\right)},\theta \right)$
3. 第2步：求${\theta }^{\left(i\right)}$使$F\left({\tilde{P}}^{\left(i+1\right)},\theta \right)$极大化
4. 重复（2）和（3），直到收敛。

3、EM算法在高斯混合模型学习中的应用

高斯混合模型

高斯混合模型$$\begin{array}{rl}& P\left(y|\theta \right)=\sum _{k=1}^K{\alpha }_k\varphi \left(y|{\theta }_k\right)\end{array}$$

其中，${\alpha }_k$是系数，${\alpha }_k\ge 0$，${\sum }_{k=1}^K{\alpha }_k=1$; $\varphi \left(y|{\theta }_k\right)$是高斯分布密度，${\theta }_k=\left({\mu }_k,{\sigma }_k^2\right)$,$$\begin{array}{rl}& \varphi \left(y|{\theta }_k\right)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{{\left(y-{\mu }_k\right)}^2}{2{\sigma }_k^2}\right)\end{array}$$

称为第$k$个分模型。

高斯混合模型参数估计的EM算法

假设观测数据$\left(y_1,y_2,\cdots ,y_N\right)$由高斯混合模型$$\begin{array}{rl}& P\left(y|\theta \right)=\sum _{k=1}^K{\alpha }_k\varphi \left(y|{\theta }_k\right)\end{array}$$

生成，其中，$\theta =\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K;{\theta }_1,{\theta }_2,\cdots ,{\theta }_K\right)$是模型参数。

1. 隐变量${\gamma }_{jk}$是0-1变量，表示观测数据$y_j$来自第$k$个分模型$$\begin{array}{rl}& \\ & {\gamma }_{jk}=\{\begin{array}{l}1,第j个观测数据来自第k个分模型\\ 0,否则\end{array}\left(j=1,2,\cdots ,N;k=1,2,\cdots ,K\right)\end{array}$$

   完全数据$$\begin{array}{rl}& \\ & \left(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jk}\right)j=1,2,\cdots ,N\end{array}$$

   完全数据似然函数$$\begin{array}{rl}& \\ & P\left(y,\gamma |\theta \right)=\prod _{j=1}^{N}P\left(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jK}|\theta \right)\\ & =\prod _{k=1}^K\prod _{j=1}^N{\left[{\alpha }_k\varphi \left(y_j|{\theta }_k\right)\right]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\varphi \left(y_j|{\theta }_k\right)\right]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{{\left(y-{\mu }_k\right)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\end{array}$$
   式中，$n_k={\sum }_{j=1}^N{\gamma }_{jk}$。

   完全数据的对数似然函数$$\begin{array}{rl}& \log P\left(y,\gamma |\theta \right)=\sum _{k=1}^K\left\{\sum _{j=1}^K{\gamma }_{jk}\log {\alpha }_k+\sum _{j=1}^K{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}{\left(y_j-{\mu }_k\right)}^2\right]\right\}\end{array}$$

2. $Q\left(\theta ,{\theta }^{\left(i\right)}\right)$函数$$\begin{array}{rl}& Q\left(\theta ,{\theta }^{\left(i\right)}\right)=E\left[\log P\left(y,\gamma |\theta \right)|y,{\theta }^{\left(i\right)}\right]\\ & =E\left\{\sum _{k=1}^K\left\{\sum _{j=1}^K{\gamma }_{jk}\log {\alpha }_k+\sum _{j=1}^K{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}{\left(y_j-{\mu }_k\right)}^2\right]\right\}\right\}\\ & =\sum _{k=1}^K\left\{\sum _{j=1}^{K}E\left({\gamma }_{jk}\right)\log {\alpha }_k+\sum _{j=1}^{K}E\left({\gamma }_{jk}\right)\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}{\left(y_j-{\mu }_k\right)}^2\right]\right\}\\ & =\sum _{k=1}^K\left\{\sum _{j=1}^K{\hat{\gamma }}_{jk}\log {\alpha }_k+\sum _{j=1}^K{\hat{\gamma }}_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}{\left(y_j-{\mu }_k\right)}^2\right]\right\}\end{array}$$

   其中，分模型$k$对观测数据$y_j$的响应度${\hat{\gamma }}_{jk}$是在当前模型参数下第$j$个观测数据来自第$k$个分模型的概率。$$\begin{array}{rl}& {\hat{\gamma }}_{jk}＝E\left({\gamma }_{jk}|y,\theta \right)=P\left({\gamma }_{jk}=1|y,\theta \right)\\ & =\frac{P\left({\gamma }_{jk}=1,y_j|\theta \right)}{{\sum }_{k=1}^{K}P\left({\gamma }_{jk}=1,y_j|\theta \right)}\\ & =\frac{{\alpha }_k\varphi \left(y|{\theta }_k\right)}{{\sum }_{k=1}^K{\alpha }_k\varphi \left(y|{\theta }_k\right)}\left(j=1,2,\cdots ,N;k=1,2,\cdots ,K\right)\end{array}$$

3. 求$Q\left(\theta ,{\theta }^{\left(i\right)}\right)$函数对$\theta$的极大值$$\begin{array}{r}{\theta }^{\left(i+1\right)}=\arg \max Q\left(\theta ,{\theta }^{\left(i\right)}\right)\end{array}$$

   得$$\begin{array}{rl}& {\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K\\ & {\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{\left(y_j-{\mu }_k\right)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K\\ & {\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,\cdots ,K\end{array}$$

高斯混合模型参数估计得EM算法：

• 输入：观测数据$y_1,y_2,\cdots ,y_N$，高斯混合模型；
• 输出：高斯混合模型参数

1. 取参数的初始值开始迭代
2. $E$步：计算分模型$k$对观测数据$y_i$的响应度$$\begin{array}{rl}& {\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi \left(y|{\theta }_k\right)}{{\sum }_{k=1}^K{\alpha }_k\varphi \left(y|{\theta }_k\right)}j=1,2,\cdots ,N;k=1,2,\cdots ,K\end{array}$$
3. $M$步：计算新迭代的模型参数$$\begin{array}{rl}& {\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K\\ & {\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{\left(y_j-{\mu }_k\right)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K\\ & {\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,\cdots ,K\end{array}$$
4. 重复2.步和3.步，直到收敛。

4、概要总结

1．EM算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法。含有隐变量的概率模型的数据表示为$\theta$ )。这里，$Y$是观测变量的数据，$Z$是隐变量的数据，$\theta$ 是模型参数。EM算法通过迭代求解观测数据的对数似然函数$L(\theta )=\log P(\mathrm{Y}|\theta )$的极大化，实现极大似然估计。每次迭代包括两步：

$E$步，求期望，即求$logP\left(Z|Y,\theta \right)$ )关于$ P\left(Z | Y, \theta^{(i)}\right)$)的期望：

$$Q\left(\theta ,{\theta }^{(i)}\right)=\sum _Z\log P(Y,Z|\theta )P\left(Z|Y,{\theta }^{(i)}\right)$$

称为$Q$函数，这里${\theta }^{(i)}$是参数的现估计值；

$M$步，求极大，即极大化$Q$函数得到参数的新估计值：

$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q\left(\theta ,{\theta }^{(i)}\right)$$

在构建具体的EM算法时，重要的是定义$Q$函数。每次迭代中，EM算法通过极大化$Q$函数来增大对数似然函数$L(\theta )$。

2．EM算法在每次迭代后均提高观测数据的似然函数值，即

$$P\left(Y|{\theta }^{(i+1)}\right)⩾P\left(Y|{\theta }^{(i)}\right)$$

在一般条件下EM算法是收敛的，但不能保证收敛到全局最优。

3．EM算法应用极其广泛，主要应用于含有隐变量的概率模型的学习。高斯混合模型的参数估计是EM算法的一个重要应用，下一章将要介绍的隐马尔可夫模型的非监督学习也是EM算法的一个重要应用。

4．EM算法还可以解释为$F$函数的极大-极大算法。EM算法有许多变形，如GEM算法。GEM算法的特点是每次迭代增加$F$函数值（并不一定是极大化$F$函数），从而增加似然函数值。