EM算法（1）—— 彻底理解EM算法推导

最近看了很多关于EM算法推导的文章，包括《统计学习方法（第二版）》中关于EM算法的内容，总感觉说得不够清楚。虽然公式都写得挺详细，但是没有说清楚为什么要这样算，怎么想到这样变形的，本文总结一下我理解的EM算法推导过程，尽量把EM的核心思想说清楚。

• 参考：

1. 《统计学习方法（第二版）》—— 第九章
2. ​​浅析：从最大似然估计（MLE）、最大后验估计（MAP）到期望最大算法（EM）​​
3. ​​为什么不把EM算法说清楚？​​
4. ​​EM算法详细推导和讲解​​
5. ​​如何感性地理解EM算法？​​
6. ​​NLP —— 图模型（零）：EM算法简述及简单示例（三硬币模型）​​
7. ​​EM算法的九重境界​​

---

文章目录

• ​​0. 补充​​
• ​​1. 引入​​

• ​​1.1 三硬币模型​​

• ​​1.1.1 建模​​
• ​​1.1.2 分析​​
• ​​1.1.2 使用EM方法解​​

• ​​1.2 问题形式化​​

• ​​2. EM算法​​

• ​​2.1 EM算法流程​​
• ​​2.2 EM算法的推导​​

• ​​2.2.1 EM核心思想：迭代优化​​
• ​​2.2.2 EM推导​​

• ​​2.2.2.1 两个难点​​
• ​​2.2.2.2 解难点2：Jensen不等式​​
• ​​2.2.2.3 解难点1：分步迭代​​

• ​​2.3 EM算法的收敛性​​

• ​​2.3.1 定理1​​
• ​​2.3.2 定理2​​
• ​​2.3.3 小结​​

• ​​3. 用EM算法解三硬币问题​​

• ​​3.1 E步​​

• ​​3.1.1 公式法直接计算 Q 函数​​
• ​​3.1.2 从期望角度计算 Q 函数​​

• ​​3.2 M步​​
• ​​3.3 对比​​
• ​​3.4 代码​​

0. 补充

1. 极大似然估计、最大后验估计：​​一文看懂 “极大似然估计” 与 “最大后验估计”​​
2. 关于Jensen不等式：​​Jensen不等式​​

1. 引入

• 根据观测性，可以把模型变量分为 ​​观测变量​​​ 和 ​​隐变量/潜在变量​​。如果模型的变量都是观测变量，可以直接使用极大似然估计或最大后验估计法估计模型参数。但当模型含有隐变量时，这些方法就不能直接使用了。
• EM算法是一种迭代算法，用于对含有隐变量的概率模型参数进行极大似然估计或最大后验概率估计。EM的适用情况和 MLE 或 MAP 一样，都是在模型确定的情况下估计模型参数，区别在于EM算法允许模型中存在隐变量

1.1 三硬币模型

• 先看一个经典的三硬币模型示例，感受变量的观测性和EM算法流程。

1.1.1 建模

• 假设有 A,B,C 三枚硬币，正面出现的概率分别为 $\pi,p,q$。进行如下实验：先掷硬币A，正面选B，反面选C；然后掷选出的硬币。我们不能观测到A的结果（隐变量），但是可以观测到选出硬币B或C的掷币结果（观测变量）。假如10次实验观测结果为
  $正，正，反，正，反，反，正，反，正，正$
  现在请估计三枚硬币掷出正面的概率，即三硬币模型的参数
• 用1表示正面，0表示反面，设随机变量 $y,z$ 分别为观测变量（B/C的掷币结果）和隐变量（A的掷币结果），观测值分别为 $Y,Z \in \{0,1\}$ ，模型参数为 $\theta = \{\pi,p,q\}$ 。注意到掷硬币是一种两点分布，则三硬币问题可以建模表示为
  $\begin{aligned} P(y|\theta) &= \sum_zP(y,z|\theta) = \sum_zP(z|\theta)P(y|z,\theta) \\ & = \pi p^y(1-p)^{1-y} + (1-\pi)q^y(1-q)^{1-y} \end{aligned}$
• 将 ​​观测数据​​​ 和 ​​未观测数据​​​ 分别表示为来自总体 $y$ 和 $z$ 的容量为 $n$ 的简单随机样本 $y = \{y_1,y_2,...,y_n\}^T,z = \{z_1,z_2,...,z_n\}^T$ ，观测值和未观测值分别为 $Y = \{Y_1,Y_2,...,Y_n\}^T,Z = \{Z_1,Z_2,...,Z_n\}^T$，考虑求模型参数 $\theta$

