【李航】统计学习方法--9. EM算法及其推广(详细推导）

EM算法是一种迭代算法,用于含有隐变量(hidden variable)的概率模型参数的极大似然估计,或极大后验概率估计。
EM算法的每次迭代由两步组成:E步,求期望(expectation);M步,求极大(maximization)。所以这一算法称为期望极大算法(expectation maximization algorithm),简称EM算法。

---

9.1 EM算法的引入

---

概率模型有时既含有观测变量( observable variable),又含有隐变量或潜在变量( (latent variable)。EM算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法。

---

9.1.1 EM算法

---

开局一个经典EM算法例子

三硬币模型 假设有 3 枚硬币, 分别记作 $\mathrm{A},\mathrm{B},\mathrm{C}$ 。这些硬币正面出现 的概率分别是 $\pi ,p$ 和 $q_{\circ }$ 进行如下掷硬币试验: 先掷硬币 $\mathrm{A}$, 根据其结果选出硬币 $\mathrm{B}$ 或硬币 $\mathrm{C}$, 正面选硬币 B, 反面选硬币 $\mathrm{C}$; 然后掷选出的硬币, 掷硬币的结果, 出现正 面记作 1 , 出现反面记作 0 ; 独立地重复 $n$ 次试验 (这里, $n=10$ ), 观测结果如下:

$$1,1,0,1,0,0,1,0,1,1$$

假设只能观测到掷硬币的结果, 不能观测掷硬币的过程。问如何估计三硬币正面出现的概率, 即三硬币模型的参数。

解 三硬币模型可以写作

$$\begin{array}{rl}P(y\mid \theta )& =\sum _{z}P(y,z\mid \theta )=\sum _{z}P(z\mid \theta )P(y\mid z,\theta )\\ & =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}$$

这里, 随机变量 $y$ 是观测变量, 表示一次试验观测的结果是 1 或 0 ; 随机变量 $z$ 是隐 变量, 表示未观测到的掷硬币 $\mathrm{A}$ 的结果; $\theta =(\pi ,p,q)$ 是模型参数。这一模型是以上数据的生成模型。注意, 随机变量 $y$ 的数据可以观测,随机变量 $z$ 的数据不可观测。

将观测数据表示为 $Y={\left(Y_1,Y_2,\cdots ,Y_n\right)}^{\mathrm{T}}$, 未观测数据表示为 $Z={\left(Z_1,Z_2,\cdots ,Z_n\right)}^{\mathrm{T}}$
则观测数据的似然函数为

$$P(Y\mid \theta )=\sum _{Z}P(Z\mid \theta )P(Y\mid Z,\theta )$$

即

$$P(Y\mid \theta )=\prod _{j=1}^n\left[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}\right]$$

考虑求模型参数 $\theta =(\pi ,p,q)$ 的极大似然估计, 即

$$\hat{\theta }=\arg \underset{\theta }{\max }\log P(Y\mid \theta )$$

下面给出针对以上问题的 EM 算法

EM 算法首先选取参数的初值, 记作 ${\theta }^{(0)}=\left({\pi }^{(0)},p^{(0)},q^{(0)}\right)$, 然后通过下面的 步骤迭代计算参数的估计值, 直至收敛为止。第 $i$ 次迭代参数的估计值为 ${\theta }^{(i)}=$ $\left({\pi }^{(i)},p^{(i)},q^{(i)}\right)。$ EM 算法的第 $i+1$ 次迭代如下。
$\mathrm{E}$ 步: 计算在模型参数 ${\pi }^{(i)},p^{(i)},q^{(i)}$ 下观测数据 $y_j$ 来自掷硬币 $\mathrm{B}$ 的概率

$${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}{\left(p^{(i)}\right)}^{y_j}{\left(1-p^{(i)}\right)}^{1-y_j}}{{\pi }^{(i)}{\left(p^{(i)}\right)}^{y_j}{\left(1-p^{(i)}\right)}^{1-y_j}+\left(1-{\pi }^{(i)}\right){\left(q^{(i)}\right)}^{y_j}{\left(1-q^{(i)}\right)}^{1-y_j}}$$

M 步: 计算模型参数的新估计值

$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}$$

$$\begin{array}{c}{p}^{(i+1)}=\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{(i+1)}}\\ q^{(i+1)}=\frac{{\sum }_{j=1}^n\left(1-{\mu }_j^{(i+1)}\right)y_j}{{\sum }_{j=1}^n\left(1-{\mu }_j^{(i+1)}\right)}\end{array}$$

EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

