狭义EM

1. 狭义EM

Reference：[统计学习方法 第二版 第九章]

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

1.1 EM 算法的引入

概率模型有时既含有观测变量（observable variable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数（此时只有固定的参数未知）。EM 算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。

例1（三硬币模型）：假设有 3 枚硬币，分别记作 A，B，C。这些硬币正面出现的概率分别是 $\pi$，$p$，$q$。进行如下掷硬币试验：先掷硬币 A，根据其结果选出硬币 B 或硬币 C，正面选硬币 B，反面选硬币 C；然后掷选出的硬币，掷硬币的结果，出现正面记作 1，出现反面记作 0；独立地重复 $n$ 次试验（这里，$n=10$），观测结果如下，如何估计三硬币模型的参数： $1,1,0,1,0,0,1,0,1,1$ 三硬币模型可以写作 P(y\mid \theta)=\sum_z P(y,z\mid\theta)=\sum_zP(z\mid\theta)P(y\mid z,\theta)=\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y} \tag{1} 这里随机变量 $y$ 是观测变量，表示一次试验观测到的结果是 1 或 0；随机变量 $z$ 是隐变量，表示未观测到的掷硬币 A 的结果；$\theta=(\pi,p,q)$ 是模型参数。则所有观测数据的似然函数为 P(Y\mid \theta)=\sum_ZP(Z\mid \theta)P(Y\mid Z, \theta)=\prod_{j=1}^n[\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]\tag{2} 考虑求模型参数 $\theta=(\pi, p, q)$ 的极大似然估计，即 \hat\theta=\mathrm{arg} \max_\theta \log P(Y\mid \theta)\tag{3} 这个问题没有解析解，因为 $\log$ 中有相加的两项，只有通过迭代的方法求解。

E 步：计算在模型参数 $\pi^{(i)},p^{(i)},q^{(i)}$ 下观测数据 $y_j$ 来自掷硬币 B 的概率（将上一步的参数作为已知值，在算法一中应该为M步） \mu_j^{(i+1)}=\frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}+(1-\pi^{(i)})(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}}\tag{4} M 步：计算模型参数的新估计值

$$\begin{aligned} & \pi^{(i+1)}=\frac{1}{n}\sum_{j=1}^n \mu_j^{(i+1)} \\ & 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 (1-\mu_j^{(i+1)})y_j}{\sum_{j=1}^n (1-\mu_j^{(i+1)})} \end{aligned}\tag{5}$$

按照上述迭代步骤直至收敛，若假设模型参数初值为 $\pi^{(0)}=0.5, p^{(0)}=0.5, q^{(0)}=0.5$，则模型参数的极大似然估计为 $\hat\pi =0.5, \hat p =0.6, \hat q= 0.6$。若假设模型参数初值为 $\pi^{(0)}=0.4, p^{(0)}=0.6, q^{(0)}=0.6$，则模型参数的极大似然估计为 $\hat\pi =0.4064, \hat p =0.5368, \hat q= 0.6432$。

一般地，用 $Y$ 表示观测随机变量的数据，$Z$ 表示隐随机变量的数据。$Y$ 和 $Z$ 连在一起称为完全数据（complete-data），观测数据 $Y$ 又称为不完全数据（incomplete-data）。

首先写出所有观测的似然函数 P(Y \mid \theta) = \prod_{j=1}^n P(y_j \mid \theta) = \prod_{j=1}^n [\pi p + (1-\pi) q]^{y_j} [\pi (1-p)+(1-\pi)(1-q)]^{(1-y_j)} \tag{6} 计算 $P(z_j =1 \mid y_j, \theta^{(i)})$： \mu_j^{(i+1)} = P(z_j = 1 \mid y_j, \theta^{(i)}) = \frac{P(y_j \mid z_j =1,\theta^{(i)})P(z_j =1 \mid \theta^{(i)})}{P(y_j \mid \theta^{(i)})}=\begin{cases} \frac{\pi^{(i)} p^{(i)}}{\pi^{(i)} p ^{(i)} + (1-\pi^{(i)})q^{(i)}} &\quad \text{if } y_j=1\\ \frac{\pi^{(i)} (1-p^{(i)})}{\pi^{(i)} (1-p^{(i)})+ (1-\pi^{(i)})(1-q^{(i)})} &\quad \text{if } y_j=0 \end{cases} \tag{7} 计算完全数据的对数似然函数的期望

