EM 算法：期望与最大化的迭代

小象实战讲义 · 机器学习算法精讲

上一节的三硬币问题停在了最关键的一步：A 决定用 B 还是 C，我们只看得到最终的正反面——$\pi, p, q$ 三个参数怎么估？这一节，讲师「不休息、把问题讲完」：先写出整个系统的似然，看清它是一个混合模型；再用一个「如果知道是哪枚硬币」的思想实验找到估计的钥匙；最后把钥匙变成可迭代的算法——EM（期望最大化）。

💡 核心导读

• 定位：三硬币问题是一个有参数模型——$\theta = (\pi, p, q)$ 三个参数决定整个系统的性能，目标是极大似然。
• 似然：引入隐含变量 $Z$（$z=1$ 掷 B、$z=0$ 掷 C）边缘化后，单次观测的概率是 $P(y|\theta) = \pi p^{y}(1-p)^{1-y} + (1-\pi)q^{y}(1-q)^{1-y}$。
• 混合模型：$\pi$ 与 $1-\pi$ 是两枚硬币的权重——课件笔记称之为伯努利混合模型（Bernoullian Mixture）。
• EM 两步：E 步用当前参数算每个结果由 B 掷出来的概率 $\mu_j$（软归属）；M 步按 $\mu_j$ 的比例「数红绿」，更新 $\pi, p, q$。
• 迭代：随机给定 $\theta_0$，$\mu \to \theta \to \mu \to \theta$ 交替进行；$\pi, p, q$ 不再变化时，就得到对整个系统的估计。

先把问题摆正：一个有参数模型

开场两个问答，把问题钉住。

第一，它是一个什么模型？——一个有参数模型（Parametric Model）。这个系统里有 $\pi, p, q$ 一共三个参数，这三个参数决定了整个系统的性能，所以把它们合起来记作 $\theta$：$\theta = (\pi, p, q)$。

第二，我们想求什么？——求给定了这些数据之后的 likelihood：说白了，就是想求「给定 $\theta$ 时，$y$ 的概率是多大」。我们观察到了很多 $y$，根据这些 $y$ 来估计 $\theta$；最终想要的计算模型，是预测下一个 $y$ 是正面还是反面的概率——那需要算出最好的参数 $\theta^*$（根据以前给定的数据训练出来）。

如果对这件事还有疑惑（讲师照顾基础稍弱的同学），回想正态分布的例子：给定一堆数据，假设它服从一个正态分布 $N(\mu, \sigma)$——训练数据就是想知道 $\mu$ 和 $\sigma$ 是什么；做极大似然，训练出最好的 $\mu^*$ 和 $\sigma^*$。知道了整个分布长什么样，下次就可以再采样：下一个数据大概就在这个位置，再下一个又能采出来一个。对未来数据的估计，实际上就是算出那个最好的参数。三硬币这里是一模一样的事情——只不过分布更复杂，是 $\pi, p, q$ 混合构成的一个分布，但最终仍然想求一个最优的参数，来预测下一个 $y$。

写出似然：隐含变量登场

那怎么写这个「复杂分布」？用刚学过的工具——引入隐含变量。

$Z$ 是控制 $Y$ 的变量：$Y$ 是硬币的正面反面，而「也许有一个变量控制着是 B 还是 C 掷出这个正反面」——很显然，$Z$ 就是我们看不到的硬币 A 的结果。$Z$ 也是二值的：

• $Z = 1$：这次掷的是硬币 B；
• $Z = 0$：这次掷的是硬币 C。

把 $Z$ 的所有可能累加进公式（上一节的边缘化），再给每一项都带上条件 $\theta$：

$$P(y|\theta) = \sum_{z \in \lbrace 0,1\rbrace} P(y, z|\theta) = \sum_{z \in \lbrace 0,1\rbrace} P(y|z,\theta)\,P(z|\theta)$$

（把 $\theta$ 都去掉，它就是熟悉的联合概率与条件概率的关系：$P(Y,Z) = P(Y|Z)P(Z)$——只是现在每一项都多了一个条件。）

拆成两半算。

