HMM（隐马尔科夫模型）这几个字，我最早是在概率论教材的目录里看到的，真正把它用起来是在一个语音端点检测的改造项目里。当时要在一段带噪声的信号里判断每一帧属于静音、说话还是过渡状态，观测特征倒是现成的，但真正驱动这些特征的状态序列是看不见的，我们只能拿观测序列反推状态。HMM就是做这件事的标准工具。这篇整理不打算把公式抄一遍完事，而是按我自己“理解、手写、调参、踩坑”的顺序，把HMM的直觉、三个经典问题和工程经验串成一份可以直接参考的笔记。

如果你是刚接触概率图模型的初学者，看完应该能独立写出前向算法和维特比算法；如果你是在准备面试或者做序列标注类项目，这篇里的数值稳定性、初值设置、状态数选择这些细节，大概率比教科书目录更有用。

1. 先建立直觉：状态看不见，但能通过观测猜出来

1.1 从一个“猜心情”的场景说起

想象你和一个不太爱说话的室友同住，每天只能通过TA带回来的外卖口味判断当天心情：心情好可能买奶茶，心情差可能直接啃面包。你看到的是奶茶、面包、炸鸡这些观测值，但室友真正处于“高兴”还是“低落”是看不到的，只能从外卖反推。

这就是HMM的核心感觉。HMM全称隐马尔科夫模型，里面有两层结构：一层是状态序列，一层是观测序列。状态是隐藏的、看不见的，观测是暴露的、我们能拿到的。语音识别里，观测是声学特征，状态是音素；词性标注里，观测是词语本身，状态是名词动词形容词；生物信息里，观测是碱基序列，状态是编码区和非编码区。

我最早犯的一个理解错误是：以为HMM是个分类器，给它一个观测就能预测一个状态。其实它建模的是整条序列的联合概率，状态之间还会随时间演化。这个“时间演化”体现在状态转移矩阵上，也是马尔可夫假设发挥作用的地方。

1.2 马尔可夫链和“隐”字到底指什么

先复习马尔可夫链。马尔可夫链说的是：当前状态只依赖于前一个状态，跟更早的历史没有关系。就像你判断今天会不会下雨，只看昨天的天气就够了，不用翻一周前的记录。这就是“齐次马尔可夫假设”。

HMM在马尔可夫链基础上加了一层“隐”：链上的状态你观察不到，你能观察到的只是状态“发射”出来的观测值。比如状态是“晴天”，你观察到的可能是“室友带了一杯冰美式”；状态是“雨天”，你观察到的可能是“室友带了一把伞”。状态本身没有直接暴露，观测是它的间接信号。

一个完整的HMM由五样东西组成，俗称五元组：

下面代码里我会用 NumPy 数组直接表示 A 和 B，行代表当前状态，列代表目标状态或观测类别，这个方向别搞反了，不然调参的时候会很痛苦。

1.3 两个核心假设：换来高效计算的关键

HMM能高效工作的基础是两条假设：

• 齐次马尔可夫假设：t时刻的状态只依赖t-1时刻的状态，不依赖更早的历史。
• 观测独立假设：t时刻的观测只依赖t时刻的状态，不依赖当前状态之外的信息。

这两条假设严格来说在真实场景里经常不成立。语音识别中，当前音素往往受到前后音素协同发音的影响，不只是看前一个音素；文本里一个词的性质也可能依赖上下文多个词。后来大家用线性链CRF、RNN、Transformer来弱化这些假设，但代价是参数更多、训练更贵、可解释性下降。

在样本量不大、需要快速上线、模型要能解释得清的场景里，HMM这两条假设反而成了优势。它把问题压缩到一个非常小的参数空间里，训练和推断都非常稳定。我个人的看法是：HMM不是被淘汰的旧工具，而是一个你手里的“最小可行模型”。

1.4 贯穿全文的例子：天气和冰淇淋

为了说话方便，我用一个经典到不能再经典的例子，来自早期概率图模型教材里的“冰淇淋问题”：

假设有两个隐藏状态：热天(Hot)和冷天(Cold)。你每天记录某人吃的冰淇淋数量，观测值有三种：1个、2个、3个。

初始概率：第一天是热天的概率0.8，冷天0.2。

状态转移矩阵 A：

发射概率矩阵 B，行是状态，列是观测1/2/3：