$$\begin{aligned} Q(\theta \mid \theta^{(i)}) = &\mathbb{E}_{P(Z \mid Y, \theta^{(i)})}[\log P(Y,Z \mid \theta)] \\ =& \mathbb{E}_{P(Z\mid Y ,\theta^{(i)})} [\sum_{j=1}^n \log P(y_j, z_j \mid \theta)]\\ =& \sum_{j=1}^n \mathbb{E}_{P(z_j \mid y_j , \theta^{(i)})}[\log P(y_j, z_j \mid \theta)]\\ =& \sum_{j=1}^n \sum_{z_j} P(z_j \mid y_j ,\theta^{(i)}) \log P(y_j, z_j \mid \theta)\\ =& \sum_{j=1}^n \left[\mu_j^{(i+1)} \log \left(\pi p^{y_j} (1-p)^{(1-y_j)}\right) + (1-\mu_j^{(i+1)}) \log \left((1-\pi)q^{y_j}(1-q)^{(1-y_j)}\right)\right] \\ \end{aligned}\tag{8}$$

对各参数求导，并令其满足一阶条件可得公式（5）（可拆解验证其凹性）。

算法1 （EM 算法）： 输入：观测变量数据 $Y$，隐变量数据 $Z$，联合分布 $P(Y,Z\mid\theta)$，条件分布 $P(Z\mid Y, \theta)$； 输出：模型参数 $\theta$。 （1）选择参数的初值 $\theta^{(0)}$，开始迭代； （2）E 步：记 $\theta^{(i)}$ 为第 $i$ 次迭代参数的估计值，在第 $i+1$ 次迭代的 E 步，计算 Q(\theta, \theta^{(i)})=E_Z[\log P(Y, Z \mid \theta) \mid Y, \theta^{(i)}]=\sum_Z P(Z\mid Y, \theta^{(i)})\log P(Y,Z\mid \theta)\tag{9} （3）M 步：求使 $Q(\theta,\theta^{(i)})$ 极大化的 $\theta$，确定第 $i+1$ 次迭代的参数的估计值 $\theta^{(i+1)}$ \theta^{(i+1)} =\mathrm{arg} \max_{\theta} Q(\theta, \theta^{(i)})\tag{10} （4）重复第（2）步和第（3）步，直到收敛。

注意：参数的初值可以任意选择，但 EM 算法对初值是敏感的。

定义1：$Q$ 函数（$Q$ function） 完全数据的对数似然函数 $\log P(Y,Z \mid \theta)$ 关于在给定观测数据 $Y$ 和当前参数 $\theta^{(i)}$ 下对未观测数据 $Z$ 的条件概率分布 $P(Z\mid Y, \theta^{(i)})$ 的期望称为 $Q$ 函数 Q(\theta, \theta^{(i)})=E_Z[\log P(Y, Z \mid \theta) \mid Y, \theta^{(i)}]\tag{11}

1.2 EM 算法的导出（其本质都是对似然 $\log{(Y\mid \theta)}$ 的最大化）

1.2.1 方法一（正向推导，只需要有进步即可）

Reference：[统计学习方法 第二版 179页] 面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据） $Y$ 关于参数 $\theta$ 的对数似然函数，即极大化 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)\tag{12} 上式极大化的主要困难在于未观测数据以及对数里的和（或者积分）。