• 用 $Y$ 表示观测随机变量的数据, $Z$ 表示隐随机变量的数据。
• $Y$ 和 $Z$ 连 在一起称为完全数据 (complete-data), 观测数据 $Y$ 又称为不完全数据 (incompletedata）。
• 假设给定观测数据 $Y$, 其概率分布是 $P(Y\mid \theta )$, 其中 $\theta$ 是需要估计的模型参数, 那么不完全数据 $Y$ 的似然函数是 $P(Y\mid \theta )$, 对数似然函数 $L(\theta )=\log P(Y\mid \theta );$ 假设 $Y$ 和 $Z$ 的联合概率分布是 $P(Y,Z\mid \theta )$, 那么完全数据的对数似然函数是 $\log P(Y,Z\mid \theta )$ 。
• $\mathrm{E}\mathrm{M}$ 算法通过迭代求 $L(\theta )=\log P(Y\mid \theta )$ 的极大似然估计。每次迭代包含两步: $\mathrm{E}$
  步, 求期望; $\mathrm{M}$ 步, 求极大化。

EM 算法 输入: 观测变量数据 $Y$, 隐变量数据 $Z$, 联合分布 $P(Y,Z\mid \theta )$, 条件分布 $P(Z\mid Y,\theta )$; 输出: 模型参数 $\theta$ 。

1. 选择参数的初值 ${\theta }^{(0)}$, 开始迭代;

2. $\mathrm{E}$ 步: 记 ${\theta }^{(i)}$ 为第 $i$ 次迭代参数 $\theta$ 的估计值, 在第 $i+1$ 次迭代的 $\mathrm{E}$ 步, 计算
   $$\begin{array}{rl}Q\left(\theta ,{\theta }^{(i)}\right)& =E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]\\ & =\sum _Z\log P(Y,Z\mid \theta )P\left(Z\mid Y,{\theta }^{(i)}\right)\end{array}$$

   这里, $P\left(Z\mid Y,{\theta }^{(i)}\right)$ 是在给定观测数据 $Y$ 和当前的参数估计 ${\theta }^{(i)}$ 下隐变量数据 $Z$ 的条 概率分布;

3. M 步: 求使 $Q\left(\theta ,{\theta }^{(i)}\right)$ 极大化的 $\theta$, 确定第 $i+1$ 次迭代的参数的估计值 ${\theta }^{(i+1)}$
   $${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q\left(\theta ,{\theta }^{(i)}\right)$$

4. 重复第(2) 步和第 $(3)$ 步, 直到收敛。

$(Q$ 函数） 完全数据的对数似然函数 $\log P(Y,Z\mid \theta )$ 关于在给定观测数 据 $Y$ 和当前参数 ${\theta }^{(i)}$ 下对未观测数据 $Z$ 的条件概率分布 $P\left(Z\mid Y,{\theta }^{(i)}\right)$ 的期望称为 $Q$ 函数，即

$$Q\left(\theta ,{\theta }^{(i)}\right)=E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]$$

下面关于 EM 算法作几点说明:
步骤 (1）参数的初值可以任意选择, 但需注意 $\mathrm{E}\mathrm{M}$ 算法对初值是敏感的。
步骤 (2) $\mathrm{E}$ 步求 $Q\left(\theta ,{\theta }^{(i)}\right)。Q$ 函数式中 $Z$ 是未观测数据, $Y$ 是观测数据。注 意, $Q\left(\theta ,{\theta }^{(i)}\right)$ 的第 1 个变元表示要极大化的参数, 第 2 个变元表示参数的当前估计 值。每次迭代实际在求 $Q$ 函数及其极大。
步骤 (3) $\mathrm{M}$ 步求 $Q\left(\theta ,{\theta }^{(i)}\right)$ 的极大化, 得到 ${\theta }^{(i+1)}$, 完成一次迭代 ${\theta }^{(i)}\to {\theta }^{(i+1)}$ 。 后面将证明每次迭代使似然函数增大或达到局部极值。
步骤 (4） 给出停止迭代的条件, 一般是对较小的正数 ${\epsilon }_1,{\epsilon }_2$, 若满足

$$\Vert {\theta }^{(i+1)}-{\theta }^{(i)}\Vert <{\epsilon }_1\text{ 或 }\Vert Q\left({\theta }^{(i+1)},{\theta }^{(i)}\right)-Q\left({\theta }^{(i)},{\theta }^{(i)}\right)\Vert <{\epsilon }_2$$

则停止迭代。备注：实际工程中可能会设置epoch迭代停止

---

9.1.2 $\mathrm{E}\mathrm{M}$ 算法的导出

---

通过近似求解观测数据的对数似然函数的极大化问题来导出 $\mathrm{E}\mathrm{M}$ 算法, 由此 可以清楚地看出 EM算法的作用。 我们面对一个含有隐变量的概率模型, 目标是极大化观测数据 (不完全数据) $Y$ 关于参数 $\theta$ 的对数似然函数, 即极大化

$$\begin{array}{rl}L(\theta )& =\log P(Y\mid \theta )=\log \sum _{Z}P(Y,Z\mid \theta )\\ & =\log \left(\sum _{Z}P(Y\mid Z,\theta )P(Z\mid \theta )\right)\end{array}$$

事实上, $\mathrm{E}\mathrm{M}$ 算法是通过迭代逐步近似极大化 $L(\theta )$ 的。假设在第 $i$ 次迭代后 $\theta$ 的 估计值是 ${\theta }^{(i)}$ 。我们希望新估计值 $\theta$ 能使 $L(\theta )$ 增加, 即 $L(\theta )>L\left({\theta }^{(i)}\right)$, 并逐步达到极 大值。为此, 考虑两者的差:

$$L(\theta )-L\left({\theta }^{(i)}\right)=\log \left(\sum _{Z}P(Y\mid Z,\theta )P(Z\mid \theta )\right)-\log P\left(Y\mid {\theta }^{(i)}\right)$$

利用 Jensen 不等式 (Jensen inequality)$\log {\sum }_j{\lambda }_j{y}_j⩾{\sum }_j{\lambda }_j\log y_j$, 其中 ${\lambda }_j⩾0,{\sum }_j{\lambda }_j=1$ 。

${\sum }_{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)=1$,

未观测数据 $Z$ 的条件概率分布

$$\begin{array}{rl}L(\theta )-L\left({\theta }^{(i)}\right)& =\log \left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}\right)-\log P\left(Y\mid {\theta }^{(i)}\right)\\ & ⩾\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}-\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\cdot \log P\left(Y\mid {\theta }^{(i)}\right)\\ & =\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)}\end{array}$$