数字随便给的，但方向符合直觉：冷天吃3个冰淇淋的现象反而更扎眼，所以发射概率高；热天吃1个也不稀奇，概率反而给得高一点。这个例子小小的，却能跑通后面所有算法。

2. 三大经典问题：算概率、找路径、学参数

2.1 概率计算：这个观测序列有多“正常”

第一个问题：给定模型参数 λ=(A,B,π)，计算观测序列 O=(o1,o2,...,oT) 出现的概率 P(O|λ)。

这能力有什么用？最典型的场景是异常检测。你给一段正常的语音序列算一个分数，再给一段可能异常的语音序列算一个分数，正常度明显低的就有问题。也可以用多个模型分别算概率，概率最大的那个模型对应的类别就是分类结果。

最朴素的做法是枚举所有可能的状态序列：假设N个状态，序列长度T，状态序列有N的T次方种，每个都要计算连乘概率再求和。T一长就是天文数字，根本不可能。

前向算法解决这个问题。它定义一个前向变量：α_t(j) 表示“到t时刻为止，观测序列的前t个值都已经看到，并且t时刻处于状态j”的概率。递推关系是：

α_1(j) = π_j × B_j(o1)

α_t(j) = B_j(ot) × Σ_i α_{t-1}(i) × A_{i,j}

最后 P(O|λ) = Σ_j α_T(j)。

这个过程其实就是动态规划。如果你写过斐波那契数列，就理解它的思路：每一步都站在上一步已经算好的结果上，避免重复计算。我习惯把它理解成“填一张 T×N 的表，每个格子只算一次”。

拿冰淇淋例子来说，观测序列是[3,3,1]，索引对应观测类别分别是2、2、0（0表示吃1个，1表示吃2个，2表示吃3个）。手动跑一下：

第1天：

• 热：0.8 × 0.3 = 0.24
• 冷：0.2 × 0.5 = 0.10

第2天，观测还是3个：

• 热：0.3 × (0.24×0.6 + 0.10×0.5) = 0.0582
• 冷：0.5 × (0.24×0.4 + 0.10×0.5) = 0.073

第3天，观测变成1个：

• 热：0.3 × (0.0582×0.6 + 0.073×0.5) = 0.021426
• 冷：0.2 × (0.0582×0.4 + 0.073×0.5) = 0.011956

整条序列概率约等于0.0334。

这里有个非常直观的观察：序列越长，概率越小，因为每一步都在乘小于1的数。到后面概率会小到超出浮点数能表达的范围，这个问题我在第4章详细说。

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

第二个问题：已知观测序列 O，求最可能的状态序列 Q。这就是解码，也叫状态标注。

词性标注就是典型应用：给你一句话里每个词的词性序列，就是HMM的解码过程。中文分词也可以用这个思路，把“词位”（词首、词中、词尾、单独成词）当状态，把每个汉字当观测。

解码用维特比算法，跟前向算法长得像，但前向做的是求和，维特比做的是取最大值。它定义两个变量：δ_t(j) 是“到t时刻为止，沿着某条最优路径走到状态j，产生前t个观测的最大概率”；ψ_t(j) 则是“这个最大值是从上一时刻哪个状态走过来的”。

递推公式：

δ_t(j) = B_j(ot) × max_i [ δ_{t-1}(i) × A_{i,j} ]

ψ_t(j) = argmax_i [ δ_{t-1}(i) × A_{i,j} ]

计算完最后一列后，从最大概率的那个状态开始，沿着 ψ 往回追溯，就是完整的最优状态序列。

用冰淇淋例子直接算序列[3,3,1]：

δ1：热0.0582，冷0.073 δ2：热0.01095，冷0.0073

注意第2天那个观测“吃3个”的时候，冷天的单点概率更高，但维特比最终回退出来的路径是：热 → 冷 → 热。

为什么中间出现冷，后面又热回去？因为状态切换有成本，而发射概率的差异会权衡。第2天冷天出3个的概率0.5远高于热天0.3，所以全路径被拉成冷；第3天观测只有1个，热天出1的概率0.3略高于冷天0.2，加上转移路径的贡献，最终又回到热。你如果只会做贪心选择，看到哪一步概率大就选哪一步，很容易选出一条比全局最优差很多的路径。维特比的价值就是：它忍得住单步的诱惑，把整条路径的乘积放到一起比较。

2.3 学习问题：从观测反推模型参数

第三个问题：只有一堆观测序列，没有状态标注，怎么把 A、B、π 学出来？