EM 算法是通过迭代逐步近似极大化 $L(\theta)$ 的，假设在第 $i$ 次迭代后 $\theta$ 的估计值是 $\theta^{(i)}$ 。我们希望新估计值 $\theta$ 能使 $L(\theta)$ 增加，即 $L(\theta)>L(\theta^{(i)})$，并逐步达到极大值。 L(\theta) -L(\theta^{(i)})=\log \left(\sum_Z P(Y\mid Z, \theta)P(Z\mid \theta)\right) - \log P(Y\mid \theta^{(i)})\tag{13} 利用 Jensen 不等式得到其下界：

$$\begin{aligned} L(\theta) -L(\theta^{(i)}) = &\log \left(\sum_Z P(Z\mid Y,\theta^{(i)}) \frac{P(Y\mid Z, \theta)P(Z\mid \theta)}{P(Z\mid Y,\theta^{(i)})} \right)-\log P(Y\mid \theta^{(i)}) \\ \geq & \sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y\mid Z, \theta)P(Z\mid \theta)}{P(Z\mid Y, \theta^{(i)})} -\log P(Y\mid \theta^{(i)}) \\ =&\sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y\mid Z, \theta)P(Z\mid \theta)}{P(Z\mid Y, \theta^{(i)})} - \sum_Z P(Z\mid Y,\theta^{(i)})\log P(Y\mid \theta^{(i)})\\ =& \sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y\mid Z,\theta)P(Z\mid \theta)}{P(Z\mid Y, \theta^{(i)})P(Y\mid \theta^{(i)})}\\ \end{aligned}\tag{14}$$

令

$$\begin{aligned} B(\theta, \theta^{(i)})=&L(\theta^{(i)})+ \sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y\mid Z,\theta)P(Z\mid \theta)}{P(Z\mid Y,\theta^{(i)})P(Y\mid \theta^{(i)})}\tag{15}\\ \end{aligned}$$

则 $L(\theta) \geq B(\theta, \theta^{(i)})$，即函数 $B(\theta, \theta^{(i)})$ 是 $L(\theta)$ 的一个下界，而且 $B(\theta^{(i)}, \theta^{(i)}) = L(\theta^{(i)})$ 。因此，任何使 $B(\theta, \theta^{(i)})$ 相较于在 $\theta^{(i)}$ 处增大的 $\theta$ 也可以使相应的 $L(\theta)$ 增大，即 $L(\theta^{(i)}) = B(\theta^{(i)}, \theta^{(i)})\leq B(\theta^{(i+1)}, \theta^{(i)})\leq L(\theta^{(i+1)})$ 。

$$\begin{aligned} \theta^{(i+1)} =& \mathrm{arg} \max_{\theta} B(\theta, \theta^{(i)})\\ =& \mathrm{arg} \max_{\theta} \left(L(\theta^{(i)}) +\sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y\mid Z, \theta) P(Z\mid \theta)}{P(Z\mid Y,\theta^{(i)})P(Y\mid \theta^{(i)})} \right)\\ =& \mathrm{arg} \max_{\theta} \left( \sum_Z P(Z\mid Y, \theta^{(i)}) \log P(Y\mid Z, \theta) P(Z\mid \theta)\right)\\ =& \mathrm{arg} \max_{\theta} Q(\theta, \theta^{(i)})\\ \end{aligned}\tag{16}$$

下图给出 EM 算法的直观解释：

1.2.2 方法二（通过条件概率公式引入隐变量）

Reference：[变分推断PPT] 对等式两边 $\log P(Y \mid \theta)=\log P(Y,Z \mid \theta)-\log P(Z \mid Y,\theta)$ 分别关于隐变量的后验分布求期望

左边得到

$$\begin{aligned} \text{Left} =& \sum_Z P(Z\mid Y,\theta^{(i)})\log P(Y\mid \theta)\\ =& \log P(Y\mid \theta)\sum_Z P(Z\mid Y,\theta^{(i)})\\ =& \log P(Y\mid \theta)\\ \end{aligned}\tag{17}$$

右边得到

