统计学习方法第九章：EM算法从浅入深（扩散模型的前身，强烈建议学习）

最新推荐文章于 2026-06-23 10:28:18 发布

原创最新推荐文章于 2026-06-23 10:28:18 发布 · 223 阅读

9 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#学习方法 #人工智能 #算法 #机器学习 #概率论

统计学习方法专栏收录该内容

12 篇文章

订阅专栏

第一部分：概率建模与隐变量的困境

在进入 EM 算法之前，我们必须先把脚下的土地踩实。这就意味着我们需要回到概率模型、极大似然估计这些最基本的概念，然后亲眼看看，为什么隐变量的出现会让原本优美的数学框架突然变得棘手。你将会发现，正是这个「和对数」的障碍，才催生了 EM 算法那巧妙的策略。

1.1 概率模型与参数

我们要学会用数学语言去讲述数据是如何产生的。这个讲述的工具，就是概率模型。

随机变量与概率分布
随机变量 $X$ 是一个取值不确定的量，但它取各个值的可能性由概率分布刻画。如果 $X$ 只能取有限个或可数个值，它就是离散型随机变量，用概率质量函数（PMF） $P(X=x)$ 描述；如果 $X$ 可以在一个区间内连续取值，它就是连续型随机变量，用概率密度函数（PDF） $f(x)$ 描述，此时 $P(a \le X \le b) = \int_a^b f(x)\,dx$ 。

参数的作用

一个概率分布的具体形状通常由一组参数 $\theta$ 控制。例如，均值和方差决定一个正态分布在哪儿、有多胖。我们把这样的分布记作 $P(X|\theta)$ 或 $f(x|\theta)$ 。参数一改变，整个概率故事就会随之改变。

独立同分布（i.i.d.）假设

当我们有一系列观测数据 $y_1, y_2, \dots, y_n$ 时，最常见的简化假设是这些数据是独立同分布的：每个数据点都来自同一个分布，且彼此独立。数学上，它们的联合概率就等于各个边缘概率的乘积：

$P(Y_1=y_1,\dots,Y_n=y_n|\theta) = \prod_{i=1}^n P(Y_i=y_i|\theta).$

实例：抛一枚不均匀硬币 假设我们有一枚可能不均匀的硬币，正面朝上的概率为 $p$ ，反面朝上的概率为 $1-p$ 。这里唯一的参数就是 $\theta = p$ 。抛一次的结果是一个伯努利随机变量：

$P(Y=1)=p, \quad P(Y=0)=1-p.$

在 i.i.d. 假设下，抛 $n$ 次得到一串结果 $y_1,\dots,y_n$ （每个 $y_i$ 是 0 或 1）的联合概率就是：

$P(y_1,\dots,y_n|p) = \prod_{i=1}^n p^{y_i} (1-p)^{1-y_i} = p^{\sum y_i} (1-p)^{n-\sum y_i}.$

你看，仅仅是换一个 p 的值，整串数据看上去的合理性就会截然不同。

一句话总结：概率模型是用参数化分布“解释”观测数据的数学故事。

1.2 极大似然估计 (MLE)

有了模型和参数，下一步就是反过来从观测数据推测最有可能的参数值。这就是极大似然估计的核心思想。

似然函数 = 联合概率

在拿到具体观测数据 $\mathbf{y} = (y_1,\dots,y_n)$ 之后，我们把模型的联合概率看作参数 $\theta$ 的函数，称为似然函数：

$L(\theta) = P(\mathbf{Y}=\mathbf{y}|\theta).$

它衡量的是：在参数 $\theta$ 下，这个已经发生的样本数据 $\mathbf{y}$ 出现的概率（或密度）有多大。

对数似然将乘法变加法 因为 $L(\theta)$ 通常是很多项的乘积，直接处理乘积在计算上很不方便（比如连乘容易导致数值下溢），所以我们几乎总是取自然对数，得到对数似然函数：

$\ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log P(Y_i=y_i|\theta).$

取对数不改变函数的单调性，所以最大化 $\ell(\theta)$ 等价于最大化 $L(\theta)$ 。

通过求导/数值优化寻找最大值

当 $\ell(\theta)$ 形式良好时，我们可以对 $\theta$ 求偏导（梯度），令其为零，然后解出参数的最优值。这就是分析解法。如果无法解析求解，就用梯度上升、牛顿法等数值优化方法。

实例：抛硬币10次，7次正面

沿用伯努利模型，设观测到 $\sum y_i = 7$ ，总次数 $n=10$ 。似然函数为

$L(p) = \binom{10}{7} p^7 (1-p)^3.$

对数似然为

$\ell(p) = \log\binom{10}{7} + 7\log p + 3\log(1-p).$

常数项 $\log\binom{10}{7}$ 不影响最优点，我们略去。求导：

$\frac{d\ell}{dp} = \frac{7}{p} - \frac{3}{1-p} = 0 \quad\Rightarrow\quad \frac{7}{p} = \frac{3}{1-p} \quad\Rightarrow\quad 7(1-p) = 3p \quad\Rightarrow\quad p = 0.7.$

直观上非常合理：10次中出现了7次正面，那么正面概率的最好估计就是 $0.7$ 。极大似然估计用最符合观测的比例来猜测参数。

一句话总结：MLE 就是寻找让已发生事件“显得最合理”的参数。

1.3 隐变量的引入

MLE 如此简洁优美，但它是建立在所有相关变量都能被观测的假设之上的。一旦某些关键变量隐藏起来，情况就急转直下。

不完全数据 vs. 完全数据

在很多实际问题中，我们只能观察到一部分变量 $\mathbf{Y}$ ，而另一部分变量 $\mathbf{Z}$ 是不可观测的，称为隐变量。如果能同时观测到 $\mathbf{Y}$ 和 $\mathbf{Z}$ ，就构成了完全数据；只有 $\mathbf{Y}$ 时，则是不完全数据。

隐变量导致似然函数中出现求和/积分在对数内

有了隐变量，观测数据的似然就必须通过边缘化来获得：

$P(\mathbf{Y}|\theta) = \sum_{\mathbf{z}} P(\mathbf{Y}, \mathbf{Z}=\mathbf{z}|\theta),$

（若 $\mathbf{Z}$ 连续，则是积分）。于是，对数似然变成

$\ell(\theta) = \log P(\mathbf{Y}|\theta) = \log \left[ \sum_{\mathbf{z}} P(\mathbf{Y}, \mathbf{Z}=\mathbf{z}|\theta) \right].$

注意，这个对数内部有一个求和（或积分），即所谓的「对数-和」结构。它无法像 1.2 节那样被拆成独立项之和，因为 $\log(\sum)$ 不是 $\sum \log$ 。

