隐马尔可夫链与 HMM 推导

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

在自然语言处理中，我们常常需要处理序列数据，例如词性标注、语音识别等任务。这些任务的核心是推断序列背后隐含的结构或状态。隐马尔可夫模型（Hidden Markov Model, HMM）是处理这类序列标注问题的经典概率图模型。本节我们将深入理解 HMM 的核心思想，并推导其核心算法——前向-后向算法，为后续学习更复杂的序列模型（如条件随机场、RNN/LSTM）打下坚实的理论基础。

💡 核心导读

• HMM 是什么：理解 HMM 作为“动态混合模型”的本质，它由隐含状态序列、观测序列、状态转移概率和发射概率构成。
• 三大核心问题：学习问题（估计参数）、解码问题（寻找最可能的状态序列）、评估问题（计算观测序列的概率）。
• 前向-后向算法：掌握如何高效计算观测序列的概率，以及如何为每个时刻的状态进行概率估计，这是 EM 算法在 HMM 中的具体实现（Baum-Welch 算法）的基础。
• 从 HMM 到现代模型：理解 HMM 的局限性，以及其核心思想（状态转移、序列建模）如何被后来的 RNN、LSTM 等神经网络模型所继承和发展。

HMM：一个动态的混合模型

隐马尔可夫模型可以看作是对混合模型（如高斯混合模型 GMM）的推广。在混合模型中，每个观测数据点独立地由一个隐含的“成分”（如高斯分布）生成。而在 HMM 中，这些隐含的“状态”不再是独立的，而是通过一个马尔可夫链连接起来，形成了一个沿着时间轴演变的动态过程。

一个 HMM 由以下五元组定义 $(\mathcal{Q}, \mathcal{V}, A, B, \pi)$：

• $\mathcal{Q} = {q_1, q_2, ..., q_N}$：所有可能的隐含状态集合，共 N 个状态。
• $\mathcal{V} = {v_1, v_2, ..., v_M}$：所有可能的观测值集合，共 M 个观测。
• $A = [a_{ij}]_{N \times N}$：状态转移概率矩阵。$a_{ij} = P(q_t = j | q_{t-1} = i)$，表示从状态 i 转移到状态 j 的概率。
• $B = [b_j(k)]_{N \times M}$：发射概率矩阵。$b_j(k) = P(o_t = v_k | q_t = j)$，表示在状态 j 下观测到值 $v_k$ 的概率。
• $\pi = [\pi_i]_{N \times 1}$：初始状态概率向量。$\pi_i = P(q_1 = i)$，表示序列起始时处于状态 i 的概率。

一个经典的例子是“天气-活动”模型：

• 隐含状态：晴天（Sunny, S）、雨天（Rainy, R）。
• 观测值：散步（Walk）、购物（Shop）、打扫（Clean）。
• 我们只能观察到一个人每天的活动（观测序列），但无法直接看到天气（状态序列）。HMM 的任务就是，给定一段活动记录，推断最可能的天气变化序列，或者估计天气变化的规律（参数 A, B）。

状态转移概率 A (天气变化规律)：
        S(t+1)   R(t+1)
S(t)   [ 0.8      0.2  ]
R(t)   [ 0.3      0.7  ]

发射概率 B (不同天气下的活动偏好)：
        Walk    Shop    Clean
S     [ 0.6     0.3     0.1  ]
R     [ 0.1     0.4     0.5  ]

初始概率 π: [P(S)=0.6, P(R)=0.4]

HMM 的图表示与联合概率

HMM 是一个典型的动态贝叶斯网络。其图结构清晰地展示了变量间的依赖关系：

时间轴: t=1        t=2        t=3       ...      t=T
隐含状态: q1 -----> q2 -----> q3 -----> ... ----> qT
            \        \        \                   \
             v        v        v                   v
观测变量:   o1        o2        o3       ...       oT

关键假设（齐次马尔可夫性 & 观测独立性）：