1. 极大似然估计：观测数据的似然函数、对数似然函数、参数 $\theta$ 估计值分别为
   $\begin{aligned} P(Y|\theta) &= \sum_ZP(Z|\theta)P(Y|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}] \\ logP(Y|\theta) &= \sum_{j=1}^nlog[\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}]\\ \end{aligned}$
   $\hat{\theta} = \argmax_\theta logP(Y|\theta)$
2. 最大后验估计：参数 $\theta$ 的后验概率、估计值分别为
   $\begin{aligned} P(\theta|Y) &= \frac{P(\theta,Y)}{P(Y)} \\ &= \frac{P(\theta)P(Y|\theta)}{P(Y)} \\ &= \frac{P(\theta)P(Z|\theta)P(Y|Z,\theta)}{P(Y)} \\ &= \frac{P(\theta)P(Z|\theta)P(Y|Z,\theta)}{\int P(Y|Z,\theta)P(Z|\theta)P(\theta)d\theta}\\ &\propto P(\theta)P(Z|\theta)P(Y|Z,\theta) \end{aligned}$
   $\hat{\theta} = \argmax_\theta P(\theta)P(Z|\theta)P(Y|Z,\theta)$
   这里需要 $\theta$ 的先验概率 $p(\theta)$，可以设成三个 $\mu = 0.5$ 的正态分布连乘，比较麻烦，所以这个例子后面就只写最大似然估计的过程了。

1.1.2 分析

• 先尝试直接极大似然估计求解。令偏导为0列方程组
  $\left\{ \begin{aligned} \frac{\partial logP(Y|\theta)}{\partial \pi} & = \sum_{j=1}^n \frac{p^{y_j}(1-p)^{1-y_j}-q^{y_j}(1-q)^{1-y_j}}{\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}} \stackrel{令}{=}0\\ \frac{\partial logP(Y|\theta)}{\partial p} & = \sum_{j=1}^n \frac{\pi p^{y_j-1}(1-p)^{-y_j}(y_j-p)}{\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}} \stackrel{令}{=}0\\ \frac{\partial logP(Y|\theta)}{\partial q} & = \sum_{j=1}^n \frac{(1-\pi) q^{y_j-1}(1-q)^{-y_j}(y_j-q)}{\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}} \stackrel{令}{=}0 \end{aligned} \right.$
  带入观测数据 $Y$，有
  $\left\{ \begin{aligned} &-10\pi q^2-10\pi p^2 + 10q^2+20\pi pq-10pq+6p-6q &= 0 \space\space\space\space (1)\\ &10\pi^2 q-10\pi^2 p-10\pi q+6\pi &= 0 \space\space\space\space(2) \\ &-10\pi^2 q+10\pi^2 p + 20\pi q -10\pi p - 10q -6\pi+6 &= 0 \space\space\space\space(3) \end{aligned} \right.$
  发现式子(2)(3)有关系：$(3) =\frac{\pi}{1-\pi} (2)$，这两个式子等价，因此有无穷解。最终结果为 $p=q= 0.6$，$\pi$。由于未观测数据 $Z$ 不可见，直接最大化观测数据 $Y$ 的似然会导致 $Z$ 及相关参数被当作不存在，使得模型实质上退化到了一个硬币（不存在隐变量）的情况

1.1.2 使用EM方法解

• 使用EM算法可迭代求出更好的唯一解，过程如下

1. 先选取参数初值，$\theta^{(0)} = \{\pi^{(0)},p^{(0)},q^{(0)}\}$
2. 反复执行以下两步，迭代计算参数 $\theta$ 估计值直到收敛为止。设第 $i$ 次迭代参数估计值为 $\theta^{(i)} = \{\pi^{(i)},p^{(i)},q^{(i)}\}$，则第 $i+1$ 次迭代如下

1. E步：计算模型在参数 $\theta^{(i)}$ 下观测数据 $y_j$ 来自掷硬币 B 的概率 $\mu_j^{(i+1)}$
   $\begin{aligned} \mu_j^{(i+1)} = P(B|y_j) &= \frac{P(B,y_j)}{P(y_j)} = \frac{P(B)P(y_j|B)}{P(B)P(y_j|B)+P(C)P(y_j|C)} \\ &= \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}} \end{aligned}$
2. M步：计算模型参数的新估计值
   $\begin{aligned} \pi^{(i+1)} &= P(B|Y) = \frac{1}{n} \sum_{i=1}^n \mu_j^{(i+1)} \\ p^{(i+1)} &= \frac{P(B|Y=1)}{P(B|Y)} = \frac{\sum_{i=1}^n \mu_j^{(i+1)}y_j}{\sum_{i=1}^n \mu_j^{(i+1)}} \\ q^{(i+1)} &= \frac{P(C|Y=1)}{P(C|Y)} = \frac{\sum_{i=1}^n (1-\mu_j^{(i+1)})y_j}{\sum_{i=1}^n (1-\mu_j^{(i+1)})} \end{aligned}$

