EM算法

Posted on

EM算法

1. 概念

概率模型有时既含有观测变量（observable variable），又含有隐变量（hidden variable）或潜在变量（latent variable），如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或者贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单的使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，是一种迭代算法。
EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。
用\(Y\)表示观测随机变量的数据，\(Z\)表示隐随机变量的数据，\(Y\)和\(Z\)连在一起称为完全数据（complete-data），观测数据\(Y\)又称为不完全数据（incomplete-data）。不完全数据的对数似然函数是\(L(\theta) = \log{P(Y| \theta)}\)，完全数据的对数似然函数为\(\log{P(Y,Z| \theta )}\)。
EM算法通过迭代求\(L(\theta)=\log{P(Y| \theta )}\)。每次迭代包含两步：E步，求期望；M步，求极大化。
EM算法：
输入：观测变量数据\(Y\)，隐变量数据\(Z\)，联合分布\(P(Y,Z| \theta )\)，条件分布\(P(Z|Y, \theta )\)；
输出：模型参数\(\theta\)。
（1）选择参数的初始值\(\theta^{(0)}\)，开始迭代；
（2）E步：记\(\theta^{(i)}\)为第\(i\)次迭代参数\(\theta\)的估计值，在第\(i+1\)次迭代的E步，计算： \[\begin{align} Q(\theta,\theta^{(i)}) &= E_Z[\log{P(Y,Z|\theta)|Y,\theta^{(i)}}] \nonumber \\ &= \sum_Z \log{P(Y,Z|\theta)P(Z|Y,\theta^{(i)})} \nonumber \end{align}\] \(P(Z|Y, \theta^{(i)})\)是在给定观测数据\(Y\)和当前的参数估计\(\theta^{(i)}\)下隐变量数据\(Z\)的条件概率分布；
（3）M步：求使\(Q(\theta,\theta^{(i)})\)极大化的\(\theta\)，确定第\(i+1\)次迭代的参数的估计值\(\theta^{(i+1)}\)： \[\theta^{(i+1)} = \arg \max_{\theta} Q(\theta, \theta^{(i)})\] （4）重复第（2）步和第（3）步，直到收敛。
Q函数：完全数据的对数似然函数\(\log{P(Y,Z|\theta)}\)关于在给定观测数据\(Y\)和当前参数\(\theta^{(i)}\)下对未观测数据\(Z\)的条件概率分布\(P(Z|Y,\theta^{(i)})\)的期望称为Q函数，即： \[Q(\theta,\theta^{(i)}) = E_Z[\log{P(Y,Z|\theta)|Y,\theta^{(i)}}]\] \(Q(\theta, \theta^{(i)})\)的第1个变元表示要极大化的参数，第2个变元表示参数当前估计值，每次迭代实际在求Q函数及其极大。
一个直观了解EM算法思路的例子就是K-Means算法：在K-Means聚类时，每个聚类簇的质心是隐变量。假设K个初始化质心，即EM算法中的E步；然后计算得到每个样本最近的质心，并把样本积累到这个质心，即EM算法中的M步；重复E步和M步，直到质心不在改变。

2. 导出

关于EM算法的导出一般有两种方式：

2.1. “通俗”方式

