EM 算法与混合模型

小象实战讲义 · 自然语言处理

在自然语言处理中，我们常常会遇到数据中存在“隐藏”信息的情况。例如，在词性标注中，我们能看到单词序列，但看不到每个单词对应的词性标签（隐变量）。如何从“不完整”的观测数据中，估计出模型的参数呢？期望最大化（EM）算法为此类问题提供了一套优雅的迭代求解框架。本节我们将通过经典的“三硬币问题”深入理解EM算法的核心思想，并掌握其在混合模型中的应用，为后续学习隐马尔可夫模型（HMM）等更复杂的隐变量模型打下坚实基础。

💡 核心导读

• 问题引入：通过“三硬币问题”理解什么是混合模型，以及为什么在存在隐变量时，直接极大似然估计会失效。
• EM算法思想：掌握EM算法的两步迭代核心：E步（期望步）计算隐变量的后验概率，M步（最大化步）基于“完整数据”更新模型参数。
• 算法推导：从概率图模型视角，形式化推导伯努利混合模型的EM更新公式，理解其数学本质。
• 收敛性理解：简要了解EM算法通过优化证据下界（ELBO）来保证似然函数单调递增的原理。
• 代码实践：通过Python实现三硬币问题的EM算法，直观观察参数如何迭代收敛。

混合模型与隐变量

想象一个场景：你面前有两个外观相同的袋子（A和B），袋A中装有大量红色硬币（正面概率为p），袋B中装有大量绿色硬币（正面概率为q）。现在，实验者先随机选择一个袋子（选择袋A的概率为π），然后从选中的袋子中取出一枚硬币，掷出并记录结果（正面H或反面T），再将硬币放回原袋。重复此过程N次，你只看到了一系列掷硬币的结果（如 H, T, H, T, T, H...），但不知道每次掷的是哪个袋子的硬币。

这就是一个典型的混合模型。我们观测到的数据（掷币结果）是由两个不同的成分（红色硬币和绿色硬币）以一定比例混合生成的。其中，每次实验选择哪个袋子，就是一个隐变量（Latent Variable），我们无法直接观测到它。

为什么直接估计参数困难？ 如果我知道每次掷的是红色还是绿色硬币，问题就很简单：

• 参数π（选择红色硬币的比例） = （红色硬币出现的次数）/ N
• 参数p（红色硬币正面概率） = （红色硬币出现正面的次数）/ （红色硬币出现的总次数）
• 参数q（绿色硬币正面概率）同理。

但问题在于，我们不知道隐变量的取值。EM算法的核心思想就是：既然不知道隐变量，我们就先“猜”一个（E步），然后用这个“猜测”去更新参数（M步），如此反复迭代，让猜测越来越准，参数也越来越接近真实值。

EM算法：期望最大化

EM算法是一个迭代算法，每次迭代包含两个步骤：

1. E步（Expectation Step）：基于当前参数估计值，计算隐变量的后验概率（即给定观测数据，隐变量取各个值的概率）。
2. M步（Maximization Step）：将E步计算得到的后验概率视为“权重”，重新估计模型参数，使得在“完整数据”（观测数据+隐变量）下的期望似然函数最大化。

我们形式化地定义问题。设观测变量为 $x$（例如，$x=1$ 代表正面，$x=0$ 代表反面），隐变量为 $z$（例如，$z=1$ 代表来自红色硬币，$z=0$ 代表来自绿色硬币）。模型参数为 $\theta = (\pi, p, q)$。

观测数据 $x$ 的（边际）概率分布是一个混合分布： $$P(x | \theta) = \sum_{z} P(x, z | \theta) = \sum_{z} P(z | \theta) P(x | z, \theta)$$ 具体到我们的伯努利混合模型： $$P(x | \theta) = \pi \cdot p^{x}(1-p)^{1-x} + (1-\pi) \cdot q^{x}(1-q)^{1-x}$$ 这个公式清晰地展示了混合：它是两个伯努利分布（一个参数为p，一个参数为q）的加权和，权重分别是π和(1-π)。

对于N个独立同分布的观测数据 $D = {x_1, x_2, ..., x_N}$，其对数似然函数为： $$\mathcal{L}(\theta; D) = \log P(D | \theta) = \sum_{i=1}^{N} \log \left[ \pi p^{x_i}(1-p)^{1-x_i} + (1-\pi) q^{x_i}(1-q)^{1-x_i} \right]$$ 由于对数内部是求和，直接对 $\theta$ 求导令其为零来最大化 $\mathcal{L}$ 非常困难（没有闭式解）。EM算法巧妙地绕开了这个难题。