$$\begin{aligned} \text{Right} =& \sum_Z P(Z\mid Y,\theta^{(i)}) \log P(Y,Z\mid \theta) -\sum_Z P(Z\mid Y,\theta^{(i)}) \log P(Z\mid Y,\theta)\\ =& Q(\theta, \theta^{(i)})- H(\theta, \theta^{(i)})\\ \end{aligned}\tag{18}$$

此处 $Q(\theta,\theta^{(i)})$ 即为 EM 算法中 M 步的优化目标，因此有 $Q(\theta^{(i+1)}, \theta^{(i)}) \geq Q(\theta^{(i)}, \theta^{(i)})$ 。

而对于 $H(\theta, \theta^{(i)})$ ，可以证明

$$\begin{aligned} &H(\theta^{(i+1)}, \theta^{(i)}) - H(\theta^{(i)}, \theta^{(i)}) \\ =& \sum_Z P(Z\mid Y, \theta^{(i)})\log P(Z\mid Y,\theta^{(i+1)}) - \sum_Z P(Z\mid Y, \theta^{(i)})\log P(Z \mid Y, \theta^{(i)})\\ =& \sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Z\mid Y, \theta^{(i+1)})}{P(Z\mid Y, \theta^{(i)})}\\ \leq & \log \sum_Z P(Z\mid Y,\theta^{(i)})\cdot \frac{P(Z\mid Y, \theta^{(i+1)})}{P(Z\mid Y, \theta^{(i)})}\\ =& 0\\ \end{aligned}\tag{19}$$

从而得到

$$\begin{aligned} & \log P(Y\mid \theta^{(i+1)})-\log P(Y\mid \theta^{(i)}) \\ =& [Q(\theta^{(i+1)}, \theta^{(i)})-H(\theta^{(i+1)}, \theta^{(i)})] - [Q(\theta^{(i)}, \theta^{(i)})-H(\theta^{(i)}, \theta^{(i)})]\\ =& [Q(\theta^{(i+1)}, \theta^{(i)}) - Q(\theta^{(i)}, \theta^{(i)})]-[H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)}, \theta^{(i)})]\\ \geq & 0 \end{aligned}\tag{20}$$

1.2.3 方法三（引入隐变量的近似分布，承接变分推断内容）

Reference：[变分推断PPT] 引入隐变量 Z 的某种分布 $q_{\phi}(Z)$

$$\begin{aligned} \log P(Y \mid \theta) =& \log P(Y,Z\mid \theta) - \log P(Z\mid Y, \theta)\\ =& \log \frac{P(Y,Z \mid \theta)}{q(Z)}-\log \frac{P(Z\mid Y, \theta)}{q(Z)}\\ \end{aligned}\tag{21}$$

对上式两边分别关于分布 $q(Z)$ 求期望，左边得到

$$\begin{aligned} \text{Left} =& \sum_Z q(Z) \log P(Y\mid \theta) \\ =&\log P(Y\mid \theta)\\ \end{aligned}\tag{22}$$

右边得到

$$\begin{aligned} \text{Right} =& \sum_Z q(Z) \log \frac{P(Y,Z\mid \theta)}{q(Z)} - \sum_Z q(Z) \log \frac{P(Z\mid Y, \theta)}{q(Z)} \end{aligned}\tag{23}$$

联立得到

$$\begin{aligned} \underbrace{\log P(Y\mid \theta)}_{\text{evidence}} =& \sum_Z q(Z) \log \frac{P(Y,Z\mid \theta)}{q(Z)} - \sum_Z q(Z) \log \frac{P(Z\mid Y, \theta)}{q(Z)}\\ =& \underbrace{\sum_Z q(Z) \log \frac{P(Y,Z\mid \theta)}{q(Z)}}_{\text{ELBO}} + \underbrace{\sum_Z q(Z) \log \frac{q(Z)}{P(Z\mid Y,\theta)}}_{\mathrm{KL}(q(Z)\mid \mid P\left(Z\mid Y, \theta\right))}\\ \end{aligned}\tag{24}$$

