如果你做过序列建模哪怕是两年经验一定绕不开隐马尔可夫模型这个名字。学的时候感觉公式全认识一上手写代码就发现到处是坑前向算法算到第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_{t1} 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的每一行之和为1B的每一行之和也为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)α_{t1}(j) [Σ_{i1}^N α_t(i) * a_ij] * b_j(o_{t1})最后得到P(O|λ) Σ_{j1}^N α_T(j)我解释一下这个递推为什么是对的。t1时刻落在状态j必须先从t时刻的某个状态i跳到j这一步的贡献是α_t(i) × a_ij把所有i的可能加起来就是在t1时刻“到达状态j且观测到o1到o_t”的累积概率。然后再乘上在状态j下观测到o_{t1}的发射概率b_j(o_{t1})就补上了当前时刻的观测。这个算法的时间复杂度是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)) # 初始化t0时各状态的概率 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(fP(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)δ_{t1}(j) max_i [δ_t(i) * a_ij] * b_j(o_{t1})同时要维护一个回溯表ψ_t(j)记录δ_{t1}(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), dtypeint) # 初始化 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_{t1} s_j | O, λ)这两个量需要同时用到前向变量α和新的后向变量β。β_t(i)定义为t时刻状态为s_i从t1到T时刻观测已知的概率也叫后向变量递推公式直接反着来。M步用这些期望去重新估计参数。转移矩阵的重估公式是a_ij_new Σ_t ξ_t(i, j) / Σ_t γ_t(i)发射矩阵的重估公式是b_j(k)new Σ{t:o_tvk} γ_t(j) / Σ_t γ_t(j)初始分布就直接取γ_1(i)。迭代到对数似然收敛就用学到的参数作为最终模型。5.2 训练代码的骨架把前向、后向和重估合在一起就是一个可训练的HMM类。我这里只写关键逻辑完整的封装长度太长拆开讲更清楚。def baum_welch(A, B, pi, observations, iterations100): 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[t1]] * beta[t1, 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 常见问题速查表现象可能原因排查方向前向概率很快变成0未做缩放序列太长加scaling检查每步归一化系数维特比路径死循环T不前进回溯表索引错误检查psi[t, path[-1]]的维度鲍姆-韦尔奇对数似然下降EM代码bug或缩放错误单独验证E步的γ和ξ求和是否为1训练不收敛参数振荡学习率过大或数据太少检查似然曲线调整初始化或换模型解码路径一直停留一个状态转移矩阵对角线过大查看A的对角元调小或重新初始化6.5 一个从实践中得来的调参心得我刚做HMM项目时最爱犯的错误是忽略状态个数的选择。N设得太多模型过拟合训练慢N设得太少表达力不够序列里真实存在的结构拟合不出。有一种相对有效的做法是画对数似然随N变化的曲线找拐点。N从2开始往上加每加1训练一次模型看对数似然提升幅度当提升明显变缓时就别再加了。这个办法是经验性的但比拍脑袋选N靠谱得多。7. 沿这个方向还能往哪走HMM本身是一个很“底层”的工具箱。掌握三大算法之后往上可以扩展出很多实用的变体高斯混合HMMGMM-HMM把发射概率从离散分布换成高斯混合语音识别里非常主流输入输出HMMIOHMM把观测也纳入状态转移的条件处理带外部输入的序列层次HMM、分段HMM处理长时依赖和状态驻留时间建模和循环神经网络结合用RNN做特征提取用HMM做序列解码取两者之长如果你主要是为了解决序列标注问题还可以对比条件随机场的方案。HMM是生成式模型优化的是联合概率CRF是判别式模型直接建模条件概率。数据量大的时候CRF通常精度更高数据量小或者需要生成式特性的时候HMM更好使。没有哪个绝对最优取决于你的数据和任务形态。我个人在实际项目里体会最深的一点是不要一上来就让HMM解决所有问题。它最适合的场景是状态转移结构相对清晰、观测噪声可控、序列长度适中。一旦序列极长或者状态依赖跨度很大HMM的效果就会明显下滑这时候就该考虑用循环网络或者状态空间模型去替代它。模型选择本身比调参重要得多。