$\log\sum$ 难以直接求导优化 假设我们对参数 $\theta$ 求导，

$\frac{\partial}{\partial\theta} \log \sum_{\mathbf{z}} P(\mathbf{Y},\mathbf{z}|\theta) = \frac{ \sum_{\mathbf{z}} \frac{\partial}{\partial\theta} P(\mathbf{Y},\mathbf{z}|\theta) }{ \sum_{\mathbf{z}} P(\mathbf{Y},\mathbf{z}|\theta) }.$

这个表达式通常极其复杂，而且没有简单的闭式解。如果隐变量有 $K$ 种可能，那么求和符号内部就混叠着各种参数，极难分离。

实例：一个“混合”硬币实验（三硬币模型的雏形）

设想这样一个实验：桌上有两枚外观相同但可能不匀的硬币 A 和 B。每次实验时，你先以概率 $\pi$ 选择硬币 A，以概率 $1-\pi$ 选择硬币 B；然后抛掷所选硬币，记录其正面（1）或反面（0）的结果。我们只能看到最终的正反面序列 $\mathbf{y} = (y_1,\dots,y_n)$ ，却看不到每次是选了哪枚硬币。这个“选择哪枚硬币”就是隐变量 $z_i \in {A,B}$ 。

记硬币 A 的正面概率为 $p$ ，硬币 B 的正面概率为 $q$ 。整个模型的参数为 $\theta = (\pi, p, q)$ 。根据生成过程，第 $i$ 次观测的似然是混合形式：

$P(y_i|\theta) = \pi \cdot p^{y_i}(1-p)^{1-y_i} \;+\; (1-\pi) \cdot q^{y_i}(1-q)^{1-y_i}.$

全部 $n$ 次观测的联合似然为乘积，再取对数得到：

$\ell(\pi,p,q) = \sum_{i=1}^n \log\Bigl[\,\pi p^{y_i}(1-p)^{1-y_i} + (1-\pi) q^{y_i}(1-q)^{1-y_i}\Bigr].$

对数内部是两项之和，因此即使对单个参数（比如 $p$ ）求导，也会得到形如

$\frac{\partial \ell}{\partial p} = \sum_{i=1}^n \frac{ \pi(\frac{y_i}{p} - \frac{1-y_i}{1-p}) p^{y_i}(1-p)^{1-y_i} }{\pi p^{y_i}(1-p)^{1-y_i} + (1-\pi) q^{y_i}(1-q)^{1-y_i}}$

的繁琐形式。令其为零根本无法得到 $p$ 的显式解。如果隐变量 $z_i$ 是已知的，我们只需把数据按硬币种类分成两组，分别用伯努利 MLE 即可轻松得到 $p$ 和 $q$ 。但现在 $z_i$ 缺失了，我们直接最大化 $\ell(\theta)$ 就陷入了对数-和的泥潭。

这就是 EM 算法登场的根本动机。

一句话总结：隐变量使得对数似然中出现“和对数”，打破了直接求导的优雅，成为 EM 算法登场的动机。

到这里，我们已经用概念搭建好了舞台，并亲眼见证了隐变量如何制造出 MLE 无法直接处理的困境。下一部分，我们将通过三硬币模型的迭代猜想，直观地感受 EM 算法如何用“猜测—再猜测”的策略从困境中开辟出一条向上的路。

第二部分：EM 算法的直觉

前面我们看到，隐变量的存在让对数似然函数变成了 $\log\sum$ 的形式，MLE 的求导之路被堵死了。但不要绝望：如果隐变量的真实值是已知的，问题就会退化为简单的分组计数。这就催生了 EM 算法的核心直觉——“如果不知道，就猜一个；猜了之后再去优化参数；优化了参数再回头修正猜测”。本节我们以三硬币模型为舞台，一步步建立这种“猜测与修正”的迭代逻辑。

2.1 三硬币模型的全概率描述

首先，我们把上一节提到的“混合硬币实验”精确化，建立完整的概率生成故事。这是整个 EM 入门最经典的 toy model。

原子：观测数据生成过程（先选硬币再抛） 假设有三枚硬币，分别记为 A、B、C。它们的正面向上的概率分别为 $\pi$ 、 $p$ 、 $q$ 。每次观测由以下两步生成：

抛掷硬币 A：
- 若 A 为正面（概率 $\pi$ ），则选择硬币 B 进行抛掷。
- 若 A 为反面（概率 $1-\pi$ ），则选择硬币 C 进行抛掷。
抛掷所选硬币（B 或 C），记录结果：正面记为 1，反面记为 0。

我们独立重复这个实验 $n$ 次，得到观测序列 $\mathbf{y} = (y_1, y_2, \dots, y_n)$ ，其中 $y_i \in {0, 1}$ 。但是，我们观察不到每次是选择了硬币 B 还是硬币 C，也无法直接看到硬币 A 的结果——真正观察不到的是“被选的硬币”这个中间状态。

原子：参数 $\pi, p, q$ 的角色 此模型的全部参数为 $\theta = (\pi, p, q)$ ，其中：

$\pi$ ：硬币 A 正面的概率，也就是选择硬币 B 的概率。
$p$ ：硬币 B 正面的概率。
$q$ ：硬币 C 正面的概率。

如果我们把“选哪枚硬币”视为隐变量 $z_i$ ，不妨约定：

$z_i = 1$ 表示选择硬币 B（即硬币 A 正面），
$z_i = 0$ 表示选择硬币 C（即硬币 A 反面）。

那么，给定隐变量 $z_i$ ，观测 $y_i$ 的条件分布是：

$P(y_i | z_i=1, \theta) = p^{y_i}(1-p)^{1-y_i}$

$P(y_i | z_i=0, \theta) = q^{y_i}(1-q)^{1-y_i}$

而隐变量本身的先验分布为：

$P(z_i=1|\theta) = \pi, \quad P(z_i=0|\theta) = 1-\pi.$

于是，观测 $y_i$ 的边缘分布就是混合形式：

$P(y_i|\theta) = \pi \cdot p^{y_i}(1-p)^{1-y_i} + (1-\pi) \cdot q^{y_i}(1-q)^{1-y_i}$ （2.1）

实例：假设 $\pi=0.4, p=0.8, q=0.3$ ，模拟生成一串数据 按照上述生成过程，我们可以编写一个随机模拟：

第 1 次：抛硬币 A，假设出现正面（概率 $\pi=0.4$ ），于是选用 B 抛掷。B 抛掷结果为 1（概率 0.8），记录 $y_1=1$ ， $z_1=1$ 。
第 2 次：抛硬币 A，出现反面（概率 0.6），选用 C 抛掷，C 抛出 0（概率 0.7），记录 $y_2=0$ ， $z_2=0$ 。
... 反复 20 次，可能得到如下观测序列（1/0串），但隐变量 $z_i$ 对我们不可见。