这是最接近“机器学习”气息的问题。教科书标准答案叫鲍姆-韦尔奇算法，本质是EM算法的一个特例。思路很朴素，一句话概括：先胡乱猜一组参数，用这组参数推断状态的后验分布，再用这些后验概率重新算参数，反复迭代直到收敛。

E步（期望步）里需要两个中间量：

• γ_t(i)：t时刻处于状态i的概率，由前向概率α_t(i)和后向概率β_t(i)乘积归一化得到。
• ξ_t(i,j)：t时刻在状态i、t+1时刻在状态j的概率。

M步（最大化步）用这两个量重新估计参数：

• π_i 取 γ_1(i)
• A_{i,j} 是 Σ_t ξ_t(i,j) 除以 Σ_t γ_t(i)
• B_j(k) 是所有观测为k的时刻的 γ_t(j) 之和，除以所有时刻的 γ_t(j) 之和

每次更新后，P(O|λ) 一定不会下降，这是EM算法保证的性质。实操中你会看到似然值每次都涨一点点，涨到某个水平就平稳了。不过它只能保证收敛到局部最优，所以初始参数怎么给很关键，这个坑我在后面单独讲。

如果你有标注好的状态序列，学习就更简单了，直接用频率统计：A_{i,j} 就是状态i后面跟状态j的次数除以状态i出现的总次数，B同理。现实情况里标注往往很贵，鲍姆-韦尔奇的价值就是能在没有标注时通过迭代挖掘出隐含结构。

3. 代码实操：用NumPy手写核心流程

3.1 准备数据和两个基础函数

先明确一下代码里的状态编号和观测编号，避免写乱：

状态：0=热，1=冷。观测：0=吃1个，1=吃2个，2=吃3个。

前向算法用数组实现非常简洁。核心一步是alpha[t-1] @ A，这个矩阵乘法的结果是一个长度为N的向量，第j个元素就是 Σ_i α_{t-1}(i) × A_{i,j}，然后逐元素乘上 B 的对应列。

import numpy as np def forward(obs, A, B, pi): T = len(obs) N = A.shape[0] alpha = np.zeros((T, N)) alpha[0] = pi * B[:, obs[0]] for t in range(1, T): alpha[t] = B[:, obs[t]] * (alpha[t - 1] @ A) return alpha

维特比算法稍微长一点，主要多了一个回溯指针。