令

$$B\left(\theta ,{\theta }^{(i)}\right)\widehat{=}L\left({\theta }^{(i)}\right)+\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)}①$$

则

$$L(\theta )⩾B\left(\theta ,{\theta }^{(i)}\right)$$

即函数 $B\left(\theta ,{\theta }^{(i)}\right)$ 是 $L(\theta )$ 的一个下界, 而且由式①可知,

$$L\left({\theta }^{(i)}\right)=B\left({\theta }^{(i)},{\theta }^{(i)}\right)$$

因此, 任何可以使 $B\left(\theta ,{\theta }^{(i)}\right)$ 增大的 $\theta$, 也可以使 $L(\theta )$ 增大。为了使 $L(\theta )$ 有尽可能大 的增长, 选择 ${\theta }^{(i+1)}$ 使 $B\left(\theta ,{\theta }^{(i)}\right)$ 达到极大, 即

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

现在求 ${\theta }^{(i+1)}$ 的表达式。省去对 $\theta$ 的极大化而訁是常数的项, 由式 ${\theta }^{(i+1)}=\arg {\max }_{\theta }B\left(\theta ,{\theta }^{(i)}\right)$ 、式 ①及式 ${\theta }^{(i+1)}=\arg {\max }_{\theta }Q\left(\theta ,{\theta }^{(i)}\right)$, 有

$$\begin{array}{rl}{\theta }^{(i+1)}& =\arg \underset{\theta }{\max }\left(L\left({\theta }^{(i)}\right)+\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)}\right)\\ & =\arg \underset{\theta }{\max }\left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log (P(Y\mid Z,\theta )P(Z\mid \theta ))\right)\\ & =\arg \underset{\theta }{\max }\left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log P(Y,Z\mid \theta )\right)\\ & =\arg \underset{\theta }{\max }Q\left(\theta ,{\theta }^{(i)}\right)\end{array}$$

上式等价于 $\mathrm{E}\mathrm{M}$ 算法的一次迭代, 即求 $Q$ 函数及其极大化。EM 算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

$\mathrm{E}\mathrm{M}$ 算法的直观解释

---

9.1.3 EM 算法在无监督学习中的应用

---

监督学习是由训练数据 { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x N , y N ) } \left\{\left(x_{1}, y_{1}\right),\left(x_{2}, y_{2}\right), \cdots,\left(x_{N}, y_{N}\right)\right\} {(x1​,y1​),(x2​,y2​),⋯,(xN​,yN​)} 学习条件概率分布$P(Y\mid X)$ 或决策函数 $Y=f(X)$ 作为模型, 用于分类、回归、标注等任务。这时训练数据中的每个样本点由输入和输出对组成。
有时训练数据只有输入没有对应的输出 { ( x 1 , ∙ ) , ( x 2 , ∙ ) , ⋯   , ( x N , ∙ ) } \left\{\left(x_{1}, \bullet\right),\left(x_{2}, \bullet\right), \cdots,\left(x_{N}, \bullet\right)\right\} {(x1​,∙),(x2​,∙),⋯,(xN​,∙)}, 从这样的数据学习模型称为无监督学习问题。

EM 算法可以用于生成模型的无监督学习。生成模型由联合概率分布 $P(X,Y)$ 表示, 可以认为无监督学习训练数据是联合概率分布产生的数据。 $X$ 为观测数据, $Y$ 为未观测数据。

---

9.2 EM算法的收敛性

---

设 $P(Y\mid \theta )$ 为观测数据的似然函数, ${\theta }^{(i)}(i=1,2,\cdots )$ 为 $\mathrm{E}\mathrm{M}$ 算法得到的参数估计序列, $P\left(Y\mid {\theta }^{(i)}\right)(i=1,2,\cdots )$ 为对应的似然函数序列, 则 $P\left(Y\mid {\theta }^{(i)}\right)$ 是单调递增的，即

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

证明 由于

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

条件联合分布的分解

条件联合分布的分解

$P(X=a,Y=b\mid Z=c)=P(X=a\mid Y=b,Z=c)P(Y=b\mid Z=c)$

取对数有

$$\log P(Y\mid \theta )=\log P(Y,Z\mid \theta )-\log P(Z\mid Y,\theta )$$

由式 $Q\left(\theta ,{\theta }^{(i)}\right)=E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]$

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

令

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

为什么？

