EM 算法

来自统计学习方法

还可以参考博客：https://blog.csdn.net/v_JULY_v/article/details/81708386

知乎 https://www.zhihu.com/question/40797593

EM算法理解的九层境界

EM 就是 E + M

EM 是一种局部下限构造

K-Means是一种Hard EM算法

从EM 到广义EM

广义EM的一个特例是VBEM

广义EM的另一个特例是WS算法

广义EM的再一个特例是Gibbs抽样算法

WS算法是VAE和GAN组合的简化版

KL距离的统一

EM 算法的引入

书中提到了三硬币模型，

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

假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数。这个模型中，不能观测的 A 的结果就是隐变量。

将观测数据表示为 $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 \max _{\theta} \log P(Y \mid \theta) \]

这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法。下面给出针对以上问题的 EM算法，其推导过程省略。

$\mathrm{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}}}\tag{9.5} \]

$\mathrm{M}$ 步：计算模型参数的新估计值

\[\pi^{(i+1)}=\frac{1}{n} \sum_{j=1}^{n} \mu_{j}^{(i+1)}\tag{9.6} \]

...

\[p^{(i+1)}=\frac{\sum_{j=1}^{n} \mu_{j}^{(i+1)} y_{j}}{\sum_{j=1}^{n} \mu_{j}^{(i+1)}} \tag{9.7} \]

...

\[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)} \tag{9.8} \]

进行数字计算。假设模型参数的初值取为

\[\pi^{(0)}=0.5, \quad p^{(0)}=0.5, \quad q^{(0)}=0.5 \]

由式 (9.5), 对 $y_{j}=1$ 与 $y_{j}=0$ 均有 $\mu_{j}^{(1)}=0.5$ 。

利用迭代公式 $(9.6) \sim$ 公式 (9.8), 得到

\[\pi^{(1)}=0.5, \quad p^{(1)}=0.6, \quad q^{(1)}=0.6 \]

由式 (9.5),

\[\mu_{j}^{(2)}=0.5, \quad j=1,2, \cdots, 10 \]

继续迭代, 得

\[\pi^{(2)}=0.5, \quad p^{(2)}=0.6, \quad q^{(2)}=0.6 \]

于是得到模型参数 $\theta$ 的极大似然估计:

\[\hat{\pi}=0.5, \quad \hat{p}=0.6, \quad \hat{q}=0.6 \]

$\pi=0.5$ 表示硬币 $\mathrm{A}$ 是均匀的, 这一结果容易理解。

如果取初值 $\pi^{(0)}=0.4, p^{(0)}=0.6, q^{(0)}=0.7$, 那么得到的模型参数的极大似然估计是 $\hat{\pi}=0.4064, \hat{p}=0.5368, \hat{q}=0.6432$ 。这就是说, $\mathrm{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)$ 。

EM算法通过迭代求 $的极大似然估计。每次迭代包含两步:E步，求期望；M步，求极大化。下面来介绍EM算法。

...

定义 $9.1(\boldsymbol{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] \tag{9.11} \]

下面关于 $\mathrm{EM}$ 算法作几点说明：

步骤 (1) 参数的初值可以任意选择，但需注意 EM算法对初值是敏感的。

步骤（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)} \rightarrow \theta^{(i+1)}$ 。后面将证明每次迭代使似然函数增大或达到局部极值。

步骤 (4) 给出停止迭代的条件, 一般是对较小的正数 $\varepsilon_{1}, \varepsilon_{2}$, 若满足

\[\left\|\theta^{(i+1)}-\theta^{(i)}\right\|<\varepsilon_{1} \quad \text { 或 }\left\|Q\left(\theta^{(i+1)}, \theta^{(i)}\right)-Q\left(\theta^{(i)}, \theta^{(i)}\right)\right\|<\varepsilon_{2} \]

则停止迭代。

Jensen 不等式

设f是定义域为实数的函数

如果对于所有的实数x，f(x)的二次导数 $f''(x) \geq 0$，那么f是凸函数。
当x是向量时，如果其hessian矩阵H是半正定的（ $\bold{H} \geq 0$），那么f是凸函数。
如果 $f''(x) > 0$ 或者 $\bold{H} > 0$，那么称f是严格凸函数。

Jensen不等式表述如下：

如果 $f$ 是凸函数，X是随机变量，那么：$E[f(X)]>=f(E[X])$，通俗的说法是函数的期望大于等于期望的函数。
特别地，如果f是严格凸函数，当且仅当 $P(X = EX) = 1$，即X是常量时，上式取等号，即 $E[f(X)] = f(EX)$。