此时，一个不知情的观测者只拥有一长串 0 和 1，他的目标是反推出这三个概率 $\pi, p, q$ 。这就是典型的“不完全数据”下的参数估计问题。

一句话总结：三硬币模型是隐变量模型的“Hello World”，观测结果背后隐藏着“选了哪枚硬币”的未观察选择。

2.2 直接极大化的死胡同

我们已经有了观测数据的似然函数2.1及其乘积形式。让我们明确地写出对数似然，正视它带来的困境。

展示对数似然 $\log P(\mathbf{Y}|\theta)$ 中无法分离的求和项 对于 $n$ 次独立观测，对数似然为：

$\ell(\theta) = \sum_{i=1}^n \log\bigl[ \pi p^{y_i}(1-p)^{1-y_i} + (1-\pi) q^{y_i}(1-q)^{1-y_i} \bigr]$ （2.2）

注意，这里的 $\log$ 函数内部是两项之和。如果我们试图对某个参数，比如 $p$ ，求偏导，会得到什么？利用链式法则：

$\frac{\partial \ell}{\partial p} = \sum_{i=1}^n \frac{ \pi \cdot \frac{\partial}{\partial p}\bigl[p^{y_i}(1-p)^{1-y_i}\bigr] }{ \pi p^{y_i}(1-p)^{1-y_i} + (1-\pi) q^{y_i}(1-q)^{1-y_i} }.$

由于分母也是两项之和，我们无法将 $\pi$ 和 $p$ 、 $q$ 解耦。令导数为零后，方程里 $p$ 同时出现在分子和分母的复杂分式中，得不到像“ $p = \frac{\sum y_i}{n}$ ”这样的显式解。

实例：具体写出序列概率表达式，指出 $\pi, p, q$ 被锁在 $\log\sum$ 中 假设我们仅仅有 5 次观测： $\mathbf{y} = (1,0,0,1,1)$ 。则对数似然为：

$\ell = \log(\pi p + (1-\pi)q) + \log(\pi(1-p) + (1-\pi)(1-q)) + \log(\pi(1-p) + (1-\pi)(1-q)) + \log(\pi p + (1-\pi)q) + \log(\pi p + (1-\pi)q).$

试图通过解 $\frac{\partial \ell}{\partial \pi}=0, \frac{\partial \ell}{\partial p}=0, \frac{\partial \ell}{\partial q}=0$ 来求极值，会发现得到的是一个超越方程组，没有解析解。当然可以用数值优化（例如梯度上升）强行求解，但那样计算复杂，且丢失了对问题结构的洞察。更致命的是，在实际更复杂的隐变量模型中，这种直接数值优化往往极其低效。

一句话总结：显含隐变量的模型使对数似然内部的和式无法拆开，逼迫我们寻找变通的优化策略。

2.3 启发式迭代

既然直接攻破“和对数”防线很困难，我们能不能曲线救国？突破口在于一个简单的观察：如果隐变量 $z_i$ 不再是隐藏的，一切就变得异常简单。

原子：若已知隐变量的状态，问题退化为简单计数 假设在某个平行宇宙里，我们不仅记录了每次抛掷的结果 $y_i$ ，还记录了当时是硬币 B 还是硬币 C 抛出的（即 $z_i$ ）。那么，参数估计就退化为三个独立的伯努利实验：

硬币 A 的正面概率 $\pi$ ：直接统计 $z_i=1$ 的比例即可， $\hat{\pi} = \frac{\sum_{i=1}^n z_i}{n}$ 。
硬币 B 的正面概率 $p$ ：把所有 $z_i=1$ 的那些次观测挑出来，统计其中 $y_i=1$ 的比例， $\hat{p} = \frac{\sum{i=1}^n z_i y_i}{\sum{i=1}^n z_i}$ 。
硬币 C 的正面概率 $q$ ：类似地， $\hat{q} = \frac{\sum{i=1}^n (1-z_i) y_i}{\sum{i=1}^n (1-z_i)}$ 。

这不过是分组之后分别做 MLE，毫无难度。可惜现实中 $z_i$ 是缺失的。但我们不禁想：能不能先猜一个 $z_i$ ？

以当前参数估计隐变量的后验期望，替代硬分配

一种幼稚的做法是“硬分配”：先随便设定一组参数，根据参数决定每个样本最可能来自哪个硬币（例如若 $\pi p > (1-\pi) q$ 就认为 $z_i=1$ ），然后按上述公式重新估计参数，再回头修正硬分配。但这样的硬分配容易造成信息丢失（尤其在模棱两可的点上），且可能不收敛或不稳定。

更合理的方法是“软分配”：不强行把每个样本划归为 B 或 C，而是计算在当前参数下，每个样本 $i$ 属于硬币 B 的后验概率（也就是“责任”）：

$\gamma_i = P(z_i=1 | y_i, \theta) = \frac{\pi p^{y_i}(1-p)^{1-y_i}}{\pi p^{y_i}(1-p)^{1-y_i} + (1-\pi) q^{y_i}(1-q)^{1-y_i}}.{2.3}$

同理，属于硬币 C 的后验概率为 $1-\gamma_i$ 。这个 $\gamma_i$ 是一个介于 0 和 1 之间的权重，表示“这个样本 $i$ 有多大概率是由硬币 B 生成的”。用这些权重替代硬分配的 0/1 标签，我们就可以“软加权”地使用上面的计数公式：

估计 $\pi$ ：不是用 $\sum z_i$ ，而是用 $\sum \gamma_i$ 作为“软计数”，那么 $\hat{\pi} = \frac{\sum \gamma_i}{n}$ 。
估计 $p$ ：用 $\gamma_i$ 加权正面的次数， $\hat{p} = \frac{\sum \gamma_i y_i}{\sum \gamma_i}$ 。
估计 $q$ ：用 $1-\gamma_i$ 加权， $\hat{q} = \frac{\sum (1-\gamma_i) y_i}{\sum (1-\gamma_i)}$ 。

这就形成了一个绝妙的循环：

给定参数 $\theta^{(t)} = (\pi^{(t)}, p^{(t)}, q^{(t)})$ ，我们可以用公式 (2.3) 算出每一个样本的“责任” $\gamma_i^{(t)}$ 。
然后，把这些责任当作“软标签”，用上述加权公式算出一组新的参数 $\theta^{(t+1)}$ 。
然后用新参数重新计算责任，如此反复。

实例：随机初始化参数，手工迭代两轮，观察参数向真实值靠近

