HMM隐马尔可夫模型Python实现：前向算法、Viterbi与Baum-Welch实战-育师

HMM（Hidden Markov Model，隐马尔可夫模型）这个名词，做 NLP、语音识别、行为序列分析的开发者应该都不陌生。我记得刚入门的时候，教材翻了好几遍，图模型和公式推导都能看懂，但一到写代码就卡住——alpha、beta、delta 这些变量到底怎么组织，矩阵每一行每一列代表什么，完全没有实感。后来我把三套核心算法各手写了一遍，才发现 HMM 的代码难点根本不在算法本身，而在索引对齐和工程细节。这篇文章我准备用一个可运行的 Python 实现，把前向算法、Viterbi 解码、Baum-Welch 训练逐行拆开讲清楚。适合刚啃完公式想动手复现的读者，也适合工作中要拿 HMM 做序列标注、行为识别、异常检测的朋友。我会先用最朴素的写法把原理讲明白，再给出能直接搬到工程里的优化版本。

1. HMM 到底是什么：先搞懂它在解决什么问题

1.1 两个序列：状态隐藏，观测可见

还是用那个最经典的"天气与行为"例子。假设有位朋友每天出门，行为只有三种：散步、购物、打扫，行为完全由当天天气决定，而天气只有两种：晴天、雨天。问题来了：我看不到天气，只能连续记录他几天的行为，我想知道他这几天分别是什么天气。

这里就有两个序列。一个序列是每天的天气，它是"隐藏"的，我观测不到；另一个序列是每天的行为，这是"可见"的观测。HMM 要解决的，就是在这两个序列之间建立概率关系。

注意 HMM 的名字里有个 Hidden——状态本身不是随机独立冒出来的，而是由一条马尔可夫链生成的，只是链条本身被藏起来了，我们只能通过观测去推测它。很多人第一次接触会问：这跟普通马尔可夫链有什么区别？普通马尔可夫链的状态是直接可见的，比如用天气预测明天下不下雨，直接看今天天气就行。而 HMM 里你根本看不到状态，只能通过观测的"投影"去推断。这个"看不见的状态"设定，恰好对应了大量真实场景：NLP 里做词性标注，词性是隐藏的、句子是可见的；语音识别里音素是隐藏的、声学特征是可见的；设备故障诊断里机器状态是隐藏的、传感器读数是可见的。理解了这一层，HMM 的适用范围就清楚了。

1.2 两个核心假设：HMM 凭什么能算

HMM 敢用这么简洁的结构，靠的是两个假设。

第一，齐次马尔可夫假设：t 时刻的状态只由 t-1 时刻的状态决定，跟更早的历史无关。就像你判断明天天气只参考今天，不看前天和昨天。第二，观测独立性假设：t 时刻的观测只由 t 时刻的状态决定，跟其他时刻的状态和观测都无关。说白了，就是今天的散步行为只由今天的天气决定，跟昨天天气没有直接关系。

这两个假设在真实世界里当然不严格成立，但它们是 HMM 能高效求解的基石。如果没有这两个假设，状态和观测之间的关系会变成指数级复杂的网络，别说训练，连推理都做不了。工程上我们常说所有模型都是错的，但有些是有用的——HMM 就是典型"刻意简化到能算"的模型，而且简化完之后在不少场景依然相当能打。所以别上来就抱怨假设太强，先想想这个假设换来了什么：把联合概率拆成了一串可连乘的小因子，这就是所有 HMM 算法能跑起来的根本原因。

1.3 五个要素：一个 HMM 的完整定义

具体到一个 HMM，由五个东西组成，代码里对应的就是几个矩阵和向量。我先把符号和后面代码里的变量对应起来，这是写代码前最重要的一步，多少 bug 都是索引对应错导致的。

数学符号	含义	代码中的变量	维度
Q / N	状态集合 / 状态数	num_states	N
V / M	观测集合 / 观测数	num_obs	M
A	状态转移矩阵，A[i][j] 表示从状态 i 转移到状态 j 的概率	self.A	N×N
B	发射矩阵，B[i][k] 表示状态 i 下观测到符号 k 的概率	self.B	N×M
pi	初始状态概率向量，pi[i] 表示初始时刻处于状态 i 的概率	self.pi	N
O	观测序列	obs_seq	T（序列长度）