通俗角度的话，求极大就是求似然函数的极大，一般都是对数似然。通常的做法是对似然函数求偏导，然后令偏导等于0，即求得参数的近似最优值，但是包含隐变量的模型不能直接进行似然函数的偏导，如果假设已经知道隐变量的值，就可以将似然函数简化，并求偏导。
对于\(m\)个样本观察数据\(x=(x^{(1)},x^{(2)},\cdots, x^{(m)})\)，找出样本的模型参数\(\theta\)，极大化模型分布的对数似然函数： \[\theta = \arg \max_{\theta} \sum_{i=1}^m \log{P(x^{(i)}| \theta)}\] 如果观察数据中包含为观测到的隐变量\(z=(z^{(1)},z^{(2)},\cdots,z^{(m)})\)，此时极大化模型分布的对数似然函数如下： \[\theta =\arg \max_{\theta} \sum_{i=1}^m \log{P(x^{(i)}| \theta) = \arg \max_{\theta} \sum_{i=1}^m} \log{\sum_{z^{(i)}}P(x^{(i)},z^{(i)}|\theta)}\] 上式无法直接求出\(\theta\)，首先利用Jensen不等式对式子进行缩放： \[\begin{align} \sum_{i=1}^m \log{\sum_{z^{(i)}} P(x^{(i)},z^{(i)}|\theta)} &= \sum_{i=1}^m \log{\sum_{z^{(i)}} [Q_i(z^{(i)}) \frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}]} \nonumber \\ &\geq \sum_{i=1}^m \sum_{z^{(i)}}[Q_i(z^{(i)}) \log{\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}}] \nonumber \end{align}\] 上式中引入一个未知的新的分布\(Q_i(z^{(i)})\)，注意这个函数并不是Q函数，其中Jensen不等式：由于对数函数是凹函数，所以有\(f(E[X]) \geq E[f(X)]\)，其中要满足Jensen不等式中的等号，那么有： \[\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})} = c\] 其中\(c\)为常数，并且\(Q_i(z^{(i)})\)是一个分布，满足： \[\sum_z Q_i(z^{(i)}) = 1\] 从上面两式可以得到： \[Q_i(z^{(i)}) = \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)\] 如果\(Q_i(z^{(i)}) = P(z^{(i)}|x^{(i)},\theta^{(i)})\)，则\(\sum_{i=1}^m \sum_{z^{(i)}}[Q_i(z^{(i)}) \log{\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}}]\)是包含隐变量的对数似然的一个下界，如果能极大化这个下界，也是在尝试极大化对数似然，即需要极大化下式： \[\arg \max \sum_{i=1}^m \sum_{z{(i)}}[Q_i(z^{(i)}) \log{\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}}]\] 去掉上式中的常数部分，需要极大化的对数似然的下界为： \[\arg \max_{\theta} \sum_{i=1}^m \sum_{z^{(i)}} [Q_i(z^{(i)}) \log{P(x^{(i)},z^{(i)}|\theta)}]\] 上式就是EM算法中的M步，而\(\sum_{z^{(i)}} [Q_i(z^{(i)}) \log{P(x^{(i)},z^{(i)}|\theta)}]\)就是Q函数所定义的式子。

2.2. ”正式“方式