第一半：$P(z|\theta)$。由谁决定？由硬币 A。A 是正面时掷 B——所以 $P(z=1|\theta) = \pi$，$P(z=0|\theta) = 1 - \pi$。「可被计算」。

第二半：$P(y|z,\theta)$。看 $z=1$ 的情形：给定掷的是硬币 B，得到正面（$y=1$）的概率多大？这就是关于硬币 B 的伯努利分布——参数是 $p$，所以是 $p^{y}(1-p)^{1-y}$。同理 $z=0$ 时是关于硬币 C 的伯努利分布，参数 $q$：$q^{y}(1-q)^{1-y}$。

合并。把两种情况加在一起，就得到单次观测的分布——$\pi$ 乘上 B 的伯努利项，加 $1-\pi$ 乘上 C 的伯努利项：

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

整个系统里那个复杂的关系，就用 $\pi, p, q$ 三个参数全部解释出来了。

别忘了我们掷了 $N$ 次：$y_1, y_2, \dots, y_N$。每一次的分布都由上面这串算出，所以似然是它们的连乘；照例取 log 变成连加：

$$\ell(\theta) = \sum_{i=1}^{N} \log\left[\pi p^{y_i}(1-p)^{1-y_i} + (1-\pi) q^{y_i}(1-q)^{1-y_i}\right]$$

问题来了：公式（系统）我们已经知道，但 $\pi, p, q$ 这三个参数是多少？并不知道。我们唯一拥有的数据就是 $y_i$——一大堆正面和反面。怎么估计？

这是一个混合模型

先看清 $P(y|\theta)$ 的结构。它的一部分是硬币 B 掷出正面的概率（一个伯努利分布），另一部分是硬币 C 的；而「有多大可能是 B」是 $\pi$、「有多大可能是 C」是 $1-\pi$——$\pi$ 和 $1-\pi$ 相当于两种情况的权重。这种「权重乘密度再相加」的结构，就叫混合模型（Mixture Model）。

为什么「$\pi$ 乘某物加 $(1-\pi)$ 乘另一物」就是混合？一个非常直观的想象（讲师原话场景）：有两个数，一个 $A$ 一个 $B$。$\pi$ 是 0 到 1 之间的概率——比如取 $\pi = 0.3$（讲师先说 0.5，马上自我纠正：「0.5 不可能，0.3 更好一点」），那么选 $A$ 的概率是 0.3、选 $B$ 的概率是 0.7。求它们之间的「均值」时，$0.3 \times A + 0.7 \times B$——相当于选了它们之间的一个中间点，更靠近 B 一点。这就是一种混合。两类可以推广到多类：$\alpha_1 P_1 + \alpha_2 P_2 + \dots + \alpha_n P_n$，只要所有 $\alpha_i$ 加起来等于 1，就代表 $n$ 个事件各按权重 $\alpha_i$ 的多模态混合。

所以三硬币的问题，很容易看成「B 和 C 这两种伯努利分布的混合」。课件笔记把它讲得更精确：这是一个 Bernoullian Mixture model（伯努利混合模型），$\pi$ 可以看作 B 和 C 的混合权重；而盒子里的东西可以换成骰子、换成任何服从某种分布的事件——把伯努利换成高斯，就是第 13 章的高斯混合模型。

如果知道是哪枚硬币：红绿标色

估计的钥匙，藏在一个思想实验里。

我们看到这些 $y$，比如 $0, 1, 1, 1, 1, 0, 0, 1, 1, 0$——一堆 0、1 的组合，但没法估计 $\pi, p, q$，关键就在于不知道每个结果是哪枚硬币掷出来的。可是一旦假设我们知道——用课件笔记的画法：把硬币 B 掷出来的结果标成红色、硬币 C 的标成绿色——这件事立刻就可以做了：

• $\pi$ 的估计：$\pi$ 代表 A 掷出正面（从而选 B）的概率，所以就是红色占全部数据的比例——红色个数比上红色加绿色；
• $p$ 的估计：红色默认全是 B 掷的，$p$ 是 B 掷出正面的概率——所以是红色的一，占所有红色的比例；
• $q$ 的估计：同理——绿色的一，占所有绿色的比例。