回到前面 2.1 节中用真实参数 $\pi=0.4, p=0.8, q=0.3$ 模拟出的 20 次观测数据。现在我们假装不知道真实值，随机初始化参数：

初始 $\pi^{(0)} = 0.5, p^{(0)} = 0.5, q^{(0)} = 0.5$ 。

第一轮迭代（E步）：对每个 $y_i$ ，计算 $\gamma_i^{(0)} = \frac{0.5 \times 0.5^{y_i}0.5^{1-y_i}}{0.5 \times 0.5^{y_i}0.5^{1-y_i} + 0.5 \times 0.5^{y_i}0.5^{1-y_i}} = 0.5$ 。

也就是说，在初始完全无信息时，每个样本的权重都是 0.5。

然后 M 步： $\pi^{(1)} = \frac{\sum 0.5}{20} = 0.5$ 。 $p^{(1)} = \frac{\sum 0.5 y_i}{\sum 0.5} = \frac{\sum y_i}{20}$ —— 假设 20 次里正面出现了 12 次，则 $p^{(1)} = 0.6$ 。 $q^{(1)} = \frac{\sum 0.5 y_i}{\sum 0.5} = 0.6$ ，两者相同，没有区分。

第一次迭代结束时参数变成了 $(0.5,0.6,0.6)$ ，区分度尚未显现。

第二轮迭代（E步）：用新参数 $\pi=0.5, p=0.6, q=0.6$ 重新计算 $\gamma_i$ 。对于一次正面（ $y=1$ ）：

$\gamma_i = \frac{0.5 \times 0.6}{0.5 \times 0.6 + 0.5 \times 0.6} = 0.5.$

对于反面（ $y=0$ ）：

$\gamma_i = \frac{0.5 \times 0.4}{0.5 \times 0.4 + 0.5 \times 0.4} = 0.5.$

还是 0.5！看来当 $p=q$ 时，模型陷入对称陷阱，权重无法区分。这是初始化敏感的一个案例。我们需要打破对称性。

现在换一个稍微不对称的初始值： $\pi^{(0)}=0.4, p^{(0)}=0.6, q^{(0)}=0.3$ 。仍然用 20 次模拟数据（真实值 0.4,0.8,0.3，正面次数可能占多）。我们只展示其中的 5 次观测： $y = (1,0,1,1,0)$ 。第一轮 E 步：

对于 $y=1$ ，

$\gamma = \frac{0.4 \times 0.6}{0.4 \times 0.6 + 0.6 \times 0.3} = \frac{0.24}{0.24+0.18} \approx 0.571.$

对于 $y=0$ ，

$\gamma = \frac{0.4 \times 0.4}{0.4 \times 0.4 + 0.6 \times 0.7} = \frac{0.16}{0.16+0.42} \approx 0.276.$

得到 $\gamma$ 序列：约 0.571, 0.276, 0.571, 0.571, 0.276。 M 步：

$\pi^{(1)} = \frac{0.571+0.276+0.571+0.571+0.276}{5} \approx 0.453$

$\sum \gamma_i \approx 2.265$ $p^{(1)} = \frac{0.571\cdot1 + 0.276\cdot0 + 0.571\cdot1 + 0.571\cdot1 + 0.276\cdot0}{2.265} = \frac{1.713}{2.265} \approx 0.756$ $q^{(1)} = \frac{(1-0.571)\cdot1 + (1-0.276)\cdot0 + \dots}{5 - 2.265} \approx \frac{0.429+0+0.429+0.429+0}{2.735} \approx 0.470$

可以看到，经过仅仅一次迭代， $p$ 从 0.6 跳升到了 0.756，已经开始往真实值 0.8 靠近； $q$ 从 0.3 升到 0.47（这是因为初始 $q$ 给得太低，数据中来自 C 的正样本将其拉高）。 $\pi$ 也做了调整。如果继续迭代，参数会逐渐收敛到真实值附近（但可能停在局部极值或鞍点）。

再来一个简单的例子：

假设观测到一个结果：1。我们随机瞎猜一组初始参数： $\pi=0.4, p=0.6, q=0.7$ 。

猜测（E步）：如果结果是 1，它来自硬币 B 的概率是多少？

来自路径 B 的概率分子是 $\pi \cdot p = 0.4 \times 0.6 = 0.24$ 。

来自路径 C 的概率分子是 $(1-\pi) \cdot q = 0.6 \times 0.7 = 0.42$ 。

那么它来自 B 的相对概率就是 $0.24 / (0.24 + 0.42) \approx 0.36$ 。

这说明，针对这个观测结果 1，硬币 A 有 36% 的可能是正面，64% 的可能是反面。
最大化（M步）：我们把这个 36% 和 64% 的“软计数”当作真实数据，回过头去更新我们的参数 $\pi, p, q$ 。
循环：拿着新的 $\pi, p, q$ ，再次重新评估每个数据点来自 B 还是 C 的概率……你会发现，经历几轮迭代后，参数奇迹般地向着真实的底层设定靠拢了。

原子：迭代流程收敛，参数不再明显改变 随着反复进行 E 步和 M 步，参数序列 $\theta^{(0)}, \theta^{(1)}, \theta^{(2)}, \dots$ 最终会稳定下来。每一步都保证对数似然 $\ell(\theta)$ 不下降（我们将在后面证明）。当参数变化极小或对数似然增量低于某个阈值（如 $10^{-6}$ ）时，我们停止迭代。

这整个流程就是 EM 算法在三硬币模型上的朴素实现。它没有直接去啃那个 $\log\sum$ 的硬骨头，而是通过引入“软标签” $\gamma_i$ ，构造了一个容易优化的替代目标（每次迭代中，伪似然是加权的伯努利似然），然后交替更新。

一句话总结：EM 算法的朴素思想就是“参数不好猜？先猜隐变量；隐变量不好猜？先猜参数”，二者交替提升。

第三部分：数学基石 —— Jensen 不等式与下界构造

第二部分结束时，我们用三硬币模型体验了 EM 算法“猜测-修正”的迭代过程。但这个软加权的把戏为什么在数学上有效？为什么对数似然能单调上升？要回答这些问题，我们必须深入 EM 算法的底层数学结构。一切的秘密都藏在一个看似不起眼的不等式中——Jensen 不等式，以及由它构造出的似然函数下界。

3.1 Jensen 不等式的解释

Jensen 不等式是凸分析中最基本也最强大的工具之一。它刻画了凸函数（或凹函数）与期望运算之间的不等关系。EM 算法的整个推导都可以看作是对它的巧妙应用。

凸函数的定义（弦在函数上方/下方）

