隐马尔可夫模型三大算法与工程实现：从原理到避坑指南

如果你做过序列建模，哪怕是两年经验，一定绕不开隐马尔可夫模型这个名字。学的时候感觉公式全认识，一上手写代码就发现到处是坑：前向算法算到第20步就下溢、维特比路径初始化搞错、鲍姆-韦尔奇训练出来一堆NaN。这篇就把HMM的知识点串一遍，从模型定义到三大经典问题，再到可复现的numpy实现，最后说说我在实际项目中踩过的坑和排查思路。适合正在学语音识别、基因序列分析、行为序列建模，或者只是面试前临时抱佛脚的读者。

1. HMM到底在解决什么问题

1.1 从一个可见序列说起

先把场景放简单。假设你在赌场里观察一个荷官连续抛硬币，记录每一轮的结果：正、正、反、正、反、反、正……如果只有一枚硬币，你完全可以假设每次抛掷独立同分布，用伯努利模型就能拟合。但问题是，这个荷官手边有两枚硬币，一枚是公平币，正面概率0.5，另一枚是作弊币，正面概率0.8。他什么时候换币你不知道，你只能看到他抛出的序列。

这种场景的本质是：你看到的观测序列背后，有一条看不见的状态链在驱动。每一时刻都有一个状态（当前用的是哪枚硬币），状态之间按概率转移（他可能偷偷换币），每个状态以一定概率发射出观测（硬币的正反面）。你要做的，就是从这个不完整的信息里反推出状态结构。

这就是HMM的全部意义：一个描述“隐藏状态”和“可见观测”之间概率关系的生成式模型。它不直接预测下一个观测，而是先建立“状态→观测”的因果关系，再在这条因果链上做推理。

1.2 隐藏的状态与可观测的发射

HMM的两个核心假设，是整个模型能成立的前提。

第一个假设是齐次马尔可夫性。通俗讲就是：t时刻的状态只依赖于t-1时刻的状态，跟更早的历史没有关系。注意这个假设很苛刻，真实序列里往往有更长程依赖，但HMM为了数学可解，必须做这个截断。

第二个假设是观测独立性。t时刻的观测只由t时刻的状态决定，和之前的观测、状态都无关。这一点在实际应用中经常被质疑，比如语音识别里相邻帧的声学特征明显存在相关性，但HMM依然能用，是因为后续还会引入差分特征、状态驻留扩展等手段来弥补。

这两个假设听起来限制很多，恰恰是它们让HMM的三个基本问题都有了高效的动态规划解法。换个复杂模型，比如带长程依赖的因果模型，计算复杂度就完全不可控了。

1.3 HMM的形式化定义

一个完整的HMM由五要素组成。

• 状态集合S = {s1, s2, ..., sN}，状态个数N，时刻t的状态记为q_t
• 观测集合V = {v1, v2, ..., vM}，观测个数M，时刻t的观测记为o_t
• 状态转移矩阵A = [a_ij]，a_ij = P(q_{t+1} = s_j | q_t = s_i)，表示从状态i跳到状态j的概率
• 发射矩阵B = [b_j(k)]，b_j(k) = P(o_t = v_k | q_t = s_j)，表示在状态j下观测到v_k的概率
• 初始状态分布π = [π_i]，π_i = P(q_1 = s_i)，第一时刻落在状态i的概率

注意A的每一行之和为1，B的每一行之和也为1，π之和为1。这是后面参数学习时做归一化的基础。

提示：很多人把HMM和马尔可夫链搞混。马尔可夫链的“状态”就是观测本身，是可见的；HMM的“状态”不可见，观测是状态的概率函数。两者差了一个隐变量层，分析方法和参数学习方法也就完全不同。

2. 三大核心问题，一个模型

2.1 评估问题：给定模型和观测，算概率

评估问题要回答的是：已知模型参数λ = (A, B, π)，给定观测序列O = (o1, o2, ..., oT)，求P(O|λ)是多少。

这个概率看起来很简单，把所有可能的状态序列Q的概率都累加起来就行：

P(O|λ) = Σ_Q P(O|Q, λ)P(Q|λ)

但状态序列数量是N^T量级，T稍微一大就指数爆炸。直接枚举在工程上完全不可行，必须用动态规划来做。

评估问题的实际用途很广。语音识别里，每个词对应一个HMM，输入一段音频特征，算它和每个词模型的匹配概率，选概率最高的作为识别结果。生物信息学里，用HMM扫描基因组序列，看它和已知基因结构的兼容程度。

2.2 解码问题：给定模型和观测，反推状态路径