• $\log P(Y\mid \theta)$ 被称为证据（evidence）
• $\sum_Z q(Z)\log \frac{P(Y,Z\mid \theta)}{q(Z)}$ 被称为证据下界（evidence lower bound，ELBO）
• $\sum_Z q(Z) \log \frac{q(Z)}{P(Z\mid Y,\theta)} = KL(q(Z)\mid \mid P(Z\mid Y,\theta))$ 是分布 $q(Z)$ 相对于分布 $P(Z\mid Y,\theta)$ 的 [[KL散度]]（Kullback-Leibler divergence）

因为 KL 散度非负，从而得到下式，当且仅当 $q(Z) = P(Z\mid Y,\theta)$ 时取等号

$$\begin{aligned} \underbrace{\log P(Y\mid \theta)}_{\text{evidence}} \geq \underbrace{\sum_Z q(Z) \log \frac{P(Y,Z\mid \theta)}{q(Z)}}_{\text{ELBO}} \end{aligned}\tag{25}$$

E 步：固定参数 $\theta^{(i)}$， 取 $q(Z) =P(Z\mid Y, \theta^{(i)})$ ，此时有（不严谨，为何此时取等号，疑为ppt错误）

$$\underbrace{\log P(Y\mid \theta)}_{\text{evidence}} = \underbrace{\sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y,Z\mid \theta)}{P(Z\mid Y,\theta^{(i)})}}_{\text{ELBO}}\tag{26}$$

M 步：ELBO 关于参数 $\theta$ 求最大，更新参数

$$\begin{aligned} \theta^{(i+1)} =& \mathrm{arg} \max_{\theta} \sum_Z P(Z\mid Y,\theta^{(i)}) \log \frac{P(Y,Z\mid \theta)}{P(Z\mid Y,\theta^{(i)})} \\ =& \mathrm{arg} \max_{\theta} \underbrace{\sum_Z P(Z\mid Y, \theta^{(i)}) \log P(Y,Z\mid \theta)}_{Q(\theta, \theta^{(i)})}\\ \end{aligned}\tag{27}$$

笔者更正：固定参数 $\theta^{(i)}$， 取 $q(Z) =P(Z\mid Y, \theta^{(i)})$ ，此时有（根据式24）

$$\log P(Y\mid \theta) = \underbrace{\sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y,Z\mid \theta)}{P(Z\mid Y,\theta^{(i)})}}_{A} + \underbrace{\sum_Z P(Z\mid Y,\theta^{(i)}) \log \frac{P(Z\mid Y,\theta^{(i)})}{P(Z\mid Y,\theta)}}_{B}\tag{28}$$

而当$\theta = \theta^{(i)}$ 时有

$$\log P(Y\mid \theta^{(i)}) = \underbrace{\sum_Z P(Z\mid Y, \theta^{(i)}) \log \frac{P(Y,Z\mid \theta^{(i)})}{P(Z\mid Y,\theta^{(i)})}}_{C} + \underbrace{\sum_Z P(Z\mid Y,\theta^{(i)}) \log \frac{P(Z\mid Y,\theta^{(i)})}{P(Z\mid Y,\theta^{(i)})}}_{D}\tag{29}$$

由KL散度性质可知 $B\geq D=0$：

$$\begin{aligned} \theta^{(i+1)} &=\arg \max_{\theta} \underbrace{\sum_Z P(Z\mid Y,\theta^{(i)}) \log P(Y,Z\mid \theta) }_{Q(\theta,\theta^{(i)})} \\ &=\arg \max_{\theta} A \\ \Rightarrow \quad \quad & A(\theta^{(i+1)}) \geq C \\ \Rightarrow \quad \quad & \log P(Y\mid \theta^{(i+1)}) \geq \log P(Y\mid \theta^{(i)}) \end{aligned}\tag{30}$$

1.3 EM 算法的收敛性