• 带入观测值计算

1. 初值选 $\theta^{(0)} = \{0.5,0.5,0.5\}$ 时，EM方法解出的估计值为 $\hat{\theta} = \{0.5,0.6,0.6\}$
2. 初值选 $\theta^{(0)} = \{0.4,0.6,0.7\}$ 时，EM方法解出的估计值为 $\hat{\theta} = \{0.4064,0.5368,0.6432\}$

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

1.2 问题形式化

• 一般地，用 $Y$ 表示观测变量的数据，用 $Z$ 表示隐变量的数据，用 $\theta$

1. ​​完全数据​​：模型中所有随机变量的数据，$Y$ 和 $Z$ 连在一起
2. ​​不完全数据​​：观测数据 $Y$ 又称不完全数据

• 当完全数据已知时，可以直接使用MLE估计或MAP估计 $\theta$；只知道不完全数据时，EM算法通过迭代最大化似然函数或后验概率。每次迭代分成两步

1. E步：求Q函数（期望）
2. M步：求Q函数（期望）极大

2. EM算法

2.1 EM算法流程

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

1. 选择合适的初始参数 $\theta^{(0)}$，开始迭代。初值可以任意选择，但需注意EM算法是对初值敏感的
2. E步：记 $\theta^{(i)}$ 为第 $i$ 轮迭代参数估计值，在第 $i+1$ 次迭代的 E 步，计算
   $\begin{aligned} Q(\theta,\theta^{(i)}) &= \mathbb{E}_Z[logP(Y,Z|\theta) | Y,\theta^{(i)}] \\ &= P(Z|Y,\theta^{(i)}) logP(Y,Z|\theta) \end{aligned}$
   这里 $P(Z|Y,\theta^{(i)})$ 是给定观测数据 $Y$ 和当前估参数估计值 $\theta^{(i)}$ 下隐变量数据 $Z$
3. M步：求使得 $Q(\theta,\theta^{(i)})$ 最大化的 $\theta$，确定第 $i+1$ 轮的参数估计值 $\theta^{(i+1)}$
   $\theta^{(i+1)} = \argmax_\theta Q(\theta,\theta^{(i)})$
4. 重复第 2/3 步，直到收敛。收敛条件一般是：对于给定小量 $\varepsilon_1,\varepsilon_2$，满足下式则收敛，停止迭代
   $|\theta^{(i+1)} - \theta^{(i)}|| < \varepsilon_1 \space\space或 \space\space ||Q(\theta^{(i+1)},\theta^{(i)}) - Q(\theta^{(i)},\theta^{(i)}) || < \varepsilon_1 \\$

• ​​Q函数​​ ：定义完全数据的对数似然函数 $logP(Y,Z|\theta)$ 关于给定观测数据 $Y$ 和当前估计参数 $\theta^{(i)}$ 下对未观测数据 $Z$ 的条件概率分布 $P(Z|Y,\theta^{(i)})$为Q函数，即
  $\begin{aligned} Q(\theta,\theta^{(i)}) &= \mathbb{E}_Z[logP(Y,Z|\theta) | Y,\theta^{(i)}] \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ &=\sum_Z P(Z|Y,\theta^{(i)})logP(Y,Z|\theta) \\ &= \sum_{j=1}^N \sum_Z P(Z|Y_j,\theta^{(i)})logP(Y_j,Z|\theta) \end{aligned}$
  EM 算法的每次迭代，实际是在求Q函数这个期望及其极大，这也是EM算法全程 “期望最大化”（Expection-Maximization）的由来

2.2 EM算法的推导

• EM算法目标：极大化参数 $\theta$ 下观测数据 $Y$
• 方法：利用Jensen不等式推出目标对数似然函数下界，迭代优化下界以优化目标，不保证全局最优

2.2.1 EM核心思想：迭代优化

• 在推导之前，有必要探讨一下为什么直接最大化似然函数没有解析解。先看极大化目标：参数 $\theta$ 下观测数据 $Y$ 的对数似然函数。即
  $\begin{aligned} L(\theta) &= logP(Y|\theta) = log\sum_Z P(Y,Z|\theta) \\ &= log(\sum_ZP(Y|Z,\theta)P(Z|\theta)) \end{aligned}$
  回顾极大似然估计 MLE 的一般形式，假设可以观测到完全数据 $X$，对数似然函数为
  $L(\theta) = log \prod_X P(X|\theta) = \sum_X log P(X|\theta)$
  要得到 $\argmax_\theta L(\theta)$ 的解析解，需要保证 $L(\theta)$ 中只有 $\theta$ 一个变量。在隐变量 $z$ 不能观测时，优化目标中同时存在估计量 $P(Z|\theta)$ 和 $\theta$