EM 算法的导出

我们面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据)Y关于参数0的对数似然函数，即极大化

\[\begin{aligned} 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{aligned} \tag{9.12} \]

注意到这一极大化的主要困难是式 (9.12) 中有末观测数据并有包含和 (或积分) 的对数。

事实上, EM 算法是通过迭代逐步近似极大化 $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 \geq \lambda_j \sum_j log y_j$，其中 $\lambda_j \geq 0,\ \sum_j \lambda_j =1$

\[\begin{aligned} 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) \\ & \geqslant \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)}-\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{aligned} \]

令：

\[B\left(\theta, \theta^{(i)}\right) \hat{=} 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)}\tag{9.13} \]

则

\[L(\theta) \geqslant B\left(\theta, \theta^{(i)}\right)\tag{9.14} \]

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

\[L\left(\theta^{(i)}\right)=B\left(\theta^{(i)}, \theta^{(i)}\right) \tag{9.15} \]

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

\[\theta^{(i+1)}=\arg \max _{\theta} B\left(\theta, \theta^{(i)}\right) \tag{9.16} \]

现在求 $\theta^{(i+1)}$ 的表达式。省去对 $\theta$ 的极大化而言是常数的项, 由式 $(9.16) 、$ 式 (9.13) 及式 $(9.10)$, 有

\[\begin{aligned} \theta^{(i+1)} &=\arg \max _{\theta}\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 \max _{\theta}\left(\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log (P(Y \mid Z, \theta) P(Z \mid \theta))\right) \\ &=\arg \max _{\theta}\left(\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log P(Y, Z \mid \theta)\right) \\ &=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right) \end{aligned} \tag{9.17} \]

式 (9.17) 等价于 EM 算法的一次迭代, 即求 $Q$ 函数及其极大化。EM 算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

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

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

监督学习是由训练数据 $\left\{\left(x_{1}, y_{1}\right),\left(x_{2}, y_{2}\right), \cdots,\left(x_{N}, y_{N}\right)\right\}$ 学习条件概率分布 $P(Y \mid X)$ 或决策函数 $Y=f(X)$ 作为模型, 用于分类、回归、标注等任务。这时训练数据中的每个样本点由输入和输出对组成。

有时训练数据只有输入没有对应的输出 $\left\{\left(x_{1}, \cdot\right),\left(x_{2}, \cdot\right), \cdots,\left(x_{N}, \cdot\right)\right\}$, 从这样的数据学习模型称为无监督学习问题。EM算法可以用于生成模型的无监督学习。生成模型由联合概率分布 $P(X, Y)$ 表示, 可以认为无监督学习训练数据是联合概率分布产生的数据。 $X$ 为观测数据, $Y$ 为末观测数据。

EM 算法的收敛性

书上给出两个定理，这些定理的证明可以看《统计学习方法，第二版》P181、182

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。EM算法的最大优点是简单性和普适性。我们很自然地要问:EM算法得到的估计序列是否收敛?如果收敛，是否收敛到全局最大值或局部极大值?下面给出关于EM算法收敛性的两个定理。

定理 9.1 设 $P(Y \mid \theta)$ 为观测数据的似然函数, $\theta^{(i)}(i=1,2, \cdots)$ 为 $\mathrm{EM}$ 算法得到的参数估计序列, $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) \geqslant P\left(Y \mid \theta^{(i)}\right) \tag{9.18} \]

定理 9.2 设 $L(\theta)=\log P(Y \mid \theta)$ 为观测数据的对数似然函数，$\theta^{(i)}(i=1,2, \cdots)$ 为 $\mathrm{EM}$ 算法得到的参数估计序列，$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^{*}$;

(2) 在函数 $Q\left(\theta, \theta^{\prime}\right)$ 与 $L(\theta)$ 满足一定条件下, 由 $\mathrm{EM}$ 算法得到的参数估计序列 $\theta^{(i)}$ 的收敛值 $\theta^{*}$ 是 $L(\theta)$ 的稳定点。

定理9.2关于函数 $Q(\theta,\theta^{(i)})$ 与 $L(\theta)$ 的条件在大多数情况下都是满足的。EM算法的收敛性包含关于对数似然函数序列 $L(\theta^{(i)})$ 的收敛性和关于参数估计序列 $\theta^{(i)}$ 的收敛性两层意思，前者并不蕴涵后者。此外，定理只能保证参数估计序列收敛到对数似然函数序列的稳定点，不能保证收敛到极大值点。所以在应用中，初值的选择变得非常重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

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