定理1：设 $L(\theta)=\log P(Y \mid \theta)$ 为观测数据的对数似然函数， $\theta^{(i)}(i=1,2, \cdots)$为 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) 在函数 $\color{red}{Q\left(\theta, \theta^{\prime}\right)}$ 与 $\color{red}{L(\theta)}$ 满足一定条件下，由 EM 算法得到的参数估计序列 $\theta^{(i)}$ 的收敛值 $\theta^*$ 是 $L(\theta)$ 的稳定点。

证明： (1) 由 $L(\theta^{(i)})=\log P\left(Y \mid \theta^{(i)}\right)$ 的单调性及 $P(Y \mid \theta)$ 的有界性得到。 (2) 证明从略，参阅文献 [1983 On the convergence properties of the EM algorithm]。

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

定义2：高斯混合模型 高斯混合模型是指具有如下形式的概率分布模型： P(y\mid \theta) =\sum_{k=1}^K \alpha_k \cdot \phi(y\mid \theta_k)\tag{31} 其中，$\alpha_k$ 是系数，$\alpha_k \geq 0$，$\sum_{k=1}^K \alpha_k = 1$；$\phi(y \mid \theta_k)$ 是高斯分布密度，$\theta_k = (\mu_k,\sigma_k^2)$， \phi(y\mid \theta_k) = \frac{1}{\sqrt{2 \pi} \sigma_k} \exp{\left(-\frac{(y-\mu_k)^2}{2\sigma_k^2}\right)}\tag{32} 称为第 $k$ 个分模型。一般混合模型可以由任意概率分布密度代替式（29）中的高斯分布密度，此处只介绍最常用的高斯混合模型。

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

假设观测数据 $y_1,y_2,\dots,y_N$ 由高斯混合模型生成， P(y\mid \theta) = \sum_{k=1}^K \alpha_k \cdot \phi(y\mid \theta_k) \tag{33} 其中，$\theta = (\alpha_1,\alpha_2,\dots,\alpha_K;\theta_1,\theta_2,\dots,\theta_K)$。

观测数据的产生过程：首先依概率 $(\alpha_1,\dots,\alpha_K)$ 选择第 $k$ 个高斯分布模型，然后依第 $k$ 个分模型的概率分布 $\phi(y\mid \theta_k)$ 生成观测数据 $y_j$。这时观测数据 $y_j$，$j=1,2,\dots,N$ 是已知的；反映观测数据 $y_j$ 来自第 $k$ 个分模型的数据是未知的，$k=1,2,\dots,K$，以隐变量 $\gamma_{jk}$ 表示，其定义如下： \gamma_{jk}=\begin{cases} 1, &\quad \text{第 } j \text{ 个观测来自第 } k \text{ 个分模型} \\ 0, &\quad \text{否则} \end{cases}\tag{34} 有了观测数据 $y_j$ 及未观测数据 $\gamma_{jk}$，那么完全数据是 (y_j,\gamma_{j1},\gamma_{j2},\dots,\gamma_{jK}), \quad j=1,2,\dots,N \tag{35} 于是可以写出完全数据的似然函数

$$\begin{aligned} P(y,\gamma\mid \theta)&=\prod_{j=1}^N P(y_j,\gamma_{j1},\gamma_{j2},\dots,\gamma_{jK} \mid \theta)\\ &=\prod_{k=1}^K \prod_{j=1}^N[\alpha_k \cdot \phi(y_j \mid \theta_k)]^{\gamma_{jk}}\\ &=\prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N [\phi(y_j\mid \theta_k)]^{\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{(y_j-\mu_k)^2}{2\sigma_k^2}\right)} \right]^{\gamma_{jk}} \end{aligned}\tag{36}$$

式中，$n_k=\sum_{j=1}^N \gamma_{jk}$，$\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\}\tag{37} 进一步计算 Q 函数

$$\begin{aligned} Q(\theta, \theta^{(i)}) & =E_{\gamma}\left[\log P(y, \gamma \mid \theta) \mid y, \theta^{(i)}\right] \\ & =E_{\gamma}\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_{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{aligned}\tag{38}$$

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

$$\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}\tag{39}$$

