从零实现HMM：前向算法、Viterbi与Baum-Welch工程指南

隐马尔可夫模型hidden Markov modelHMM这个名字听起来像是概率论教材里最不好惹的那一章可如果你写过语音识别、做过中文分词、跑过基因序列标注回头一看会发现早就跟它打过照面了。它的核心设定朴素得有点反直觉系统在背后按一条你根本看不见的状态链悄悄演化每个时刻只吐给你一个观测值你要靠这串观测值去反推状态是什么、参数该怎么估。这种看得见结果、看不见原因的结构在现实里比比皆是——用户点没点广告你能看到他心里想买什么你看不到设备报出来的温度你能读到内部是正常运转还是已经开始劣化你不知道。我写这篇的动机很简单。市面上讲 HMM 的材料要么停在公式推导层面三个问题列出来、公式贴上去就结束了读者抄完公式不知道代码该长什么样要么直接甩一句调 hmmlearn 就行等到真出问题——概率全部下溢成 0、模型收敛到一个明显不合理的结果、状态数怎么选——完全无从下手。我打算把这两头接上从为什么这个算法必须这么设计讲到我的实际代码长什么样中间那些只有踩过坑才知道的细节一个都不省。适合谁看如果你已经能看懂基本的条件概率和矩阵乘法想真正把 HMM 用起来这篇是给你写的。如果你只是想快速回顾一下公式也可以只看第 1 节和第 7 节一个帮你建立直觉一个帮你判断这东西今天还值不值得用。1. 从看不见的状态说起HMM 要解决的到底是什么问题1.1 两个骰子的赌局先建立对隐藏状态的物理直觉抛开公式先想象一个场景。桌上有两个骰子一个是正常骰子六个面概率均等另一个被做过手脚掷出 6 的概率是 0.5其余五个面平分剩下的 0.5。现在有人每次随机挑一个骰子掷一次记下点数然后换不换骰子你也不知道——但你只能看到一长串点数6 6 1 6 6 6 2 6 ...。问题来了这串点数背后每一轮到底用的是哪个骰子你脑子里会自然冒出两个判断。第一连续看到好几个 6那大概率是那个动了手脚的骰子。第二骰子的切换本身也有惯性一个人不太可能每一轮都换着掷上一轮用了哪个这一轮继续用的概率应该更高。这两个判断恰好就是 HMM 的两根支柱发射概率某个状态下观测值出现的可能性和转移概率状态与状态之间的切换倾向。再换个更贴近应用的例子。假设你有个朋友住在很远的地方每天只能在社交平台上看到他发了什么——散步购物打扫房间。你想知道他那边天气怎么样但天气预报你拿不到。这里晴/雨是隐藏状态散步/购物/打扫是你能看到的观测。晴天的日子里他更可能去散步雨天的日子里他更可能待在家打扫。你手上有一串他的活动记录想反推这几天的天气序列。这就是 HMM 最典型的用法观测序列已知状态序列未知。我说这两个例子的意思是想强调一点HMM 不是先有数学再有应用而是先有隐藏原因驱动可见结果这类问题然后才被人用概率语言写下来的。很多人学 HMM 卡住就是卡在直接看公式脑子里没有一个具体场景在跑。你先把骰子或者天气这个画面刻进脑子后面所有的 α、β、γ 都只是给这个画面配上的符号而已。1.2 五个要素各自管什么π、A、B 的分工不能混一个标准的 HMM 由五样东西确定通常记成一个元组 λ (N, M, π, A, B)。N 是隐藏状态的个数M 是观测值的个数离散情形下就是词典大小。剩下三个才是真正需要估的参数也是新手最容易混的地方。参数名称形状物理含义约束π初始状态分布长度 N 的向量序列第一个时刻处于各状态的概率非负和为 1A状态转移矩阵N × N 矩阵从状态 i 跳到状态 j 的概率 a_ij每行非负行和为 1B发射概率矩阵N × M 矩阵状态 i 下观测到符号 k 的概率 b_i(k)每行非负行和为 1这张表我建议你对着代码多看几遍因为九成的实现 bug 都出在这三个矩阵的行列方向上。A 是行是当前状态、列是下一状态也就是说A[i][j]表示 i → jB 是行是状态、列是观测符号B[i][k]表示状态 i 生成符号 k。一旦你把 A 写成了转置前向算法照样能跑出数字但那个数字是错的而且不会报错这才是最要命的。另外注意行和为 1 这个约束。π 是一个概率分布A 的每一行、B 的每一行都是一个概率分布。你在做参数更新比如 Baum-Welch或者手工调参的时候如果破坏了行归一化整个模型就失去了概率解释算出来的似然再也无法比较大小。我在代码里习惯在每次参数更新后强制加一步归一化宁可多一次除法也不留下脏参数。状态数 N 的选择是个独立话题第 5 节会专门讲。这里先记住一个原则N 是你对系统内部有几种模式的先验假设它不是一个能靠数据自动确定的量只能靠实验去试。1.3 三大问题与算法对应别把它们记成三个孤立知识点教科书里 HMM 一定会列出三个基本问题很多人把它们当作三个独立的考点背下来其实它们是同一条流水线上的三道工序。问题一句话描述经典算法时间复杂度评估问题给定 λ 和观测序列 O算 P(O | λ)前向算法或前向-后向O(N²T)解码问题给定 λ 和 O找最可能的状态序列 Q*Viterbi 算法O(N²T)学习问题只给 O反推最合适的 λBaum-WelchEM 的一种特例O(N²T) 每轮它们的关系是这样的学习问题内部要反复调用评估问题来算似然和新旧参数而解码问题在训练完之后才用得上。所以真正动手实现的时候顺序一定是先写前向、再写 Viterbi、最后写 Baum-Welch。反过来写会非常痛苦因为 Baum-Welch 的正确性完全依赖前向和后向的中间结果而这两个东西写错了训练出来的参数虽然能跑但毫无意义。还有个容易被忽略的点这三个问题的复杂度都是 O(N²T)而不是 O(N^T)。这个差距是决定性的也是第 2 节要展开的核心。简单说如果复杂度是 N 的 T 次方那 HMM 根本不可能被用在真实数据上一条一百帧的语音就足以让任何计算机算到天荒地老。2. 前向算法为什么不能暴力枚举一次复杂度上的生死抉择2.1 先算算暴力枚举到底有多离谱要算 P(O | λ)最直觉的做法是把所有可能的状态序列都列出来对每一条序列算这条序列产生这串观测的概率最后加总。公式写出来很干净P(O | λ) Σ_Q π_{q1} b_{q1}(o_1) a_{q1q2} b_{q2}(o_2) ... a_{q_{T-1}q_T} b_{qT}(o_T)问题出在那个 Σ_Q 上。Q 是长度为 T 的状态序列每个位置有 N 种取值所以一共有 N^T 条序列。每条序列要算 2T 次乘法。总的计算量是 O(T · N^T)。我们代入具体数字感受一下。假设 N 3一个很小的模型T 100一百帧语音这在语音识别里算极短的。3 的 100 次方大约是 5×10^47。就算你的机器每秒能算十亿次乘法也需要 10^31 秒。宇宙年龄才 4×10^17 秒。这个差距不是靠优化常数或者换硬件能弥补的它是量级上的不可能。提示这个指数爆炸是所有序列模型都要面对的第一道墙。前向算法、Viterbi、CTC 的前向后向本质上都是在用动态规划把这堵墙推倒把指数复杂度压成多项式复杂度。2.2 前向变量的递推结构把公共子问题缓存下来关键洞察是不同路径之间共享了大量前缀。比如两条状态序列1→2→1→3和1→2→1→2它们前三个位置完全一样产生前三个观测的概率也一样只是最后一个位置分岔了。暴力枚举会把这段前缀算两遍而动态规划只算一遍。于是定义前向变量α_t(i) P(o_1, o_2, ..., o_t, q_t i | λ)读法是到第 t 个时刻为止观测序列是 o_1 到 o_t并且第 t 个时刻恰好处于状态 i 的联合概率。注意这是联合概率不是条件概率这一点很多资料讲得不清楚导致后面跟 β 相乘的时候有人想不通为什么可以直接乘。递推关系来自一个很自然的分解α_t(j) [ Σ_i α_{t-1}(i) · a_ij ] · b_j(o_t)翻译成人话要想到达第 t 刻在状态 j 且看到了 o_t这个局面可以从上一刻的任意状态 i 出发跳过来概率是 α_{t-1}(i) · a_ij把所有 i 的可能性加起来再乘上在状态 j 下吐出观测 o_t的概率。初始条件是 α_1(i) π_i · b_i(o_1)终止时 P(O | λ) Σ_i α_T(i)。复杂度一下就降到了 O(N²T)。还是 N 3、T 100计算量大概是 900 次乘法眨眼就完事了。从 10^47 到 900这就是动态规划在序列建模里的分量。2.3 数值下溢浮点数会在你毫无察觉的时候把结果抹成 0上面那段推导数学上完全正确但直接照着写代码跑到 T 稍微大一点就会出事而且是那种结果直接变成 0.0的静默错误。原因是你一直在做连乘。α_1(i) 是个小于 1 的数每往前推一步就要乘一堆小于 1 的概率数值单调递减。double 类型能表示的最小正规格化数大约在 1e-308 量级一旦乘积掉到这个线以下就成了非规格化数再掉下去直接变成 0。粗略估算如果每一步平均乘一个 0.1 量级的因子那么大约 T 308 步之后α 就全线归零了。后果不只是结果不准而是整个后续计算全部失效。似然变成 0 意味着对数似然变成负无穷Baum-Welch 里 α · β 全是 0γ 做归一化的时候是 0/0得到 NaNNaN 一旦进入参数矩阵后面所有迭代都是 NaN。我见过不止一个项目在这上面浪费了一整天以为是数据有问题实际就是没做数值处理。注意只要你的序列长度可能超过几十就一定要处理下溢。这不是优化项是必选项。2.4 缩放因子与 log 域两条路线的取舍主流解法有两条。第一条是缩放因子法scaling也是教科书上最常推的做法。思路是每一步算完 α_t(i) 之后除以当前时刻所有状态的和 c_tc_t 1 / Σ_i α_t(i)α̂_t(i) α_t(i) · c_t因为做了归一化α̂_t(i) 的数值始终在 0 到 1 之间而且每时刻所有状态的和恰好是 1永远不会下溢。同时这个过程可逆原始的似然满足 P(O | λ) 1 / Π_t c_t取对数就是 log P(O | λ) -Σ_t log c_t。所以缩放不但解决了下溢还免费给了你一个数值稳定的对数似然。第二条是全程对数域log domain。把所有概率取对数乘法变加法加法用 log-sum-exp 处理。这条路的好处是不会下溢负数域可以到 -700 以下坏处是要自己实现 log-sum-exp而且不能再用 numpy 的矩阵乘法直接操作——log(exp(a) exp(b))每次都要手工展开代码会变啰嗦。方案优点缺点适用场景缩放因子法可复用矩阵乘法代码简洁顺带得到对数似然概率本身接近 0 时缩放因子会很大绝大多数离散 HMM 实现对数域数值最稳无下溢风险需要 log-sum-exp向量化困难序列超长、或要做混合精度直接计算代码最少必下溢只适合 T 20 的演示我个人的习惯是缩放法为主。Viterbi 那部分则直接全程对数域因为它只涉及加法和 max不需要 log-sum-exp天然适合对数域。这两条路线在同一份代码里并存其实非常自然后文第 6 节会给出完整实现。3. Viterbi 解码从观测序列反推最可能的状态轨迹3.1 动态规划表里每一项代表什么评估问题问的是这串观测整体出现的概率有多大解码问题问的是最可能是哪条状态链产生了它。这两个问题的目标不同算法结构却惊人地相似只是把求和换成了取最大。Viterbi 定义了一个新变量δ_t(j) max over q_1..q_{t-1} of P(q_1, ..., q_{t-1}, q_t j, o_1, ..., o_t | λ)注意这里取的是路径上的最大值而不是所有路径的总和。递推式写成δ_t(j) [ max_i δ_{t-1}(i) · a_ij ] · b_j(o_t)直觉上很好理解我要找一条最优路径走到第 t 刻在状态 j那么它的前缀一定是走到第 t-1 刻某个状态 i的最优路径。这是最优子结构也是动态规划能成立的前提。这里有一个必须点破的细节Viterbi 求出来的路径不等于每个时刻单独取最可能状态拼起来的那条路径。前者是所有时刻联合起来最可能的整条路径后者是逐点最大化边缘概率。这两个结果经常会不一样。当年我第一次接触这个区别时也觉得别扭但仔细想想就通了逐点最优拼起来的那条路径很可能包含了一个概率极低的转移比如从状态 A 一步跳到状态 D而这种跳转在实际模型里几乎不会发生。联合最优则会绕开它。3.2 回溯指针别忘了记下从哪来Viterbi 的实现必须维护两张表一张存 δ 值一张存指针 ψ。ψ_t(j) 记录的是在第 t 刻到达状态 j 的那条最优路径上一刻是从哪个状态来的。为什么必须存指针因为动态规划是自底向上填表的填到最后一刻你只知道终点在哪不知道路径怎么走。没有指针的话你得从头重新推一遍那就退化成指数搜索了。存了指针最后从 δ_T 的最大值位置开始一步步顺着 ψ 往回跳就能在 O(T) 时间里还原完整路径。我见过有人写 Viterbi 只在最后对 δ_T 取了个 argmax 就当成答案输出那其实只输出了最后一个时刻的状态前面全是错的。这个 bug 特别隐蔽因为输出长度对不上通常会被立刻发现但如果他恰好在最后加了个 padding 或者只关心最后一个状态就完全看不出来了。3.3 一段可以直接跑的 Viterbi 实现import numpy as np def viterbi(obs, pi, A, B): obs: 观测序列元素为 0..M-1 的整数 pi : (N,) 初始状态分布 A : (N, N) 转移矩阵A[i, j] P(i - j) B : (N, M) 发射矩阵B[i, k] P(观测 k | 状态 i) 返回: (最优状态路径, 该路径的对数概率) T len(obs) N len(pi) EPS 1e-300 log_pi np.log(pi EPS) log_A np.log(A EPS) log_B np.log(B EPS) delta np.full((T, N), -np.inf) psi np.zeros((T, N), dtypeint) # 初始化 delta[0] log_pi log_B[:, obs[0]] # 递推 for t in range(1, T): for j in range(N): # 上一刻所有状态转移到 j 的候选值 candidates delta[t - 1] log_A[:, j] best_i int(np.argmax(candidates)) psi[t, j] best_i delta[t, j] candidates[best_i] log_B[j, obs[t]] # 回溯 path np.zeros(T, dtypeint) path[T - 1] int(np.argmax(delta[T - 1])) for t in range(T - 2, -1, -1): path[t] psi[t 1, path[t 1]] return path, delta[T - 1].max()代码里有几个地方值得单独说一下。第一加EPS再取对数。如果 A 或 B 里有精确的 0np.log(0) 会返回 -inf 并抛出运行时警告。加上一个极小的正数就把 -inf 变成了一个很大的负数约 -690这在后续做加法时行为正确不会污染结果。这也是工程实践和教科书推导的差别数学上 0 就是 0工程上你得给它一个地板值。第二delta初始化为 -inf 而不是 0。如果初始化成 0那些理论上不可能的状态在第一天就会有个 0 的对数概率比真实可能状态的负数值还大argmax 会直接选错。用一个明确的负无穷表示不可达语义才干净。第三内层循环里的delta[t-1] log_A[:, j]是向量化写法比再套一层 Python 循环快得多。别看 N 小的时候差别不明显当 N 到几十、T 到几千的时候这个改写能带来十倍以上的速度差异。3.4 Viterbi 与前向算法只差一个运算符把两个算法并排放在一起看你会发现它们的骨架几乎完全一致。对比项前向算法Viterbi 算法中间量语义前缀联合概率之和前缀路径的最大概率递推核心Σ_i α_{t-1}(i)·a_ijmax_i δ_{t-1}(i)·a_ij是否需要指针不需要需要 ψ 用于回溯数值处理缩放因子或对数域天然适合对数域输出标量 P(O|λ)状态序列 路径得分这个对比给我的启发是很多看起来不同的算法骨子里是同一个动态规划框架只是聚合算子不同——求和得到的是整体可能性取最大得到的是最优解。同样的思路你在 CTC 的前向后向、条件随机场的维特比解码里都会再遇到一次。理解了这一层再看那些算法就不会觉得是一个个孤立的黑盒了。顺带提一个工程上的细节Viterbi 的输出得分不是概率而是对数域上的路径得分。你没法拿它和 P(O|λ) 直接比大小也不应该对它做 exp 之后再解释成这条路径的概率。它是最优路径的得分仅此而已。我见过有人拿这个分数去做阈值判断结果阈值怎么调都不对就是因为混淆了这两个量。4. Baum-Welch 训练没标注数据时参数是怎么学出来的4.1 EM 的套路先猜状态再更新参数来回拉锯前面两节都假设 λ (π, A, B) 已知。可现实里最常见的情况是你只有一堆观测序列一个标注都没有。这时候要估参数靠的就是 Baum-Welch 算法。它的本质是 EM期望最大化在 HMM 上的具体形式。EM 的思路可以粗暴地概括成一句话如果我大概知道状态序列参数就好估如果我大概知道参数状态序列就好猜。那就先随便给个参数用参数猜状态再用猜出来的状态更新参数反复来回。具体到 HMM这个猜是通过计算两个后验量来完成的。γ_t(i) P(q_t i | O, λ)第 t 刻处于状态 i 的后验概率ξ_t(i, j) P(q_t i, q_{t1} j | O, λ)第 t 刻在 i、下一刻在 j 的后验概率这两个量都可以从前向变量和后向变量直接算出来。后向变量 β_t(i) P(o_{t1}, ..., o_T | q_t i, λ)含义是给定第 t 刻在状态 i后面那串观测出现的概率。有了它俩γ_t(i) α_t(i)·β_t(i) / P(O|λ) ξ_t(i, j) α_t(i)·a_ij·b_j(o_{t1})·β_{t1}(j) / P(O|λ)这里就能看出为什么前向和后向必须都用缩放后的版本分子分母里的缩放因子会精确地相互抵消而 P(O|λ) 用缩放因子算出来的对数似然对应上整套推导是一致的。4.2 重估公式背后其实是很朴素的数数拿到 γ 和 ξ 之后参数更新公式长这样π̂_i γ_1(i)â_ij Σ_{t1}^{T-1} ξ_t(i, j) / Σ_{t1}^{T-1} γ_t(i)b̂_i(k) Σ_{t: o_t k} γ_t(i) / Σ_{t1}^{T} γ_t(i)别看符号多翻译成人话非常直观。â_ij 的分子是整个序列里从 i 跳到 j 的期望次数分母是整个序列里待在 i 的期望次数两者一除就是待在 i 的时候往 j 走的比例。这和统计学里最原始的频率估计完全是一个路子只不过因为状态看不见我们用的是期望次数而不是实际次数。b̂_i(k) 同理分子是期望在状态 i 并且观测到符号 k 的次数分母是期望待在状态 i 的总次数。π̂_i 直接取第一时刻的状态后验。我在纸上推第一遍的时候总觉得用估计出来的 γ 去更新参数这件事有点循环论证的味道。后来想明白了EM 每一步都保证似然不下降这是它有收敛性保证的原因所以这个循环是朝上走的不会越走越差。代价是它可能停在局部最优这就引出了下一小节。4.3 初始化决定你能不能跳出局部最优Baum-Welch 是局部优化算法初始值给得不好结果就会很差。几种常见做法随机初始化最省事但结果不稳定。同一个模型跑两次可能一次收敛到似然 -320另一次 -410。如果你在做实验对比这会让你怀疑人生。我的做法是固定随机种子并且至少跑三到五次取最好的那个。均匀初始化即 π 和 A 全填 1/N、B 全填 1/M然后在上面加一点小扰动。不加扰动的话所有状态完全对称梯度方向可能也一样导致多个状态收敛到同一个解等于白设了那么多个状态。这是个很容易忽略的坑你以为设了 5 个状态训练完发现 5 行 B 几乎一模一样。基于先验知识的初始化在我做过的项目里效果最好。比如做设备状态监测我知道大概有正常、亚健康、异常三类那我就可以让 π 和 A 偏向正常状态持续时间长、异常状态一旦进入不太容易出来这种结构。把领域知识编码进初始值往往比多跑几十轮迭代更有用。还有一个技巧值得单独提参数绑定与结构化约束。如果你的应用里能确信某些转移概率应该为 0比如状态只能沿着 1→2→3 单向演进那就把对应的初始化值设成 0 并加地板值保护。这相当于给模型加了个强先验能显著减少要估的参数数量在小数据上尤其有用。4.4 收敛判据什么时候该停手标准的停止条件是似然增量小于阈值比如连续两轮的对数似然变化小于 1e-4 就停。但我实际用下来发现光看似然不够还得盯两个东西。第一是迭代次数上限。EM 在接近收敛的时候会变得非常慢一轮只涨 0.0001 的似然也能磨几十轮。设一个比如 100 轮的上限能省下大量时间。第二是过拟合。如果训练数据很少而状态数很多似然会一直涨但那是在记忆训练数据。判断方法很简单留一小部分数据做验证集看验证集似然什么时候开始掉。掉头的位置才是真正该停的地方。注意对数似然单调不降是 EM 在完全数据下的理论性质。如果你发现似然在下降那说明实现里有 bug最常见的是缩放因子没处理好、γ 的归一化用了错误的轴、或者前后向的初始化边界写错了。5. 把模型落到工程里观测表示、平滑与状态数5.1 连续观测离散化的代价前面所有讨论都假设观测是离散符号。但现实数据大多是连续的语音的 MFCC 特征、传感器的温度读数、金融的收益率。最省事的做法是先做向量量化VQ用一个聚类算法把连续向量映射成若干个离散码字然后套用离散 HMM。这么做的问题很明显量化误差是硬损失。两个相差很小的向量可能被分到不同码字相差很大的也可能因为落在聚类边界上被分到同一个。而且码本大小 M 又成了一个要调的超参数M 太小丢信息M 太大则 B 矩阵参数量暴涨N × M在小数据上直接过拟合。我一般的判断标准是如果特征维度低比如 1 到 3 维、数据量中等VQ 离散 HMM 够用且实现简单如果特征维度高十几维以上或者对精度有要求就别省这个事直接上连续观测模型。5.2 GMM-HMM语音识别时代的经典配置连续观测下最常见的选择是每个状态配一个高斯混合模型也就是所谓的 GMM-HMM。此时发射概率不再是查表而是b_i(o) Σ_{m1}^{K} w_{im} · N(o | μ_{im}, Σ_{im})每个状态 i 有 K 个高斯分量每个分量有自己的权重、均值和协方差。参数量比离散情形大得多Baum-Welch 的更新公式也要相应扩展到每个观测点对每个高斯分量各有一个后验责任度。这套结构在上世纪九十年代到 2010 年代初是语音识别的主流直到被深度神经网络取代。如果你想找一份完整、严谨的 GMM-HMM 推导那部分内容量足够单独写一篇这里只提两个我踩过的坑一是协方差矩阵的对角化假设。为了减少参数量实践中几乎都用对角协方差。这意味着你假设特征各维之间独立。如果特征之间有强相关比如你直接塞了原始 MFCC 而没有做 DCT 去相关这个假设会显著伤害性能。二是方差地板。某些高斯分量的方差会在训练中越缩越小最后接近 0导致那个分量的似然在遇到稍微偏离均值的点时就变成天文数字把整个后验压成独热向量。做法是给方差设一个下限比如特征全局方差的 1% 或者一个固定的小常数。5.3 零概率与平滑地板值该给多少离散 HMM 里如果某个观测符号在训练数据中从未在某个状态下出现过那 b_i(k) 的估计值就是 0。这会带来两个后果一是取对数变 -inf二是如果测试时真出现了这个符号整条路径的得分直接变成负无穷模型完全无法处理。标准解法是拉普拉斯平滑分子加 1、分母加 Mb̂_i(k) (count_i(k) 1) / (Σ_k count_i(k) M)如果 M 很大比如码本有几千个加 1 已经不够了更常用的是加一个更小的伪计数比如 0.01 或者 0.1。这个数怎么选我的经验是先给 0.1 跑一轮看未登录符号的处理是否正常再根据验证集表现微调。对 A 矩阵其实也需要类似处理只是问题没那么严重因为状态数通常比观测数少A 里出现零的概率低得多。但结构化约束人为把某些 a_ij 设为 0的时候就要小心了这些人为的 0 不应该被平滑掉需要在更新后重新置零。5.4 状态数怎么选一个没有标准答案的问题这是 HMM 实践里最玄学的一环。状态数不是从数据里学出来的是你自己定的。我的做法一般是这样数据规模建议状态数范围理由几百帧2 - 4参数少避免过拟合几千帧4 - 10能建模一定复杂度几万帧以上10 - 30数据足以支撑参数估计有领域知识按物理含义定比如正常/预警/故障就是 3核心判断方法是画似然-状态数曲线。状态数从 2 加到 10看训练集似然和验证集似然各怎么变。如果验证集似然在某个点之后不再上升甚至下降那个点就是上限。如果训练集似然一直涨、验证集不涨说明已经过拟合了该选更小的状态数。还有个更直观的诊断方法训练完之后把 B 矩阵每一行打印出来或者画成热力图看状态之间有没有明显的差异。如果好几个状态的发射分布几乎一样说明状态数给多了或者初始化没扰动好那几个状态是冗余的。这个检查我基本每次都会做比看数字直观得多。6. 实战从零写一个能用的 HMM 工具类6.1 数据结构与接口设计前面把三块算法分别讲完了这里把它们拼成一个完整的类。接口设计上我坚持两点一是参数全部用 numpy 数组传形状固定为 π(N,)、A(N,N)、B(N,M)二是所有对外方法都返回对数域的数值避免调用者拿到一个已经下溢成 0 的概率还不知道。import numpy as np EPS 1e-300 class DiscreteHMM: def __init__(self, n_states, n_obs, seed0): rng np.random.default_rng(seed) self.N n_states self.M n_obs # 均匀初始化 小扰动避免状态对称 self.pi np.full(n_states, 1.0 / n_states) self.A np.full((n_states, n_states), 1.0 / n_states) self.B np.full((n_states, n_obs), 1.0 / n_obs) self.A rng.normal(0, 0.01, self.A.shape) self.B rng.normal(0, 0.01, self.B.shape) self._normalize() def _normalize(self): self.pi np.clip(self.pi, EPS, None) self.pi / self.pi.sum() self.A np.clip(self.A, EPS, None) self.A / self.A.sum(axis1, keepdimsTrue) self.B np.clip(self.B, EPS, None) self.B / self.B.sum(axis1, keepdimsTrue)_normalize这一步是整套代码的安全网。不管前面的更新公式写了什么只要最后过一遍这个函数参数就一定满足概率约束。我建议在每个对外修改参数的方法末尾都调它一次成本极低收益极高。6.2 前向、后向与 Baum-Welch 的实现要点def forward(self, obs): T len(obs) alpha np.zeros((T, self.N)) c np.zeros(T) alpha[0] self.pi * self.B[:, obs[0]] c[0] 1.0 / (alpha[0].sum() EPS) alpha[0] * c[0] for t in range(1, T): # alpha[t-1] 与 A 的转置做乘法等价于 sum_i alpha[t-1,i] * A[i,j] alpha[t] (alpha[t - 1] self.A) * self.B[:, obs[t]] c[t] 1.0 / (alpha[t].sum() EPS) alpha[t] * c[t] log_prob -np.sum(np.log(c EPS)) return alpha, c, log_prob def backward(self, obs, c): T len(obs) beta np.zeros((T, self.N)) beta[T - 1] c[T - 1] for t in range(T - 2, -1, -1): beta[t] (self.A (self.B[:, obs[t 1]] * beta[t 1])) * c[t] return beta前向里的alpha[t-1] self.A是关键。alpha[t-1]是长度 N 的行向量self.A是 N×N矩阵乘法的结果第 j 个分量正好是 Σ_i alpha[t-1,i]·A[i,j]一步到位把内层循环消掉了。这个改写带来的速度提升非常明显值得记住。后向里的beta[T-1] c[T-1]容易写错。标准推导里如果前向做了缩放后向的初始化也必须用同样的缩放因子才能保证 α_t · β_t 的乘积在缩放意义下自洽。很多人这里直接初始化成 1结果 γ 的归一化之后数值看着差不多对但训练出来的参数总有点偏。def baum_welch_step(self, obs): T len(obs) alpha, c, log_prob self.forward(obs) beta self.backward(obs, c) gamma alpha * beta gamma / (gamma.sum(axis1, keepdimsTrue) EPS) # xi: (T-1, N, 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[t 1]][None, :] * beta[t 1][None, :]) xi[t] / (xi[t].sum() EPS) new_pi gamma[0].copy() denom_A gamma[:-1].sum(axis0)[:, None] EPS new_A xi.sum(axis0) / denom_A denom_B gamma.sum(axis0)[:, None] EPS new_B np.zeros_like(self.B) obs_arr np.asarray(obs) for k in range(self.M): mask (obs_arr k) if mask.any(): new_B[:, k] gamma[mask].sum(axis0) new_B / denom_B return new_pi, new_A, new_B, log_probxi那段用的是广播把 N×N 的中间结果一次算出来。如果写成两层循环代码可读性其实更高但速度慢很多。我建议先用循环版本验证正确性确认无误后再换成广播版本两个版本的结果应该逐元素一致这也是一种很好的自测手段。6.3 用天气-活动这组玩具数据做端到端验证现在造一组数据来跑通。隐藏状态是天气观测是活动。import numpy as np rng np.random.default_rng(42) # 真实参数只用于生成数据模型不知道 pi_true np.array([0.6, 0.4]) # 晴、雨 A_true np.array([[0.7, 0.3], [0.4, 0.6]]) B_true np.array([[0.6, 0.3, 0.1], # 晴: 散步/购物/打扫 [0.1, 0.4, 0.5]]) # 雨: 散步/购物/打扫 def sample_seq(T, pi, A, B): states, obs [], [] s rng.choice(len(pi), ppi) for _ in range(T): states.append(s) obs.append(rng.choice(B.shape[1], pB[s])) s rng.choice(A.shape[0], pA[s]) return np.array(states), np.array(obs) _, obs sample_seq(3000, pi_true, A_true, B_true) model DiscreteHMM(n_states2, n_obs3, seed7) prev None for it in range(50): new_pi, new_A, new_B, logp model.baum_welch_step(obs) model.pi, model.A, model.B new_pi, new_A, new_B model._normalize() if prev is not None and abs(logp - prev) 1e-4: print(f第 {it} 轮收敛, log P {logp:.2f}) break prev logp print(学到的 A:\n, np.round(model.A, 3)) print(学到的 B:\n, np.round(model.B, 3))跑出来的 A 和 B 不会和真实值完全一样但结构应该对得上一个状态的发射分布偏向散步另一个偏向打扫。这是无监督学习的固有特性——状态的编号是不可辨识的学出来的状态 0 可能对应真实的晴天也可能对应雨天取决于初始化。所以比较结果的时候不要直接比矩阵要看两个状态的发射分布是否有明确区分。如果你希望状态编号稳定可以在训练完之后按某个规则重排比如按 B 的第一列大小排序。这个操作在需要跨实验对比的时候很有用。6.4 和现成库对拍确认自己没写错从零实现最大的风险是跑出来了但数字不对。我强烈建议用一个成熟库做交叉验证。Python 生态里常用的是 hmmlearn。这里有个很多人踩过的坑hmmlearn 在较新的版本里改了接口语义。老的MultinomialHMM在 0.3.0 之后表示的是多项分布输入是计数对分类观测要用CategoricalHMM。如果你按老教程写MultinomialHMM会发现怎么都对不上甚至直接报错。from hmmlearn.hmm import CategoricalHMM ref CategoricalHMM(n_components2, n_iter50, random_state7) ref.fit(obs.reshape(-1, 1)) print(参考实现 log 似然:, ref.score(obs.reshape(-1, 1))) print(参考实现的发射矩阵:\n, np.round(np.exp(ref.emissionprob_), 3))对拍的时候注意两点。一是似然值不能直接比大小因为随机初始化和迭代轮数不同两个实现会收敛到不同的局部最优似然差几个点很正常。二是要看结构把两个实现学出来的 B 矩阵按行排序后比较如果分布形态一致比如都是一个状态偏散步、一个偏打扫说明你的实现是对的。真正需要警惕的是量级误差。如果你的实现算出来似然是 -20000参考实现是 -300那基本可以确定你的缩放或者归一化写错了这时候就该回去逐行检查forward和backward里 c 的用法。7. HMM 的边界在哪里它输给了谁又在哪里还在赢7.1 观测独立假设是它的硬伤HMM 有两个核心假设马尔可夫性当前状态只依赖上一状态和观测独立性当前观测只依赖当前状态。第二个假设在现实中经常站不住。举个具体的例子。做中文词性标注的时候一个词的词性显然受到它前后若干个词的影响而 HMM 的生成式结构让它在建模这种观测之间的相互依赖时非常吃力。更麻烦的是HMM 是生成式模型它要对 P(O) 建模也就是要建模任意观测序列出现的概率。对于一个有几十万词的词表这件事既没必要也做不好。判别式模型比如条件随机场 CRF直接对 P(Q|O) 建模绕开了对 P(O) 的建模还天然支持任意丰富的观测特征——前后词、词的后缀、是否包含数字等等。这就是为什么在序列标注任务上CRF 和后来的神经网络模型能全面压过 HMM。7.2 与主流序列模型的定位对比模型建模对象观测特征训练方式典型强项HMMP(O, Q) 生成式仅当前观测无监督 EM 或有监督计数数据少、缺标注、需要可解释状态CRFP(Q | O) 判别式任意上下文特征有监督序列标注精度CTCP(Q | O) 边缘化对齐神经网络特征有监督无需帧级对齐端到端语音识别序列到序列模型P(Q | O)神经网络特征有监督长距离依赖、生成任务这张表不是要评判谁高谁低而是想说清楚一件事它们解决的约束条件不同。CRF 和神经网络模型精度更高但它们都需要大量标注数据。如果你的场景里根本拿不到标注HMM 这种能靠 EM 从无标注数据里学出东西的能力就是别人替代不了的。7.3 今天还值得用 HMM 的几类场景说了这么多它的短板但我在实际工作里并没有把它彻底淘汰掉反而有几类场景还在持续使用。第一类是数据量小又有强结构先验的任务。比如工业设备的退化阶段识别我大概知道健康状态是缓慢变化的、故障状态一旦进入就不容易退出这种结构可以直接写进 A 矩阵的初始值里。几百条序列就够训练出一个可用的模型换成神经网络早就过拟合了。第二类是需要可解释性的场合。HMM 学出来的状态是有具体含义的——你打开 B 矩阵一看就知道状态 2 是个高异常倾向的状态。而在很多神经网络模型里你要额外做一大堆可解释性工作还未必看得清楚。在需要给非技术同事解释为什么系统判断这个设备异常的时候HMM 的优势非常实在。第三类是作为更强模型的组件。比如用 HMM 做初步的序列切分把长序列分段之后交给下游模型处理或者用它生成伪标签再用这些伪标签去训练神经网络。这种用小模型撬动大模型的用法我在数据标注资源紧张的项目里用过几次效果比直接硬上神经网络好。第四类是快速原型验证。想验证一个序列建模的方向有没有戏用 HMM 搭个基线只要几十行代码跑几分钟就出结果。对比之下搭一个神经网络基线可能要一两天。先用最便宜的手段确认问题本身可解再决定要不要投入更多资源这个顺序我建议一直保持。最后分享一个我在实际使用中的体会HMM 的调参重点从来不在算法本身而在状态数的选择和初始值的设置。算法是确定的代码写对了就是对的但状态数和初始值这两个东西直接决定了你能不能用它解决问题。我见过太多人在前向后向公式上反复纠结却从来没打印过 B 矩阵看看到底学出了什么。多做几次训完就打开参数看一眼的动作比多读几遍推导有用得多。