1. 齐次马尔可夫性：当前状态 $q_t$ 只依赖于前一个状态 $q_{t-1}$，与更早的状态和观测无关。即 $P(q_t | q_{1:t-1}, o_{1:t-1}) = P(q_t | q_{t-1})$。
2. 观测独立性：当前观测 $o_t$ 只依赖于当前状态 $q_t$，与其他时刻的状态和观测无关。即 $P(o_t | q_{1:t}, o_{1:t-1}) = P(o_t | q_t)$。

基于这两个假设，一个长度为 T 的观测序列 $O = (o_1, o_2, ..., o_T)$ 和状态序列 $Q = (q_1, q_2, ..., q_T)$ 的联合概率可以分解为：

$$P(O, Q) = P(q_1) \cdot P(o_1 | q_1) \cdot \prod_{t=2}^{T} \left[ P(q_t | q_{t-1}) \cdot P(o_t | q_t) \right]$$

用模型参数表示为：

$$P(O, Q) = \pi_{q_1} \cdot b_{q_1}(o_1) \cdot \prod_{t=2}^{T} \left[ a_{q_{t-1}q_t} \cdot b_{q_t}(o_t) \right]$$

这个公式是理解所有 HMM 算法的基础。它表明，生成一个序列的过程是：先根据 $\pi$ 选择一个初始状态，然后根据 $B$ 生成一个观测；接着根据 $A$ 转移到下一个状态，再根据 $B$ 生成观测，如此反复。

前向算法：高效计算 P(O|λ)

给定模型参数 $\lambda = (A, B, \pi)$，计算观测序列 $O$ 出现的概率 $P(O|\lambda)$，即评估问题。最直接的方法是穷举所有可能的状态序列 $Q$，计算其联合概率 $P(O, Q|\lambda)$ 并求和：

$$P(O|\lambda) = \sum_{Q} P(O, Q|\lambda) = \sum_{q_1, q_2, ..., q_T} \pi_{q_1} b_{q_1}(o_1) \prod_{t=2}^{T} a_{q_{t-1}q_t} b_{q_t}(o_t)$$

这种方法的计算复杂度是 $O(T N^T)$，是指数级的，不可行。前向算法通过动态规划，将复杂度降低到 $O(T N^2)$。

定义前向概率 $\alpha_t(i)$：在时刻 t，观测到前 t 个观测 $o_1, o_2, ..., o_t$，且此时状态为 $q_t = i$ 的联合概率。

$$\alpha_t(i) = P(o_1, o_2, ..., o_t, q_t = i | \lambda)$$

递推步骤：

1. 初始化（t=1）： $$\alpha_1(i) = \pi_i \cdot b_i(o_1), \quad i = 1, 2, ..., N$$ 这表示序列刚开始时，处于状态 i 并观测到 $o_1$ 的概率。
2. 递推（t=2,3,…,T）： $$\alpha_t(j) = \left[ \sum_{i=1}^{N} \alpha_{t-1}(i) \cdot a_{ij} \right] \cdot b_j(o_t), \quad j = 1, 2, ..., N$$ 要到达时刻 t 的状态 j，必须从时刻 t-1 的某个状态 i 转移过来。$\sum_{i} \alpha_{t-1}(i) a_{ij}$ 汇总了所有可能到达状态 j 的路径的概率，再乘以在状态 j 下观测到 $o_t$ 的概率 $b_j(o_t)$，就得到了 $\alpha_t(j)$。
3. 终止： $$P(O|\lambda) = \sum_{i=1}^{N} \alpha_T(i)$$ 在最后一个时刻 T，所有可能的状态下观测到整个序列 O 的概率之和，就是我们要的 $P(O|\lambda)$。

# 前向算法手算示例（伪代码逻辑，非真实运行）
import numpy as np

# 定义模型参数 (N=2, M=3)
pi = np.array([0.6, 0.4])  # [S, R]
A = np.array([[0.8, 0.2],  # S->S, S->R
              [0.3, 0.7]]) # R->S, R->R
# 观测索引: Walk:0, Shop:1, Clean:2
B = np.array([[0.6, 0.3, 0.1],  # S 下观测概率
              [0.1, 0.4, 0.5]]) # R 下观测概率