由于 : $\log P(Y\mid \theta )=\log P(Y,Z\mid \theta )-\log P(Z\mid Y,\theta )$ 两边对Z求积分 :
​$\begin{array}{rl}\log P(Y\mid \theta )& amp;=\sum _{z}P\left(Z\mid Y,{\theta }^{(i)}\right)\cdot \log P(Y\mid \theta )\\ & amp;=\log P(Y\mid \theta )\sum _{z}P\left(Z\mid Y,{\theta }^{(i)}\right)\\ & amp;=\log P(Y\mid \theta )\end{array}$
对左边公式求积分，由于 $\log P(Y\mid \theta )$ 无关Z的积分，后项求积分等于1
​$\mathrm{Righ}t=\log P(Y,Z\mid \theta )-\log P(Z\mid Y,\theta )$
​$={\sum }_{z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log P(Y,Z\mid \theta )-{\sum }_{z}P(Z\mid Y,{\theta }^{(i)}\log P(Z\mid Y,\theta )$

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

$$\log P(Y\mid \theta )=Q\left(\theta ,{\theta }^{(i)}\right)-H\left(\theta ,{\theta }^{(i)}\right)$$

在上式中分别取 $\theta$ 为 ${\theta }^{(i)}$ 和 ${\theta }^{(i+1)}$ 并相减, 有

$$\begin{array}{rl}& \log P\left(Y\mid {\theta }^{(i+1)}\right)-\log P\left(Y\mid {\theta }^{(i)}\right)\\ & =\left[Q\left({\theta }^{(i+1)},{\theta }^{(i)}\right)-Q\left({\theta }^{(i)},{\theta }^{(i)}\right)\right]-\left[H\left({\theta }^{(i+1)},{\theta }^{(i)}\right)-H\left({\theta }^{(i)},{\theta }^{(i)}\right)\right]\end{array}①$$

为证式 $P\left(Y\mid {\theta }^{(i+1)}\right)⩾P\left(Y\mid {\theta }^{(i)}\right)$, 只需证上式右端是非负的。上式右端的第 1 项, 由于 ${\theta }^{(i+1)}$ 使 $Q\left(\theta ,{\theta }^{(i)}\right)$ 达到极大, 所以有

$$Q\left({\theta }^{(i+1)},{\theta }^{(i)}\right)-Q\left({\theta }^{(i)},{\theta }^{(i)}\right)⩾0$$

其第 2 项, 由式 $H\left(\theta ,{\theta }^{(i)}\right)={\sum }_Z\log P(Z\mid Y,\theta )P\left(Z\mid Y,{\theta }^{(i)}\right)$ 可得:

$$\begin{array}{rl}H\left({\theta }^{(i+1)},{\theta }^{(i)}\right)-H\left({\theta }^{(i)},{\theta }^{(i)}\right)& =\sum _Z\left(\log \frac{P\left(Z\mid Y,{\theta }^{(i+1)}\right)}{P\left(Z\mid Y,{\theta }^{(i)}\right)}\right)P\left(Z\mid Y,{\theta }^{(i)}\right)\\ & ⩽\log \left(\sum _Z\frac{P\left(Z\mid Y,{\theta }^{(i+1)}\right)}{P\left(Z\mid Y,{\theta }^{(i)}\right)}P\left(Z\mid Y,{\theta }^{(i)}\right)\right)\\ & =\log \left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i+1)}\right)\right)=0\end{array}②$$

这里的不等号由 Jensen 不等式得到。
由式 $Q\left({\theta }^{(i+1)},{\theta }^{(i)}\right)-Q\left({\theta }^{(i)},{\theta }^{(i)}\right)⩾0$ 和式②即知式①右端是非负的。
定理 设 $L(\theta )=\log P(Y\mid \theta )$ 为观测数据的对数似然函数, ${\theta }^{(i)}(i=1,2,\cdots )$
为 $\mathrm{E}\mathrm{M}$ 算法得到的参数估计序列, $L\left({\theta }^{(i)}\right)(i=1,2,\cdots )$ 为对应的对数似然函数序列。

( 1 ）如果 $P(Y\mid \theta )$ 有上界, 则 $L\left({\theta }^{(i)}\right)=\log P\left(Y\mid {\theta }^{(i)}\right)$ 收敛到某一值 $L^{\ast }$;
( 2 ) 在函数 $Q\left(\theta ,{\theta }^{\prime }\right)$ 与 $L(\theta )$ 满足一定条件下, 由 $\mathrm{E}\mathrm{M}$ 算法得到的参数估计序列 ${\theta }^{(i)}$ 的收敘值 ${\theta }^{\ast }$ 是 $L(\theta )$ 的稳定点。

证明 \quad (1) 由 $L(\theta )=\log P\left(Y\mid {\theta }^{(i)}\right)$ 的单调性及 $P(Y\mid \theta )$的有界性立即得到。

---

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

---

9.3.1 高斯混合模型

---

高斯混合模型是指具有如下形式的概率分布模型:

$$P(y\mid \theta )=\sum _{k=1}^K{\alpha }_k\varphi \left(y\mid {\theta }_k\right)$$

其中, ${\alpha }_k$ 是系数, ${\alpha }_k⩾0,{\sum }_{k=1}^K{\alpha }_k=1;\varphi \left(y\mid {\theta }_k\right)$ 是高斯分布密度, ${\theta }_k=\left({\mu }_k,{\sigma }_k^2\right)$,

$$\varphi \left(y\mid {\theta }_k\right)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{{\left(y-{\mu }_k\right)}^2}{2{\sigma }_k^2}\right)①$$