E步：计算隐变量的后验概率

给定当前参数估计 $\theta^{(t)} = (\pi^{(t)}, p^{(t)}, q^{(t)})$，对于第 $i$ 个观测 $x_i$，我们计算隐变量 $z_i=1$（即来自红色硬币）的后验概率。根据贝叶斯公式： $$\mu_i^{(t+1)} = P(z_i=1 | x_i, \theta^{(t)}) = \frac{P(z_i=1 | \theta^{(t)}) P(x_i | z_i=1, \theta^{(t)})}{P(x_i | \theta^{(t)})}$$ 代入具体分布： $$\mu_i^{(t+1)} = \frac{\pi^{(t)} \cdot (p^{(t)})^{x_i}(1-p^{(t)})^{1-x_i}}{\pi^{(t)} \cdot (p^{(t)})^{x_i}(1-p^{(t)})^{1-x_i} + (1-\pi^{(t)}) \cdot (q^{(t)})^{x_i}(1-q^{(t)})^{1-x_i}}$$ $\mu_i$ 可以理解为，在现有参数下，观测到 $x_i$ 这个结果时，它“应该”有多大比例归因于红色硬币。同理，$1 - \mu_i$ 就是归因于绿色硬币的比例。

M步：重新估计参数

现在，我们有了每个数据点属于红色硬币的“软分配”权重 $\mu_i$。我们可以利用这些权重，像处理“完整数据”一样来更新参数：

• 更新π：π是红色硬币被选中的先验概率。我们用所有数据点属于红色硬币的权重之和来估计： $$\pi^{(t+1)} = \frac{1}{N} \sum_{i=1}^{N} \mu_i^{(t+1)}$$
• 更新p：p是红色硬币出正面的概率。我们用所有“属于”红色硬币的数据点中，正面所占的加权比例来估计。注意，$x_i$ 为1时才是正面。 $$p^{(t+1)} = \frac{\sum_{i=1}^{N} \mu_i^{(t+1)} \cdot x_i}{\sum_{i=1}^{N} \mu_i^{(t+1)}}$$
• 更新q：q是绿色硬币出正面的概率。我们用所有“属于”绿色硬币的数据点（权重为 $1-\mu_i$）中，正面所占的加权比例来估计。 $$q^{(t+1)} = \frac{\sum_{i=1}^{N} (1 - \mu_i^{(t+1)}) \cdot x_i}{\sum_{i=1}^{N} (1 - \mu_i^{(t+1)})}$$

E步和M步不断迭代，直到参数的变化小于某个阈值，或似然函数不再显著增加，算法收敛。

算法实现与演示

下面我们用Python实现三硬币问题的EM算法。为了清晰和自包含，我们将所有步骤写在一个代码块内。

import numpy as np

# 设置随机种子以保证结果可复现
np.random.seed(42)

# 1. 生成模拟数据（真实参数：π=0.6, p=0.7, q=0.3）
true_pi = 0.6
true_p = 0.7
true_q = 0.3
N = 1000  # 数据点数量

# 生成隐变量z（0:绿币，1:红币）
z_true = np.random.binomial(1, true_pi, N)
# 根据z生成观测数据x
x = np.zeros(N, dtype=int)
for i in range(N):
    if z_true[i] == 1:  # 红币
        x[i] = np.random.binomial(1, true_p)
    else:  # 绿币
        x[i] = np.random.binomial(1, true_q)

print(f"生成 {N} 个观测数据。")
print(f"正面(H)出现次数: {np.sum(x)}， 反面(T)出现次数: {N - np.sum(x)}")
print(f"真实参数: π={true_pi:.2f}, p={true_p:.2f}, q={true_q:.2f}")

生成 1000 个观测数据。
正面(H)出现次数: 556， 反面(T)出现次数: 444
真实参数: π=0.60, p=0.70, q=0.30

现在，我们假装不知道真实参数和隐变量，仅凭观测数据 x 来运行EM算法。

import numpy as np

# 重新构建上文用到的微型数据（本节各代码块相互独立，可单独运行）
np.random.seed(42)
true_pi = 0.6
true_p = 0.7
true_q = 0.3
N = 1000  # 数据点数量
z_true = np.random.binomial(1, true_pi, N)
x = np.zeros(N, dtype=int)
for i in range(N):
    if z_true[i] == 1:  # 红币
        x[i] = np.random.binomial(1, true_p)
    else:  # 绿币
        x[i] = np.random.binomial(1, true_q)