顺手记住一个约束：A 的每一行求和等于 1，B 的每一行求和等于 1，pi 求和等于 1，因为它们是概率分布。写代码初始化的时候这个约束经常被忘掉，后面训练出来的参数看起来会非常离谱。还有一点，状态和观测的编号都是从 0 开始的下标，跟 Python 的数组索引天然对齐，这也是我推荐用 numpy 实现的原因之一。

2. 三大经典问题：HMM 的所有代码都在回答这三个问题

2.1 概率计算问题：这个序列出现的概率是多少

第一个问题是：给定模型参数 λ=(A, B, pi)，一个观测序列 O 出现的概率 P(O|λ) 是多少。

最朴素的做法是把所有可能的状态序列都枚举出来，把每条状态序列对应的转移概率和发射概率全部乘起来再求和。状态序列一共有 N^T 种，T 稍微大一点就爆炸了。比如 10 个状态、100 个时刻，N^T = 10^100，这个数字比宇宙中的原子数量还大。所以必须用动态规划。

前向算法的核心思想是定义一个前向概率 alpha[t][j]：在 t 时刻处于状态 j，并且已经观测到 O[0..t] 的概率。从 t=0 开始往后递推，每一步只复用上一步的结果。时间复杂度 O(T·N²)，空间复杂度 O(TN)。这背后就是动态规划里"记住中间结果、避免重复计算"的思路，跟斐波那契数列迭代优化是同一个套路。掌握这个算法，后面所有 HMM 变体——半隐马尔可夫、分层 HMM、隐半马尔可夫——都逃不出这个框架。

2.2 解码问题：最可能的隐藏状态序列是什么

第二个问题是：已知观测序列 O，求最可能产生它的状态序列 I*。注意这个问题和第一个问题的区别。第一个问题是"这串观测总体上有多大概率出现"；第二个问题是"如果它出现了，背后最有可能是哪条状态路径"。

Viterbi 算法就是专门干这个的动态规划。它定义 delta[t][j]：在 t 时刻到达状态 j，并且走的是当前所有路径里概率最大的那条路径的路径概率。同时用 psi[t][j] 记录"上一个时刻是从哪个状态转移过来的"，方便最后回溯。复杂度同样是 O(TN²)。

Viterbi 和前向算法有个微妙差别：前向用求和累加所有路径的概率，Viterbi 用取最大值只保留最有可能的那条路径。我第一次写的时候就把这里的 max 写成了 sum，结果输出一个"平均化"的伪路径，看起来每个状态都像，但整条路线完全不符合真实转移结构。这个区别是理解两个算法的关键。求和对应"边缘概率"，取 max 对应"最大后验路径"，语义上完全不同。

2.3 学习问题：模型参数从哪来

第三个问题是：给定一批观测序列，怎么估计 A、B、pi 这三个参数。

如果训练数据里每个观测序列都附带对应的状态序列（有标注），那很简单：直接统计状态转移频次、发射频次、初始状态频次，再归一化就行，本质上就是极大似然估计。

但大多数场景下没有标注，只有观测序列。这时候要用 Baum-Welch 算法，它是 EM 算法的一个特例。思路是：先随便给一组参数，用这组参数去算"在给定观测下，每个时刻处于某个状态的概率"（E 步），然后用这些软标签重新统计参数（M 步），再拿新参数去算，反复迭代直到收敛。这个"先猜再修正"的过程，和 K-Means 的迭代思路很像，只不过 K-Means 是硬分配，Baum-Welch 是软分配——每个时刻的状态不是非此即彼，而是带着一个概率权重分配到所有状态上。

三个问题对应三块代码，下面两章我把它们逐一写出来。

3. 手写代码：前向算法与 Viterbi 解码逐行拆解

3.1 先搭一个 HMM 的类骨架

我用 numpy 实现，因为矩阵运算是 HMM 的天然表达方式。先定义模型骨架，参数先用均匀分布初始化。