EM算法的一个重要应用是高斯混合模型的参数估计。高斯混合模型应用广泛，在许多情况下，EM算法是学习高斯混合模型(Gaussian mixture model)的有效方法。

高斯混合模型

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

\[P(y \mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right) \tag{9.24} \]

其中, $\alpha_{k}$ 是系数, $\alpha_{k} \geqslant 0, \sum_{k=1}^{K} \alpha_{k}=1 ; \phi\left(y \mid \theta_{k}\right)$ 是高斯分布密度, $\theta_{k}=\left(\mu_{k}, \sigma_{k}^{2}\right)$,

\[\phi\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)\tag{9.25} \]

称为第 $k$ 个分模型。

一般混合模型可以由任意概率分布密度代替式 (9.25) 中的高斯分布密度, 我们只介绍最常用的高斯混合模型。

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

假设观测数据 $y_{1}, y_{2}, \cdots, y_{N}$ 由高斯混合模型生成,

\[P(y \mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right) \tag{9.26} \]

其中, $\theta=\left(\alpha_{1}, \alpha_{2}, \cdots, \alpha_{K} ; \theta_{1}, \theta_{2}, \cdots, \theta_{K}\right)$ 。我们用 $\mathrm{EM}$ 算法估计高斯混合模型的参数 $\theta$ 。

1. 明确隐变量, 写出完全数据的对数似然函数

可以设想观测数据 $y_{j}, j=1,2, \cdots, N$, 是这样产生的: 首先依概率 $\alpha_{k}$ 选择第 $k$ 个高斯分布分模型 $\phi\left(y \mid \theta_{k}\right)$, 然后依第 $k$ 个分模型的概率分布 $\phi\left(y \mid \theta_{k}\right)$ 生成观测数据 $y_{j}$ 。这时观测数据 $y_{j}, j=1,2, \cdots, N$, 是已知的; 反映观测数据 $y_{j}$ 来自第 $k$ 个分模型的数据是末知的, $k=1,2, \cdots, K$, 以隐变量 $\gamma_{j k}$ 表示, 其定义如下:

\[\begin{aligned} &\gamma_{j k}= \begin{cases}1, & \text { 第 } j \text { 个观测来自第 } k \\ 0, & \text { 否则 }\end{cases} \\ &j=1,2, \cdots, N ; \quad k=1,2, \cdots, K \end{aligned} \tag{9.27} \]

$\gamma_{j k}$ 是 0-1 随机变量。

有了观测数据 $y_{j}$ 及末观测数据 $\gamma_{j k}$, 那么完全数据是

\[\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K}\right), \quad j=1,2, \cdots, N \]

于是, 可以写出完全数据的似然函数:

\[\begin{aligned} P(y, \gamma \mid \theta) &=\prod_{j=1}^{N} P\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K} \mid \theta\right) \\ &=\prod_{k=1}^{K} \prod_{j=1}^{N}\left[\alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)\right]^{\gamma_{j k}} \\ &=\prod_{k=1}^{K} \alpha_{k}^{n_{k}} \prod_{j=1}^{N}\left[\phi\left(y_{j} \mid \theta_{k}\right)\right]^{\gamma_{j k}} \\ &=\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_{j k}} \end{aligned} \]

式中, $n_{k}=\sum_{j=1}^{N} \gamma_{j k}, \sum_{k=1}^{K} n_{k}=N$ 。

那么, 完全数据的对数似然函数为