• 事实上，未观测数据 $Z$ 和参数 $\theta$。还是以三硬币模型为例，如果我们知道 $Z$，就能知道哪些观测结果来自硬币B，哪些来自硬币C，进而能分别使用极大似然估计法得到 $p,q$。在只有不完全数据 $Y$ 时，为了区分观测结果的来源以估计 $p,q$，我们需要先估计每一个观测数据 $Y_i$ 对应的 $Z_i$ 取值（或取值的分布），这个估计又是依赖于 $p,q$ 的，如图所示
  
  对于这种含有两个未知变量的优化问题，可以使用迭代的方法。选取一个合适的初始值开始迭代，每轮迭代过程中，先固定一个变量优化另一个，再固定另一个变量估计优化这个，直到二者逐步收敛为止。EM 算法在此思想基础上更进一步，在估计未观测数据时，不止估计其值，而是估计其分布，进而把优化目标转换完全数据似然函数为对未观测数据分布的期望形式。此处推荐参考一篇文章：​​如何感性地理解EM算法？​​
• 强化学习中 model-based 方法 policy-iteration 使用了这种分步迭代思想，不妨把它的示意图贴过来
  
  policy-iteration 的迭代分为两个阶段，evaluation步骤对当前策略下的状态价值进行评估，improvement步骤利用策略价值贪心地优化策略。可以把每个流程看成二维空间中的一条线，代表对于其目标的解决方案。每个流程都把价值 $V$ 或 策略 $\pi$ 推向其中一条线，由于这两条线不是正交的，所以两个目标间会产生相互作用。evaluation步会导致策略逐渐偏离最优，improvement步会导致价值偏离新策略，直接冲向某个目标会导致某种程度上偏离另一个目标，然而整体流程的结果最终会更接近最优的总目标。EM算法中的E步可以类比为evaluation步骤，M步可以类比为improvement步骤。

2.2.2 EM推导

2.2.2.1 两个难点

• 极大化目标：参数 $\theta$ 下观测数据 $Y$ 的对数似然函数。即
  $\begin{aligned} L(\theta) &= logP(Y|\theta) = log\sum_Z P(Y,Z|\theta) \\ &= log(\sum_ZP(Y|Z,\theta)P(Z|\theta)) \\ \end{aligned}$
• 注意到 $\hat{\theta} = \argmax_\theta L(\theta)$ 这一极大化的主要困难有两个
1. 需要同时完成对 $P(Z|\theta)$ 的估计和 $\theta$
2. 求和符号 $\sum$ 在 $log$

2.2.2.2 解难点2：Jensen不等式

• 首先处理第二个难点，我们希望能把 $\sum$ 符号拿到 $log$ 外面，这时就要用到 jensen 不等式，即
  $\sum_{i=1}^M \lambda_if(x_i) \geq f(\sum_{i=1}^M\lambda_ix_i)$
  此式要求 $f$ 是凸函数，$\lambda_i\geq0$ 且 $\sum_{i=1}^M \lambda_i = 1$，当随机变量 $X$ 是常数时等号成立。由于 $log$ 是凹函数，所以需要应用 jensen 不等式的凹函数形式，即把 $\geq$ 换成 $\leq$。
• 为了应用 Jensen 不等式，还要拼凑一个 $\lambda$。注意到 $log$ 里的和是关于隐变量 $Z$ 的和，而 $\lambda$ 需要满足 $\lambda_i\geq0$ 且 $\sum_{i=1}^M \lambda_i = 1$，不妨设 $\lambda$ 为隐变量 $z$ 的分布函数 $K(Z|\theta)$，于是有
  $\begin{aligned} L(\theta) &= logP(Y|\theta) \\ &=log(\sum_Z P(Y,Z|\theta)) \\ &= log(\sum_Z K(Z|\theta) \frac{P(Y,Z|\theta)}{K(Z|\theta)}) \\ &\geq \sum_Z K(Z|\theta)log\frac{P(Y,Z|\theta)}{K(Z|\theta)} \end{aligned}$
  这样我们就找到了似然函数的下界形式，并成功把 $\sum$ 拿到了 $log$ 外边，接下来我们希望通过最大化下界来间接地最大化 $L(\theta)$。为了保证优化效果，我们要求 $\sum_Z K(Z|\theta)log(\frac{P(Y,Z|\theta)}{K(Z|\theta)})$，也就是至少有一点使不等式中等号成立。因为如果这个下界特别小，离似然函数太远，优化它也没有什么用