# 观测序列: Walk, Shop, Clean (索引: 0, 1, 2)
O = [0, 1, 2]
T = len(O)
N = A.shape[0]

# 初始化 alpha
alpha = np.zeros((T, N))
alpha[0, :] = pi * B[:, O[0]]
print(f"alpha_1: {alpha[0]}")

# 递推
for t in range(1, T):
    for j in range(N):
        alpha[t, j] = np.sum(alpha[t-1, :] * A[:, j]) * B[j, O[t]]
    print(f"alpha_{t+1}: {alpha[t]}")

# 终止
P_O = np.sum(alpha[T-1, :])
print(f"P(O|λ) = {P_O:.6f}")

alpha_1: [0.36 0.04]
alpha_2: [0.09 0.04]
alpha_3: [0.0084 0.023 ]
P(O|λ) = 0.031400

后向算法与状态概率估计

为了估计在给定整个观测序列 O 的条件下，某一时刻 t 处于状态 i 的概率 $P(q_t = i | O, \lambda)$，我们需要同时利用过去和未来的信息。这引出了后向算法。

定义后向概率 $\beta_t(i)$：在时刻 t 状态为 $q_t = i$ 的条件下，从 t+1 到 T 的观测序列 $o_{t+1}, o_{t+2}, ..., o_T$ 出现的概率。

$$\beta_t(i) = P(o_{t+1}, o_{t+2}, ..., o_T | q_t = i, \lambda)$$

递推步骤：

1. 初始化（t=T）： $$\beta_T(i) = 1, \quad i = 1, 2, ..., N$$ 定义在序列末尾，未来没有观测，概率为 1。
2. 递推（t=T-1, T-2, …, 1）： $$\beta_t(i) = \sum_{j=1}^{N} a_{ij} \cdot b_j(o_{t+1}) \cdot \beta_{t+1}(j), \quad i = 1, 2, ..., N$$ 从当前状态 i 出发，可以转移到下一个时刻的任意状态 j。转移概率是 $a_{ij}$，在状态 j 下观测到 $o_{t+1}$ 的概率是 $b_j(o_{t+1})$，而从状态 j 开始能生成剩余观测序列的概率是 $\beta_{t+1}(j)$。对所有可能的 j 求和，即得到 $\beta_t(i)$。

结合前向和后向概率，我们可以计算两个重要的概率：

1. 单个状态的概率（用于状态“软”标注）： 定义 $\gamma_t(i) = P(q_t = i | O, \lambda)$。根据贝叶斯公式和前向-后向概率： $$\gamma_t(i) = \frac{P(O, q_t=i | \lambda)}{P(O|\lambda)} = \frac{\alpha_t(i) \beta_t(i)}{\sum_{j=1}^{N} \alpha_t(j) \beta_t(j)}$$ 分子 $\alpha_t(i) \beta_t(i)$ 恰好等于 $P(O, q_t=i | \lambda)$，因为它结合了到时刻 t 为止的观测和状态 i（前向部分），以及从状态 i 出发生成剩余观测的概率（后向部分）。

2. 两个连续状态的概率（用于参数估计）： 定义 $\xi_t(i, j) = P(q_t = i, q_{t+1} = j | O, \lambda)$。它表示在给定整个观测序列的条件下，时刻 t 处于状态 i 且时刻 t+1 处于状态 j 的概率。 $$\xi_t(i, j) = \frac{\alpha_t(i) \cdot a_{ij} \cdot b_j(o_{t+1}) \cdot \beta_{t+1}(j)}{P(O|\lambda)}$$ 分子是联合概率 $P(O, q_t=i, q_{t+1}=j | \lambda)$ 的分解。

# —— 重新构建上文用到的微型数据（本节各代码块相互独立，可单独运行）
import numpy as np
pi = np.array([0.6, 0.4])  # [S, R]
A = np.array([[0.8, 0.2],  # S->S, S->R
              [0.3, 0.7]]) # R->S, R->R