import numpy as np class HMM: def __init__(self, num_states, num_obs): self.N = num_states # 状态数 self.M = num_obs # 观测数 # 均匀初始化，保证每行和为 1 self.A = np.full((self.N, self.N), 1.0 / self.N) self.B = np.full((self.N, self.M), 1.0 / self.M) self.pi = np.full(self.N, 1.0 / self.N)

为什么初始化为均匀分布而不是全零？因为概率矩阵要求每行求和为 1，全零矩阵不满足概率约束；另外训练阶段如果初始值全零，梯度信号根本传不动，模型永远不会更新。第 4 章讲 Baum-Welch 时你会看到，初始化对收敛结果的影响有多大，这里先记住"均匀分布 + 一点噪声"是最稳妥的起点。

3.2 前向算法：把索引对齐再说

很多教材里 alpha 递推公式长这样：

alpha[t][j] = (sum_i alpha[t-1][i] * A[i][j]) * B[j][obs[t]]

把这个公式翻译成代码之前，先把维度对齐：alpha[t-1] 是一个长度为 N 的向量，乘上 A 的第 j 列得到 N 个数，求和后是"从任意状态转移到 j 的总概率"，再乘以 B[j][obs[t]]，就得到"到达 j 且发射出当前观测"的概率。

def forward(self, obs_seq): T = len(obs_seq) alpha = np.zeros((T, self.N)) # 初始时刻：先乘初始概率，再乘发射概率 alpha[0] = self.pi * self.B[:, obs_seq[0]] # 递推 for t in range(1, T): for j in range(self.N): alpha[t, j] = np.sum(alpha[t-1] * self.A[:, j]) * self.B[j, obs_seq[t]] return alpha def obs_prob(self, obs_seq): alpha = self.forward(obs_seq) return np.sum(alpha[-1])

最后返回的 P(O|λ)，就是把最后一个时刻的所有状态概率加起来。因为到了 T-1 时刻，不管落在哪个状态，只要走完了整个观测序列，都要计入总概率。

我拿一组具体参数验证这段代码。假设天气例子的参数为：

A = np.array([[0.7, 0.3], [0.4, 0.6]]) B = np.array([[0.6, 0.3, 0.1], [0.1, 0.4, 0.5]]) pi = np.array([0.6, 0.4])

观测序列取 [0, 1, 2]，即"散步、购物、打扫"。手算前向过程：

t=0：alpha[0] = [0.6×0.6, 0.4×0.1] = [0.36, 0.04]
t=1：alpha[1][0] = (0.36×0.7 + 0.04×0.4) × 0.3 = (0.252 + 0.016) × 0.3 = 0.0804
alpha[1][1] = (0.36×0.3 + 0.04×0.6) × 0.4 = (0.108 + 0.024) × 0.4 = 0.0528
t=2：alpha[2][0] = (0.0804×0.7 + 0.0528×0.4) × 0.1 = 0.0774 × 0.1 = 0.00774
alpha[2][1] = (0.0804×0.3 + 0.0528×0.6) × 0.5 = 0.0558 × 0.5 = 0.0279
P(O) = 0.00774 + 0.0279 = 0.03564

跑一遍代码，如果输出 0.03564，说明实现正确。这种手算校验是排查算法 bug 最有效的手段，比对着调试器干瞪眼快得多。索引有没有写反、矩阵方向有没有弄错，手算一次立刻暴露。

3.3 Viterbi 解码：动态规划加回溯指针

Viterbi 的递推公式是：

delta[t][j] = max_i (delta[t-1][i] * A[i][j]) * B[j][obs[t]]
psi[t][j] = argmax_i (delta[t-1][i] * A[i][j])

def viterbi(self, obs_seq): T = len(obs_seq) delta = np.zeros((T, self.N)) psi = np.zeros((T, self.N), dtype=int) delta[0] = self.pi * self.B[:, obs_seq[0]] for t in range(1, T): for j in range(self.N): probs = delta[t-1] * self.A[:, j] psi[t, j] = np.argmax(probs) delta[t, j] = probs[psi[t, j]] * self.B[j, obs_seq[t]] # 从最后时刻往前回溯，得到完整状态路径 path = np.zeros(T, dtype=int) path[-1] = np.argmax(delta[-1]) for t in range(T-2, -1, -1): path[t] = psi[t+1, path[t+1]] return path