先看凸函数。直观地说，一个函数 $f(x)$ 如果“向上弯”，它就是凸的。形式化定义：对于区间上任意两点 $x_1, x_2$ 以及权重 $\lambda \in [0,1]$ ，有

$f(\lambda x_1 + (1-\lambda) x_2) \le \lambda f(x_1) + (1-\lambda) f(x_2).$

这意味着连接曲线上任意两点的弦位于函数图像的上方（或重合）。典型的凸函数有 $f(x)=x^2$ 、 $f(x)=e^x$ 。反过来，如果不等号方向相反，即 $f(\lambda x_1 + (1-\lambda) x_2) \ge \lambda f(x_1) + (1-\lambda) f(x_2)$ ，则函数是凹函数，如 $f(x)=\log x$ 、 $f(x)=\sqrt{x}$ 。

概率权重下的期望视角

如果将 $\lambda$ 理解为概率：取 $x_1$ 的概率为 $\lambda$ ， $x_2$ 的概率为 $1-\lambda$ ，那么 $\lambda x_1+(1-\lambda)x_2$ 就是随机变量 $X$ 的期望 $\mathbb{E}[X]$ 。 $\lambda f(x_1)+(1-\lambda)f(x_2)$ 就是 $f(X)$ 的期望 $\mathbb{E}[f(X)]$ 。于是凸函数的定义可以推广到任意概率分布：

$f(\mathbb{E}[X]) \le \mathbb{E}[f(X)] \quad .$

对于凹函数，不等号反向：

$f(\mathbb{E}[X]) \ge \mathbb{E}[f(X)] \quad$

这就是 Jensen 不等式。它的精髓在于：凸函数把期望“拉高”了，凹函数把期望“压低”了。

实例： $f(x) = -\log x$ 在两个点上的数值验证 取凹函数 $\log x$ ，则 $-\log x$ 是凸函数。我们来验证凸函数的 Jensen 不等式。设 $x_1=1, x_2=3$ ，概率权重各 0.5。那么

$\mathbb{E}[X] = 0.5\times 1 + 0.5\times 3 = 2.$

$f(\mathbb{E}[X]) = -\log 2 \approx -0.693.$

$\mathbb{E}[f(X)] = 0.5\times(-\log 1) + 0.5\times(-\log 3) = 0 + 0.5\times(-1.099) = -0.5495.$

显然 $-0.693 \le -0.5495$ ，满足 $f(\mathbb{E}[X]) \le \mathbb{E}[f(X)]$ 。同时我们也观察到，对于凹函数 $\log x$

$\log(\mathbb{E}[X]) = \log 2 \approx 0.693, \quad \mathbb{E}[\log X] = 0.5\log 1 + 0.5\log 3 \approx 0.5495,$

所以 $\log(\mathbb{E}[X]) \ge \mathbb{E}[\log X]$ ，即凹函数的期望不大于函数在期望处的值。这正是我们后面要用到的形式。

一句话总结：Jensen 不等式规定，凸函数在期望处的值不大于函数值的期望，这是构造似然下界的“扳手”。

3.2 构建似然下界 (ELBO)

现在，我们手持 Jensen 不等式这把扳手，来拆解那个令人头疼的 $\log\sum$ 结构。目标是给对数似然 $\ell(\theta) = \log P(\mathbf{Y}|\theta)$ 找一个容易优化的“替身”——一个下界函数。

引入分布 $Q(Z)$

隐变量 $Z$ 可以取 $K$ 种离散值（或者连续，我们以离散为例）。我们引入一个关于 $Z$ 的任意概率分布 $Q(Z)$ ，满足：

$Q(Z) \ge 0, \quad \sum_{Z} Q(Z) = 1.$

$Q(Z)$ 可以是任何合法的概率质量函数，它只是我们暂时借用的“工具分布”。

对数似然的恒等变换

观测数据的对数似然为：

$\ell(\theta) = \log P(\mathbf{Y}|\theta) = \log \sum_{Z} P(\mathbf{Y}, Z|\theta).$

现在，我们往求和符号里塞进 $Q(Z)$ ，用一种“乘以1”的技巧：

$\ell(\theta) = \log \sum_{Z} Q(Z) \frac{P(\mathbf{Y}, Z|\theta)}{Q(Z)}$

为什么这么做？因为 $\sum_Z Q(Z) \cdot \frac{P(\mathbf{Y},Z|\theta)}{Q(Z)}$ 可以看作是对随机变量 $R = \frac{P(\mathbf{Y},Z|\theta)}{Q(Z)}$ 在分布 $Q(Z)$ 下的期望：

$\ell(\theta) = \log \mathbb{E}_{Z \sim Q}\left[ \frac{P(\mathbf{Y}, Z|\theta)}{Q(Z)} \right].$

原子：应用 Jensen 不等式（凹函数版本） $\log$ 是凹函数，根据 Jensen 不等式，对于凹函数 $f$ 有 $f(\mathbb{E}[R]) \ge \mathbb{E}[f(R)]$ 。代入 $f = \log$ ，我们立刻得到：

$\log \mathbb{E}_{Q}\left[ \frac{P(\mathbf{Y}, Z|\theta)}{Q(Z)} \right] \ge \mathbb{E}_{Q}\left[ \log \frac{P(\mathbf{Y}, Z|\theta)}{Q(Z)} \right]$

右边这一坨，我们把它展开：

$\mathbb{E}_{Q}\left[ \log \frac{P(\mathbf{Y}, Z|\theta)}{Q(Z)} \right] = \sum_{Z} Q(Z) \log P(\mathbf{Y}, Z|\theta) - \sum_{Z} Q(Z) \log Q(Z).$

这就是一个关于 $\theta$ 和 $Q$ 的函数，我们记为 $\mathcal{L}(Q, \theta)$ ：

$\mathcal{L}(Q, \theta) = \sum_{Z} Q(Z) \log P(\mathbf{Y}, Z|\theta) - \sum_{Z} Q(Z) \log Q(Z)$

注意，第二项就是 $Q$ 分布的熵 $H(Q)$ 的负数（在信息论中熵是 $-\sum Q\log Q$ ），因此 $\mathcal{L}$ 也可以写成：

$\mathcal{L}(Q, \theta) = \mathbb{E}_{Q}[\log P(\mathbf{Y}, Z|\theta)] + H(Q).$

对于任意的分布 $Q$ 和参数 $\theta$ ，都有：

$\ell(\theta) \ge \mathcal{L}(Q, \theta).$

$\mathcal{L}(Q, \theta)$ 就是对数似然 $\ell(\theta)$ 的一个下界（lower bound）。在机器学习文献中，它经常被称为 ELBO（Evidence Lower BOund，证据下界）。