B = np.array([[0.6, 0.3, 0.1],  # S 下观测概率
              [0.1, 0.4, 0.5]]) # R 下观测概率
O = [0, 1, 2]
T = len(O)
N = A.shape[0]
# 前向递推：alpha_t(i) = P(o_1..o_t, q_t = i)
alpha = np.zeros((T, N))
alpha[0, :] = pi * B[:, O[0]]
for t in range(1, T):
    alpha[t, :] = (alpha[t-1, :] @ A) * B[:, O[t]]
for t in range(T):
    print(f"alpha_{t+1}: {alpha[t]}")

# 后向算法：beta_t(i) = P(o_{t+1}..o_T | q_t = i)
# 初始化 beta
beta = np.zeros((T, N))
beta[T-1, :] = 1.0  # 注意 Python 索引从0开始，T-1对应时刻T
print(f"beta_{T}: {beta[T-1]}")

# 递推 beta
for t in range(T-2, -1, -1): # 从 T-2 到 0
    for i in range(N):
        beta[t, i] = np.sum(A[i, :] * B[:, O[t+1]] * beta[t+1, :])
    print(f"beta_{t+1}: {beta[t]}") # t是索引，对应时刻 t+1

# gamma_t(i) = alpha_t(i) * beta_t(i) / P(O|lambda)
P_O = np.sum(alpha[T-1, :])
gamma = np.zeros((T, N))
for t in range(T):
    gamma[t, :] = (alpha[t, :] * beta[t, :]) / P_O
    print(f"gamma_{t+1}: {gamma[t]} (状态概率)")

alpha_1: [0.36 0.04]
alpha_2: [0.09 0.04]
alpha_3: [0.0084 0.023 ]
beta_3: [1. 1.]
beta_2: [0.18 0.38]
beta_1: [0.0736 0.1226]
gamma_1: [0.84382166 0.15617834] (状态概率)
gamma_2: [0.51592357 0.48407643] (状态概率)
gamma_3: [0.26751592 0.73248408] (状态概率)

结果显示，给定观测序列 [Walk, Shop, Clean]，第一天处于晴天（状态 0）的概率为 0.8438，第二天 0.5159，第三天 0.2675；对应的雨天概率依次为 0.1562、0.4841、0.7325。最后一天雨天概率反而最高，是因为第三步观测是 Clean，雨天发射该观测的概率 0.5 远高于晴天的 0.1。三行概率各自相加均为 1，且 $\alpha_3$ 两项之和正是 $P(O|\lambda)=0.0314$，与后向递推的回代结果一致。

Baum-Welch 算法：无监督参数学习

当模型参数 $\lambda = (A, B, \pi)$ 未知时，我们需要从观测序列 O 中学习它们，即学习问题。Baum-Welch 算法是 EM 算法在 HMM 中的具体应用。

E步：在当前参数估计 $\lambda^{(old)}$ 下，利用前向-后向算法计算 $\gamma_t(i)$ 和 $\xi_t(i, j)$。

M步：利用 E 步计算出的期望统计量，重新估计参数 $\lambda^{(new)}$。

1. 初始状态概率： $$\hat{\pi}_i = \gamma_1(i)$$
2. 状态转移概率： $$\hat{a}_{ij} = \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{\sum_{t=1}^{T-1} \gamma_t(i)}$$ 分子是从状态 i 转移到 j 的期望次数，分母是处于状态 i 的期望次数。
3. 发射概率： $$\hat{b}_j(k) = \frac{\sum_{t=1}^{T} \gamma_t(j) \cdot \mathbb{I}(o_t = v_k)}{\sum_{t=1}^{T} \gamma_t(j)}$$ 分子是在状态 j 下观测到值 $v_k$ 的期望次数，分母是处于状态 j 的期望次数。$\mathbb{I}$ 是指示函数。

重复 E 步和 M 步直到参数收敛，得到的是观测似然的一个局部极大值（或驻点）附近的参数估计；Baum-Welch 不保证找到全局极大似然解，实践中常多次初始化并保留似然最高的结果。