def viterbi(obs, A, B, pi): T, N = len(obs), A.shape[0] delta = np.zeros((T, N)) psi = np.zeros((T, N), dtype=int) delta[0] = pi * B[:, obs[0]] psi[0] = -1 for t in range(1, T): for j in range(N): tmp = delta[t - 1] * A[:, j] psi[t, j] = np.argmax(tmp) delta[t, j] = np.max(tmp) * B[j, obs[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 初始化为 -1，因为它只在 t≥1 时才会用到。用回溯生成路径时，从最后一个时刻开始倒推，每一步查上一时刻的指针。

3.2 跑一遍前向和维特比，验证手工计算结果

基于前面的冰淇淋参数，我直接跑：

pi = np.array([0.8, 0.2]) A = np.array([ [0.6, 0.4], [0.5, 0.5] ]) B = np.array([ [0.3, 0.4, 0.3], [0.2, 0.3, 0.5] ]) # 观测依次是：吃3个、吃3个、吃1个 obs = np.array([2, 2, 0]) alpha = forward(obs, A, B, pi) print("最后时刻alpha:", alpha[-1]) print("观测序列概率:", alpha[-1].sum()) print("维特比路径:", viterbi(obs, A, B, pi))

我这边跑出来的最后时刻 alpha 向量约是 [0.0214, 0.0120]，求和得到0.0334左右，跟第2章手动推导结果一致。维特比输出 [0, 1, 0]，也就是热 → 冷 → 热。

拿到 alpha 之后还能做一件有价值的事：把每一时刻的 alpha 归一化，就能看到“在观测已知的情况下，每一时刻最可能是哪个状态”的概率分布。这跟维特比不同，它给出的是逐点概率最大的状态，不保证整体最优。两个结果放在一起对比，能帮你判断模型是“模棱两可”还是“非常确定”。

3.3 补上后向算法和一个极简Baum-Welch

前向算法只往前看，Baum-Welch还需要后向概率：β_t(i) 表示“已知t时刻处于状态i，从t+1到T的观测序列”的概率。后向递推反过来跑，从序列末尾往前算。

def backward(obs, A, B, pi): T, N = len(obs), A.shape[0] beta = np.zeros((T, N)) beta[-1] = 1.0 for t in range(T - 2, -1, -1): beta[t] = (A * B[:, obs[t + 1]]) @ beta[t + 1] return beta

有了前向和后向，就能写一个最短的Baum-Welch迭代。我用一个单序列版本做示例，重点看更新公式是怎么落地的：

def baum_welch(obs, A, B, pi, max_iter=50): T, N = len(obs), A.shape[0] K = B.shape[1] for _ in range(max_iter): alpha = forward(obs, A, B, pi) beta = backward(obs, A, B, pi) gamma = np.zeros((T, N)) xi = np.zeros((T - 1, N, N)) for t in range(T - 1): denom_xi = (alpha[t][:, None] * A * B[:, obs[t + 1]] * beta[t + 1][None, :]).sum() xi[t] = alpha[t][:, None] * A * B[:, obs[t + 1]] * beta[t + 1][None, :] / denom_xi denom_gamma = (alpha[t] * beta[t]).sum() gamma[t] = alpha[t] * beta[t] / denom_gamma denom_gamma_last = (alpha[-1] * beta[-1]).sum() gamma[-1] = alpha[-1] * beta[-1] / denom_gamma_last pi = gamma[0] A = xi.sum(axis=0) / xi.sum(axis=0).sum(axis=1, keepdims=True) B.fill(0) for k in range(K): B[:, k] = gamma[obs == k].sum(axis=0) B = B / gamma.sum(axis=0)[:, None] return A, B, pi

初值我一般不会全随机，而是先用一个粗略的均匀分布，或者用k-means对观测做硬聚类，再用聚类结果数一下状态之间的频率。无监督学习这玩意儿，初值不好，最后收敛出来的模型可能完全没有业务含义。代码跑通是一回事，跑到有意义的模型是另一回事。

4. 工程化避坑：数值稳定性、初值与状态数选择

4.1 数值下溢：长序列概率会小到“消失”

在前向算法里，每一步都在乘小于1的概率，序列稍微长一点，alpha 的值就会小到浮点数无法表示，出现下溢。假设每步乘0.3，100步之后就变成10的负53次方量级，活活把 double 的精度跑干。

工程上有两种主流处理手段。

第一种是对数域计算：把前向递推里所有乘法变成加法，log α_t(j) = log B_j(ot) + log(Σ_i α_{t-1}(i) × A_{i,j})。注意里面还有个求和，不能直接把log项相加，需要用 log-sum-exp 技巧：

log(Σ exp(x_i)) = max(x_i) + log(Σ exp(x_i - max(x_i)))

这个操作在代码里非常常见，自己写一个小函数封装一下就好。

第二种是缩放归一化：每算完一列 alpha，就除以这一列的和，把概率压回0到1的范围内。这样做会额外记录一个缩放因子向量，最后整个观测序列的似然等于所有缩放因子的乘积乘以最后归一化列的和。前向后向算法里用了缩放之后，Baum-Welch 的 gamma、xi 计算也会稳定得多。

我自己的经验是：短序列（几十步）用原始算法没问题；序列长度超过几百就一定要做数值稳定化处理。维特比算法也同理，它取的是最大值路径，直接把概率换成 log 概率，加和变成取最大，根本不会出现下溢。

4.2 初值敏感：无监督HMM最容易被坑的地方

鲍姆-韦尔奇算法保证的是局部收敛，不是全局收敛。如果 B 的初值给得离谱，比如某个状态的发射概率全是0，那么迭代过程中 gamma 对应列永远是0，这个状态就“死”掉了，后面怎么更新都救不回来。

所以初始化 B 的时候有两个原则：

• 所有概率必须严格为正，别给0。给个0.01这种小地板值，等于留了一条活路。
• 发射概率的初值最好能体现状态之间的区分度，而不是全都均值化。全均值化会让迭代前期非常漫长，而且容易掉进对称解里。

我常用的初始化套路是：先用 k-means 对观测序列做聚类，把聚类中心映射成状态的初始发射概率；再用聚类后的硬标签统计状态转移频率，得到 A 的初值。这样等于先用无监督方法找了一个“半醒”的起点，再交给 EM 去微调。多跑几个随机种子，选出最终似然值最高的那个，也是常规操作。

4.3 状态数K怎么定：不是越大越好

HMM里的N（状态数）是要人为指定的。业务上有明确语义的状态，比如“静音/说话/过渡”，直接按语义定。没有语义时，就得用数据说话。

一个简单有效的方法是试多个K值，分别训练HMM，然后看模型在验证集上的平均对数似然。K太小，模型拟合不足，似然很低；K太大，模型过拟合，训练似然很高但验证集似然反而下降。取验证集似然曲线的“肘部”位置，或者用BIC这类带惩罚项的准则。

我踩过最典型的坑是：为了追求训练集上的拟合效果，把K调得很大，结果每个状态都退化成只发射一种观测的死板单元，完全丧失泛化能力。序列建模问题里，状态过多带来的过拟合比状态过少带来的欠拟合更难排查，因为训练日志看起来非常漂亮，上线就露馅。

4.4 多序列训练和不等长序列处理

现实里你手里的样本不是一条长序列，而是很多条长度不等的短序列。Baum-Welch 迭代时要对每条序列分别跑前向后向，然后把所有序列的 gamma、xi 对应相加，最后用累加结果统一更新参数。

这里有个容易踩的细节：不要把每条序列单独拿去做一轮EM更新。正确的做法是全部序列跑完一轮“E步”，再统一做“M步”，不然参数会被最后看到的序列带偏。

不等长序列做 batch 训练时，一般要 padding 到统一样式，同时用 mask 把无效位置盖住。我在代码里会在 forward 函数里传一个 mask 数组，把 padding 部分的 beta 设成1、alpha 设成0，这样它们既不影响递推，也不参与参数更新。这个小细节，很多开源库里隐藏得很深，但一旦你自己实现就会发现它是多序列训练能否收敛的关键。

5. 常见问题排查与调参速查

5.1 问题速查表

5.2 我踩过的三个印象深的坑

第一个坑是发射概率矩阵里出现了0。当时有个业务场景，某个观测类别在训练集里一直没出现，B对应列直接是0。训练阶段看着没事，一上线新数据出现了这个观测，前向概率立刻变0，整个链路崩掉。后来我把所有概率初始化加了一个0.01的地板值，并且规定每次M步更新后都做一次同样处理，再没出过问题。

第二个坑是Baum-Welch迭代时我把 π 和 A 的更新顺序写错了。EM的原则是：用旧参数跑完所有统计量，再一次性更新所有参数。我一开始图省事，π 更新完立刻用新 π 跑后续公式，结果迭代几次后似然反而往下掉，查了大半天才醒悟过来。这个错误在代码层面非常隐蔽，因为单步看起来都合理。

第三个坑跟状态数选择有关。某次做序列聚类，业务方觉得“模式越细越好”，K给了20。模型训练完，一半状态只有几条序列支撑，随便一点扰动就会让状态归属大翻转。后面我改成先画不同K下的验证集对数似然曲线，发现K在8左右就到平台期，再往上纯粹是过拟合。从那时候起，K的选择在我这儿就是一门“先看曲线再拍脑袋”的学问。

5.3 一个很实用的诊断手段

模型训练完别急着看指标。把某个长度为几百的观测序列跑一遍前向算法，把每一时刻 alpha 的归一化向量画成热力图。横轴是时间，纵轴是状态，格子颜色代表概率高低。这张图能非常直观地告诉你：模型是不是真的把某些时序模式识别出来了，还是只是把所有状态糊成一团。

我还喜欢把维特比路径叠加在原始观测值上面画。状态切换的位置是否对应观测值明显变化的拐点，往往一眼就能判断这个HMM是学到了结构化规律，还是在硬背训练集的细节。这个方法我推荐给每个刚开始写HMM的人，比盯着似然值可靠得多。

最后说一点个人体会。HMM不是那种能给你惊喜的模型，它不会突然把准确率抬到天上，但它会清清楚楚地把序列里的转移结构和发射偏好摆在你面前。第一次自己写完前向后向维特比和Baum-Welch那一整套，你再看CRF、再看RNN，会突然明白很多东西是打哪儿来的。这几年我见过不少项目一上来就想用大模型暴力解序列问题，结果样本量只有几千条，训出来还没HMM稳。遇到序列数据，先问自己三个问题：状态大概有几种、转移是否有明确含义、标注成本高不高。只要这三个问题回答得上，HMM就还是一个值得先试的方案。