解码问题要回答的是：观测序列已知，最可能的状态序列是什么？

公式上看就是求argmax_Q P(Q|O, λ)。注意这里和评估问题不同，不是把所有路径的概率加起来，而是找出那一条概率最大的路径。

这个问题的典型场景是分词和词性标注。把每个词的词性看成隐藏状态，把词本身看成观测，解码得到的路径就是最可能的词性序列。命名实体识别也常用这个思路。

2.3 学习问题：只有观测，反推模型参数

学习问题最硬核：只知道一组或几组观测序列，不知道状态序列，也不知道任何参数，要求出λ = (A, B, π)，使得P(O|λ)最大。

这比前两个问题难得多，因为中间隔着一层隐变量。直接求导基本不可能，标准解法是EM算法的一个特例——鲍姆-韦尔奇算法。后面第5章详细展开。

第三个问题不是可有可无的附加题，它是HMM能落地的前提。大部分真实场景你手里只有标注好的观测序列，根本没有状态序列标注。能自己从数据里学出参数，模型才算真正能用。

3. 前向算法：把概率算到不爆炸

3.1 前向变量的递推思路

评估问题的动态规划解法就是前向算法，核心是定义一个前向变量：

α_t(i) = P(o1, o2, ..., o_t, q_t = s_i | λ)

这个变量的含义是：到t时刻为止的全部观测都已经出现，并且t时刻状态是s_i的联合概率。

它的递推关系特别漂亮：

α_1(i) = π_i * b_i(o1)

α_{t+1}(j) = [Σ_{i=1}^N α_t(i) * a_ij] * b_j(o_{t+1})

最后得到P(O|λ) = Σ_{j=1}^N α_T(j)

我解释一下这个递推为什么是对的。t+1时刻落在状态j，必须先从t时刻的某个状态i跳到j，这一步的贡献是α_t(i) × a_ij；把所有i的可能加起来，就是在t+1时刻“到达状态j且观测到o1到o_t”的累积概率。然后再乘上在状态j下观测到o_{t+1}的发射概率b_j(o_{t+1})，就补上了当前时刻的观测。

这个算法的时间复杂度是O(T×N²)，比指数级的暴力枚举降了好几个数量级。

3.2 前向算法的numpy实现

我用numpy写一个最小可运行版本，并用一个“两枚硬币”的例子做验证。

import numpy as np def forward(A, B, pi, observations): N = A.shape[0] T = len(observations) alpha = np.zeros((T, N)) # 初始化：t=0时各状态的概率 alpha[0] = pi * B[:, observations[0]] # 递推 for t in range(1, T): for j in range(N): alpha[t, j] = np.sum(alpha[t-1] * A[:, j]) * B[j, observations[t]] return alpha, np.sum(alpha[T-1])

测试数据设置：状态0是公平币（正面概率0.5），状态1是偏币（正面概率0.8）。转移矩阵设为A = [[0.9, 0.1], [0.1, 0.9]]，初始分布π = [0.5, 0.5]。观测序列用0表示反面，1表示正面。

A = np.array([[0.9, 0.1], [0.1, 0.9]]) B = np.array([[0.5, 0.5], [0.2, 0.8]]) pi = np.array([0.5, 0.5]) obs = [1, 1, 0, 1, 1, 0, 0, 1] alpha, prob = forward(A, B, pi, obs) print(f"P(O|lambda) = {prob:.6f}")

跑一遍得到的结果是一个很接近1的小数，因为模型参数都是归一化的，观测序列的各种组合概率加在一起是1，具体那个序列占的比例就是它。用这个方法算出来的概率，能直接拿来做两个模型的横向对比。

注意：前向算法算的是联合概率，不是条件概率P(O|λ)。所以在比较不同模型时，只用相对大小就够了，绝对值偶尔会小到让人觉得“是不是算错了”，其实数值范围正常。

4. 维特比算法：找出那条最可能的路径

4.1 动态规划、回溯和路径重构

解码问题的思路和前向算法类似，但不是求和而是求最大，这里用到维特比算法。定义一个变量：

δ_t(i) = max P(q1, q2, ..., q_t = s_i, o1, o2, ..., o_t | λ)

含义是：到t时刻状态是s_i，且观测序列完整出现的前提下，前t-1步状态路径能取到的最优概率。

递推关系：

δ_1(i) = π_i * b_i(o1)

δ_{t+1}(j) = max_i [δ_t(i) * a_ij] * b_j(o_{t+1})

同时要维护一个回溯表ψ_t(j)，记录δ_{t+1}(j)取到最大值时，上一时刻是哪个状态。这样算完最后一层之后，从δ_T里挑最大的那个状态作为终点，再顺着ψ表一路回溯，就得到完整的最优状态路径。