Baum-Welch 算法流程（伪代码）：
1. 初始化模型参数 λ = (A, B, π)，可以随机或根据先验知识设定。
2. 循环直到收敛：
   a. (E步) 运行前向-后向算法，计算 γ_t(i) 和 ξ_t(i, j)。
   b. (M步) 根据上述公式，用 γ 和 ξ 重新估计 A, B, π。
   c. 检查参数变化是否小于阈值，或似然值 P(O|λ) 增长是否小于阈值。
3. 输出学习到的参数 λ。

维特比算法：解码最可能状态序列

给定模型参数 $\lambda$ 和观测序列 O，找出最可能（最大概率）的隐含状态序列 $Q^*$，即解码问题。这需要用维特比算法（Viterbi Algorithm），它也是一种动态规划算法。

定义 $\delta_t(i)$：在时刻 t，沿着一条路径到达状态 $q_t = i$，且产生观测 $o_1, o_2, ..., o_t$ 的最大概率。

$$\delta_t(i) = \max_{q_1, q_2, ..., q_{t-1}} P(q_1, q_2, ..., q_{t-1}, q_t = i, o_1, o_2, ..., o_t | \lambda)$$

同时，用 $\psi_t(i)$ 记录使 $\delta_t(i)$ 最大的前一个状态（回溯指针）。

递推步骤：

1. 初始化： $$\delta_1(i) = \pi_i b_i(o_1), \quad \psi_1(i) = 0$$
2. 递推： $$\delta_t(j) = \max_{1 \leq i \leq N} [\delta_{t-1}(i) a_{ij}] \cdot b_j(o_t)$$ $$\psi_t(j) = \arg\max_{1 \leq i \leq N} [\delta_{t-1}(i) a_{ij}]$$
3. 终止： $$P^* = \max_{1 \leq i \leq N} \delta_T(i)$$ $$q_T^* = \arg\max_{1 \leq i \leq N} \delta_T(i)$$
4. 路径回溯（t=T-1, …, 1）： $$q_t^* = \psi_{t+1}(q_{t+1}^*)$$

维特比算法与前向算法结构相似，但将求和 $\sum$ 换成了取最大值 $\max$，并增加了回溯步骤以记录路径。

# —— 重新构建上文用到的微型数据（本节各代码块相互独立，可单独运行）
import numpy as np
pi = np.array([0.6, 0.4])  # [S, R]
A = np.array([[0.8, 0.2],  # S->S, S->R
              [0.3, 0.7]]) # R->S, R->R
B = np.array([[0.6, 0.3, 0.1],  # S 下观测概率
              [0.1, 0.4, 0.5]]) # R 下观测概率
O = [0, 1, 2]
T = len(O)
N = A.shape[0]
alpha = np.zeros((T, N))
beta = np.zeros((T, N))

# 维特比算法示例（接续之前参数）
delta = np.zeros((T, N))
psi = np.zeros((T, N), dtype=int)

# 初始化
delta[0, :] = pi * B[:, O[0]]
psi[0, :] = 0
print(f"delta_1: {delta[0]}")

# 递推
for t in range(1, T):
    for j in range(N):
        # 计算所有可能的前驱状态 i 的 delta[t-1, i] * A[i, j]
        prob = delta[t-1, :] * A[:, j]
        psi[t, j] = np.argmax(prob)
        delta[t, j] = prob[psi[t, j]] * B[j, O[t]]
    print(f"delta_{t+1}: {delta[t]}")

# 终止
P_star = np.max(delta[T-1, :])
qT_star = np.argmax(delta[T-1, :])
print(f"最优路径概率 P* = {P_star:.6f}")
print(f"最终状态 q{T}_* = {qT_star} (0:S, 1:R)")

# 回溯
q_star = np.zeros(T, dtype=int)
q_star[T-1] = qT_star
for t in range(T-2, -1, -1):
    q_star[t] = psi[t+1, q_star[t+1]]

print(f"最可能状态序列: {q_star} (对应: {['S' if s==0 else 'R' for s in q_star]})")