实例：在三硬币模型中代入一个任意 Q，展示下界低于对数似然

回到三硬币模型。假设我们只有一次观测 $y=1$ ，参数 $\pi=0.4, p=0.8, q=0.3$ 。那么真实对数似然为：

$\ell = \log( \pi p + (1-\pi) q ) = \log(0.4\times0.8 + 0.6\times0.3) = \log(0.32 + 0.18) = \log 0.5 \approx -0.6931.$

现在，我们随便取一个 Q分布。隐变量 Z 表示选择硬币，

取 $Q(Z=B)=0.2, Q(Z=C)=0.8$

计算下界：

$\mathcal{L} = 0.2 \log P(y=1, Z=B|\theta) + 0.8 \log P(y=1, Z=C|\theta) - [0.2\log 0.2 + 0.8\log 0.8].$

其中， $P(y=1, Z=B|\theta) = \pi p = 0.32$ ， $P(y=1, Z=C|\theta) = (1-\pi) q = 0.18$ 。

代入：

$\mathcal{L} = 0.2 \log 0.32 + 0.8 \log 0.18 - (0.2\log 0.2 + 0.8\log 0.8).$

计算：

$0.2 \log 0.32 = 0.2 \times (-1.1394) = -0.2279$ , $0.8 \log 0.18 = 0.8 \times (-1.7148) = -1.3718$ , 求和 $= -1.5997$ 。

第二项（负熵）： $-(0.2\log 0.2 + 0.8\log 0.8) \approx -(0.2\times -1.6094 + 0.8\times -0.2231) = -(-0.3219 - 0.1785) = 0.5004$ 所以 $\mathcal{L} \approx -1.5997 + 0.5004 = -1.0993$ 。显然 $-1.0993 \le -0.6931$ ，下界确实在对数似然的下方。如果调整 $Q$ ，这个下界会上下浮动，但永远不会超过 $\ell(\theta)$ 。

一句话总结：通过对隐变量引入任意分布 Q，Jensen 不等式将讨厌的 $\log\sum$ 转换为期望 $\log$ ，生成一个容易优化的下界函数。

3.3 下界收紧条件： $Q(Z) = P(Z|\mathbf{Y}, \theta)$

$\mathcal{L}(Q,\theta)$ 只是一个下界，它可能离真实的 $\ell(\theta)$ 很远。如果我们要用优化下界来间接爬升真实似然，自然希望每一步下界至少在当前参数处是“紧贴”真实似然的。否则，下界提升了，真实似然未必跟着升。那么，什么时候这个下界能够取到等号，紧紧地贴上对数似然？

Jensen 不等式取等条件

Jensen 不等式 $\log \mathbb{E}[R] \ge \mathbb{E}[\log R]$ 取等的条件是：随机变量 $R$ 在其分布下几乎处处为常数（即方差为零）。换句话说， $\frac{P(\mathbf{Y}, Z|\theta)}{Q(Z)}$ 必须是一个与 $Z$ 无关的常数 $c$ 。

假设存在常数 $c$ ，使得对所有可能的 $Z$ ，有

$\frac{P(\mathbf{Y}, Z|\theta)}{Q(Z)} = c.$

那么 $Q(Z) = \frac{1}{c} P(\mathbf{Y}, Z|\theta)$ 。又因为 $\sum_Z Q(Z) = 1$ ，可得 $c = \sum_Z P(\mathbf{Y}, Z|\theta) = P(\mathbf{Y}|\theta)$ 。代入回去：

$Q(Z) = \frac{P(\mathbf{Y}, Z|\theta)}{P(\mathbf{Y}|\theta)} = P(Z|\mathbf{Y}, \theta).$

这就是我们熟悉的后验分布！结论非常明确：当且仅当 $Q(Z)$ 等于在当前参数 $\theta$ 下隐变量的后验概率分布时，Jensen 不等式取等号，即 $\mathcal{L}(Q,\theta) = \ell(\theta)$ 。

后验概率 $P(Z|\mathbf{Y},\theta)$ 是使下界紧贴似然的唯一 $Q$ 因此，如果我们固定参数 $\theta^{(t)}$ ，然后选择 $Q^{(t+1)}(Z) = P(Z|\mathbf{Y}, \theta^{(t)})$ ，那么在点 $\theta^{(t)}$ 处，下界函数会与对数似然函数相切（函数值相等，且梯度也往往一致，虽然这需要额外条件，但至少是紧贴的）。这就提供了一条清晰的思路：

在当前的参数 $\theta^{(t)}$ ，用后验 $P(Z|\mathbf{Y}, \theta^{(t)})$ 更新 $Q$ ，使下界上升到与 $\ell(\theta)$ 平齐；
然后固定这个 $Q$ ，优化 $\theta$ 使下界 $\mathcal{L}(Q,\theta)$ 最大化。因为 $\mathcal{L}$ 在 $\theta^{(t)}$ 处已经等于 $\ell(\theta^{(t)})$ ，任何对 $\mathcal{L}$ 的提升都会使 $\ell(\theta)$ 至少获得同等的提升（因为 $\ell \ge \mathcal{L}$ ）。

实例：在三硬币模型中，给定参数计算后验，验证下界等于对数似然

回到 3.2 节的单样本例子 y=1，参数 $\pi=0.4, p=0.8, q=0.3$ 。真实后验概率为：

$P(Z=B|y=1, \theta) = \frac{\pi p}{\pi p + (1-\pi) q} = \frac{0.32}{0.5} = 0.64,$

$P(Z=C|y=1, \theta) = 1 - 0.64 = 0.36.$

现在我们取 $Q(B)=0.64, Q(C)=0.36$ ，计算下界：

$\mathcal{L} = 0.64 \log 0.32 + 0.36 \log 0.18 - (0.64\log 0.64 + 0.36\log 0.36).$

各项数值： $\mathcal{L} = -1.3465 + 0.6534 = -0.6931$ 。正是 $\log 0.5$ 的值！下界完美地贴合了对数似然。

如果我们在这一点上增大下界（比如调整参数），由于下界从与真实似然相切的位置出发，真实似然至少也会有同等的增长——这就是 EM 算法收敛性背后的直观几何：先贴上去，再向上爬。

一句话总结：若取 Q 为当前参数下隐变量的后验分布，则下界在当前位置与对数似然函数相切，提供局部精确替代。

第三部分到这里结束。我们已经打下最核心的数学地基：Jensen 不等式如何从凸/凹性推导出来，又如何通过引入 $Q$ 分布将带隐变量的对数似然转化为一个可处理的下界 $\mathcal{L}(Q,\theta)$ ，并且明确指出这个下界在 Q 取后验分布时是紧的。这正是 EM 算法两步走——E步构造下界，M步抬高下界——的数学原理。下一部分，我们将把这些洞见凝固为 EM 算法的完整通用推导，让公式说话，给出严谨的 E 步和 M 步定义。