psi 数组就是回溯指针。从最后一个时刻最大概率对应的状态开始，向前倒推每一步是从哪个状态来的。一开始我经常忘掉回溯，只返回 delta 的最大值，结果只拿到一个"最后时刻的状态"而没有整条路径。Viterbi 的完整输出一定包含回溯过程，这点要格外注意，面试里也经常拿这个点考人。

还是用前面那组参数，观测序列 [0, 1, 2]。手算：

t=0：delta[0] = [0.36, 0.04]
t=1：j=0，probs = [0.36×0.7, 0.04×0.4] = [0.252, 0.016]，psi=0，delta=0.252×0.3=0.0756
t=1：j=1，probs = [0.36×0.3, 0.04×0.6] = [0.108, 0.024]，psi=0，delta=0.108×0.4=0.0432
t=2：j=0，probs = [0.0756×0.7, 0.0432×0.4] = [0.05292, 0.01728]，psi=0，delta=0.005292
t=2：j=1，probs = [0.0756×0.3, 0.0432×0.6] = [0.02268, 0.02592]，psi=1，delta=0.01296

回溯：最后最大的是状态 1（雨天），psi[2][1]=1，所以 t=1 是状态 1；psi[1][1]=0，所以 t=0 是状态 0（晴天）。最优路径为 [晴, 雨, 雨]。

有意思的是，前向算法给出整体概率 0.03564，而 Viterbi 告诉我们最可能的天气路径是"晴→雨→雨"。实际项目中这两个结果经常一起看：一个回答"这个序列正常吗"，一个回答"如果正常，背后发生了什么"。

3.4 把两层循环改成矩阵运算

上面的嵌套循环直观，但效率低。numpy 的向量化改写很简单：

def forward_vec(self, obs_seq): T = len(obs_seq) alpha = np.zeros((T, self.N)) alpha[0] = self.pi * self.B[:, obs_seq[0]] for t in range(1, T): alpha[t] = (alpha[t-1] @ self.A) * self.B[:, obs_seq[t]] return alpha

alpha[t-1] @ self.A 的结果是一个长度为 N 的向量，第 j 个元素就是 sum_i(alpha[t-1][i] * A[i][j])。这一步替换掉了最内层的 j 循环。同理，Viterbi 的 delta 也可以用向量化写，但 psi 记录 argmax 的时候需要一点技巧：先广播算出所有 probs，再用 argmax 沿行方向取值。

什么时候用循环版本，什么时候用向量化版本？我的建议是：学习阶段先写循环版，逻辑清楚、好调试；工程阶段换向量化版。两个版本跑相同输入比对结果，还能顺便验证实现有没有写错。这一点值得反复强调，我见过很多人在一个函数上直接上向量化，结果出了 bug 半天定位不到。

3.5 工程必做：在 log 域里算

前面所有代码都在概率域（0~1）里乘来乘去。观测序列一长，比如语音识别里几十上百帧，alpha 的值会指数级缩小，很快就下溢成 0.0。一旦所有概率都变成 0，模型就完全失去区分能力，训练时 log-likelihood 变成负无穷，什么信息都学不到。

解决办法是在对数域计算，乘法变成加法。Viterbi 里本来就全是乘法和 max，log 化非常自然：

def viterbi_log(self, obs_seq): T = len(obs_seq) logA = np.log(self.A) logB = np.log(self.B) log_pi = np.log(self.pi) delta = np.zeros((T, self.N)) psi = np.zeros((T, self.N), dtype=int) delta[0] = log_pi + logB[:, obs_seq[0]] for t in range(1, T): for j in range(self.N): probs = delta[t-1] + logA[:, j] psi[t, j] = np.argmax(probs) delta[t, j] = probs[psi[t, j]] + logB[j, obs_seq[t]] path = np.zeros(T, dtype=int) path[-1] = np.argmax(delta[-1]) for t in range(T-2, -1, -1): path[t] = psi[t+1, path[t+1]] return path