维特比和前向算法长得几乎一样，只有两个符号的区别：一个是“求和”，一个是“取最大值”。理解了这点，两个算法就变成一套思路了。

4.2 维特比的代码实现

def viterbi(A, B, pi, observations): N = A.shape[0] T = len(observations) delta = np.zeros((T, N)) psi = np.zeros((T, N), dtype=int) # 初始化 delta[0] = pi * B[:, observations[0]] psi[0] = 0 # 递推 for t in range(1, T): for j in range(N): temp = delta[t-1] + np.log(A[:, j]) psi[t, j] = np.argmax(temp) delta[t, j] = np.max(temp) + np.log(B[j, observations[t]]) # 回溯路径 path = [np.argmax(delta[T-1])] for t in range(T-1, 0, -1): path.append(psi[t, path[-1]]) path.reverse() return np.array(path)

这里有个细节：我用log概率来避免连乘下溢。在前向算法里我们直接连乘，因为最后求和时数值还在可接受范围；但维特比要比较不同路径的概率大小，连乘一旦下溢，结果可能全是0，比较就失去了意义。换成log空间，max操作等价于在原始概率空间做max，结果不受影响。

验证一下，还是用上面那个硬币例子：

path = viterbi(A, B, pi, obs) print(path)

输出的状态序列会清楚显示哪几步模型认为在“公平币状态”，哪几步在“偏币状态”。这个路径就是给定观测下，最有解释力的状态走法。

4.3 一个容易忽视的问题：多条最优路径

理论上，可能存在多条路径得到完全相同的最大概率。比如状态个数少、观测序列又短的时候，A和B的参数碰巧让两条路径的打分一模一样。维特比算法只会返回其中一条，至于返回哪一条，取决于argmax的判定顺序。

这在绝大多数应用里无所谓。但如果你依赖解码结果做后续统计，比如用状态路径去统计每个状态的频次，就要注意这个偏差，必要时可以加上一个微小的随机扰动打破平局，或改用“前向后向算法”计算每个时刻每个状态的概率期望。

5. 鲍姆-韦尔奇算法：从观测里自学参数

5.1 EM思想的落地

学习问题没有状态序列标注，相当于一个带隐变量的极大似然估计，直接用MLE是无解的。鲍姆-韦尔奇算法本质上就是EM算法。

EM的思路分两步走。

E步：用当前参数λ，计算隐变量的后验期望。具体到HMM，就是算两个量：

γ_t(i) = P(q_t = s_i | O, λ)

ξ_t(i, j) = P(q_t = s_i, q_{t+1} = s_j | O, λ)

这两个量需要同时用到前向变量α和新的后向变量β。β_t(i)定义为：t时刻状态为s_i，从t+1到T时刻观测已知的概率，也叫后向变量，递推公式直接反着来。

M步：用这些期望去重新估计参数。转移矩阵的重估公式是：

a_ij_new = Σ_t ξ_t(i, j) / Σ_t γ_t(i)

发射矩阵的重估公式是：

b_j(k)new = Σ{t:o_t=vk} γ_t(j) / Σ_t γ_t(j)

初始分布就直接取γ_1(i)。

迭代到对数似然收敛，就用学到的参数作为最终模型。

5.2 训练代码的骨架

把前向、后向和重估合在一起，就是一个可训练的HMM类。我这里只写关键逻辑，完整的封装长度太长，拆开讲更清楚。

def baum_welch(A, B, pi, observations, iterations=100): N = A.shape[0] T = len(observations) for _ in range(iterations): # E步 alpha, _ = forward(A, B, pi, observations) beta, _ = backward(A, B, pi, observations) gamma = np.zeros((T, N)) xi = np.zeros((T-1, N, N)) for t in range(T-1): denom = np.sum(alpha[t] * beta[t]) for i in range(N): gamma[t, i] = alpha[t, i] * beta[t, i] / denom for j in range(N): xi[t, i, j] = (alpha[t, i] * A[i, j] * B[j, observations[t+1]] * beta[t+1, j]) / denom # M步 pi = gamma[0] / np.sum(gamma[0]) for i in range(N): for j in range(N): A[i, j] = np.sum(xi[:, i, j]) / np.sum(gamma[:, i]) for j in range(N): for k in range(B.shape[1]): mask = np.array(observations) == k B[j, k] = np.sum(gamma[mask, j]) / np.sum(gamma[:, j]) return A, B, pi