数数而已。所以讲师在这里点破：隐含变量的关键，就在于我们并不知道它是哪个硬币掷出来的；但一旦知道，这些事都解决了。我们要做的，就是对每一次结果「它由 B 还是 C 掷出来」产生一个估计——只要得到它，剩下的事情都可以做。

用概率标色：μ 与 E 步

可是开头我们并没有红绿——那就近似地估计：给一个新的参数 $\mu_j$，代表「第 $j$ 次的结果，有多大比例是由 B 掷出来的」。既然 $P(y|\theta)$ 是个混合模型，把 B 的那一半劈出来放在分子、全部可能性放在分母，就是 $\mu_j$（贝叶斯公式的又一次现身）：

$$\mu_j = \frac{\pi \cdot p^{y_j}(1-p)^{1-y_j}}{\pi \cdot p^{y_j}(1-p)^{1-y_j} + (1-\pi) \cdot q^{y_j}(1-q)^{1-y_j}}$$

上面的部分是 B 发生的概率，底下是 B 加上 C 发生的概率——所以在给定当前 $\pi, p, q$ 时，它就是第 $j$ 次结果由 B 掷出来的概率。

注意它跟红绿标色的关系：我们不可能真的知道每个结果到底是红色还是绿色掷出来的——红色和绿色，是期望得到的结果。我们能做到的，是给每次结果一个「属于红色」的概率：它也许属于红色是 0.7、它是 0.2、它是 0.9、它是 0.4、0.8……「带着比例的红色」——把 0.7 个红色、0.2 个红色这样加在一起，就是所有的红色。用给定的这些 $y$，对红绿给出一个更好的估计，这个估计就是 $\mu$。

这一步，就是 EM 里的 E step（Expectation，求期望）：根据给定的参数，求每一个掷出来的结果有多大概率是由 B 掷出来的。

M 步：按比例数数

有了 $\mu_j$，就可以按红绿标色的方式「数数」了——只是现在每一个数都带着比例。

$\pi$ 的更新：$\pi$ 代表「掷了 $n$ 次里面，平均有多少次用了 B」。掷 100 次、$\pi = 0.3$，就代表约 30 次的结果是 B 掷出来的。所以把每一次「属于 B 的份额」$\mu_j$ 合在一起求均值：

$$\pi^{new} = \frac{1}{N}\sum_{j=1}^{N} \mu_j$$

$p$ 的更新：红色里有多少红色的一。红色的总数是 $\sum_j \mu_j$（把带着比例的红色加在一起）；「红色的一」——$y_j$ 只有 0 或 1：$y_j = 1$ 时，$\mu_j$ 这一份红色算进去；$y_j = 0$ 时，乘上 0 直接归零。所以分子是 $\sum_j \mu_j y_j$：

$$p^{new} = \frac{\sum_{j=1}^{N} \mu_j y_j}{\sum_{j=1}^{N} \mu_j}$$

$q$ 的更新：同理，绿色是 $1 - \mu_j$，绿色的一是 $(1-\mu_j) y_j$：

$$q^{new} = \frac{\sum_{j=1}^{N} (1-\mu_j)\, y_j}{\sum_{j=1}^{N} (1-\mu_j)}$$

这一步，就是 M step（Maximization，最大化）：根据 $\mu$ 的估计，对 $\pi, p, q$ 做最大化的估计。

完整算法：互相迭代

把两步接起来，算法的全貌是：

1. 随机给定初始值 $\theta_0 = (\pi_0, p_0, q_0)$——「随便比如给 0.3、0.4、0.7 都可以」；
2. 用第 $t$ 时刻的参数 $\pi^t, p^t, q^t$ 代入 $\mu$ 的公式，算出第 $t+1$ 时刻的 $\mu^{t+1}$（E 步）；
3. 用 $\mu^{t+1}$ 数数，更新出 $\pi^{t+1}, p^{t+1}, q^{t+1}$（M 步）；
4. 重复 2–3——「就是一个互相迭代的一个算法」。当 $\pi, p, q$ 不再发生任何变化时，我们就得到了对整个系统的估计。