第四部分：EM 算法的严密推导

第三部分我们用 Jensen 不等式打造了一把“下界”钥匙，并揭示了它紧贴真实似然的条件。现在，我们要将这些几何直觉凝固为一套严谨的、可机械执行的算法框架。这一部分将给出 EM 算法的标准符号体系、E 步和 M 步的正式定义，并在三硬币模型上展示从抽象公式到具体迭代公式的完整推导链路。

4.1 通用框架与符号约定

在进入公式推导前，我们先统一语言。这能确保接下来的公式既精确又通用，适用于任何带隐变量的概率模型。

观测数据 $\mathbf{Y}$ ，隐变量 $\mathbf{Z}$ ，完全数据 $\mathbf{X}=(\mathbf{Y}, \mathbf{Z})$ 假设我们有观测数据（不完全数据） $\mathbf{Y}$ 。同时存在未观测的隐变量 $\mathbf{Z}$ ，它可以是离散的或连续的（连续时以下求和换为积分即可）。若我们能同时观测到 $\mathbf{Y}$ 和 $\mathbf{Z}$ ，则称 $(\mathbf{Y}, \mathbf{Z})$ 为完全数据。完全数据的联合概率分布为 $P(\mathbf{Y}, \mathbf{Z} | \theta)$ ，其中 $\theta$ 是模型参数。

目标最大化对数似然 $\log P(\mathbf{Y}|\theta)$

我们的终极目标不变：寻找参数 $\theta$ ，极大化观测数据的对数似然（或称证据下界对应的证据本身）：

$\ell(\theta) = \log P(\mathbf{Y}|\theta) = \log \sum_{\mathbf{Z}} P(\mathbf{Y}, \mathbf{Z}|\theta)$

由于 $\log\sum$ 的直接求导障碍，我们转而使用 EM 算法的两步迭代策略。

一句话总结：EM 算法用完全数据的对数似然期望作为辅助函数，间接攀爬不完全数据的似然山。

4.2 E 步：求 $Q$ 函数

E 步的全称是 Expectation step（期望步骤）。它的任务是利用当前参数 $\theta^{(t)}$ ，计算出隐变量的条件后验分布，并据此构造一个关于参数 $\theta$ 的“替身函数”，即 $Q$ 函数。

$Q(\theta|\theta^{(t)}) = \mathbb{E}_{\mathbf{Z}|\mathbf{Y},\theta^{(t)}} [\log P(\mathbf{Y}, \mathbf{Z}|\theta)]$

给定 $\theta^{(t)}$ ，我们先求出隐变量在当前观测下的后验分布 $P(\mathbf{Z}|\mathbf{Y}, \theta^{(t)})$ 。然后，我们定义 $Q$ 函数为：在这个后验分布下，完全数据对数似然 $\log P(\mathbf{Y}, \mathbf{Z}|\theta)$ 的期望：

$Q(\theta | \theta^{(t)}) = \sum_{\mathbf{Z}} P(\mathbf{Z}|\mathbf{Y}, \theta^{(t)}) \log P(\mathbf{Y}, \mathbf{Z}|\theta)$

注意，这里的 $\log$ 内部是联合概率 $P(\mathbf{Y}, \mathbf{Z}|\theta)$ ，它依赖于待优化的新参数 $\theta$ ；而加权求和的权重 $P(\mathbf{Z}|\mathbf{Y}, \theta^{(t)})$ 则基于旧参数 $\theta^{(t)}$ 计算，在 E 步中被视为常数。

本质是计算隐变量后验下的加权平均

如果隐变量 $\mathbf{Z}$ 是离散的， $Q$ 函数就是对所有可能的 $\mathbf{Z}$ 的求和；若连续则为积分。它本质上是给每个可能的隐变量状态“分发责任”，用这些责任为权重，计算完全数据对数似然的加权平均值。

实例：对三硬币模型写出 $Q$ 函数的闭式表达式

回顾三硬币模型。观测 $\mathbf{Y} = {y_i}{i=1}^n$ ，隐变量 $\mathbf{Z} = {z_i}{i=1}^n$ ，其中 $z_i \in {0,1}$ 表示第 $i$ 次是否选用硬币 B。参数 $\theta = (\pi, p, q)$ 。

完全数据的联合似然为：

$P(\mathbf{Y}, \mathbf{Z}|\theta) = \prod_{i=1}^n \bigl[ \pi p^{y_i}(1-p)^{1-y_i} \bigr]^{z_i} \bigl[ (1-\pi) q^{y_i}(1-q)^{1-y_i} \bigr]^{1-z_i}.$

取对数：

$\log P(\mathbf{Y}, \mathbf{Z}|\theta) = \sum_{i=1}^n \Bigl[ z_i \bigl( \log\pi + y_i\log p + (1-y_i)\log(1-p) \bigr) + (1-z_i) \bigl( \log(1-\pi) + y_i\log q + (1-y_i)\log(1-q) \bigr) \Bigr].$

在 E 步，我们需要在当前参数 $\theta^{(t)}$ 下，计算隐变量 $z_i$ 的后验期望。根据公式 (2.3)，每个样本的“责任”为：

$\gamma_i^{(t)} = P(z_i=1 | y_i, \theta^{(t)}) = \frac{\pi^{(t)} (p^{(t)})^{y_i} (1-p^{(t)})^{1-y_i}}{\pi^{(t)} (p^{(t)})^{y_i} (1-p^{(t)})^{1-y_i} + (1-\pi^{(t)}) (q^{(t)})^{y_i} (1-q^{(t)})^{1-y_i}}.$

同时 $P(z_i=0|y_i, \theta^{(t)}) = 1-\gamma_i^{(t)}$ 。

现在，我们将 $\log P(\mathbf{Y}, \mathbf{Z}|\theta)$ 中的 $z_i$ 和 $1-z_i$ 替换为它们的期望 $\gamma_i^{(t)}$ 和 $1-\gamma_i^{(t)}$ ，得到 Q 函数：

$Q(\theta|\theta^{(t)}) = \sum_{i=1}^n \Bigl[ \gamma_i^{(t)} \bigl( \log\pi + y_i\log p + (1-y_i)\log(1-p) \bigr) + (1-\gamma_i^{(t)}) \bigl( \log(1-\pi) + y_i\log q + (1-y_i)\log(1-q) \bigr) \Bigr]$