• 现在来讨论 $K(Z|\theta)$ 的形式，为了得到紧下界，回头看 Jensen 不等式中使得等号成立的条件，需要满足
  $\frac{P(Y,Z|\theta)}{K(Z|\theta)} = C$
  因为 $K(Z|\theta)$ 是 $Z$ 的分布函数，有
  $\begin{aligned} &\because \sum_ZK(Z|\theta) = \sum_Z \frac{P(Y,Z|\theta)}{C} = 1 \\ &\therefore C = \sum_Z P(Y,Z|\theta) = P(Y|\theta) \\ &\therefore K(Z|\theta) = \frac{P(Y,Z|\theta)}{C} = \frac{P(Y,Z|\theta)}{P(Y|\theta)} = P(Z|Y,\theta) \end{aligned}$
• 这样我们就得到了紧下界的表示形式，由于始终满足等号成立条件 $\frac{P(Y,Z|\theta)}{P(Z|Y,\theta)} = P(Y|\theta) \equiv C$，因此不等号转换为恒等号，可以证明如下
  $\begin{aligned} \sum_Z P(Z|Y,\theta)log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta)} &= \sum_Z P(Z|Y,\theta)logP(Y|\theta) \\ &= logP(Y|\theta) \\ &= L(\theta) \end{aligned}$

2.2.2.3 解难点1：分步迭代

• 难点1的本质就是没法同时估计 $P(Z|Y,\theta)$ 和优化 $\theta$ ，我们上面的一堆分析并没有解决这个问题。牢记 2.2.1 节的分析，解决这种问题可以使用分步迭代法，每轮迭代先控制一个变量优化另一个变量，再控制另一个变量优化这个变量。具体到EM算法中，第 $i$

  1. 控制 $\theta$ 估计 $Z$：固定模型参数 $\theta^{(i)}$，估计 $Z$ 条件概率分布 $P(Z|Y,\theta^{(i)})$
  3. 控制 $Z$ 优化 $\theta$：利用估计的 $P(Z|Y,\theta^{(i)})$ 重新构造紧下界函数，最大化紧下界函数得到新的参数值 $\theta^{(i+1)}$。由于
     $\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})} 不恒等于 C$
     此时等号恒成立条件不再满足，还原到不等号
     $L(\theta) =log(\sum_ZP(Y,Z|\theta)) \geq \sum_Z P(Z|Y,\theta^{(i)}) log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})}$
     不妨设紧下界为 $l(\theta,\theta^{(i)})$，即
     $l(\theta,\theta^{(i)}) = \sum_Z P(Z|Y,\theta^{(i)}) log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})} \\ L(\theta) \geq l(\theta,\theta^{(i)})$
     注意到当 $\theta = \theta^{(i)}$ 时， $L(\theta^{(i)}) =l(\theta^{(i)},\theta^{(i)})$ 等号成立

  

  4. 如图所示，求 $\theta^{(i+1)} = \argmax_\theta l(\theta,\theta^{(i)})$ ，就可以通过增大紧下界间接地增大 $L(\theta)$

• 忽略掉对 $l(\theta,\theta^{(i)})$ 中对极大化而言无关的项，得到 Q函数如下
  $\begin{aligned} Q(\theta,\theta^{(i)}) &= \sum_Z P(Z|Y,\theta^{(i)}) logP(Y,Z|\theta) \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \end{aligned}$

• 大功告成！至此我们已经对EM算法流程有了清晰的理解，并完整地推导了EM的核心函数 $Q$。不妨比较一下 $Q$ 函数和原始的优化目标 $L(\theta)$
  $\begin{aligned} Q(\theta,\theta^{(i)}) &= \sum_Z P(Z|Y,\theta^{(i)}) logP(Y,Z|\theta) \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}log(P(Y|Z,\theta)P(Z|\theta)) \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ L(\theta) &= log(\sum_ZP(Y,Z|\theta)) \\ &= log(\sum_ZP(Y|Z,\theta)P(Z|\theta))\\ &= log(\mathbb{E}_{Z|\theta}P(Y|Z,\theta)) \end{aligned}$
  主要的变化有两个

  1. 使用 Jensen 不等式把 $\sum$ 拿到 $log$
  2. 引入后验概率 $P(Z|Y,\theta^{(i)})$，使迭代求解成为可能。 这里可能会有疑惑，$Q(\theta,\theta^{(i)})$ 里有个 $P(Y,Z|\theta) = P(Y|Z,\theta)P(Z|\theta)$，这不还是要用 $\theta$ 估计未观测数据分布 $P(Z|\theta)$ 吗。我个人理解是：式子中存在 $P(Z|\theta)$ 不是关键，关键要看优化目标是关于什么的期望，或者说优化目标这个随机变量服从的分布是什么。

2.3 EM算法的收敛性

• EM算法提供了一种近似计算含有隐变量概率模型的极大似然估计的方法。其最大优点是简单性和普适性。
• EM算法得到的估计序列是否收敛，如果收敛，是否收敛到全局最大值或局部最大值？下面给出关于EM算法收敛性的两个定理。证明见《统计学习方法（第二版）》第九章

