企业官网建设流程全解析

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, seed=42): 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(axis=1, keepdims=True) self.B = rng.random((n_states, n_obs)) self.B /= self.B.sum(axis=1, keepdims=True) self.pi = rng.random(n_states) self.pi /= self.pi.sum() # 这几行断言建议留在代码里，以后改参数时能救你 assert np.allclose(self.A.sum(axis=1), 1.0) assert np.allclose(self.B.sum(axis=1), 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 的所有应用场景，最后都会落到三个经典问题上。搞懂这三个问题，代码的组织结构就清楚了。

用代码的话来说，评估就是一个返回标量的函数，解码就是一个返回整数数组的函数，学习就是一个同时更新 pi、A、B 的函数。每次看到"评估、解码、学习"这三个词，就在脑子里把它们映射回这三个函数签名，很多书架感就消失了。

2.2 各自对应哪个算法，复杂度是多少

评估问题用前向算法或后向算法，解码问题用维特比算法，学习问题用 Baum-Welch 算法。前两者的时间复杂度都是 O(T·N²)，学习问题每轮迭代也是 O(T·N²)，差别在于迭代次数。

为什么不直接枚举？因为状态序列的可能数量是 N 的 T 次方。假设 N=3、T=100，这个数字比宇宙中的原子数还大，暴力枚举完全不可行。动态规划的核心思想是复用中间结果：前向算法把"到 t 时刻为止的联合概率"存下来，维特比把"到 t 时刻为止的最优路径分数"存下来，这才把指数级问题降到了多项式级。

2.3 为什么"参数学习"是最容易写错的一环

评估和解码都是给定参数做推理，公式单一，实现一次基本不会错。学习问题不一样，它先要同时算前向和后向，再用它们组合出两个期望统计量 gamma 和 xi，最后才更新参数。任何一个地方的下标错一位，结果看起来都挺合理，但训练出来的模型就是不对。

如果你写完了整套代码，我强烈建议先单独验证前向和维特比，确认没问题之后再写 Baum-Welch。下面两节我会按这个顺序来，因为学习问题复用前向、后向的结果，前面的地基打不牢，后面根本没法调试。

3. 前向算法代码：评估序列概率，以及数值下溢这个老熟人

3.1 前向变量和递归公式的直觉

前向算法定义了一个变量 α_t(i)，表示"到 t 时刻为止，已经观测到 o_1 到 o_t，并且 t 时刻处于状态 i"的联合概率。它有一个非常漂亮的递推关系：

α_{t+1}(j) = ( Σ_i α_t(i) · A[i, j] ) · B[j, o_{t+1}]

用大白话解释：所有状态 i 在 t 时刻的概率，先通过转移矩阵把"接力棒"传给状态 j，累加起来，然后乘以状态 j 发出观测 o_{t+1} 的概率。每一轮递推都是"先转移、再发射"两步操作。最终序列的总概率就是把最后一列的 α_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 一长就会出事。我第一版跑在 T=500 左右的数据上时，alpha 里的数值已经变成 1e-300 级别，再往下直接下溢成 0，总概率算出来是 0，log 更是直接变 -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 如何验证前向算法算对了

写算法最怕"看起来对但实际上不对"。前向算法有两个很便宜的验证方法。第一个是取 N=1，也就是只有一个状态的退化情况。这时 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), repeat=T): 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 一个 N=2 的完整手算例子

光讲公式不够直观，我拿一个具体的小例子走一遍。假设 N=2，M=2，观测序列 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]

为了手算方便我用自然对数，先在 t=0 时刻初始化：

δ_0(0) = log(0.5) + log(0.6) ≈ -1.204 δ_0(1) = log(0.5) + log(0.2) ≈ -2.303

t=1 时，计算每个状态的最优前驱：

t=2 时继续递推：

所以最优路径的最后一刻是状态 0。回溯时看 psi[2, 0]=0，说明 t=1 时刻最优状态是 0；再看 psi[1, 0]=0，说明 t=0 时刻最优状态是 0。整条最优状态序列是 [0, 0, 0]。如果只看发射概率会觉得 t=1 时刻观测为 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), dtype=int) 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, dtype=int) states[-1] = delta[-1].argmax() for t in range(T-2, -1, -1): states[t] = psi[t+1, states[t+1]] return states, delta

回溯一共就两行核心代码，但这里是最容易写错的地方。psi[t, j] 记录的是"为了在 t 时刻到达状态 j，t-1 时刻应当选哪个状态"。所以在回溯时，要先用当前已经确定好的 t+1 时刻状态去查 psi[t+1, states[t+1]]，得到 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 且 t+1 时刻处于状态 j 的后验概率。前者是"单点"的后验，后者是"相邻两点"的后验。

它们的计算都需要后向变量 beta。后向变量 β_t(i) 表示"给定 t 时刻处于状态 i，从 t+1 到 T 时刻观测到 o_{t+1} 到 o_T 的概率"。注意后向递推的方向是从 T 往 1 走，而且要复用前向算法里存下来的缩放系数 c：

def 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[t+1]] * beta[t+1]) beta[t] /= c[t+1] return beta

缩放后的前向 alpha 和缩放后的后向 beta 逐元素相乘，再进行一行归一化，就得到 gamma。xi 的计算稍复杂，它等于 alpha_t(i) × A[i,j] × B[j,o_{t+1}] × beta_{t+1}(j)，同样要除以 c[t+1] 保持缩放一致。

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_iter=100, tol=1e-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(axis=1, keepdims=True) xi = np.zeros((T-1, N, N)) for t in range(T-1): xi[t] = (alpha[t][:, None] * self.A) * \ (self.B[:, O[t+1]] * beta[t+1])[None, :] xi[t] /= c[t+1] # M步：更新参数 new_pi = gamma[0] new_A = xi.sum(axis=0) / gamma[:-1].sum(axis=0)[:, 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(axis=0) / gamma.sum(axis=0) # 处理训练集中没出现过的观测符号，避免log(0) new_B[new_B == 0] = 1e-6 new_B /= new_B.sum(axis=1, keepdims=True) 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, seed=None): rng = np.random.default_rng(seed) states = np.zeros(T, dtype=int) obs = np.zeros(T, dtype=int) states[0] = rng.choice(self.n_states, p=self.pi) obs[0] = rng.choice(self.n_obs, p=self.B[states[0]]) for t in range(1, T): states[t] = rng.choice(self.n_states, p=self.A[states[t-1]]) obs[t] = rng.choice(self.n_obs, p=self.B[states[t]]) return obs, states

然后用 baum_welch 训练，最后比较学到的 A、B 和真实参数。这里有个坑：HMM 的状态编号是有任意性的，你学到"状态 0"不一定对应真实的"状态 0"，可能整体换了个顺序。所以在验证时要枚举所有状态排列，找到和真实参数最接近的那个排列再比较。N=3 时只有 6 种排列，直接暴力枚举就行；N 大了就随机抽一些排列做近似比较。

我就是靠这个方法抓到一个隐藏很深的问题：我的 alpha、beta、gamma 单独看都合法，但 xi 的缩放系数少除了一个 c[t+1]，导致学习得到的转移矩阵系统性偏差。用真实参数生成数据后一对比，偏差立刻暴露出来。

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 变化曲线画出来，保存下来，以后每次改动代码都能一眼看出收敛行为是否正常。一个健康的训练过程，那条曲线应该平滑上升，最后稳稳地停住；如果它震荡或者下降，赶紧回头查缩放因子，十有八九问题就出在那里。