（板书花絮：讲师写到迭代轮次时，在上角标与下角标之间纠结了一阵——「哎呀，这个 notation 有点乱了……不然还是这么写吧」——最终把迭代轮次放在上角标：$\pi^t$ 指第 $t$ 轮迭代的 $\pi$，别和幂次混淆。）

这就是 EM algorithm，EM 算法——Expectation Maximization，期望最大化：分两个部分，E step 求期望（每个结果由 B 掷出来的概率），M step 求最大化（给定 $\mu$，更新 $\pi, p, q$），不断迭代，最后得到对 $\pi, p, q$ 最好的估计。「这就是我给大家所讲述的 EM 算法」——本节视频到此收尾。

亲手跑一遍：三硬币的 EM

下面把整个算法落成代码（纯标准库）。先用已知的「真实参数」模拟生成一串观测——然后把这串观测交给 EM，看它能不能把参数恢复出来：

import random

def em_coin(observations, init_params, max_iters=20, tol=1e-6):
    """三硬币 EM: E 步算 mu_j(第 j 次由 B 掷出的概率), M 步按比例数红绿"""
    pi, p, q = init_params
    N = len(observations)
    history = [(0, pi, p, q)]
    for it in range(1, max_iters + 1):
        # E-step: 逐个计算 mu_j
        mu = []
        for y in observations:
            prob_b = pi * (p ** y) * ((1 - p) ** (1 - y))    # B 的那一半
            prob_c = (1 - pi) * (q ** y) * ((1 - q) ** (1 - y))  # C 的那一半
            mu.append(prob_b / (prob_b + prob_c))           # B / (B + C)
        # M-step: 按比例数数
        sum_mu = sum(mu)                                      # 红色总数
        pi_new = sum_mu / N                                   # 红色占比
        p_new = sum(m * y for m, y in zip(mu, observations)) / sum_mu   # 红色的一 / 红色
        q_new = sum((1 - m) * y for m, y in zip(mu, observations)) / (N - sum_mu)
        if max(abs(pi_new - pi), abs(p_new - p), abs(q_new - q)) < tol:
            pi, p, q = pi_new, p_new, q_new
            history.append((it, pi, p, q))
            break
        pi, p, q = pi_new, p_new, q_new
        history.append((it, pi, p, q))
    return (pi, p, q), history

# 模拟生成观测: 真实参数 pi=0.6, p=0.7, q=0.3
random.seed(42)
true_pi, true_p, true_q = 0.6, 0.7, 0.3
obs = []
for _ in range(10):
    if random.random() < true_pi:            # 掷 A: 正面 -> 用 B
        obs.append(1 if random.random() < true_p else 0)
    else:                                    # A 反面 -> 用 C
        obs.append(1 if random.random() < true_q else 0)

print("观测序列:", obs)
print("序列中1的比例:", sum(obs) / len(obs))

# EM: 初始参数故意偏离真实值
params, history = em_coin(obs, init_params=(0.4, 0.6, 0.3))
print(f"\n真实参数: pi={true_pi:.4f}, p={true_p:.4f}, q={true_q:.4f}")
print(f"EM估计:   pi={params[0]:.4f}, p={params[1]:.4f}, q={params[2]:.4f}")
print("\n迭代过程:")
print("  迭代      pi       p       q")
for it, pi_v, p_v, q_v in history:
    print(f"{it:5d}   {pi_v:.4f}  {p_v:.4f}  {q_v:.4f}")

观测序列: [1, 1, 0, 1, 1, 1, 1, 0, 1, 1]
序列中1的比例: 0.8

真实参数: pi=0.6000, p=0.7000, q=0.3000
EM估计:   pi=0.5123, p=0.8923, q=0.7030

迭代过程:
  迭代      pi       p       q
    0   0.4000  0.6000  0.3000
    1   0.5123  0.8923  0.7030
    2   0.5123  0.8923  0.7030