# 2. EM 算法实现（固定初始化，保证结果可复现）
def em_bernoulli_mixture(x, pi0=0.5, p0=0.7, q0=0.3, max_iter=50, tol=1e-8):
    """
    伯努利混合模型的EM算法。
    参数:
        x: 一维数组，观测数据（0或1）
        pi0, p0, q0: 固定初始化参数
        max_iter: 最大迭代次数
        tol: 收敛阈值（参数变化小于此值则停止）
    返回:
        pi_est, p_est, q_est: 估计的参数
        history: 记录每次迭代的参数值
    """
    N = len(x)
    pi_est, p_est, q_est = pi0, p0, q0
    history = [(pi_est, p_est, q_est)]

    for it in range(max_iter):
        # --- E步：后验概率 mu_i = P(z_i=1 | x_i, theta) ---
        prob_red = pi_est * (p_est ** x) * ((1 - p_est) ** (1 - x))
        prob_green = (1 - pi_est) * (q_est ** x) * ((1 - q_est) ** (1 - x))
        mu = prob_red / (prob_red + prob_green)  # 后验概率向量，长度N

        # --- M步：更新参数 ---
        pi_new = np.mean(mu)
        sum_mu = np.sum(mu)
        p_new = np.sum(mu * x) / sum_mu if sum_mu > 0 else p_est
        sum_one_minus_mu = np.sum(1 - mu)
        q_new = np.sum((1 - mu) * x) / sum_one_minus_mu if sum_one_minus_mu > 0 else q_est

        history.append((pi_new, p_new, q_new))
        delta = np.abs(pi_new - pi_est) + np.abs(p_new - p_est) + np.abs(q_new - q_est)
        pi_est, p_est, q_est = pi_new, p_new, q_new
        if delta < tol:
            print(f"算法在第 {it+1} 次迭代收敛。")
            break
    else:
        print(f"达到最大迭代次数 {max_iter}。")

    return pi_est, p_est, q_est, history

# 运行EM算法
pi_final, p_final, q_final, history = em_bernoulli_mixture(x)
print(f"\nEM算法估计结果:")
print(f"π_est = {pi_final:.4f}")
print(f"p_est = {p_final:.4f}")
print(f"q_est = {q_final:.4f}")
print(f"\n与真实值的绝对误差:")
print(f"|π - π_est| = {abs(true_pi - pi_final):.4f}")
print(f"|p - p_est| = {abs(true_p - p_final):.4f}")
print(f"|q - q_est| = {abs(true_q - q_final):.4f}")

# 3. 查看迭代过程
print("\n迭代历史（前6次）：")
print("迭代 |     π    |     p    |     q    ")
print("-" * 40)
for i, (pi_i, p_i, q_i) in enumerate(history[:7]):
    print(f"{i:4d} | {pi_i:.6f} | {p_i:.6f} | {q_i:.6f}")

算法在第 2 次迭代收敛。

EM算法估计结果:
π_est = 0.5224
p_est = 0.7450
q_est = 0.3492

与真实值的绝对误差:
|π - π_est| = 0.0776
|p - p_est| = 0.0450
|q - q_est| = 0.0492

迭代历史（前6次）：
迭代 |     π    |     p    |     q
----------------------------------------
   0 | 0.500000 | 0.700000 | 0.300000
   1 | 0.522400 | 0.745023 | 0.349246
   2 | 0.522400 | 0.745023 | 0.349246

有趣的是，在这个例子中，算法从给定的初始值出发，一轮 M 步后参数就不再变化（这组数据上 E 步责任度在下一轮几乎不变），之后便稳定下来；这是当前示例和初始化的结果，并不代表一般情形下 EM 一轮即收敛。

但这里必须先泼一盆冷水：这三个参数在这个数据集上根本不可识别（unidentifiable）。每个数据点只掷了一次硬币、只产生一个 0/1，于是我们能观测到的唯一统计量是样本均值：

$$r = \Pr(x=1) = \pi p + (1-\pi) q$$

一个方程、三个未知数——解是一整条曲线，而不是一个点。样本均值是 0.556，所以任何满足 πp + (1−π)q ≈ 0.556 的三元组都是同等好的解，似然完全相同。举两个例子：(π=0.5, p=0.90, q=0.212) 与 (π=0.3, p=0.80, q=0.4514) 都给出 r≈0.556，与程序输出的 (0.5224, 0.7450, 0.3492) 一样好。

因此上面那张”与真实值的绝对误差”表没有诊断价值：它比较的是”算法落在曲线上的哪一点”，而不是”有没有估准”。EM 在这种退化设定下只会停在初始化附近的那一点，并沿曲线做极小幅度的漂移，看起来像”收敛到真值附近”，其实是数据本身没有提供足够信息。