注意 log 版本返回的是路径。如果想要概率值，得把 delta 的最后结果做指数还原并归一化，但工程里通常只需要路径，因为 log 域的数值大小关系跟原始概率完全一致，max 的结果不变。

前向算法在 log 域里麻烦一点，因为里面有求和。不能直接给整体取 log，需要用 logsumexp 技巧。scipy 提供了现成函数，自己实现也就三行：

def logsumexp(x): m = np.max(x) return m + np.log(np.sum(np.exp(x - m)))

原理是把最大的数提出来，让剩余部分落在安全的指数范围内，避免 exp 溢出。没有这一步，长序列上照样会出 -inf。

4. 参数训练：Baum-Welch 的完整实现

4.1 有标注数据：直接统计就行

先讲最简单的。如果我有 100 条带状态标注的序列，每条序列是 (观测串, 对应的状态串)，那么参数估计就是数数然后归一化：

pi[i] = 以状态 i 开头的序列数 / 总序列数
A[i][j] = 从 i 转移到 j 的次数 / 从 i 出发的总转移次数
B[i][k] = 状态 i 下观测到 k 的次数 / 状态 i 出现的总次数

这种极大似然估计一次搞定，不需要迭代。在数据量够大、标注可靠时效果最好。代码就是几个 Counter 的事，我就不贴了。实际项目里如果有少量标注数据，我通常会先统计出初始参数，再用 Baum-Welch 在更大的无标注数据上微调，这样比纯随机初始化收敛快得多，结果也更稳。

4.2 没有标注：先猜一组参数，再逐步修正

现实场景里标注很贵，更多时候手上只有一堆观测序列。这时候 Baum-Welch 出场。

算法的每一步需要两组中间量。第一组是 gamma[t][i]，表示在给定整个观测序列和当前参数下，t 时刻处于状态 i 的概率。第二组是 xi[t][i][j]，表示 t 时刻在状态 i、t+1 时刻在状态 j 的联合概率。这两个量都要靠前向算法和后向算法配合计算。

后向算法跟前向算法是对称的。定义 beta[t][i]：在 t 时刻处于状态 i 的条件下，从 t+1 到最后生成剩余观测的概率。初始化时 beta[T-1] = 1，然后从后往前递推。

def backward(self, obs_seq): T = len(obs_seq) beta = np.zeros((T, self.N)) beta[-1] = 1.0 # 最后一个时刻向后看，没有观测要生成 for t in range(T-2, -1, -1): for i in range(self.N): beta[t, i] = np.sum(self.A[i, :] * self.B[:, obs_seq[t+1]] * beta[t+1]) return beta

为什么 beta[-1] 初始化为 1 而不是 0？因为 T-1 时刻之后没有观测了，"生成空序列"的概率天然是 1。如果你把它初始化为 0，整个递推就全变成 0。这是新手最容易踩的坑之一，还有一版经常写反方向，把 beta 从前往后推，结果语义全乱，排查半天。

有了 alpha 和 beta，gamma 和 xi 就能算出来：

def compute_gamma_xi(self, obs_seq): T = len(obs_seq) alpha = self.forward(obs_seq) beta = self.backward(obs_seq) prob = np.sum(alpha[-1]) # P(O|lambda) gamma = alpha * beta / prob # T×N 矩阵 xi = np.zeros((T-1, self.N, self.N)) for t in range(T-1): xi[t] = alpha[t][:, None] * self.A * self.B[:, obs_seq[t+1]] * beta[t+1] xi[t] /= prob return gamma, xi, prob

gamma[t][i] 的理解很直观：t 时刻处于状态 i 的概率。xi[t][i][j] 更细一层：t 时刻在 i、t+1 时刻在 j 的联合概率。这两个量就是 E 步的"软计数"，所有参数更新都建立在它们上面。

4.3 M 步：把期望统计量变成新参数

更新公式长这样：

pi[i] = gamma[0][i]
A[i][j] = sum_t xi[t][i][j] / sum_t gamma[t][i]
B[i][k] = sum_{t: obs[t]=k} gamma[t][i] / sum_t gamma[t][i]

