机器学习数学推导（十五）：隐马尔可夫模型
从一个原理推出 HMM 的三大经典算法：把联合分布按时间因子化，再用动态规划复用跨时间的子计算。覆盖前向后向的边缘与平滑、Viterbi 的 MAP 解码，以及 Baum-Welch（EM）的参数学习。
雾里有人在你身后走过。看不见人，只听见脚步——短促、轻、急。从节奏和音色，你能猜出对方是在走、跑，还是跛着腿吗？如果听到一整段声音呢？哪条步态序列最可能产生它？又或者，在你对"走路"建立的模型下，这段声音本身有多大概率出现？
这就是 HMM 的三个基本问题。令人惊讶的是，三者归根结底是同一个把戏：把联合概率 $P(\mathbf{O}, \mathbf{I})$ 按时间因子化，再用动态规划跨时间复用子计算。暴力枚举要算 $O(N^T)$ 次，前向后向、Viterbi、Baum-Welch 全部只要 $O(N^2 T)$。指数塌缩成多项式，靠的就是马尔可夫假设——给定现在，未来与过去条件独立。
你将学到什么
- HMM 的联合分布与两个条件独立性假设
- 前向 / 后向算法：边缘概率 $P(\mathbf{O}\mid\lambda)$ 与平滑后验 $\gamma_t(i), \xi_t(i,j)$
- Viterbi：把
sum换成max，DP 框架原封不动地解 MAP - Baum-Welch：EM 在 HMM 上的特化；为什么 M 步公式就是"期望计数"
- 数值陷阱（下溢、缩放）与现代后裔（CRF、RNN、CTC）
前置知识
- 马尔可夫链与随机矩阵
- 动态规划 / 记忆化
- EM 算法与 ELBO（参见第十三篇 ）
1. 模型本身
HMM 是序列建模中最简单又非平凡的隐变量模型。两条平行的链沿时间展开：
- 隐状态链 $\mathbf{I}=(i_1,\dots,i_T)$，$i_t\in\{1,\dots,N\}$；
- 观测链 $\mathbf{O}=(o_1,\dots,o_T)$，$o_t\in\{v_1,\dots,v_M\}$。
模型由三组参数 $\lambda=(\boldsymbol{\pi}, \mathbf{A}, \mathbf{B})$ 完整描述：
| 模块 | 符号 | 含义 |
|---|---|---|
| 初始分布 | $\pi_i = P(i_1 = i)$ | 链在 $t=1$ 处怎么"出生" |
| 转移矩阵 | $a_{ij} = P(i_{t+1}=j \mid i_t=i)$ | 状态如何演化 |
| 发射矩阵 | $b_j(k) = P(o_t = v_k \mid i_t = j)$ | 状态如何辐射证据 |
两个条件独立性假设让一切变得可计算：
- 状态的一阶马尔可夫性：$P(i_{t+1}\mid i_{1:t}) = P(i_{t+1}\mid i_t)$。
- 观测独立性：$o_t$ 只依赖 $i_t$，与其它状态、其它观测无关。
两条假设合在一起，让我们能把整张联合分布按时间因子化：
$$P(\mathbf{O}, \mathbf{I} \mid \lambda) = \pi_{i_1}\,b_{i_1}(o_1) \prod_{t=2}^{T} a_{i_{t-1},i_t}\,b_{i_t}(o_t).$$后文所有算法，都是在不枚举 $N^T$ 条隐路径的前提下，对这个乘积做求和或求最大的精巧办法。
三状态天气玩具模型
全文用一个三状态天气模型——晴 / 雨 / 多云——做演示，转移概率如图。“晴"具有黏性（$a_{\text{晴晴}}=0.70$），“多云"则是连接其它状态的枢纽。规模小到可以手算，结构又足以暴露算法所有细节。
三个问题，一个联合分布
给定 $\lambda$，我们能问的问题穷举如下：
| # | 问题 | 输入 | 输出 | 算法 |
|---|---|---|---|---|
| 1 | 概率计算 | $\lambda, \mathbf{O}$ | $P(\mathbf{O}\mid\lambda)$ | 前向 / 后向 |
| 2 | 解码 | $\lambda, \mathbf{O}$ | $\arg\max_{\mathbf{I}} P(\mathbf{I}\mid\mathbf{O},\lambda)$ | Viterbi |
| 3 | 学习 | $\mathbf{O}$（与 $N$） | $\hat\lambda = \arg\max_\lambda P(\mathbf{O}\mid\lambda)$ | Baum-Welch（EM） |
对（1）若直接枚举所有隐路径，$N{=}50, T{=}100$ 时大约要算 $10^{170}$ 次。下面三节把全部三个问题压到多项式时间。
2. 前向算法：从左向右扫一遍
定义前向变量
$$\alpha_t(i) = P(o_1, o_2, \dots, o_t,\; i_t = i \mid \lambda),$$它表示"已经生成了到时刻 $t$ 为止的观测，并且此刻处在状态 $i$“的联合概率。
初始化。 在 $t=1$ 时还没有发生转移，只有初始分布与第一次发射：
$$\alpha_1(i) = \pi_i\, b_i(o_1).$$递推。 把 $\alpha$ 从 $t-1$ 推到 $t$：对前一时刻的状态做条件、再边缘化：
$$\begin{aligned} \alpha_t(j) &= P(o_1,\dots,o_t,\, i_t=j) \\ &= \sum_{i=1}^N P(o_1,\dots,o_{t-1},\, i_{t-1}=i)\,P(i_t=j\mid i_{t-1}=i)\,P(o_t\mid i_t=j)\\ &= \left[\sum_{i=1}^N \alpha_{t-1}(i)\, a_{ij}\right] b_j(o_t). \end{aligned}$$方括号里那个求和，是唯一跨越时间边界的计算；其余都是本地的。这就是动态规划最干净的样子。
终止。 把最后一时刻的状态求和掉：
$$P(\mathbf{O}\mid\lambda) = \sum_{i=1}^N \alpha_T(i).$$复杂度。 共 $T$ 个时刻，每步 $N$ 个目标状态，每个目标状态需要对 $N$ 个前驱求和，总计 $O(N^2 T)$。当 $N{=}50, T{=}100$ 时约 $2.5\times 10^{5}$ 次操作，比暴力法快约 $10^{165}$ 倍。
下溢问题。 概率几何级数地相乘，$\alpha_t$ 早晚会溢出。两种标准做法：
- 对数空间 +
logsumexp：$\log\alpha_t(j) = \mathrm{lse}_i\!\big(\log\alpha_{t-1}(i) + \log a_{ij}\big) + \log b_j(o_t)$。 - 缩放前向：每一步把 $\alpha_t$ 除以 $c_t = \sum_j \alpha_t(j)$，并累加 $\log P(\mathbf{O}) = \sum_t \log c_t$。
3. 后向算法与后验平滑
只有前向，能算 $P(\mathbf{O})$；要拿到每一时刻隐状态的后验分布，还得从右往左再扫一遍。
定义后向变量
$$\beta_t(i) = P(o_{t+1}, o_{t+2}, \dots, o_T \mid i_t = i, \lambda).$$注意条件方向：$\beta$ 是在已知当前状态下、关于未来观测的条件概率；$\alpha$ 则是关于过去观测与状态的联合概率。两者结构对称。
边界。 $T$ 之后没有观测了，所以对所有 $i$ 都有 $\beta_T(i) = 1$。
递推。 从 $t+1$ 反推到 $t$，对下一个状态求和：
$$\beta_t(i) = \sum_{j=1}^N a_{ij}\, b_j(o_{t+1})\, \beta_{t+1}(j).$$自洽性检查。 两次扫描合起来必须能复原边缘概率：
$$P(\mathbf{O}\mid\lambda) = \sum_i \pi_i\,b_i(o_1)\,\beta_1(i) = \sum_i \alpha_T(i).$$驱动学习的两个后验
把 $\alpha$、$\beta$ 都缓存好之后，下面两个量都能用 $O(1)$ 额外开销算出来：
单点状态后验（平滑分布）：
$$\gamma_t(i) = P(i_t=i \mid \mathbf{O},\lambda) = \frac{\alpha_t(i)\,\beta_t(i)}{P(\mathbf{O}\mid\lambda)}.$$相邻对后验（转移责任度）：
$$\xi_t(i,j) = P(i_t=i, i_{t+1}=j \mid \mathbf{O},\lambda) = \frac{\alpha_t(i)\,a_{ij}\,b_j(o_{t+1})\,\beta_{t+1}(j)}{P(\mathbf{O}\mid\lambda)}.$$它们正是 Baum-Welch 在 E 步要算的两个统计量。
4. Viterbi：把求和换成求最大
解码问的是最可能的状态序列，而不是总概率：
$$\mathbf{I}^* = \arg\max_{\mathbf{I}}\, P(\mathbf{I}\mid\mathbf{O},\lambda) = \arg\max_{\mathbf{I}}\, P(\mathbf{O},\mathbf{I}\mid\lambda).$$（分母 $P(\mathbf{O})$ 不依赖 $\mathbf{I}$，可以扔掉。）
定义
$$\delta_t(j) = \max_{i_1,\dots,i_{t-1}} P(i_1,\dots,i_{t-1}, i_t=j,\, o_1,\dots,o_t\mid \lambda),$$也就是"长度为 $t$ 且终点在 $j$“的所有路径中概率最大值。递推与前向一模一样，只把一个算子换掉：
$$\boxed{\;\delta_t(j) = \max_{i}\big[\delta_{t-1}(i)\, a_{ij}\big]\, b_j(o_t).\;}$$为了能复原路径本身，还要存回溯指针：
$$\psi_t(j) = \arg\max_i\big[\delta_{t-1}(i)\, a_{ij}\big].$$前向扫完后取 $i_T^* = \arg\max_i \delta_T(i)$，再回溯：$i_t^* = \psi_{t+1}(i_{t+1}^*)$。
为什么把 sum 换成 max 还能成立？ 因为两个算子在按时间因子化的联合分布上都满足分配律——“max-times” 和 “sum-times” 都构成可交换的半环。动态规划的骨架完全一致，仅本地约简方式不同。
数值版本。 实际中一律在对数空间跑：
$$\log\delta_t(j) = \max_i\big[\log\delta_{t-1}(i) + \log a_{ij}\big] + \log b_j(o_t).$$全部变成加法，下溢从根本上不会发生。
5. Baum-Welch：HMM 上的 EM
当 $\lambda$ 未知时，目标是最大化 $\log P(\mathbf{O}\mid\lambda)$。隐状态 $\mathbf{I}$ 是隐变量，自然要用 EM。这就是 Baum-Welch 算法——它甚至比通用的 EM 框架还早几年（Baum & Petrie, 1966），是 EM 应用中最优雅的范例之一。
E 步：后验统计量
给定当前 $\lambda^{(k)}$，跑一次前向后向，得到 $\gamma_t(i)$ 与 $\xi_t(i,j)$。它们正是指示变量 $\mathbb{1}[i_t=i]$ 和 $\mathbb{1}[i_t=i, i_{t+1}=j]$ 在当前模型下的期望值。
完全数据对数似然按时间因子化：
$$\log P(\mathbf{O},\mathbf{I}\mid\lambda) = \log\pi_{i_1} + \sum_{t=1}^{T-1}\log a_{i_t i_{t+1}} + \sum_{t=1}^{T}\log b_{i_t}(o_t).$$在后验 $P(\mathbf{I}\mid \mathbf{O},\lambda^{(k)})$ 下取期望，每个 $\mathbb{1}[\cdot]$ 都被它的后验概率替换掉：
$$Q(\lambda;\lambda^{(k)}) = \sum_i \gamma_1(i)\log\pi_i + \sum_{t=1}^{T-1}\sum_{i,j}\xi_t(i,j)\log a_{ij} + \sum_{t=1}^{T}\sum_{j,k}\gamma_t(j)\,\mathbb{1}[o_t=v_k]\log b_j(k).$$三项之间的参数 $\boldsymbol{\pi}, \mathbf{A}, \mathbf{B}$ 完全解耦——M 步因此化为三个相互独立的带约束最大化。
M 步：三组闭式更新
每组参数都带有归一化约束（每行之和为 1）。引入拉格朗日乘子后，每条更新公式都化作：
$$\text{参数} = \frac{\text{该事件的期望计数}}{\text{条件上下文的期望计数}}.$$具体地：
$$\boxed{\; \hat\pi_i = \gamma_1(i),\qquad \hat a_{ij} = \frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)},\qquad \hat b_j(k) = \frac{\sum_{t:\,o_t=v_k}\gamma_t(j)}{\sum_{t=1}^{T}\gamma_t(j)}. \;}$$把每个比值读成”$i\to j$ 的期望转移次数 / 离开 $i$ 的期望总次数”，发射也是同理。形式与可观测数据下的 MLE 一致，只是把"实际观测到的次数"换成"后验下的期望次数”。
收敛性
EM 保证 $P(\mathbf{O}\mid\lambda^{(k+1)}) \geq P(\mathbf{O}\mid\lambda^{(k)})$——上图的曲线在理论上必然单调上升。陷阱在于似然面是多峰的，Baum-Welch 只能找到局部最优。常用对策是随机重启、用 k-means / segmental k-means 初始化，或加入信息性先验。
6. 应用：词性标注
词性标注（POS tagging）是 HMM 的教科书应用。隐状态是词性（PRON、VERB、ADJ、NOUN……），观测是词。转移矩阵编码语法规则（限定词后多接形容词或名词）；发射矩阵编码词典（love 多为动词，偶为名词）。
对句子 “I love natural language processing”，Viterbi 路径点亮为 PRON / VERB / ADJ / NOUN / NOUN——即便 processing 在词典上有歧义，结果依然正确，因为 ADJ -> NOUN 转移概率很高。
同一套引擎也活跃在语音识别（状态 = 音素，观测 = MFCC 帧）、生物信息的 profile HMM（状态 = 匹配 / 插入 / 删除列）、手势识别等任何"离散隐过程产生噪声观测流"的场景中。
深入问答
Q1：前向 vs Viterbi——换个算子为什么这么关键？ 前向给出边缘 $P(\mathbf{O}\mid\lambda) = \sum_{\mathbf{I}} P(\mathbf{O},\mathbf{I})$；Viterbi 给出 $\max_{\mathbf{I}} P(\mathbf{O},\mathbf{I})$。同一套 DP 骨架，不同的半环（sum-product vs max-product）。前向回答"这段证据有多可信”，Viterbi 回答"用什么故事最能解释它"。
Q2：Viterbi 为什么直接最大化联合而非后验？ 因为 $P(\mathbf{I}\mid\mathbf{O}) = P(\mathbf{O},\mathbf{I})/P(\mathbf{O})$，分母对 $\mathbf{I}$ 是常数。最大化联合等价于最大化后验，还省下一次归一化。
Q3：Baum-Welch 在什么情况下会翻车？ 三种典型失败模式：(a) 初始化不好——掉进平坦或退化的局部极值；(b) 标签置换 ——状态只能在置换意义下被识别；(c) 观测坍缩——某个状态的发射只覆盖训练里见过的符号，未见符号概率为零。对 (c) 的修复：加一个小 Dirichlet 先验做平滑。
Q4：序列标注为什么常用 CRF 而不用 HMM？ CRF 是判别式模型，直接建模 $P(\mathbf{I}\mid\mathbf{O})$，可以用 $\mathbf{O}$ 的全局重叠特征（大小写、后缀模式、上下文窗口），不必受条件独立的束缚。但当训练数据少、或者需要生成序列时，HMM 仍然是更自然的选择。
Q5：深度学习时代 HMM 还有用吗？ 作为端到端独立建模的方案，基本被 RNN / Transformer 替代。但推断算法仍在：现代语音里的 CTC 解码本质上就是在对齐网格上的前向算法；序列级蒸馏用的是 Viterbi；结构化输出 Transformer 又借鉴了 CRF，而 CRF 又借鉴了 HMM。
Q6：状态数 $N$ 怎么选？ 先小后大，再用信息准则（$\text{AIC} = -2\log L + 2|\lambda|$，$\text{BIC} = -2\log L + |\lambda|\log T$）或留出似然来选；贝叶斯非参数（用层次 Dirichlet 过程作先验的 iHMM）则把 $N$ 也看作随机变量交给数据决定。
Q7：连续观测怎么处理？ 把发射矩阵换成密度即可：单个高斯、混合高斯（语音里经典的 GMM-HMM）、或神经密度（混合密度网络 -> “neural HMM”）。前向后向递推一字不变，只是 $b_j(o_t)$ 变成一次似然评估。
习题
E1. 取 $N=2$，$\boldsymbol{\pi}=(0.6, 0.4)$，$\mathbf{B} = \begin{pmatrix}0.5 & 0.5 \\ 0.4 & 0.6\end{pmatrix}$，$o_1 = v_1$，求 $\alpha_1$。
解答
$\alpha_1(1) = 0.6 \cdot 0.5 = 0.30$，$\alpha_1(2) = 0.4 \cdot 0.4 = 0.16$。E2. 当 $N{=}50, T{=}100$ 时，比较前向算法与暴力枚举的复杂度。
解答
前向：$N^2 T = 2.5\times 10^{5}$ 次乘法。暴力：$N^T \approx 10^{170}$ 条路径——不可计算。E3. 解释 $\hat a_{ij} = \frac{\sum_t \xi_t(i,j)}{\sum_t \gamma_t(i)}$ 的物理意义。
解答
$i\to j$ 的期望转移次数 / 从状态 $i$ 离开的期望总次数。是把硬性 MLE 计数推广到软性后验计数。E4. 证明：先跑前向再跑后向，得到的 $P(\mathbf{O}\mid\lambda)$ 与单独跑前向一致。
解答
$\sum_i \pi_i b_i(o_1)\beta_1(i) = \sum_i \alpha_1(i)\beta_1(i)/1 = \sum_i \alpha_t(i)\beta_t(i)$ 对任意 $t$ 成立（链式法则）；当 $t = T$ 时 $\beta_T \equiv 1$，于是结果就是 $\sum_i \alpha_T(i)$。E5. 解释为什么 Viterbi 时间复杂度是 $O(N^2 T)$，而回溯指针的空间是 $O(NT)$。
解答
正向网格扫一遍要 $O(N^2 T)$（与前向同级），每个 (状态, 时刻) 单元再额外存一个整数回溯指针——共 $NT$ 个槽位——回溯阶段线性时间消耗。参考文献
- Rabiner, L. R. (1989). A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proc. IEEE 77(2), 257-286.
- Baum, L. E., & Petrie, T. (1966). Statistical Inference for Probabilistic Functions of Finite State Markov Chains. Annals of Math. Statist. 37(6).
- Viterbi, A. J. (1967). Error Bounds for Convolutional Codes and an Asymptotically Optimum Decoding Algorithm. IEEE Trans. IT 13(2).
- Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective, 第 17 章. MIT Press.
- Bishop, C. M. (2006). Pattern Recognition and Machine Learning, 第 13 章. Springer.
- Jelinek, F. (1997). Statistical Methods for Speech Recognition. MIT Press.
- Lafferty, J., McCallum, A., & Pereira, F. (2001). Conditional Random Fields. ICML.
- Eddy, S. R. (1998). Profile Hidden Markov Models. Bioinformatics 14(9).