2.3.1 定理1

• 设 $P(Y|\theta)$ 为观测数据的似然函数，$\theta^{(i)}(i=1,2,...)$ 为EM算法得到的参数估计序列，$P(Y|\theta^{(i)})(i=1,2,...)$ 为对应的似然函数序列，则 $P(Y|\theta^{(i)})$ 是单调递增的，即
  $P(Y|\theta^{(i+1)}) \geq P(Y|\theta^{(i)})$

2.3.2 定理2

• 设 $P(Y|\theta)$ 为观测数据的似然函数，$\theta^{(i)}(i=1,2,...)$ 为EM算法得到的参数估计序列，$P(Y|\theta^{(i)})(i=1,2,...)$ 为对应的似然函数序列，$L(\theta^{(i)})$

  1. 如果 $P(Y|\theta)$ 有上界，则 $L(\theta^{(i)}) = logP(Y|\theta^{(i)})$ 收敛到某一值 $L^*$
  2. 在函数 $Q(\theta,\theta')$ 与 $L(\theta)$ 满足一定条件下（大多数情况下都满足），由 EM 算法得到的参数估计序列 $\theta^{(i)}$ 的收敛值 $\theta^*$ 是 $L(\theta)$

• 参考文献：Wu C . On the Convergence Properties of the EM Algorithm[J]. Annals of Statistics, 1983, 11(1):95-103.

2.3.3 小结

• EM算法的收敛性包含两层意思，前者不蕴含后者
1. 对数似然函数 $L(\theta^{(i)})$
2. 参数估计序列 $\theta^{(i)}$
• 定理仅保证参数估计序列 $\theta^{(i)}$ 收敛到对数似然函数序列 $L(\theta^{(i)})$。如果 $L(\theta^{(i)})$
• 在应用中，初值的选择变得非常重要，常用的办法是选取几个不同的初值进行迭代，对各个初始值加以比较，从中选择最好的。

3. 用EM算法解三硬币问题

• 了解EM算法原理之后，回看 1.1节 最后的三硬币问题求解过程，分析M步和E步的计算

• 把三硬币问题的 “观测数据” 和 “未观测数据” 分别表示为来自总体 $y$ 和 $z$ 的容量为 $n$ 的简单随机样本 $y = \{y_1,y_2,...,y_n\}^T,z = \{z_1,z_2,...,z_n\}^T$ ，观测值和未观测值分别为 $Y = \{Y_1,Y_2,...,Y_n\}^T,Z = \{Z_1,Z_2,...,Z_n\}^T$，考虑求模型参数 $\theta$

• 注意：$y_i,Y_i$ 和 $z_i,Z_i$

  1. 在带入观测值前（计算Q函数前），概率相关计算一律使用 $y_i,z_i$；带入观测值后（计算Q函数时）一律使用 $Y_i,Z_i$
  2. 求和符号 $\sum_{j=1}^N$ 是对每个随机变量$y_j$的取值或观测值$Y_j$求和；$\sum_z$

  第1、2节符号可能有不规范的，都按这个理解

3.1 E步

• 先求 $P(z|y_j,\theta^{(i)})$
  $P(z|y_j,\theta^{(i)}) = \left\{ \begin{aligned} &\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}} & ,if \space z = 1\\ &1-\mu_j^{(i+1)} &,if \space z=0 \end{aligned} \right.$