说白了：新参数 = 期望下的计数 / 期望下的归一化分母。E 步算期望，M 步更新参数，反复迭代直到参数变化量小于阈值。

def baum_welch(self, obs_seqs, max_iter=100, tol=1e-6): for it in range(max_iter): A_num = np.zeros((self.N, self.N)) A_den = np.zeros(self.N) B_num = np.zeros((self.N, self.M)) B_den = np.zeros(self.N) pi_num = np.zeros(self.N) for obs in obs_seqs: gamma, xi, prob = self.compute_gamma_xi(obs) pi_num += gamma[0] for t in range(len(obs) - 1): A_num += xi[t] A_den += gamma[t] for t in range(len(obs)): B_num[:, obs[t]] += gamma[t] B_den += gamma[t] new_A = A_num / A_den[:, None] new_B = B_num / B_den[:, None] new_pi = pi_num / len(obs_seqs) diff = np.max(np.abs(new_A - self.A)) self.A, self.B, self.pi = new_A, new_B, new_pi if diff < tol: print(f"第 {it+1} 轮迭代后收敛") break

这段代码有几个容易出错的地方，我说一下。

第一，pi_num 只累加每条序列的 gamma[0]，然后除以序列条数。因为 gamma[0] 本身已经归一化过，所有状态加起来等于 1，所以平均下来就是初始状态分布的期望。第二，A_den 累加的是第 0 到 T-2 时刻的 gamma，不是全部时刻。因为转移发生在相邻时刻之间，最后一个时刻没有向后的转移。如果你把所有时刻的 gamma 都加进去，分母会比分子大，更新出来的 A 会系统性偏小。第三，B_num 累加的是 gamma 在观测等于 k 的各个时刻的值，也就是把软计数按观测符号分类汇总。第四，收敛判断我只看了 A 的最大变化量，实际项目中建议同时看 A 和 B，把两者的最大绝对值差一并拿来判断，更稳妥。

4.4 初始化对 Baum-Welch 的影响

Baum-Welch 是 EM 算法，EM 保证收敛到局部最优，不保证全局最优。初始化参数不同，最后可能收敛到完全不同的局部最优，这是这类算法的宿命，不是 bug。我在实际项目里常用的策略有两种：第一种是随机初始化多次，比如随机跑 20 次，每次用不同随机种子，最后取似然值最高的那组参数；第二种是先用均匀分布初始化跑几个轮次，然后把接近 0 的概率加一个极小值，防止模型卡死在 0 概率状态上。

还有一个小技巧：如果观测序列很长，可以把每条序列切成几段来训练，降低单次计算量。HMM 的时间复杂度对序列长度是线性的，切分很安全，不会破坏模型的本质。另外要注意，如果初始转移矩阵里有 0，训练过程中这个 0 永远不可能被更新成正数，因为统计量是从当前参数算出来的，0 概率对应的 xi 期望永远是 0。所以初始化时用均匀分布加一点噪声，比用稀疏矩阵稳妥得多。

5. 实际工程里的坑与排查实录

5.1 概率下溢：算着算着全变成 0

前面反复提过，前向算法在长序列下会概率下溢。我排查过不少线上异常，现象就是 obs_prob 输出 0.0，或者训练时 log-likelihood 变成 -inf。解决办法就是 log 域 + logsumexp。Viterbi 换成 log 版本就行，前向算法需要写 logsumexp 版本：

def forward_log(self, obs_seq): T = len(obs_seq) logA = np.log(self.A) logB = np.log(self.B) log_pi = np.log(self.pi) log_alpha = np.zeros((T, self.N)) log_alpha[0] = log_pi + logB[:, obs_seq[0]] for t in range(1, T): for j in range(self.N): log_alpha[t, j] = logsumexp(log_alpha[t-1] + logA[:, j]) + logB[j, obs_seq[t]] return logsumexp(log_alpha[-1])

注意 logA、logB、log_pi 在迭代前算好缓存，不要在每步重复取对数。如果你对性能敏感，这条优化能省不少时间。同样的做法也适用于 Viterbi 的 log 版本。

5.2 零概率问题：平滑是必须的