$\hat{\gamma}_{jk}$ 是在当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分模型的概率，称为分模型 $k$ 对观测数据 $y_j$ 的响应度。将$\hat{\gamma}_{jk}=E[{\gamma_{jk}}]$ 及 $\hat{n}_k = \sum_{j=1}^N E[\gamma_{jk}]$ 代入式（38），即得 Q\left(\theta, \theta^{(i)}\right)=\sum_{k=1}^K\left\{\hat{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{40} 迭代的 M 步是求函数 $Q(\theta, \theta^{(i)})$ 对 $\theta$ 的极大值，即求新一轮迭代的模型参数 \theta^{(i+1)} = \arg \max_{\theta} Q(\theta, \theta^{(i)})\tag{41} 用 $\hat{\mu}_k$，$\hat{\sigma}_{k}^2$ 及 $\hat{\alpha}_k$ ，$k=1,2,\dots, K$，表示 $\theta^{(i+1)}$ 的各参数。求 $\hat{\mu}_k$，$\hat{\sigma}_{k}^2$ 只需将式（40）分别对 $\hat{\mu}_k$，$\hat{\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}_{jk} \cdot y_j}{\sum_{j=1}^N \hat{\gamma}_{jk}}, \quad k=1,2,\dots,K \tag{42} \hat{\sigma}_{k}^2 =\frac{\sum_{j=1}^N \hat{\gamma}_{jk} (y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}}, \quad k=1,2,\dots,K \tag{43} \hat{\alpha}_k = \frac{\hat{n}_k}{N} =\frac{\sum_{j=1}^N \hat{\gamma}_{jk}}{N}, \quad k=1,2,\dots, K \tag{44} 重复以上计算，直到对数似然函数值不再有明显的变化为止。

1.5 EM 算法的推广（可参考变分推断PPT更简单易理解）

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

定义3：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}) \tag{45} 称为 $F$ 函数，式中 $H(\tilde{P})=-E_{\tilde{P}}\log \tilde{P}(Z)$ 是分布 $\tilde{P}(Z)$ 的熵。

在定义3中，通常假设 $P(Y,Z\mid \theta)$ 是 $\theta$ 的连续函数，因而 $F(\tilde{P},\theta)$ 是 $\tilde{P}$ 和 $\theta$ 的连续函数。函数 $F(\tilde{P}, \theta)$ 还有以下重要性质。 引理1： 对于固定的 $\theta$ ，存在唯一的分布 $\tilde{P}_{\theta}$ 极大化 $F(\tilde{P},\theta)$，这时 $\tilde{P}_{\theta}$ 由下式给出 \tilde{P}_{\theta}(Z)=P(Z\mid Y,\theta) \tag{46} 并且 $\tilde{P}_{\theta}$ 随 $\theta$ 连续变化。

证明： 对于固定的 $\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)\tag{47} 将其对 $\tilde{P}$ 求偏导数（针对特定 $Z$） \frac{\partial L}{\partial \tilde{P}(Z)}=\log P(Y,Z\mid \theta)-\log \tilde{P}(Z) -1-\lambda \tag{48} 令偏导数等于 0，得出 \lambda=\log P(Y,Z\mid \theta)-\log \tilde{P}_{\theta}(Z)-1 \tag{49} 由此推出 $\tilde{P}_{\theta}(Z)$ 与 $P(Y,Z\mid \theta)$ 成比例 \frac{P(Y,Z\mid \theta)}{\tilde{P}_{\theta}(Z)}=\exp{(1+\lambda)} \tag{50} 再从约束条件 $\sum_Z \tilde{P}_{\theta}(Z)=1$ 得到式（46）。 由假设 $P(Y,Z\mid\theta)$ 是 $\theta$ 的连续函数，得到 $\tilde{P}_{\theta}$ 是 $\theta$ 的连续函数。

引理2： 若 $\tilde{P}_{\theta}(Z)=P(Z\mid Y,\theta)$，则 F(\tilde{P},\theta)=\log P(Y\mid \theta) \tag{51} 证明：