怎么才能让它可识别？把实验设计改掉，而不是换初始化：让每个数据点连续掷多次（比如每轮掷 10 次、只记录正面次数），此时观测到的是二项分布而不是单次伯努利，πp + (1−π)q 这条约束就被”多枚硬币的组合分布”替换掉，三个参数才可能被唯一确定。这也是经典 EM 硬币示例（Do & Batzoglou, 2008）里”每轮连掷 5 次”的原因。

EM算法的收敛性

为什么EM算法能保证收敛？其背后有一个优美的数学解释。EM算法可以看作是在最大化证据下界（Evidence Lower BOund, ELBO）。

我们的目标是最大化观测数据的对数似然 $\log P(X | \theta)$。对于任意一个关于隐变量 $Z$ 的分布 $Q(Z)$，根据詹森不等式（Jensen‘s inequality），有以下不等式成立： $$\log P(X | \theta) = \log \sum_{Z} P(X, Z | \theta) \ge \sum_{Z} Q(Z) \log \frac{P(X, Z | \theta)}{Q(Z)} \triangleq \mathcal{L}(Q, \theta)$$ $\mathcal{L}(Q, \theta)$ 就是对数似然的一个下界（ELBO）。EM算法的每一步都在提升这个下界：

• E步：固定 $\theta$，找到使下界 $\mathcal{L}(Q, \theta)$ 最大的 $Q$。可以证明，最优的 $Q$ 恰好是隐变量的后验分布 $Q(Z) = P(Z | X, \theta)$。
• M步：固定 $Q$，找到使下界 $\mathcal{L}(Q, \theta)$ 最大的 $\theta$。这等价于最大化期望完整数据对数似然 $\mathbb{E}_{Q}[\log P(X, Z | \theta)]$。

由于每一步都提升了ELBO，而ELBO是对数似然的下界，因此对数似然本身在每次迭代中也是非递减的：$\log P(X | \theta^{(t+1)}) \ge \log P(X | \theta^{(t)})$。这个单调性只保证观测对数似然不再下降，并不直接保证参数序列收敛，也不保证极限点一定是局部极大值。

需要注意的是，EM 的极限通常是似然的一个驻点（stationary point）；在可识别且目标函数性质良好的情形下，它常表现为局部极大值，但理论上不排除鞍点，最终落点仍依赖初始参数的选择。实践中，常采用多次随机初始化并选择似然最高的结果。

📝 动手练一练

1. 概念辨析：在伯努利混合模型的M步中，更新公式 $p^{(new)} = \frac{\sum \mu_i x_i}{\sum \mu_i}$。如果某个数据点 $x_j = 0$（反面），它对分子 $\sum \mu_i x_i$ 的贡献是多少？这对更新p有影响吗？请解释原因。

2. 代码实验：修改上面EM算法实现中的参数初始化部分，尝试几组不同的初始值（例如 pi=0.1, p=0.9, q=0.1 或 pi=0.9, p=0.1, q=0.9）。观察算法是否总能收敛到相同的解？收敛速度有差异吗？

本章小结

本节我们深入探讨了期望最大化（EM）算法，它是处理含有隐变量模型参数估计的基石性方法。

• 核心思想：通过迭代执行E步（计算隐变量后验期望）和M步（最大化期望完整数据似然），来间接优化难以直接求解的观测数据似然函数。
• 算法流程：我们以伯努利混合模型（三硬币问题）为例，一步步推导了E步和M步的具体计算公式，并实现了完整的算法。
• 收敛性：E 步把下界抬到”等于当前对数似然”的紧处，M 步再抬下界，因此对数似然逐步非减。但”单调非减”只保证收敛到驻点（通常是一个局部极大），并不保证收敛到全局最优，也不意味着每步都真的在改善 ELBO 的最优值——初始值仍然决定落到哪个驻点，多初始化取似然最大的做法正是为此。
• 应用广泛：EM算法是许多重要模型的基础，包括接下来要学习的高斯混合模型（GMM）、隐马尔可夫模型（HMM） 的参数学习（Baum-Welch算法），以及概率主题模型（如PLSA）等。

📋 行动清单

• 理解混合模型中“隐变量”和“混合系数”的概念，能用自己的话解释三硬币问题。
• 掌握EM算法E步和M步的目标，并能写出伯努利混合模型对应的更新公式。
• 运行提供的代码，尝试修改数据生成的真实参数，观察EM算法估计的准确性。
• （选做）思考：如果袋子不是两个而是K个，即K个伯努利分布的混合，EM算法的公式应该如何推广？

—— 小象教研组