• 再求 $P(z,y_j|\theta) = P(y_j|z,\theta)P(z|\theta)$
  $P(z,y_j|\theta) = \left\{ \begin{aligned} &\pi p^{y_j} (1-p)^{1-y_j} & ,if \space z = 1 \\ &(1-\pi) q^{y_j} (1-q)^{1-y_j} & ,if \space z = 0 \end{aligned} \right.$

3.1.1 公式法直接计算 Q 函数

• 直接带入公式计算 $Q(\theta,\theta^{(i)})$
  $\begin{aligned} Q(\theta,\theta^{(i)}) &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ &= \sum_{j=1}^N \mathbb{E}_{Z|Y_j,\theta^{(i)}}logP(Y_j,Z|\theta) \\ &= \sum_{j=1}^N\{\sum_ZP(Z|Y_j,\theta^{(i)})logP(Y_j,Z|\theta)\} \\ &= \sum_{j=1}^N\{P(Z=1|Y_j,\theta^{(i)})log(Y_j,Z=1|\theta) + P(Z=0|Y_j,\theta^{(i)})log(Y_j,Z=0|\theta)\} \\ &= \sum_{j=1}^N\{\mu_j^{(i+1)}log[\pi p^{Y_j} (1-p)^{1-Y_j}] + (1-\mu_j^{(i+1)})log[(1-\pi) q^{Y_j} (1-q)^{1-Y_j}]\} \end{aligned}$

3.1.2 从期望角度计算 Q 函数

• 从期望角度计算
  $\begin{aligned} \because P(y_j,z|\theta) &= \left\{ \begin{aligned} &\pi p^{y_j} (1-p)^{1-y_j} & ,if \space z= 1 \\ &(1-\pi) q^{y_j} (1-q)^{1-y_j} & ,if \space z= 0 \end{aligned} \right. \\ \therefore logP(y_j,z|\theta) &= \left\{ \begin{aligned} &log[\pi p^{y_j} (1-p)^{1-y_j}] & ,if \space z= 1 \\ &log[(1-\pi) q^{y_j} (1-q)^{1-y_j}] & ,if \space z= 0 \end{aligned} \right. \\ \therefore \mathbb{E}_{z|y_i,\theta^{(i)}}logP(y_j,z|\theta) &= \left\{ \begin{aligned} & \mathbb{E}_{z=1|y_i,\theta^{(i)}}logP(y_j,z|\theta) & ,if \space z= 1 \\ & \mathbb{E}_{z=0|y_i,\theta^{(i)}}logP(y_j,z|\theta) & ,if \space z= 0 \end{aligned} \right. \\ &= \left\{ \begin{aligned} & P(z=1|y_i,\theta^{(i)}) log[\pi p^{y_j} (1-p)^{1-y_j}] & ,if \space z= 1 \\ & P(z=0|y_i,\theta^{(i)}) log[(1-\pi) q^{y_j} (1-q)^{1-y_j}] & ,if \space z= 0 \end{aligned} \right. \\ &= \left\{ \begin{aligned} & \mu_j^{(i+1)} log[\pi p^{y_j} (1-p)^{1-y_j}] & ,if \space z= 1 \\ & (1-\mu_j^{(i+1)}) log[(1-\pi) q^{y_j} (1-q)^{1-y_j}] & ,if \space z= 0 \end{aligned} \right. \\ \therefore Q(\theta,\theta^{(i)}) &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ &= \sum_{j=1}^N\{\mathbb{E}_{Z=0|Y_j,\theta^{(i)}}logP(Y_j,Z|\theta)+ \mathbb{E}_{Z=1|Y_j,\theta^{(i)}}logP(Y_j,Z|\theta)\}\\ &= \sum_{j=1}^N\{\mu_j^{(i+1)}log[\pi p^{Y_j} (1-p)^{1-Y_j}] + (1-\mu_j^{(i+1)})log[(1-\pi) q^{Y_j} (1-q)^{1-Y_j}]\} \end{aligned}$

3.2 M步

• 对 $Q(\theta,\theta^{(i)})$ 求偏导
  $\begin{aligned} \frac{\partial Q(\theta,\theta^{(i)})}{\partial \pi} &= \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{p^{Y_j}(1-p)^{1-Y_j} }{\pi p^{Y_j}(1-p)^{1-Y_j}} + (1-\mu_j^{(i+1)}) \frac{-q^{Y_j}(1-q)^{1-Y_j}}{(1-\pi)q^{Y_j}(1-q)^{1-Y_j}}\} \\ &= \sum_{j=1}^N\{\frac{\mu_j^{(i+1)}-\pi}{\pi(1-\pi)}\} \\ &=\frac{(\sum_{j=1}^N\mu_j^{(i+1)})-N\pi}{\pi(1-\pi)}\\ \frac{\partial Q(\theta,\theta^{(i)})}{\partial p} &= \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{\pi p^{Y_j-1}(1-p)^{-Y_j}(Y_j-p) }{\pi p^{Y_j}(1-p)^{1-Y_j}}\} \\ &= \sum_{j=1}^N\{\frac{\mu_j^{(i+1)}(Y_j-p)}{p(1-p)}\} \\ &=\frac{(\sum_{j=1}^N\mu_j^{(i+1)}Y_j)-(p\sum_{j=1}^N\mu_j^{(i+1)})}{p(1-p)}\\ \end{aligned}$
• 令偏导为 0 解 $\theta^{(i+1)}$
  $\begin{aligned} &\frac{\partial Q(\theta,\theta^{(i)})}{\partial \pi} \stackrel{令}{=}0 \Rightarrow \pi^{(i+1)} = \frac{1}{n} \sum_{i=1}^n \mu_j^{(i+1)} \\ &\frac{\partial Q(\theta,\theta^{(i)})}{\partial p}\stackrel{令}{=}0 \Rightarrow \left\{ \begin{aligned} &p^{(i+1)} = \frac{\sum_{i=1}^n \mu_j^{(i+1)}Y_j}{\sum_{i=1}^n \mu_j^{(i+1)}} \\ &q^{(i+1)} = \frac{\sum_{i=1}^n (1-\mu_j^{(i+1)})Y_j}{\sum_{i=1}^n (1-\mu_j^{(i+1)})} （对称性） \end{aligned} \right. \end{aligned}$

3.3 对比

• 对比 EM 和直接最大似然估计的优化目标
  $\begin{aligned} &Q(\theta,\theta^{(i)}) = \sum_{j=1}^N\{\mu_j^{(i+1)}log[\pi p^{Y_j} (1-p)^{1-Y_j}] + (1-\mu_j^{(i+1)})log[(1-\pi) q^{Y_j} (1-q)^{1-Y_j}]\} \\ &logP(Y|\theta) = \sum_{j=1}^Nlog[\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}]\\ \end{aligned}$
• 对比 EM 和直接最大似然估计的偏导数方程组（因为最后要令等于0，重点对比分子即可）
  $\begin{aligned} &\left\{ \begin{aligned} &\frac{\partial Q(\theta,\theta^{(i)})}{\partial \pi} = \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{p^{Y_j}(1-p)^{1-Y_j} }{\pi p^{Y_j}(1-p)^{1-Y_j}} + (1-\mu_j^{(i+1)}) \frac{-q^{Y_j}(1-q)^{1-Y_j}}{(1-\pi)q^{Y_j}(1-q)^{1-Y_j}}\} \\ &\frac{\partial Q(\theta,\theta^{(i)})}{\partial p} = \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{\pi p^{Y_j-1}(1-p)^{-Y_j}(Y_j-p) }{\pi p^{Y_j}(1-p)^{1-Y_j}}\} \end{aligned} \right.\\ &\left\{ \begin{aligned} \frac{\partial logP(Y|\theta)}{\partial \pi} & = \sum_{j=1}^n \frac{p^{Y_j}(1-p)^{1-Y_j}-q^{Y_j}(1-q)^{1-Y_j}}{\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}}\\ \frac{\partial logP(Y|\theta)}{\partial p} & = \sum_{j=1}^n \frac{\pi p^{Y_j-1}(1-p)^{-Y_j}(Y_j-p)}{\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}} \end{aligned} \right. \end{aligned}$

3.4 代码

• 最后给出三硬币问题的代码。网上搜到的基本都是同一份代码的复制粘贴，其实那份代码是适用于 jupyter notebook 那种交互式环境的，用了比较少见的 ​​yield​​ 关键字，可读性也比较差

• 我按照 1.1.2 节给出的流程重新写了一遍，如下

  import numpy as np
  import math
  
  class EM:
      def __init__(self, theta, y):
          self.theta = theta
          self.y = y
          self.mu = [0]*len(y)
  
          print('init pi={:.3f}; p={:.3f}; q={:.3f}'.format(theta[0],theta[1],theta[2]))
      
      # E_step
      def EStep(self):
          for j in range(len(self.y)):
              yj = self.y[j]
              pi,p,q = self.theta
              P_B_yj = pi * math.pow(p, yj) * math.pow((1-p), 1-yj)
              P_C_yj = (1-pi) * math.pow(q, yj) * math.pow((1-q), 1-yj)
              
              self.mu[j] = P_B_yj/(P_B_yj+P_C_yj)
  
      # m_step
      def MStep(self):
          sum_mu,sum_muy,sum_y = 0,0,0
          N = len(self.y)
          for j in range(N):
              sum_mu += self.mu[j]
              sum_muy += self.mu[j]*self.y[j]
              sum_y += self.y[j]
  
          self.theta[0] = sum_mu/N                    # pi
          self.theta[1] = sum_muy/sum_mu              # p
          self.theta[2] = (sum_y-sum_muy)/(N-sum_mu)  # q
  
      def em(self):
          num = 1
          while True:
              delta = self.theta.copy()
              self.EStep()
              self.MStep()
              print('ite{} pi={:.3f}; p={:.3f}; q={:.3f}'.format(num,self.theta[0],self.theta[1],self.theta[2]))
              num += 1
  
              # 退出条件：theta变化量二范数小于阈值
              for i in range(3):
                  delta[i] -= self.theta[i]
              if np.linalg.norm(delta) < 0.005:
                  print('')
                  break
                  
  data=[1,1,0,1,0,0,1,0,1,1]
  
  test1 = EM([0.5, 0.5, 0.5],data)
  test1.em()
  
  test2 = EM([0.4, 0.6, 0.7],data)
  test2.em()
  
  '''
  init pi=0.500; p=0.500; q=0.500
  ite1 pi=0.500; p=0.600; q=0.600
  ite2 pi=0.500; p=0.600; q=0.600
  
  init pi=0.400; p=0.600; q=0.700
  ite1 pi=0.406; p=0.537; q=0.643
  ite2 pi=0.406; p=0.537; q=0.643
  '''