$$\begin{aligned} F(\tilde{P},\theta)=&E_{\tilde{P}}[\log P(Y,Z\mid \theta)]+H(\tilde{P})\\ =&E_{\tilde{P}}[\log \frac{P(Y,Z\mid\theta)}{P(Z\mid Y,\theta)}]\\ =&\sum_Z P(Z\mid Y,\theta)\cdot \log \frac{P(Y,Z\mid\theta)}{P(Z\mid Y,\theta)} \\ =&\sum_Z P(Z\mid Y,\theta) \cdot \log P(Y\mid \theta) \\ =&\log P(Y\mid \theta) \end{aligned}\tag{52}$$

由以上引理，可以得到关于 EM 算法用 $F$ 函数的极大-极大算法的解释。 定理2： 设 $L(\theta)=\log P(Y\mid \theta)$ 为观测数据的对数似然函数，$\theta^{(i)}, i=1,2,\dots$，为 EM 算法得到的参数估计序列，函数 $F(\tilde{P},\theta)$ 由式（45）定义。如果 $F(\tilde{P},\theta)$ 在 $\tilde{P}^*$ 和 $\theta^*$ 有局部极大值，那么 $L(\theta)$ 也在 $\theta^*$ 有局部极大值。类似地，如果 $F(\tilde{P},\theta)$ 在 $\tilde{P}^*$ 和 $\theta^*$ 达到全局最大值，那么 $L(\theta)$ 也在 $\theta^*$ 达到全局最大值。 证明： 由引理1和引理2可知，$L(\theta)=\log P(Y\mid \theta)=F(\tilde{P}_{\theta},\theta)$ 对任意 $\color{red}{\theta}$ 成立。特别地，对于使 $F(\tilde{P},\theta)$ 达到极大的参数 $\theta^*$，有 L(\theta^*)=F(\tilde{P}_{\theta^*}, \theta^*)=F(\tilde{P}^*, \theta^*)\tag{53} 为了证明 $\theta^*$ 是 $L(\theta)$ 的极大点，需要证明不存在接近 $\theta^*$ 的点 $\theta^{**}$，使 $L(\theta^{**})>L(\theta^*)$。假如存在这样的点 $\theta^{**}$，那么应有 $F(\tilde{P}^{**},\theta^{**})>F(\tilde{P}^*,\theta^*)$，这里 $\tilde{P}^{**}=\tilde{P}_{\theta^{**}}$。但因 $\tilde{P}_{\theta}$ 是随 $\theta$ 连续变化的，$\tilde{P}^{**}$ 应接近 $\tilde{P}^*$，这与 $\tilde{P}^*$ 和 $\theta^*$ 是 $F(\tilde{P},\theta)$ 的局部极大点的假设矛盾。类似可以证明关于全局最大值的讨论。

定理3： 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(Z \mid Y, \theta^{(i)}) \tag{54}$$

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

$$\begin{aligned} F(\tilde{P}^{(i+1)}, \theta) & =E_{\tilde{P}^{(i+1)}}[\log P(Y, Z \mid \theta)]+H(\tilde{P}^{(i+1)}) \\ & =\sum_Z \log P(Y, Z \mid \theta) P(Z \mid Y, \theta^{(i)})+H(\tilde{P}^{(i+1)}) \end{aligned}\tag{55}$$

由 $Q\left(\theta, \theta^{(i)}\right)$ 的定义式 11 有

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

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

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

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

1.5.2 GEM 算法

算法2： 输入：观测数据， $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)$ 。

算法3： 输入：观测数据, $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}\tag{58}$$

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

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

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

当参数 $\theta$ 的维数为 $d(d \geqslant 2)$ 时，可采用一种特殊的 GEM 算法，它将 EM 算法的 M 步分解为 $d$ 次条件极大化，每次只改变参数向量的一个分量，其余分量不改变。

算法4： 输入：观测数据，$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}\tag{60}$$

(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)\tag{61}$$

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