这里注意，denom计算时用的是t时刻的α和β乘积和，而不是全局归一化。因为α_t(i)×β_t(i)恰好等于P(O, q_t = s_i | λ)，对所有i求和就是P(O|λ)，这就是γ的归一化常数。这个等价关系是前向后向算法最精巧的一点。

5.3 训练时的收敛判据

实操中不能傻跑固定迭代次数，通常监控每一轮迭代的对数似然log P(O|λ)，当它的增长小于某个阈值（比如1e-4）时停止。注意对数似然一定是单调不减的，这是EM算法的理论保证。如果发现对数似然突然下降，几乎可以肯定是代码bug或者数值溢出。

注意：鲍姆-韦尔奇只能收敛到局部最优，不能保证全局最优。多初始化几组参数、分别训练、选对数似然最高的模型，是常规操作。

6. 实操中常见的坑和排查思路

6.1 数值下溢：最普遍的杀手

前向算法里连乘几百个小于1的概率，很快就变成0。维特比因为用了log所以安全，前向和鲍姆-韦尔奇就没这么幸运。

标准解法是缩放（scaling）。在每一层递推之后，用该时刻全部α值的和去归一化，记录下缩放系数，最后把所有系数乘回来。具体做法：

c = np.zeros(T) for t in range(T): c[t] = np.sum(alpha[t]) alpha[t] /= c[t]

对应的β也用同样的系数缩放。到最后计算对数似然时：

log P(O|λ) = Σ_t log c_t

因为每一步的缩放因子连乘，正好把前面损失的概率量级补回来了。这个技巧是所有工程实现HMM的标配，我见过不少学习者直接拿未缩放的版本做长序列，结果第100步全部变成0，怎么查都查不出问题，其实就是这里。

6.2 初始化参数不能太随意

A和B的初始值如果设成全0或全相等，EM算法很可能直接收敛到一个收敛点，但模型完全没区分度。常见做法是：

• 转移矩阵：每个状态转移概率设一个中等数值（0.5到0.9），再留一点转移概率给其他状态
• 发射矩阵：用均匀分布，但加一点轻微扰动
• 初始分布：设为均匀，或根据先验知识设成某个状态更可能

还有一个更重要的初始化技巧：如果计算资源允许，做10次随机初始化，选对数似然最高的模型。这比调一大堆超参数实在多了。

6.3 局部最优和过拟合

鲍姆-韦尔奇对初始值敏感，这是EM的通病。我实操时遇到过几次：同一个数据集，从不同初始值出发，最后收敛出的模型结构天差地别，其中一个模型的对数似然明显更高。

处理思路有几个层面。一是多随机初始化，二是降低自由度，比如约束发射矩阵的结构，三是加入正则化，比如对转移概率做Dirichlet先验的MAP估计。对于小数据集，后两点尤其有效，因为纯MLE会在参数空间里乱跑，把噪声也拟合进去。

6.4 常见问题速查表

6.5 一个从实践中得来的调参心得

我刚做HMM项目时，最爱犯的错误是忽略状态个数的选择。N设得太多，模型过拟合，训练慢；N设得太少，表达力不够，序列里真实存在的结构拟合不出。有一种相对有效的做法是：画对数似然随N变化的曲线，找拐点。N从2开始往上加，每加1训练一次模型，看对数似然提升幅度，当提升明显变缓时就别再加了。这个办法是经验性的，但比拍脑袋选N靠谱得多。

7. 沿这个方向还能往哪走

HMM本身是一个很“底层”的工具箱。掌握三大算法之后，往上可以扩展出很多实用的变体：

• 高斯混合HMM（GMM-HMM）：把发射概率从离散分布换成高斯混合，语音识别里非常主流
• 输入输出HMM（IOHMM）：把观测也纳入状态转移的条件，处理带外部输入的序列
• 层次HMM、分段HMM：处理长时依赖和状态驻留时间建模
• 和循环神经网络结合：用RNN做特征提取，用HMM做序列解码，取两者之长

如果你主要是为了解决序列标注问题，还可以对比条件随机场的方案。HMM是生成式模型，优化的是联合概率；CRF是判别式模型，直接建模条件概率。数据量大的时候，CRF通常精度更高；数据量小或者需要生成式特性的时候，HMM更好使。没有哪个绝对最优，取决于你的数据和任务形态。

我个人在实际项目里体会最深的一点是：不要一上来就让HMM解决所有问题。它最适合的场景是状态转移结构相对清晰、观测噪声可控、序列长度适中。一旦序列极长或者状态依赖跨度很大，HMM的效果就会明显下滑，这时候就该考虑用循环网络或者状态空间模型去替代它。模型选择本身，比调参重要得多。