称为第 $k$ 个分模型。
一般混合模型可以由任意概率分布密度代替式①中的高斯分布密度。

---

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

---

观测数据 $y_1,y_2,\cdots ,y_N$ 由高斯混合模型生成,

$$P(y\mid \theta )=\sum _{k=1}^K{\alpha }_k\varphi \left(y\mid {\theta }_k\right)$$

其中, $\theta =\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K;{\theta }_1,{\theta }_2,\cdots ,{\theta }_K\right)$ 。用 EM算法估计高斯混合模型的参数 $\theta$ 。

1. 明确隐变量, 写出完全数据的对数似然函数
   可以设想观测数据 $y_j,j=1,2,\cdots ,N$, 是这样产生的: 首先依概率 ${\alpha }_k$ 选择第 $k$ 个高斯分布分模型 $\varphi \left(y\mid {\theta }_k\right)$, 然后依第 $k$ 个分模型的概率分布 $\varphi \left(y\mid {\theta }_k\right)$ 生成观测数据 $y_j$ 。 这时观测数据 $y_j,j=1,2,\cdots ,N$, 是已知的; 反映观测数据 $y_j$ 来自第 $k$ 个分模 型的数据是未知的, $k=1,2,\cdots ,K$, 以隐变量 ${\gamma }_{jk}$ 表示, 其定义如下:
   $$\begin{array}{rlr}& {\gamma }_{jk}=\{\begin{array}{ll}1,& \text{ 第 }j\text{ 个观测来自第 }k\text{ 个分模型 }\\ 0,& \text{ 否则 }\end{array}& \\ & j=1,2,\cdots ,& N;k=1,2,\cdots ,K\end{array}$$

   ${\gamma }_{jk}$ 是 0-1 随机变量。
   有了观测数据 $y_j$ 及未观测数据 ${\gamma }_{jk}$, 那么完全数据是
   $$\left(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jK}\right),j=1,2,\cdots ,N$$

   完全数据的似然函数:
   ​$\begin{array}{rl}P(y,\gamma \mid \theta )& =\prod _{j=1}^{N}P\left(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jK}\mid \theta \right)\\ & =\prod _{k=1}^K\prod _{j=1}^N{\left[{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)\right]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\varphi \left(y_j\mid {\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_j-{\mu }_k\right)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\\ \text{ 式中, }{n}_k=\sum _{j=1}^N{\gamma }_{jk},\sum _{k=1}^K{n}_k=N& \end{array}$
   那么, 完全数据的对数似然函数为
   ​$\log P(y,\gamma \mid \theta )={\sum }_{k=1}^K\left\{n_k\log {\alpha }_k+{\sum }_{j=1}^N{\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\}$

2. EM 算法的 $\mathrm{E}$ 步: 确定 $Q$ 函数
   ​$\begin{array}{rl}Q\left(\theta ,{\theta }^{(i)}\right)& =E\left[\log P(y,\gamma \mid \theta )\mid y,{\theta }^{(i)}\right]\\ & =E\left\{\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\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}^N\left(E{\gamma }_{jk}\right)\log {\alpha }_k+\sum _{j=1}^N\left(E{\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\}\end{array}③$
   ​$E\left({\gamma }_{jk}\mid y,\theta \right)$, 记为 ${\hat{\gamma }}_{ik}$
   $$\begin{array}{rl}{\hat{\gamma }}_{jk}=& E\left({\gamma }_{jk}\mid y,\theta \right)=P\left({\gamma }_{jk}=1\mid y,\theta \right)\\ & =\frac{P\left({\gamma }_{jk}=1,y_j\mid \theta \right)}{{\sum }_{k=1}^{K}P\left({\gamma }_{jk}=1,y_j\mid \theta \right)}\\ =& \frac{P\left(y_j\mid {\gamma }_{jk}=1,\theta \right)P\left({\gamma }_{jk}=1\mid \theta \right)}{{\sum }_{k=1}^{K}P\left(y_j\mid {\gamma }_{jk}=1,\theta \right)P\left({\gamma }_{jk}=1\mid \theta \right)}\\ =& \frac{{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)}{{\sum }_{k=1}^K{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)},j=1,2,\cdots ,N;k=1,2,\cdots ,K\end{array}$$

   ${\hat{\gamma }}_{jk}$ 是在当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分模型的概率, 称为分模型 $k$ 对 观测数据 $y_j$ 的响应度。
   将 ${\hat{\gamma }}_{jk}=E{\gamma }_{jk}$ 及 $n_k={\sum }_{j=1}^{N}E{\gamma }_{jk}$ 代入式 $③$, 即得
   $$Q\left(\theta ,{\theta }^{(i)}\right)=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\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\}④$$

3. 确定 EM 算法的 M 步
   迭代的 $\mathrm{M}$ 步是求函数 $Q\left(\theta ,{\theta }^{(i)}\right)$ 对 $\theta$ 的极大值, 即求新一轮迭代的模型参数:
   $${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q\left(\theta ,{\theta }^{(i)}\right)$$

   用 ${\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$ 只需将
   式④分别对 ${\mu }_k,{\sigma }_k^2$ 求偏导数并令其为 0, 即可得到; 求 ${\hat{\alpha }}_k$ 是在 ${\sum }_{k=1}^K{\alpha }_k=1$ 条件 下求偏导数并令其为 0 得到的。结果如下:
   $$\begin{array}{r}{\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{n_k}{N}=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,\cdots ,K\end{array}$$

   备注：对${\sigma }_k^2$求导采用换元法，令${\sigma }_k^2$=$t$，最后求得结果
   用式$Q\left(\theta ,\theta (i)\right)={\sum }_{k=1}^K\left\{n_k\log {\alpha }_k+{\sum }_{j=1}^N{\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\}+\lambda ({\sum }_{k=1}^K{\alpha }_k-1)$对${\alpha }_k$求偏导，=>$\frac{n_k}{{\alpha }_k}+\lambda =0$=>${\alpha }_k=-\frac{n_k}{\lambda }$=>${\sum }_1^K{\alpha }_k={\sum }_1^K-\frac{n_k}{\lambda }=1$=>${\sum }_1^K{n}_k=-\lambda$∵${\sum }_{k=1}^K{n}_k=N$，∴$\lambda =-N$

   重复以上计算, 直到对数似然函数值不再有明显的变化为止。

(高斯混合模型参数估计的EM算法)
输入: 观测数据 $y_1,y_2,\cdots ,y_N$, 高斯混合模型; 输出：高斯混合模型参数。

1. 取参数的初始值开始迭代;

2. $\mathrm{E}$ 步: 依据当前模型参数, 计算分模型 $k$ 对观测数据 $y_j$ 的响应度
   
   $${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)}{{\sum }_{k=1}^K{\alpha }_k\varphi \left(y_j\mid {\theta }_k\right)},j=1,2,\cdots ,N;k=1,2,\cdots ,K$$

3. $\mathrm{M}$ 步: 计算新一轮迭代的模型参数
   $${\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$$

   $$\begin{array}{r}{\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$步, 直到收敛。

---

9.4 EM算法的推广

---

9.4.1 F函数的极大-极大算法

---

$\mathbf{F}$ 函数 假设隐变量数据 $Z$ 的概率分布为 $\tilde{P}(Z)$, 定义分布 $\tilde{P}$ 与参 数 $\theta$ 的函数 $F(\tilde{P},\theta )$ 如下:

$$F(\tilde{P},\theta )=E_{\tilde{P}}[\log P(Y,Z\mid \theta )]+H(\tilde{P})=F(\tilde{P},\theta )=\sum _Z[\tilde{P}(Z)\log P(Y,Z\mid \theta )]-\sum _Z[\tilde{P}(Z)\log \tilde{P}(Z)]$$

称为 $F$ 函数。式中 $H(\tilde{P})=-E_{\tilde{P}}\log \tilde{P}(Z)$ 是分布 $\tilde{P}(Z)$ 的熵。

引理$1$ 对于固定的 $\theta$, 存在唯一的分布 ${\tilde{P}}_{\theta }$ 极大化 $F(\tilde{P},\theta )$, 这时 ${\tilde{P}}_{\theta }$ 由下式给出:

$${\tilde{P}}_{\theta }(Z)=P(Z\mid Y,\theta )$$

并且 ${\tilde{P}}_{\theta }$ 随 $\theta$ 连续变化。

证明 \quad 对于固定的 $\theta$, 可以求得使 $F(\tilde{P},\theta )$ 达到极大的分布 ${\tilde{P}}_{\theta }(Z)$ 。为此, 引进拉 格朗日乘子 $\lambda$, 拉格朗日函数为

$$L=E_{\tilde{P}}\log P(Y,Z\mid \theta )-E_{\tilde{P}}\log \tilde{P}(Z)+\lambda \left(1-\sum _Z\tilde{P}(Z)\right)$$

将其对 $\widetilde{P_Z}$ 求偏导数:

$$\frac{\partial L}{\partial \tilde{P}(Z)}=\log P(Y,Z\mid \theta )-\log \tilde{P}(Z)-1-\lambda$$

令偏导数等于 0 , 得出

$$\lambda =\log P(Y,Z\mid \theta )-\log {\tilde{P}}_{\theta }(Z)-1$$

由此推出 ${\tilde{P}}_{\theta }(Z)$ 与 $P(Y,Z\mid \theta )$ 成比例

$$\frac{P(Y,Z\mid \theta )}{{\tilde{P}}_{\theta }(Z)}={\mathrm{e}}^{1+\lambda }$$

再从约束条件 ${\sum }_Z{\tilde{P}}_{\theta }(Z)=1$ ,KaTeX parse error: Expected 'EOF', got '&' at position 50: …e}^{1+\lambda}=&̲gt;P(Y \mid \th…,$P(Z)=\frac{P(Y,Z\mid \theta )}{{\mathrm{e}}^{1+\lambda }}=\frac{P(Y,Z\mid \theta )}{P(Y\mid \theta )}=\frac{P(Z\mid Y,\theta )\cdot P(Y\mid \theta )}{P(Y\mid \theta )}=P(Z\mid Y,\theta )$

引理 $2$ 若 ${\tilde{P}}_{\theta }(Z)=P(Z\mid Y,\theta )$, 则$F(\tilde{P},\theta )=logP(Y\mid \theta )$

$$\begin{array}{rl}F(\tilde{P},\theta )& =E_{\tilde{p}}[logP(Z|Y,\theta )\cdot P(Y|\theta )]-E_{\tilde{p}}log\tilde{P}(Z)\\ & =E_{\tilde{p}}[logP(Z|Y,\theta )+logP(Y|\theta )]-E_{\tilde{p}}log\tilde{P}(Z)\\ & =E_{\tilde{p}}logP(Z|Y,\theta )+E_{\tilde{p}}logP(Y|\theta )-E_{\tilde{p}}log\tilde{P}(Z)\\ & =E_{\tilde{p}}log\tilde{P}(Z)+E_{\tilde{p}}logP(Y|\theta )-E_{\tilde{p}}log\tilde{P}(Z)\\ & =E_{\tilde{p}}logP(Y|\theta )\\ & =\sum _{Z}P(Z)logP(Y|\theta )\\ & =logP(Y|\theta )\end{array}$$

定理 3设 $L(\theta )=\log P(Y\mid \theta )$ 为观测数据的对数似然函数, ${\theta }^{(i)},i=1,2,\cdots$,为 $\mathrm{E}\mathrm{M}$ 算法得到的参数估计序列, 函数 $F(\tilde{P},\theta )$ 由式 $(9.33)$ 定义。如果 $F(\tilde{P},\theta )$ 在 ${\tilde{P}}^{\ast }$ 和 ${\theta }^{\ast }$ 有局部极大值, 那么 $L(\theta )$ 也在 ${\theta }^{\ast }$ 有局部极大值。类似地, 如果 $F(\tilde{P},\theta )$ 在 ${\tilde{P}}^{\ast }$和 ${\theta }^{\ast }$ 达到全局最大值, 那么 $L(\theta )$ 也在 ${\theta }^{\ast }$ 达到全局最大值。

证明 由引理 $1$ 和引理 $2$ 可知, $L(\theta )=\log P(Y\mid \theta )=F\left({\tilde{P}}_{\theta },\theta \right)$ 对任意 $\theta$ 成立。
特别地, 对于使 $F(\tilde{P},\theta )$ 达到极大的参数 ${\theta }^{\ast }$, 有

$$L\left({\theta }^{\ast }\right)=F\left({\tilde{P}}_{\theta }\cdot ,{\theta }^{\ast }\right)=F\left({\tilde{P}}^{\ast },{\theta }^{\ast }\right)$$

为了证明 ${\theta }^{\ast }$ 是 $L(\theta )$ 的极大点, 需要证明不存在接近 ${\theta }^{\ast }$ 的点 ${\theta }^{\ast \ast }$, 使 $L\left({\theta }^{\ast \ast }\right)>L\left({\theta }^{\ast }\right)$ 。 假如存在这样的点 ${\theta }^{\ast \ast }$, 那么应有 $F\left({\tilde{P}}^{\ast \ast },{\theta }^{\ast \ast }\right)>F\left({\tilde{P}}^{\ast },{\theta }^{\ast }\right)$, 这里 ${\tilde{P}}^{\ast \ast }={\tilde{P}}_{\theta }\cdots$ 但因 ${\tilde{P}}_{\theta }$是随 $\theta$ 连续变化的, ${\tilde{P}}^{\ast \ast }$ 应接近 ${\tilde{P}}^{\ast }$, 这与 ${\tilde{P}}^{\ast }$ 和 ${\theta }^{\ast }$ 是 $F(\tilde{P},\theta )$ 的局部极大点的假设 矛盾。

定理 4 EM 算法的一次迭代可由 $F$ 函数的极大-极大算法实现. 设 ${\theta }^{(i)}$ 为第 $i$ 次迭代参数 $\theta$ 的估计, ${\tilde{P}}^{(i)}$ 为第 $i$ 次迭代函数 $\tilde{P}$ 的估计。在第 $i+1$

次迭代的两步为:(1) 对固定的 ${\theta }^{(i)}$, 求 ${\tilde{P}}^{(i+1)}$ 使 $F\left(\tilde{P},{\theta }^{(i)}\right)$ 极大化;( 2 ) 对固定的 ${\tilde{P}}^{(i+1)}$, 求 ${\theta }^{(i+1)}$ 使 $F\left({\tilde{P}}^{(i+1)},\theta \right)$ 极大化.
证明 $(1)$ 由引理 $1$, 对于固定的 ${\theta }^{(i)}$,

$${\tilde{P}}^{(i+1)}(Z)={\tilde{P}}_{{\theta }^{(i)}}(Z)=P\left(Z\mid Y,{\theta }^{(i)}\right)$$

使 $F\left(\tilde{P},{\theta }^{(i)}\right)$ 极大化。此时,

$$\begin{array}{rl}F\left({\tilde{P}}^{(i+1)},\theta \right)& =E_{{\tilde{P}}^{(i+1)}}[\log P(Y,Z\mid \theta )]+H\left({\tilde{P}}^{(i+1)}\right)\\ & =\sum _Z\log P(Y,Z\mid \theta )P\left(Z\mid Y,{\theta }^{(i)}\right)+H\left({\tilde{P}}^{(i+1)}\right)\end{array}$$

由 $Q\left(\theta ,{\theta }^{(i)}\right)$ 的定义式 $Q\left(\theta ,{\theta }^{(i)}\right)=E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]$ 有

$$F\left({\tilde{P}}^{(i+1)},\theta \right)=Q\left(\theta ,{\theta }^{(i)}\right)+H\left({\tilde{P}}^{(i+1)}\right)$$

(2) 固定 ${\tilde{P}}^{(i+1)}$, 求 ${\theta }^{(i+1)}$ 使 $F\left({\tilde{P}}^{(i+1)},\theta \right)$ 极大化。得到

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

---

9.4.2 GEM 算法

---

GEM 算法 1
输入: 观测数据, $F$ 函数; 输出：模型参数。
(1) 初始化参数 ${\theta }^{(0)}$, 开始迭代:
(2) 第 $i+1$ 次迭代, 第 1 步: 记 ${\theta }^{(i)}$ 为参数 $\theta$ 的估计值, ${\tilde{P}}^{(i)}$ 为函数 $\tilde{P}$ 的估计,
求 ${\tilde{P}}^{(i+1)}$ 使 $\tilde{P}$ 极大化 $F\left(\tilde{P},{\theta }^{(i)}\right);$
(3)第 2 步: 求 ${\theta }^{(i+1)}$ 使 $F\left({\tilde{P}}^{(i+1)},\theta \right)$ 极大化;
(4) 重复
(2) 和 $(3)$, 直到收敛。
在 GEM算法 1 中, 有时求 $Q\left(\theta ,{\theta }^{(i)}\right)$ 的极大化是很困难的。下面介绍的 $\mathrm{G}\mathrm{E}\mathrm{M}$ 算
法 2 和 GEM算法 3 并不是直接求 ${\theta }^{(i+1)}$ 使 $Q\left(\theta ,{\theta }^{(i)}\right)$ 达到极大的 $\theta$, 而是找一个 ${\theta }^{(i+1)}$ 使得 $Q\left({\theta }^{(i+1)},{\theta }^{(i)}\right)>Q\left({\theta }^{(i)},{\theta }^{(i)}\right).$
GEM 算法 2
输入: 观测数据, $Q$ 函数; 输出: 模型参数。
(1) 初始化参数 ${\theta }^{(0)}$, 开始迭代:
(2) 第 $i+1$ 次迭代, 第 1 步: 记 ${\theta }^{(i)}$ 为参数 $\theta$ 的估计值, 计算

$$\begin{array}{rl}Q\left(\theta ,{\theta }^{(i)}\right)& =E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]\\ & =\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\log P(Y,Z\mid \theta )\end{array}$$

(3) 第 2 步: 求 ${\theta }^{(i+1)}$ 使

$$Q\left({\theta }^{(i+1)},{\theta }^{(i)}\right)>Q\left({\theta }^{(i)},{\theta }^{(i)}\right)$$

(2) 和 $(3)$, 直到收敛。
(4) 重复
当参数 $\theta$ 的维数为 $d(d⩾2)$ 时, 可采用一种特殊的 GEM算法, 它将 EM算法的
M 步分解为 $d$ 次条件极大化, 每次只改变参数向量的一个分量, 其余分量不改变。

GEM 算法 3
输入: 观测数据, $Q$ 函数;
输出: 模型参数。
(1)初始化参数 ${\theta }^{(0)}=\left({\theta }_1^{(0)},{\theta }_2^{(0)},\cdots ,{\theta }_d^{(0)}\right)$, 开始迭代;
(2) 第 $i+1$ 次迭代, 第 1 步: 记 ${\theta }^{(i)}=\left({\theta }_1^{(i)},{\theta }_2^{(i)},\cdots ,{\theta }_d^{(i)}\right)$ 为参数 $\theta =\left({\theta }_1,{\theta }_2,\cdots ,{\theta }_d\right)$
的估计值, 计算

$$\begin{array}{rl}Q\left(\theta ,{\theta }^{(i)}\right)& =E_Z\left[\log P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]\\ & =\sum _{Z}P\left(Z\mid y,{\theta }^{(i)}\right)\log P(Y,Z\mid \theta )\end{array}$$

(3) 第 2 步: 进行 $d$ 次条件极大化:
首先, 在 ${\theta }_2^{(i)},\cdots ,{\theta }_d^{(i)}$ 保持不变的条件下求使 $Q\left(\theta ,{\theta }^{(i)}\right)$ 达到极大的 ${\theta }_1^{(i+1)}$;
然后, 在 ${\theta }_1={\theta }_1^{(i+1)},{\theta }_j={\theta }_j^{(i)},j=3,4,\cdots ,d$ 的条件下求使 $Q\left(\theta ,{\theta }^{(i)}\right)$ 达到极大
的 ${\theta }_2^{(i+1)};$
如此继续, 经过 $d$ 次条件极大化, 得到 ${\theta }^{(i+1)}=\left({\theta }_1^{(i+1)},{\theta }_2^{(i+1)},\cdots ,{\theta }_d^{(i+1)}\right)$ 使得

$$Q\left({\theta }^{(i+1)},{\theta }^{(i)}\right)>Q\left({\theta }^{(i)},{\theta }^{(i)}\right)$$

(4)重复
(2) 和 ( 3 ), 直到收敛。