\[\log P(y, \gamma \mid \theta)=\sum_{k=1}^{K}\left\{n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\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 算法的 E 步: 确定 Q 函数

\[\begin{aligned} 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_{j k}\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_{j k}\right) \log \alpha_{k}+\sum_{j=1}^{N}\left(E \gamma_{j k}\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{aligned} \tag{9.28} \]

这里需要计算 $E\left(\gamma_{j k} \mid y, \theta\right)$, 记为 $\hat{\gamma}_{j k}$ 。

\[\begin{aligned} \hat{\gamma}_{j k} &=E\left(\gamma_{j k} \mid y, \theta\right)=P\left(\gamma_{j k}=1 \mid y, \theta\right) \\ &=\frac{P\left(\gamma_{j k}=1, y_{j} \mid \theta\right)}{\sum_{k=1}^{K} P\left(\gamma_{j k}=1, y_{j} \mid \theta\right)} \\ &=\frac{P\left(y_{j} \mid \gamma_{j k}=1, \theta\right) P\left(\gamma_{j k}=1 \mid \theta\right)}{\sum_{k=1}^{K} P\left(y_{j} \mid \gamma_{j k}=1, \theta\right) P\left(\gamma_{j k}=1 \mid \theta\right)} \\ &=\frac{\alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K \end{aligned} \]

$\hat{\gamma}_{j k}$ 是在当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分模型的概言, 称为分模型 $k$ 对观测数据 $y_{j}$ 的响应度。

将 $\hat{\gamma}_{j k}=E \gamma_{j k}$ 及 $n_{k}=\sum_{j=1}^{N} E \gamma_{j k}$ 代入式 (9.28), 即得

\[Q\left(\theta, \theta^{(i)}\right)=\sum_{k=1}^{K}\left\{n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \hat{\gamma}_{j k}\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\} \tag{9.29} \]

3. 确定 EM 算法的 M 步

迭代的 $\mathrm{M}$ 步是求函数 $Q\left(\theta, \theta^{(i)}\right)$ 对 $\theta$ 的极大值, 即求新一轮迭代的模型参数:

\[\theta^{(i+1)}=\arg \max _{\theta} 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}$ 只需将式 (9.29) 分别对 $\mu_{k}, \sigma_{k}^{2}$ 求偏导数并令其为 0 , 即可得到; 求 $\hat{\alpha}_{k}$ 是在 $\sum_{k=1}^{K} \alpha_{k}=1$ 条件下求偏导数并令其为 0 得到的。结果如下:

\[\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K \tag{9.30} \]

\[\hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K \tag{9.31} \]

\[\hat{\alpha}_{k}=\frac{n_{k}}{N}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K \tag{9.32} \]

重复以上计算, 直到对数似然函数值不再有明显的变化为止。
现将什计高斯混合模型参数的 EM 算法总结如下。

算法 9.2 (高斯混合模型参数估计的EM算法)
输入: 观测数据 $y_{1}, y_{2}, \cdots, y_{N}$, 高斯混合模型;
输出: 高斯混合模型参数。
(1) 取参数的初始值开始迭代:
(2) $\mathrm{E}$ 步: 依据当前模型参数, 计算分模型 $k$ 对观测数据 $y_{j}$ 的响应度

\[\hat{\gamma}_{j k}=\frac{\alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K \]

(3) $\mathrm{M}$ 步: 计算新一轮迭代的模型参数

\[\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K \]

\[\hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K \]

\[\hat{\alpha}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K \]

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

EM 算法的推广

$\mathrm{EM}$ 算法还可以解释为 $F$ 函数（F function）的极大-极大算法 (maximizationmaximization algorithm），基于这个解释有若干变形与推广, 如广义期望极大 (generalized expectation maximization, GEM) 算法。下面予以介绍。

F函数的极大-极大算法

首先引进 $F$ 函数并讨论其性质。

定义 9.3 ( $\boldsymbol{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$ 函数。式中 $H(\tilde{P})=-E_{\tilde{P}} \log \tilde{P}(Z)$ 是分布 $\tilde{P}(Z)$ 的熵。

在定义 9.3 中, 通常假设 $P(Y, Z \mid \theta)$ 是 $\theta$ 的连续函数, 因而 $F(\tilde{P}, \theta)$ 是 $\tilde{P}$ 和 $\theta$ 的连续函数。函数 $F(\tilde{P}, \theta)$ 还有以下重要性质。

引理 9.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$ 连续变化。

证明，见《统计信息方法第二版》P188

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

\[F(\tilde{P}, \theta)=\log P(Y \mid \theta) \tag{9.36} \]

《统计信息方法第二版》让自己证

由以上引理, 可以得到关于 EM算法用 $F$ 函数的极大-极大算法的解释。

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

证明，见《统计信息方法第二版》P188、189

定理 9.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)$ 极大化。

通过以上两步完成了EM算法的一次迭代。由此可知，由EM算法与F函数的极大-极大算法得到的参数估计序列 $\theta^{(i)}$，i= 1,2,…·，是一致的。

这样，就有EM算法的推广。

GEM 算法

算法 9.3 (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)$ 的极大化是很困难的。下面介绍的 GEM算. 法 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)$ 。

算法 9.4 (GEM 算法 2)

输入: 观测数据, $Q$ 函数;
输出: 模型参数。

(1) 初始化参数 $\theta^{(0)}$, 开始迭代;
(2) 第 $i+1$ 次迭代, 第 1 步: 记 $\theta^{(i)}$ 为参数 $\theta$ 的估计值, 计算

\[\begin{aligned} 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{aligned} \]

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

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

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

算法 9.5 (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{aligned} 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{aligned} \]

(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)$ 使得