训练数据有限时，发射矩阵里很容易出现 0。比如某个状态在训练集中从来没发射过某个观测符号，B[i][k] = 0。推理时一旦遇到这个观测，整条路径概率直接变成 0，前向概率也跟着变成 0，这种一刀切的结果非常不稳健。

常用方案是拉普拉斯平滑：统计计数时给每个格子加一个平滑项 δ。比如 B[i][k] = (计数[i][k] + δ) / (sum_k 计数[i][k] + δM)。δ 一般取 0.1 或 1。在 log 域里平滑还有一个额外好处：避免对 0 取 log 得到 -inf。对数域里 -inf 一参与加法，结果直接变 -inf，回溯时整条路径断掉，这个坑我踩过不止一次。所以不管训练数据多充分，我都会在参数里保留一点平滑项，宁可让模型稍微"糊"一点，也不要出现刚性 0。

5.3 收敛到局部最优：换个种子试试

Baum-Welch 对初始化敏感，这是实践里最容易被忽略的点。有时候训练日志显示似然值停在某个水平上不去，换一组随机初始参数，似然又能跳高一截，跟 K-Means 的处境一模一样。我一般会做随机重启：用 10 到 20 组不同随机种子初始化，各训练到收敛，然后拿验证集上的似然或错误率挑最好的模型。代价是训练时间乘以一个系数，但换来的是稳定性和最终效果，这笔账划算。

另外提醒一点：观察似然值曲线时，如果发现前期涨得很快后期剧烈震荡，可能是学习率（这里的步长）没有衰减机制。Baum-Welch 本身没有学习率概念，但如果观测数据量级差异很大，可以在更新时对分母做 floor，强制最小值，防止某个状态被某一条特别长的序列主导。

5.4 状态数 N 怎么选：模型选择的实际问题

状态数 N 是 HMM 里最让人头疼的超参数。N 太小，模型表达力不够，拟合不了数据里的复杂模式；N 太大，参数数量急速膨胀，训练数据不够时严重过拟合。参数总数是 N² + NM + N，对 N 是平方级增长。

我的经验做法是：先按业务含义定一个合理的初值，再用信息准则辅助判断。比如在给定 N 下训练完模型后，计算 BIC 或 AIC，然后对 N 取几个候选值画折线图选拐点。还有个更实用的办法：直接看验证集上的似然值。当 N+1 的模型在验证集上不再显著优于 N，说明 N 已经够用。过拟合有个典型特征：训练集似然很高，验证集似然反而下降。看到这个现象，优先怀疑 N 太大。

5.5 HMM 与 CRF：什么时候别用 HMM

我经常被问：做序列标注，什么时候用 HMM，什么时候用条件随机场（CRF）？

HMM 是生成式模型，建模联合概率 P(状态, 观测)，再通过贝叶斯做推断。CRF 是判别式模型，直接建模条件概率 P(状态|观测)。生成式模型的好处是可以用无标注数据训练，Baum-Welch 正好发挥这个优势；判别式模型的优势是特征工程灵活，CRF 里可以随便加前后词、词性、位置、模板等特征，而 HMM 的发射概率本质上只能依赖当前观测这一个维度。

所以选择原则很简单：如果训练数据有大量无标注样本可以配合一小批标注样本，或者数据本身的生成过程比较清晰，比如语音信号建模，优先考虑 HMM；如果任务以标注数据为主，关键又在于特征工程，比如中文分词、命名实体识别，直接上 CRF 或基于条件随机场的深度模型，别纠结。模型没有绝对优劣，只有场景适配问题。

下面是几个高频问题的速查表，方便出问题时对照：

现象	可能原因	排查方向
前向概率一直是 0	概率下溢 / B 中有 0	换 log 域；加平滑
Viterbi 路径不合理	发射矩阵没学好 / 初始参数差	检查 B 的训练统计量；换初始化
Baum-Welch 不收敛	似然函数震荡 / 局部最优	加大 max_iter；随机重启
训练似然很高但预测很差	过拟合	减小 N；加大数据量；加平滑
参数出现 NaN	分母为 0 / 除零	检查平滑项；检查序列是否为空