通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法。当面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据）\(Y\)关于参数\(\theta\)的对数似然函数，即极大化 \[\begin{align} L(\theta) &= \log{P(Y|\theta)} = \log{\sum_Z P(Y,Z|\theta)} \nonumber \\ &= \log{\{\sum_Z P(Y|Z,\theta)P(Z|\theta)\}} \nonumber \end{align}\] 极大化的主要困难是上式中有未观测数据并包含和的对数。
EM算法通过迭代逐步近似极大化\(L(\theta)\)的，假设在第\(i\)次迭代后\(\theta\)的估计值是\(\theta^{(i)}\)，则希望新估计值\(\theta\)能使\(L(\theta)\)增加，即\(L(\theta) \geq L(\theta^{(i)})\)，并逐步达到极大值，为此考虑两者的差： \[L(\theta) - L(\theta^{(i)}) = \log{\{\sum_Z P(Y|Z,\theta)P(Z|\theta)\}} - \log{P(Y|\theta^{(i)})}\] 利用Jensen不等式得到其下界： \[\begin{align} L(\theta) - L(\theta^{(i)}) &= \log{\sum_Z P(Z|Y,\theta^{(i)}) \frac{P(Y|Z,\theta^{(i)})P(Z|\theta)}{P(Z|Y,\theta^{(i)})}} - \log{P(Y|\theta^{(i)})} \nonumber \\ &\geq \sum_Z{P(Z|Y,\theta^{(i)}) \log{\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}}} - \log{P(Y|\theta^{(i)})} \nonumber \\ &= \sum_Z P(Z|Y,\theta^{(i)}) \log{\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}} \nonumber \end{align}\] 令 \[B(\theta,\theta^{(i)}) \hat{=} L(\theta^{(i)})+\sum_Z P(Z|Y,\theta^{(i)}) \log{\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}}\] 则 \[L(\theta) \geq B(\theta,\theta^{(i)})\] 函数\(B(\theta,\theta^{(i)})\)是\(L(\theta)\)的一个下界，并且满足： \[L(\theta^{(i)}) = B(\theta^{(i)},\theta^{(i)})\] 因此，任何可以使\(B(\theta,\theta^{(i)})\)增大的\(\theta\)，也可以使\(L(\theta)\)增大，为了使\(L(\theta)\)有尽可能大的增长，选择\(\theta^{(i+1)}\)使\(B(\theta,\theta^{(i)})\)达到极大，即 \[\theta^{(i+1)}=\arg \max_{\theta}{B(\theta,\theta^{(i)})}\] 现在求\(\theta^{(i+1)}\)的表达式，省去对\(\theta\)的极大化而言的常数项，有 \[\begin{align} \theta^{(i+1)} &= \arg \max_{\theta}\{L(\theta^{(i)})+\sum_Z P(Z|Y,\theta^{(i)}) \log{\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}}\} \nonumber \\ &= \arg \max_{\theta} \{\sum_Z P(Z|Y,\theta^{(i)}) \log{P(Y|Z,\theta)P(Z|\theta)}\} \nonumber \\ &= \arg \max_{\theta} \{\sum_Z P(Z|Y,\theta^{(i)}) \log{P(Y,Z|\theta)}\} \nonumber \\ &= \arg \max_{\theta} Q(\theta, \theta^{(i)}) \nonumber \end{align}\] 上式等价于EM算法的一次迭代，即求Q函数及其极大化，EM算法通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
EM算法可以用于生成模型的非监督学习。
下图给出EM算法的直观解释：

3. 收敛性

定理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)}) \geq P(Y|\theta^{(i)})\] 证明：由于 \[P(Y|\theta) = \frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}\] 取对数有 \[\log{P(Y|\theta)} = \log{P(Y,Z|\theta)} - \log{P(Z|Y,\theta)}\] Q函数定义为： \[Q(\theta, \theta^{(i)}) = \sum_Z \log{P(Y,Z|\theta)P(Z|Y,\theta^{(i)})}\] 令 \[H(\theta,\theta^{(i)}) = \sum_Z \log{P(Z|Y,\theta)P(Z,|Y,\theta^{(i)})}\] 于是对数似然函数可以写成 \[\log{P(Y|\theta)} = Q(\theta,\theta^{(i)}) - H(\theta,\theta^{(i)})\] 将上式中分别取\(\theta\)为\(\theta^{(i)}\)和\(\theta^{(i+1)}\)并相减，有 \[\log{P(Y|\theta^{(i+1)})} - \log{P(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)})]\] 为证明定理中的不等式成立，只需证明上式中右端的式子非负，对于第1项，由于\(\theta^{(i+1)}\)使\(Q(\theta,\theta^{(i)})\)达到极大，所以有 \[Q(\theta^{(i+1)},\theta^{(i)}) - Q(\theta^{(i)},\theta^{(i)}) \geq 0\] 对于第2项，可得 \[\begin{align} H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)}) &= \sum_Z \{\log{\frac{P(Z|Y,\theta^{(i+1)})}{P(Z,Y,\theta^{(i)})}}\}P(Z|Y,\theta^{(i)}) \nonumber \\ &\leq \log{\{\sum_Z \frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}P(Z|Y,\theta^{(i)})\}} \nonumber \\ &= log{\{\sum_Z P(Z|Y,\theta^{(i+1)})\}} = 0 \nonumber \end{align}\] 有此可知右端式子非负。
定理2：设\(L(\theta) = \log{P(Y|\theta)}\)为观测数据的对数似然函数，\(\theta^{(i)}(i=1,2,\cdots)\)为EM算法得到的参数估计序列，\(L(\theta^{(i)})(i=1,2,\cdots)\)为对应的似然函数序列：

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

