HMM隐马尔可夫模型的代码我前前后后手写过五遍以上最近一次是在某个序列分析任务里现成工具包对观测分布做定制化建模时不够灵活只能从零开始写。回头看这个模型的核心代码其实不长真正拦人的是三个地方数值下溢、维特比回溯的下标、以及 Baum-Welch 里 gamma 和 xi 的推导。这篇文章就把一整套能跑的代码拆开带着你从数据结构一路看到参数训练。这篇适合已经听说过隐马尔可夫模型、知道状态转移矩阵大概是什么但还没有完整实现过一遍的读者。如果你只是想调现成库这篇可能不是最快的路径但如果你想弄懂评估、解码、学习三个核心问题在代码层面分别长什么样这篇文章应该能帮你把最后那层窗户纸捅破。所有代码基于 Python 和 NumPy我尽量用最直接的方式写不引入奇技淫巧。1. 先确定要处理的数据长什么样状态、观测和三条假设1.1 隐状态不是隐藏的变量而是你观测不到的那一层很多人一开始会被隐状态这个名字绕晕。拿一个经典例子来说你掷骰子能看到的是一串点数比如 1、3、6、2、5这是观测序列但你看不到每次用的是公平骰子还是灌铅骰子这个每次用哪种骰子的序列就是隐状态序列。HMM 的假设是观测是状态生成的状态之间按概率转移而状态本身不可见。放到代码里观测序列 O 就是一个长度为 T 的整数数组数组里每个元素是观测符号的索引隐状态序列 Q 同样是一个长度为 T 的整数数组但它是未知的是我们要求解的东西。很多时候我们在意的不是概率本身而是这条状态序列——比如把一段语音切分成音素把一段操作日志归纳成几个阶段。1.2 模型参数在代码里就是四个 numpy 数组一个离散观测的 HMM参数可以用四个数组装完pi形状 (N,)初始状态分布pi[i] 表示第一个时刻处于状态 i 的概率。A形状 (N, N)状态转移矩阵A[i, j] 表示从状态 i 转移到状态 j 的概率。B形状 (N, M)发射概率矩阵B[i, k] 表示在状态 i 下观测到符号 k 的概率。O形状 (T,)观测序列。其中 N 是隐状态数量M 是观测符号种类数T 是序列长度。注意 A 和 B 都是行随机矩阵也就是每一行概率之和必须等于 1。初始化的时候最容易犯的错就是忘了归一化后面所有 log 概率、期望计数全部跑偏。import numpy as np class HMM: def __init__(self, n_states, n_obs, seed42): self.n_states n_states self.n_obs n_obs rng np.random.default_rng(seed) # 随机初始化但保证每行归一化 self.A rng.random((n_states, n_states)) self.A / self.A.sum(axis1, keepdimsTrue) self.B rng.random((n_states, n_obs)) self.B / self.B.sum(axis1, keepdimsTrue) self.pi rng.random(n_states) self.pi / self.pi.sum() # 这几行断言建议留在代码里以后改参数时能救你 assert np.allclose(self.A.sum(axis1), 1.0) assert np.allclose(self.B.sum(axis1), 1.0) assert np.allclose(self.pi.sum(), 1.0)我在实际写代码时会把断言放在init末尾。概率矩阵谁是行归一化、谁是列归一化非常容易记混一旦 B 的某一行出现全零后面 Viterbi 的 log 计算直接炸出 -inf排查起来很痛苦。1.3 三个基本假设决定了代码可以怎么简化HMM 能写成短短几十行靠的是三个很强硬的假设。第一个是齐次马尔可夫假设当前状态只依赖前一个状态和更早的历史无关所以只需要一张 A 矩阵就够了。第二个是观测独立假设当前观测只依赖当前状态和别的时刻的观测无关所以只需要一张 B 矩阵。第三个是初始分布假设第一个状态从 pi 里采样。这三个假设听起来简单但它们是所有递推公式的根基。如果哪天你看到有人讨论二阶 HMM 或者输入输出 HMM代码复杂度会立刻上一个台阶因为状态只依赖上一状态这个简化没了。你现在写的这套代码本质上就是这三个假设的直接翻译。2. HMM 要解决的三个问题其实就是三个函数2.1 三个经典问题的定义和输入输出HMM 的所有应用场景最后都会落到三个经典问题上。搞懂这三个问题代码的组织结构就清楚了。问题你在回答什么已知信息输出评估问题这条观测序列出现的概率有多大模型参数 λ 和观测 OP(O|λ)解码问题最可能的隐状态序列是哪一条模型参数 λ 和观测 O状态序列 Q*学习问题模型参数应该定成什么只有观测 O参数未知模型参数 λ用代码的话来说评估就是一个返回标量的函数解码就是一个返回整数数组的函数学习就是一个同时更新 pi、A、B 的函数。每次看到评估、解码、学习这三个词就在脑子里把它们映射回这三个函数签名很多书架感就消失了。2.2 各自对应哪个算法复杂度是多少评估问题用前向算法或后向算法解码问题用维特比算法学习问题用 Baum-Welch 算法。前两者的时间复杂度都是 O(T·N²)学习问题每轮迭代也是 O(T·N²)差别在于迭代次数。为什么不直接枚举因为状态序列的可能数量是 N 的 T 次方。假设 N3、T100这个数字比宇宙中的原子数还大暴力枚举完全不可行。动态规划的核心思想是复用中间结果前向算法把到 t 时刻为止的联合概率存下来维特比把到 t 时刻为止的最优路径分数存下来这才把指数级问题降到了多项式级。2.3 为什么参数学习是最容易写错的一环评估和解码都是给定参数做推理公式单一实现一次基本不会错。学习问题不一样它先要同时算前向和后向再用它们组合出两个期望统计量 gamma 和 xi最后才更新参数。任何一个地方的下标错一位结果看起来都挺合理但训练出来的模型就是不对。如果你写完了整套代码我强烈建议先单独验证前向和维特比确认没问题之后再写 Baum-Welch。下面两节我会按这个顺序来因为学习问题复用前向、后向的结果前面的地基打不牢后面根本没法调试。3. 前向算法代码评估序列概率以及数值下溢这个老熟人3.1 前向变量和递归公式的直觉前向算法定义了一个变量 α_t(i)表示到 t 时刻为止已经观测到 o_1 到 o_t并且 t 时刻处于状态 i的联合概率。它有一个非常漂亮的递推关系α_{t1}(j) ( Σ_i α_t(i) · A[i, j] ) · B[j, o_{t1}]用大白话解释所有状态 i 在 t 时刻的概率先通过转移矩阵把接力棒传给状态 j累加起来然后乘以状态 j 发出观测 o_{t1} 的概率。每一轮递推都是先转移、再发射两步操作。最终序列的总概率就是把最后一列的 α_T(i) 全部加起来。3.2 第一版代码先跑通再谈优化先把最直白的版本写出来def forward_naive(self, O): T len(O) alpha np.zeros((T, self.n_states)) alpha[0] self.pi * self.B[:, O[0]] for t in range(1, T): alpha[t] (alpha[t-1] self.A) * self.B[:, O[t]] total_prob alpha[-1].sum() return alpha, total_prob这段代码最核心的一行就是alpha[t] (alpha[t-1] self.A) * self.B[:, O[t]]。alpha[t-1] self.A是在做矩阵乘法得到每个状态 j 在转移之前的累计概率后面的* self.B[:, O[t]]是逐元素相乘相当于给每个状态 j 乘上它发出当前观测符号的概率。这里用到一个小技巧用self.B[:, O[t]]一次性取出所有状态对当前观测的发射概率避免了写一个 j 循环。这个版本在 T 很小的时候没问题但 T 一长就会出事。我第一版跑在 T500 左右的数据上时alpha 里的数值已经变成 1e-300 级别再往下直接下溢成 0总概率算出来是 0log 更是直接变 -inf。这就是数值下溢问题。3.3 缩放系数与对数似然为什么直接乘会爆掉连续乘大量小于 1 的概率数字会指数级缩小最终超出 float64 的表示范围。解决办法是在每一时刻算完 α_t 后用一个缩放系数 c_t 把 α_t 归一化让这一时刻的概率总和变成 1。这个 c_t 本身会被记录下来因为最终对数似然就是所有 log(c_t) 的和。def forward_scaled(self, O): T len(O) alpha np.zeros((T, self.n_states)) c np.zeros(T) alpha[0] self.pi * self.B[:, O[0]] c[0] alpha[0].sum() alpha[0] / c[0] for t in range(1, T): alpha[t] (alpha[t-1] self.A) * self.B[:, O[t]] c[t] alpha[t].sum() alpha[t] / c[t] loglik np.log(c).sum() return alpha, c, loglik注意缩放之后的 alpha 已经不是原始的联合概率了它变成了给定前 t 个观测的条件下t 时刻各状态的条件概率的某种比例形式。但这反而是后面 Baum-Welch 想要的因为缩放后的 alpha 和同样缩放后的 beta 相乘可以直接得到状态的后验概率。这个 c 数组要保留好后向算法和参数更新都要用到它。我第一次写的时候把 c 只当成临时变量丢掉了结果后向算法怎么都对不上后来才意识到缩放因子是前后向共享的桥梁。3.4 如何验证前向算法算对了写算法最怕看起来对但实际上不对。前向算法有两个很便宜的验证方法。第一个是取 N1也就是只有一个状态的退化情况。这时 P(O|λ) 应该等于所有发射概率的连乘def forward_check_single_state(hmm, O): p 1.0 for t in range(len(O)): p * hmm.B[0, O[t]] return p第二个方法是暴力枚举。当 N 和 T 都很小的时候把所有可能的状态序列都列出来每条序列算一遍联合概率再求和结果应该和 forward 一致import itertools def forward_brute_force(hmm, O): T len(O) total 0.0 for states in itertools.product(range(hmm.n_states), repeatT): p hmm.pi[states[0]] * hmm.B[states[0], O[0]] for t in range(1, T): p * hmm.A[states[t-1], states[t]] * hmm.B[states[t], O[t]] total p return total暴力版本是 O(N^T)只适合 T≤4、N≤3 的小玩具。但如果你的 forward 在这个玩具上能和暴力枚举对上那递推逻辑基本就是对的。我习惯把这两个检查封装成测试函数每次改动代码后跑一遍。4. 维特比算法代码解码隐藏状态回溯步骤最容易错4.1 动态规划表格delta 和 psi 两个矩阵评估问题是求所有状态序列的概率之和解码问题是求哪条状态序列的概率最大。维特比算法和前向算法长得非常像唯一区别是把求和换成了取最大值。它维护两个矩阵。delta[t, j] 表示到 t 时刻为止所有能够到达状态 j 的路径里概率最大的那条路径的分数psi[t, j] 表示达到这个最优分数时t-1 时刻处于哪个状态。psi 就是回溯时要用到的指针。4.2 一个 N2 的完整手算例子光讲公式不够直观我拿一个具体的小例子走一遍。假设 N2M2观测序列 O[0, 1, 0]参数如下A [[0.7, 0.3], [0.4, 0.6]] B [[0.6, 0.4], [0.2, 0.8]] pi [0.5, 0.5]为了手算方便我用自然对数先在 t0 时刻初始化δ_0(0) log(0.5) log(0.6) ≈ -1.204 δ_0(1) log(0.5) log(0.2) ≈ -2.303t1 时计算每个状态的最优前驱目标状态 j来自状态 0 的分数来自状态 1 的分数最优前驱δ_1(j)0-1.204 log(0.7) ≈ -1.561-2.303 log(0.4) ≈ -3.2190-1.561 log(0.4) ≈ -2.4771-1.204 log(0.3) ≈ -2.408-2.303 log(0.6) ≈ -2.8130-2.408 log(0.8) ≈ -2.631t2 时继续递推目标状态 j来自状态 0 的分数来自状态 1 的分数最优前驱δ_2(j)0-2.477 log(0.7) ≈ -2.834-2.631 log(0.4) ≈ -3.5470-2.834 log(0.6) ≈ -3.3441-2.477 log(0.3) ≈ -3.681-2.631 log(0.6) ≈ -3.1421-3.142 log(0.2) ≈ -4.751所以最优路径的最后一刻是状态 0。回溯时看 psi[2, 0]0说明 t1 时刻最优状态是 0再看 psi[1, 0]0说明 t0 时刻最优状态是 0。整条最优状态序列是 [0, 0, 0]。如果只看发射概率会觉得 t1 时刻观测为 1 时状态 1 更合理但综合转移概率之后状态 0 的整体路径分数更高——这就是维特比在局部直觉和全局最优之间做取舍的体现。4.3 回溯那两行代码为什么顺序不能反下面是完整的维特比实现我已经默认在 log 域里计算def viterbi(self, O): T len(O) N self.n_states logA np.log(self.A) logB np.log(self.B) logpi np.log(self.pi) delta np.zeros((T, N)) psi np.zeros((T, N), dtypeint) delta[0] logpi logB[:, O[0]] for t in range(1, T): for j in range(N): candidate delta[t-1] logA[:, j] psi[t, j] candidate.argmax() delta[t, j] candidate.max() logB[j, O[t]] # 回溯从最后一个时刻的最优状态倒推回去 states np.zeros(T, dtypeint) states[-1] delta[-1].argmax() for t in range(T-2, -1, -1): states[t] psi[t1, states[t1]] return states, delta回溯一共就两行核心代码但这里是最容易写错的地方。psi[t, j] 记录的是为了在 t 时刻到达状态 jt-1 时刻应当选哪个状态。所以在回溯时要先用当前已经确定好的 t1 时刻状态去查 psi[t1, states[t1]]得到 t 时刻的状态。常见错误是写成states[t] psi[t, states[t]]这等于拿还没确定的东西去查表结果完全不对。我第一次写 Viterbi 就在这个下标上栽了跟头输出的状态序列总是比预期慢了一拍。后来养成了一个习惯在纸上把 psi 矩阵画出来标出回溯的箭头方向再对照代码检查这类错就会好很多。4.4 用 log 概率改写的理由维特比在前向算法的基础上更进一步全程乘积在 log 域里变成了求和。log 是单调递增函数不管原来的概率谁大谁小取完 log 之后的大小关系不变所以 argmax 结果完全一致。好处有两个数值不会下溢而且加法比乘法快一点。如果你选择不用 log 域直接在原始概率上做乘法遇到长序列同样会面临下溢问题。所以我一律建议Viterbi 直接用 log 域前向算法则用缩放方案。两种处理方式的动机相同但实现细节不一样别混在一起。5. Baum-Welch 参数学习从期望计数到参数更新5.1 前向后向变量如何合体成 gamma 和 xi学习和前两个问题最大的区别是我们手里只有观测序列不知道状态序列所以要用期望来代替真实计数。Baum-Welch 的本质是 EM 算法E 步计算状态的后验概率M 步用这些后验概率作为权重重新估计参数。需要两个统计量。gamma_t(i) 表示在给定整条观测序列的条件下t 时刻处于状态 i 的后验概率xi_t(i, j) 表示 t 时刻处于状态 i 且 t1 时刻处于状态 j 的后验概率。前者是单点的后验后者是相邻两点的后验。它们的计算都需要后向变量 beta。后向变量 β_t(i) 表示给定 t 时刻处于状态 i从 t1 到 T 时刻观测到 o_{t1} 到 o_T 的概率。注意后向递推的方向是从 T 往 1 走而且要复用前向算法里存下来的缩放系数 cdef backward_scaled(self, O, c): T len(O) beta np.zeros((T, self.n_states)) beta[-1] 1.0 for t in range(T-2, -1, -1): beta[t] self.A (self.B[:, O[t1]] * beta[t1]) beta[t] / c[t1] return beta缩放后的前向 alpha 和缩放后的后向 beta 逐元素相乘再进行一行归一化就得到 gamma。xi 的计算稍复杂它等于 alpha_t(i) × A[i,j] × B[j,o_{t1}] × beta_{t1}(j)同样要除以 c[t1] 保持缩放一致。5.2 三组更新公式的代码落地得到 gamma 和 xi 之后参数更新就非常机械了全是分子除以分母的计数过程新 pi 就是 gamma[0]也就是第一个时刻各状态的后验概率。新 A[i,j] 等于所有时刻 xi_t(i,j) 之和除以所有时刻 gamma_t(i) 之和。新 B[j,k] 等于所有处于状态 j 且观测为 k的 gamma 之和除以所有处于状态 j的 gamma 之和。我把完整实现写出来def baum_welch(self, O, max_iter100, tol1e-6): T len(O) N self.n_states loglik_prev -np.inf for iteration in range(max_iter): alpha, c, loglik self.forward_scaled(O) beta self.backward_scaled(O, c) # E步计算后验概率 gamma alpha * beta gamma / gamma.sum(axis1, keepdimsTrue) xi np.zeros((T-1, N, N)) for t in range(T-1): xi[t] (alpha[t][:, None] * self.A) * \ (self.B[:, O[t1]] * beta[t1])[None, :] xi[t] / c[t1] # M步更新参数 new_pi gamma[0] new_A xi.sum(axis0) / gamma[:-1].sum(axis0)[:, None] new_B np.zeros_like(self.B) for k in range(self.n_obs): mask (O k) if mask.sum() 0: new_B[:, k] gamma[mask].sum(axis0) / gamma.sum(axis0) # 处理训练集中没出现过的观测符号避免log(0) new_B[new_B 0] 1e-6 new_B / new_B.sum(axis1, keepdimsTrue) self.pi new_pi self.A new_A self.B new_B # 收敛判断对数似然增量小于阈值就停止 if abs(loglik - loglik_prev) tol: break loglik_prev loglik return loglik实现中有几个地方要特别说明。gamma alpha * beta之后理论上每一行已经近似归一化了但浮点误差会累积所以再显式除以行和保证后续计数干净。xi[t]那行的向量化写法本质上就是在用广播机制同时算出 N×N 个组合比三重 for 循环快很多也都容易读懂。5.3 收敛判断、初值敏感和伪计数EM 算法保证对数似然单调不减所以每一轮迭代的 loglik 如果出现下降第一反应应该是代码有 bug而不是参数问题。我调试的时候会把每次迭代的 loglik 打出来正常情况是一串严格递增的数字最后趋于平稳。初值敏感是 EM 绕不开的坑。不同的随机初始化可能收敛到不同的局部最优解所以实际训练时我习惯跑 5 到 10 个不同的随机种子选对数似然最高的那组参数。如果某些状态在训练中完全没被用到gamma 的某一列会接近零这时候就需要给 B 加一点伪计数或者干脆减少状态数 N。另外用留出集来判断状态数 N 是否合适。训练集上 N 越大对数似然越高但留出集上会先升后降那个拐点对应的 N 才是比较合理的选择。这一步几乎没人写在基础教程里但实际项目里非常管用。6. 我的调试经验怎么确认你写的 HMM 是对的6.1 用已知参数生成数据再反推回去调试 HMM 最可靠的方法是先用一组已知参数生成人工数据再用你的训练代码去学习参数看能不能学回来。先生成数据def generate(self, T, seedNone): rng np.random.default_rng(seed) states np.zeros(T, dtypeint) obs np.zeros(T, dtypeint) states[0] rng.choice(self.n_states, pself.pi) obs[0] rng.choice(self.n_obs, pself.B[states[0]]) for t in range(1, T): states[t] rng.choice(self.n_states, pself.A[states[t-1]]) obs[t] rng.choice(self.n_obs, pself.B[states[t]]) return obs, states然后用 baum_welch 训练最后比较学到的 A、B 和真实参数。这里有个坑HMM 的状态编号是有任意性的你学到状态 0不一定对应真实的状态 0可能整体换了个顺序。所以在验证时要枚举所有状态排列找到和真实参数最接近的那个排列再比较。N3 时只有 6 种排列直接暴力枚举就行N 大了就随机抽一些排列做近似比较。我就是靠这个方法抓到一个隐藏很深的问题我的 alpha、beta、gamma 单独看都合法但 xi 的缩放系数少除了一个 c[t1]导致学习得到的转移矩阵系统性偏差。用真实参数生成数据后一对比偏差立刻暴露出来。6.2 状态坍缩和零概率的排查训练 HMM 时最常见的现象是状态坍缩某个状态在 gamma 里占据了几乎所有概率其他状态形同虚设学出来的模型退化成只有一个有效状态。这种时候先检查初始化纯随机初始化很容易让 EM 走进这个局部最优。解决办法是给 B 设置一个下限避免某些状态一开始就完全不可能发出常见观测或者减小 N从小的状态数开始训练再逐步增加。还有一个隐蔽的坑是 B 里出现严格为零的概率。训练数据里某个观测符号在某个状态下从没出现过新 B 对应位置就会变成 0。这正是我在 baum_welch 里加new_B[new_B 0] 1e-6的原因——否则下一次 Viterbi 在 log 域里会遇到 -inf整个状态序列直接废掉。不过加完下限后记得重新归一化行这个顺序不能反。6.3 什么情况下别硬上 HMM每次写完这套代码我都会提醒自己一句话HMM 不是万能的。如果观测是高维连续值比如原始语音帧离散的 B 矩阵装不下需要换成高斯发射概率或者 GMM-HMM如果数据存在明显的长程依赖一阶马尔可夫假设根本不够用CRF 或循环神经网络是更好的选择如果你的数据根本没有时序概念那 HMM 压根就不该出场。另外如果你只是要做简单的序列分类HMM 还要配上分类器才能工作不如直接用带时序结构的分类模型省事。我见过很多项目把问题包装成 HMM结果模型训出来效果还不如一个简单的逻辑回归加滑窗特征。HMM 的优势在于状态本身有解释性以及能优雅地处理不等长序列这两点满足的时候它才是那个正确的工具。写这套代码的过程中我最大的体会是HMM 的数学符号看起来很重但一旦翻译成 NumPy 的矩阵操作每一步都变得非常机械。遇到 bug 不要急着怀疑公式先从前向算法开始逐层验证用玩具数据把每个模块都焊死再往上搭 Baum-Welch。这套自底向上验证的习惯帮我省下的调试时间远远超过写代码本身的时间。最后再分享一个小技巧把训练过程中的 loglik 变化曲线画出来保存下来以后每次改动代码都能一眼看出收敛行为是否正常。一个健康的训练过程那条曲线应该平滑上升最后稳稳地停住如果它震荡或者下降赶紧回头查缩放因子十有八九问题就出在那里。