读一读这个结果：观测里 1 占了八成，EM 的解读是——大约一半的试验用了 B，而两枚硬币都被判成偏正面（$p = 0.89$、$q = 0.70$）。只有 10 个观测点，EM 找到的是「能解释这串数据」的一组参数，它跟生成数据的真实参数不必相同；换成更长的序列，估计才有机会贴近真值。另外注意：这个初始值下第 1 轮迭代就到达了不动点（第 2 轮参数完全不变）——小数据上 EM 收敛得很快。

📝 动手练一练

1. 极端情形：假设初始参数恰好让每个 $\mu_j$ 算出来都等于 0 或 1，那么 E 步和 M 步会退化成什么样？对照上一节「红绿标色」的思想实验想一想：软归属（$\mu_j$ 取 0 到 1 之间的分数）的价值在哪里？

2. 换个序列跑一跑：把观测换成偏反面的序列 $[1,0,0,1,0,0,1,0,0,1]$（1 只占四成），初始参数不变。EM 会给出什么样的「解读」？再验证一个漂亮的恒等式：模型算出的「1 的比例」是否恰好等于序列里 1 的频率？

# 自包含精简版: 与正文 em_coin 同一套 E/M 公式, 跑到参数不动为止
def em_coin_flat(observations, init_params, max_iters=50, tol=1e-6):
    pi, p, q = init_params
    for _ in range(max_iters):
        mu = [pi * (p ** y) * ((1 - p) ** (1 - y))
              / (pi * (p ** y) * ((1 - p) ** (1 - y))
                 + (1 - pi) * (q ** y) * ((1 - q) ** (1 - y)))
              for y in observations]
        pi_n = sum(mu) / len(mu)
        p_n = sum(m * y for m, y in zip(mu, observations)) / sum(mu)
        q_n = sum((1 - m) * y for m, y in zip(mu, observations)) / (len(mu) - sum(mu))
        if max(abs(pi_n - pi), abs(p_n - p), abs(q_n - q)) < tol:
            return pi_n, p_n, q_n
        pi, p, q = pi_n, p_n, q_n
    return pi, p, q

obs_new = [1, 0, 0, 1, 0, 0, 1, 0, 0, 1]
params = em_coin_flat(obs_new, init_params=(0.4, 0.6, 0.3))
print(f"EM估计: pi={params[0]:.4f}, p={params[1]:.4f}, q={params[2]:.4f}")

EM估计: pi=0.3941, p=0.5800, q=0.2829

本章小结

• 问题：三硬币是一个有参数模型，$\theta = (\pi, p, q)$；目标是从观测序列做极大似然、预测下一个 $y$。
• 似然：边缘化隐含变量 $Z$（A 的结果）后，$P(y|\theta) = \pi p^{y}(1-p)^{1-y} + (1-\pi)q^{y}(1-q)^{1-y}$——一个以 $\pi$、$1-\pi$ 为权重的伯努利混合模型。
• 钥匙：如果知道每个结果来自哪枚硬币（红绿标色），参数估计就是三次数数——难点全在「不知道」。
• E 步：用当前参数算 $\mu_j$——第 $j$ 次结果由 B 掷出来的概率（B 的概率 / B+C 的概率），红色绿色的「期望版」。
• M 步：按 $\mu_j$ 的比例数红绿，更新 $\pi, p, q$ 三个分式。
• 迭代：随机给定 $\theta_0$，E、M 交替，$\pi, p, q$ 不再变化即得到估计。这套框架里，把伯努利换成别的分布，就是第 13 章的混合模型家族。

📋 行动清单

学完本节，你可以立即：

• 手推一遍 $\mu_j$：从贝叶斯公式出发，写出「B 的概率 / B+C 的概率」，确认分子分母的每一项来自混合模型的哪一半。
• 换数据做实验：修改正文代码的观测序列与初始参数，观察 EM 给出的「解读」怎么变——并检查 $\pi p + (1-\pi)q$ 是否总是等于序列的 1 频率。
• 红绿笔手算：取 5 个观测，手算一轮 E 步的 $\mu_j$ 与 M 步的 $\pi^{new}, p^{new}, q^{new}$——每个分式的分子分母都对应「哪种颜色的什么」。
• 预习下一节《琴生不等式》：EM 为什么「每次迭代都不亏」——那把数学钥匙，下一节配齐。

—— 小象教研组