定理2关于函数\(Q(\theta,\theta')\)与\(L(\theta)\)的条件大多数情况下是满足的。EM算法的收敛性包含关于对数似然函数序列\(L(\theta^{(i)})\)的收敛性和关于参数估计序列\(\theta^{(i)}\)的收敛性两层意思，前者并不蕴含后者。
定理只能保证参数估计序列收敛到对数似然函数序列的稳定点，不能保证收敛到极大值点，所以初值的选择非常重要。

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

EM算法是学习高斯混合模型（Gaussian misture model）的有效方法。
### 4.1. 高斯混合模型 高斯混合模型：高斯混合模型是指具有如下形式的概率分布模型： \[P(y|\theta) = \sum_{k=1}^K \alpha_k \phi(y|\theta_k)\] 其中，\(\alpha_k\)是系数，\(\alpha_k \geq 0\)，\(\sum_{k=1}^K \alpha_k = 1\)；\(\phi(y|\theta_k)\)是高斯分布密度，\(\theta_k = (\mu_k,\sigma_k^2)\)， \[\phi(y|\theta_k) = \frac{1}{\sqrt{2 \pi}\sigma_k} \exp(-\frac{(y-y_k)^2}{2\sigma_k^2})\] 称为第k个分模型。
一般混合模型可以由任意概率分布密度代替上式中的高斯分布密度，这里只介绍最常用的高斯混合模型。
### 4.2. 高斯混合模型参数估计的EM算法 假设观测数据\(y_1,y_2,\cdots,y_N\)由高斯混合模型生成， \[P(y|\theta) = \sum_{k=1}^K \alpha_k \phi(y|\theta_k)\] 其中，\(\theta=(\alpha_1,\alpha_2,\cdots,\alpha_K;\theta_1,\theta_2,\cdots,\theta_K)\)。这里用EM算法估计高斯混合模型的参数\(\theta\)。
观测数据\(y_j,j=1,2,\cdots,N\)的产生：依概率\(\alpha_k\)选择第k个高斯分布分模型\(\phi(y|\theta_k)\)；依第k个分模型的概率分布\(\phi(y|\theta_k)\)生成观测数据\(y_j,j=1,2,\cdots,N\)。观测数据已知，高斯分布分模型未知，\(k=1,2,\cdots,K\)，用隐变量\(\gamma_{jk}\)表示，定义如下： \[\gamma_{jk} = \begin{cases} 1, & 第j个观测来自第k个分模型 \\ 0, & 否则 \end{cases}\] \(\gamma_{jk}\)是0-1随机变量。
\(y_j\)是观测数据，\(\gamma_{jk}\)是未观测数据，完全数据是 \[(y_j,\gamma_1,\gamma_2,\cdots,\gamma_{jk}), j=1,2,\cdots,N\] 完全数据的似然函数： \[\begin{align} P(y,\gamma|\theta) &= \prod_{j=1}^N P(y_j,\gamma_1,\gamma_2,\cdots,\gamma_{jk}|\theta) \nonumber \\ &= \prod_{k=1}^K \prod_{j=1}^N[\alpha_k \phi(y_j|\theta_k)]^{\gamma_{jk}} \nonumber \\ &= \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N[\phi(y_j|\theta)]^{\gamma_{jk}} \nonumber \\ &= \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N[\frac{1}{\sqrt{2\pi}\sigma_k} \exp(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2})]^{\gamma_{jk}} \nonumber \end{align}\] 其中，\(n_k = \sum_{j=1}^N \gamma_{jk}\)，\(\sum_{k=1}^K n_k = N\)。
那么，完全数据的对数似然函数为 \[\log{P(y,\gamma|\theta)} = \sum_{k=1}^K\{n_k \log{\alpha_k}+\sum_{j=1}^N \gamma_{jk}[\log{(\frac{1}{\sqrt{2\pi}})}-\log{\sigma_k- \frac{1}{2\sigma_k^2}(y_j-\mu_k)^2}]\}\] E步：确定Q函数
\[\begin{align} Q(\theta,\theta^{(i)}) &= E[\log{P(y,\gamma|\theta)}|y,\theta^{(i)}] \nonumber \\ &= E\{\sum_{k=1}^K\{n_k \log{\alpha_k}+\sum_{j=1}^N \gamma_{jk}[\log{(\frac{1}{\sqrt{2\pi}})}-\log{\sigma_k- \frac{1}{2\sigma_k^2}(y_j-\mu_k)^2}]\}\} \nonumber \\ &= \sum_{k=1}^K\{\sum_{j=1}^N(E\gamma_{jk}) \log{\alpha_k}+\sum_{j=1}^N (E\gamma_{jk})[\log{(\frac{1}{\sqrt{2\pi}})}-\log{\sigma_k- \frac{1}{2\sigma_k^2}(y_j-\mu_k)^2}]\} \nonumber \end{align}\] 这里需要计算\(E(\gamma_{jk}|y,\theta)\)，记为\(\hat{\gamma}_{jk}\)。 \[\begin{align} \hat{\gamma}_{jk} &= E(\gamma_{jk}|y,\theta) = P(\gamma_{jk}=1|y,\theta) \nonumber \\ &= \frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^K P(\gamma_{jk}=1,y_j|\theta)} \nonumber \\ &= \frac{P(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^KP(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)} \nonumber \\ &= \frac{\alpha_k \phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k \phi(y_j|\theta_k)} \nonumber , & j=1,2,\cdots,N; k=1,2,\cdots,K \end{align}\] \(\hat{\gamma}_{jk}\)是在当前模型参数下第\(j\)个观测数据来自第\(k\)个分模型的概率，称为分模型\(k\)对观测数据\(y_j\)的响应度。
将\(\hat{\gamma}_{jk} = E\gamma_{jk}\)及\(n_k = \sum_{j=1}^N E\gamma_{jk}\)带入Q函数中得 \[Q(\theta,\theta^{(i)}) = \sum_{k=1}^K\{n_k \log{\alpha_k}+\sum_{j=1}^N \hat{\gamma}_{jk}[\log{(\frac{1}{\sqrt{2\pi}})}-\log{\sigma_k- \frac{1}{2\sigma_k^2}(y_j-\mu_k)^2}]\}\] M步：求极大值
求Q函数对\(\theta\)的极大值： \[\theta^{(i+1)} = \arg \max_{\theta} Q(\theta,\theta^{(i)})\] 用\(\hat{\mu}_k\)，\(\hat{\sigma}_{k}^2\)及\(\hat{\alpha}_k\)，\(k=1,2,\cdots,K\)，表示\(\theta^{(i+1)}\)的各参数，求\(\hat{\mu}_k\)，\(\hat{\sigma}_{k}^2\)只需分别对Q函数求偏导并令其为0；求\(\hat{\alpha}_k\)是在\(\sum_{k=1}^K \alpha_k=1\)条件下求偏导并令其为0.结果如下： \[\begin{equation} \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 \nonumber \\ \hat{\sigma}_k^2 = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}}，k=1,2,\cdots,K \nonumber \\ \hat{\alpha}_k = \frac{n_k}{N} = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}}{N}，k=1,2,\cdots,K \nonumber \end{equation}\] 重复以上计算，直到对数似然函数值不再有明显变化。