delta_1: [0.36 0.04]
delta_2: [0.0864 0.0288]
delta_3: [0.006912 0.01008 ]
最优路径概率 P* = 0.010080
最终状态 q3_* = 1 (0:S, 1:R)
最可能状态序列: [0 1 1] (对应: ['S', 'R', 'R'])

根据维特比算法，对于观测序列 [Walk, Shop, Clean]，最可能的天气序列是 Sunny → Rainy → Rainy（输出中的 [0 1 1]），对应路径概率 $P^{*} = 0.010080$；而 Sunny → Sunny → Sunny 的概率只有 0.006912，不是最优。这也正是维特比与 gamma 的分野：逐时刻取多个 gamma 的最大值，得到的是 Sunny → Sunny → Rainy（0.8438、0.5159、0.7325），而维特比给出 Sunny → Rainy → Rainy，第 2、3 天的选择不同。维特比在整条路径上做最大化，因此它给出的序列未必是逐时刻边缘最优拼起来的那个。

📝 动手练一练

1. 概念辨析：前向算法中的 $\alpha_t(i)$ 和维特比算法中的 $\delta_t(i)$ 都表示在时刻 t 到达状态 i 的概率，它们之间最本质的区别是什么？这个区别导致了算法流程和最终用途有何不同？

2. 手算练习：考虑一个更简单的 HMM，只有两个状态 (A, B)，观测也只有两种 (X, Y)。参数如下： $\pi = [0.8, 0.2]$, $A = [[0.7, 0.3], [0.4, 0.6]]$, $B = [[0.9, 0.1], [0.2, 0.8]]$ (行对应状态 A, B，列对应观测 X, Y)。 观测序列为 O = [X, Y, X]。 a) 请手动计算前向概率 $\alpha_2(B)$（它是“到时刻 2 观测为 $o_1,o_2$ 且处于状态 B”的联合概率，不是时刻 2 处于状态 B 的后验概率）。 b) 这个序列的整体概率 $P(O|\lambda)$ 是多少？

本章小结

本节我们系统学习了隐马尔可夫模型（HMM），它是处理序列标注问题的基石。我们从其作为“动态混合模型”的直观理解出发，形式化定义了 HMM 的五元组。核心在于掌握其三大问题及对应算法：

1. 评估问题（计算 $P(O|\lambda)$）：通过前向算法（或后向算法）高效求解，复杂度 $O(TN^2)$。
2. 解码问题（寻找 $Q^*$）：通过维特比算法求解，同样基于动态规划，但使用 $\max$ 操作和路径回溯。
3. 学习问题（估计 $\lambda$）：通过 Baum-Welch 算法（EM 算法在 HMM 中的实现）求解，迭代地进行 E 步（计算 $\gamma_t, \xi_t$）和 M 步（重估参数）。

此外，我们推导了后向算法，并结合前向概率定义了状态后验概率 $\gamma_t(i)$，这为软标注和参数学习提供了基础。

HMM 的强大在于它用简洁的概率图模型刻画了状态转移和观测生成的过程。然而，它也有明显的局限性：其输出独立性假设和马尔可夫性假设在复杂的自然语言场景中往往过于强。例如，一个词的词性可能依赖于更远的上下文，而不仅仅是前一个词性。尽管如此，HMM 的思想深远地影响了后续模型的发展。循环神经网络（RNN、LSTM）可以看作是用神经网络单元替代了 HMM 中的概率表（A, B），从而能够学习更复杂、更长距离的依赖关系。注意力机制和 Transformer 则进一步打破了序列的线性递推结构，实现了全局依赖建模。

📋 行动清单

• 理解并能够手动推导前向、后向算法的递推公式，明确每个变量的物理意义。
• 针对一个简单的 HMM（如 2 状态 2 观测），能够手工执行前向算法、维特比算法各一步，并验证结果。
• 思考 HMM 的独立性假设在词性标注任务中可能造成什么问题。尝试举一个中文或英文的例子，其中当前词的词性依赖于前两个词，而 HMM 无法很好地建模。

—— 小象教研组