这个 Q函数非常友好：它彻底拆分了 $\log\sum$ ， $\pi, p, q$ 分别出现在独立的加法项中。这正是 E 步的精髓：用旧参数的后验期望，把复杂的对数边缘似然“偷换”为一个加权版本的完全数据对数似然。

一句话总结：E 步用旧参数算出隐变量的后验分布，据此构建完全数据对数似然的“期望替身”。

4.3 M 步：极大化 Q 函数

M 步的全称是 Maximization step（最大化步骤）。任务非常纯粹：把 E 步给出的 $Q$ 函数当作目标函数，寻找新的参数 $\theta$ 使其最大化。

$\theta^{(t+1)} = \arg\max_\theta Q(\theta|\theta^{(t)})$ 由于 Q 函数的结构通常比原始对数似然简单得多（不再包含 $\log\sum$ ），我们可以用常规的优化方法求解其最大值。如果 Q 函数是凸的或具有闭式解，我们可以直接求导并令导数为零。

无约束或有约束（如概率和为 1）的求导/拉格朗日法

在概率模型中，参数通常受限于概率的规范性约束（例如 $\pi + (1-\pi) = 1$ 自然满足，但有时参数如混合系数需显式约束和为 1）。此时，我们用拉格朗日乘子法。对于三硬币模型， $\pi, p, q$ 各自在 [0,1] 内，且似然函数里的 $\log\pi$ 和 $\log(1-\pi)$ 自然要求它们非零，但无需显式约束和为一的拉格朗日项（因为 $\pi$ 和 $1-\pi$ 已自动配对）。我们直接求偏导。

实例：对三硬币的 Q 函数求偏导，解出 $\pi, p, q$ 的迭代更新公式

我们将 Q 函数按参数分解：

$Q(\theta|\theta^{(t)}) = \sum_i \gamma_i^{(t)} \log\pi \;+\; \sum_i (1-\gamma_i^{(t)}) \log(1-\pi)$

$\;+\; \sum_i \gamma_i^{(t)} \bigl[ y_i\log p + (1-y_i)\log(1-p) \bigr]$

$\;+\; \sum_i (1-\gamma_i^{(t)}) \bigl[ y_i\log q + (1-y_i)\log(1-q) \bigr].$

求解 $\pi$ ：涉及 $\pi$ 的项为 $\bigl(\sum_i \gamma_i^{(t)}\bigr) \log\pi + \bigl(\sum_i (1-\gamma_i^{(t)})\bigr) \log(1-\pi)$ 。对 $\pi$ 求导并令为零：

$\frac{\partial Q}{\partial \pi} = \frac{\sum_i \gamma_i^{(t)}}{\pi} - \frac{\sum_i (1-\gamma_i^{(t)})}{1-\pi} = 0.$

解得：

$\pi^{(t+1)} = \frac{1}{n} \sum_{i=1}^n \gamma_i^{(t)}$

求解 $p$ ：涉及 $p$ 的项为 $\sum_i \gamma_i^{(t)} [ y_i\log p + (1-y_i)\log(1-p) ]$ 。对 $p$ 求导：

$\frac{\partial Q}{\partial p} = \sum_i \gamma_i^{(t)} \left( \frac{y_i}{p} - \frac{1-y_i}{1-p} \right) = 0.$

整理得：

$p^{(t+1)} = \frac{\sum_{i=1}^n \gamma_i^{(t)} y_i}{\sum_{i=1}^n \gamma_i^{(t)}}.$

求解 $q$ ：对称地，

$q^{(t+1)} = \frac{\sum_{i=1}^n (1-\gamma_i^{(t)}) y_i}{\sum_{i=1}^n (1-\gamma_i^{(t)})}. \$

这三条更新公式与我们在 2.3 节朴素推导中得到的软加权公式完全一致。但这一次，它们是从严格定义的 Q函数出发，通过最大化推导出来的，具有坚实的数学合法性。每一轮 M 步，我们用上轮 E 步得到的责任 $\gamma_i^{(t)}$ 刷新参数，使 Q 函数（也就是下界 ELBO 中去掉常数熵项的部分）达到最大。

一句话总结：M 步优化这个“期望替身”，由于替身是似然的下界且相切，参数更新必然推动真实似然上升。

4.4 算法总步与停止准则

现在我们把 E 步和 M 步串成完整的算法流程，并给出实用中判断收敛的准则。

原子：初始化 → 迭代 E/M 步 → 参数变化极小或似然收敛 EM 算法的执行步骤：

初始化：选择一个初始参数 $\theta^{(0)}$ 。初始值往往随机选取（或利用先验知识），因为 EM 可能收敛到局部极大值。
迭代：对于 $t = 0, 1, 2, \dots$ ：
- E 步：计算后验分布 $P(\mathbf{Z}|\mathbf{Y}, \theta^{(t)})$ ，并构造 $Q(\theta|\theta^{(t)})$ 。
- M 步：求 $\theta^{(t+1)} = \arg\max_\theta Q(\theta|\theta^{(t)})$ 。
- 检查收敛：若满足停止条件，则终止。

常见的停止准则有两种：

参数变化量很小： $|\theta^{(t+1)} - \theta^{(t)}| < \varepsilon$ （例如 $\varepsilon = 10^{-6}$ ）。
对数似然增量很小： $\ell(\theta^{(t+1)}) - \ell(\theta^{(t)}) < \delta$ （例如 $\delta = 10^{-6}$ ）。

在实际编码中，通常监控参数变化或对数似然的相对增量，因为计算 $\ell(\theta)$ 本身在 $\log\sum$ 模型中可能有一定开销（虽然比直接优化它容易，但每步检查仍可只在必要时计算）。

输入: 观测序列 y = [y_1, y_2, ..., y_n], 停机阈值 epsilon
输出: 参数 pi, p, q 的估计

1. 初始化: 
    pi = 0.4  (或随机)
    p  = 0.6
    q  = 0.3
    prev_pi, prev_p, prev_q = pi, p, q

2. 循环:
    // E-step
    for i = 1 to n:
        计算 gamma_i:
        num = pi * (p^y_i) * ((1-p)^(1-y_i))
        den = num + (1-pi) * (q^y_i) * ((1-q)^(1-y_i))
        gamma[i] = num / den
    
    // M-step
    sum_gamma = sum(gamma)
    pi_new  = sum_gamma / n
    p_new   = sum(gamma[i] * y_i for i) / sum_gamma
    q_new   = sum((1-gamma[i]) * y_i for i) / (n - sum_gamma)

    // 收敛判断
    if max(|pi_new - pi|, |p_new - p|, |q_new - q|) < epsilon:
        break

    // 更新参数
    pi, p, q = pi_new, p_new, q_new

3. 返回 pi, p, q