HMM隐马尔可夫模型Python实现：前向、维特比与Baum-Welch算法详解

如果你已经被 HMM（Hidden Markov Model）的公式绕得有点晕，却又想看到它究竟如何用代码落地，这篇笔记应该正好适合你。我会绕过一上来就铺满的数学推导，直接用一套完整可运行的 Python 实现，把 HMM 最核心的几件事讲清楚：给定模型怎么计算观测序列概率、怎么解码最优隐藏状态路径、怎么在没有标注的情况下从数据里估计参数。文章更适合刚啃完理论、想动手复现的读者，也适合做序列建模、但被矩阵维度绕晕的算法工程师。

很多公开教程在讲隐马尔可夫模型时都会从概率图、状态转移、发射概率一路推导，推导当然重要，可不少人卡在"公式我懂了，代码就是不会写"。我这次把思路反过来：直接把代码摊开，一段一段解释每个数组在干什么、每个矩阵乘法在算什么，再回头去对照公式，你会发现很多之前想不通的地方一下就顺了。整篇文章只需要 NumPy，不依赖任何深度学习框架。

为了便于验证，我使用一个很经典的天气与活动例子：隐藏状态是"晴天/雨天"两种，观测行为是"散步/宅家/购物"三种。状态少、符号少，每一步都可以手算核对，定位错误时也不会被复杂数据干扰。如果你在做语音、文本或生物序列相关任务，把这里的 K 和 D 换成自己的状态数和符号数就行，核心逻辑完全一致。

1. HMM建模：先把 pi、A、B 三个矩阵在代码里摆正

1.1 三个参数两句话讲完

HMM 本质上是一台双层的随机机器：底层是一条看不见的马尔可夫链，链上的每个状态会按照一个概率分布向外发射一个可观测的符号。我们要建模的就是三组概率参数。pi 是初始状态概率向量，长度 K，表示序列第一个时刻处于每个隐藏状态的概率；A 是状态转移矩阵，形状 K×K，元素 A[i][j] 表示当前时刻在状态 i、下一时刻跳到状态 j 的概率；B 是发射概率矩阵，形状 K×D，元素 B[i][o] 表示在状态 i 下观测到符号 o 的概率。

代码里我用一个 HMM 类把它们串起来。初始化时用随机数并归一化，确保每一行都是合法的概率分布：

import numpy as np class HMM: def __init__(self, K, D, random_seed=None): self.K = K self.D = D if random_seed is not None: np.random.seed(random_seed) self.pi = np.random.rand(K) self.pi /= self.pi.sum() self.A = np.random.rand(K, K) self.A /= self.A.sum(axis=1, keepdims=True) self.B = np.random.rand(K, D) self.B /= self.B.sum(axis=1, keepdims=True)

这段代码在干什么很简单：rand 生成 0 到 1 的均匀随机数，除以行和是为了让概率和为 1。random_seed 只是为了复现结果，真正项目里更推荐用可复现的随机源，不然每次训练结果都不一样，定位问题会很痛苦。注意A.sum(axis=1, keepdims=True)的 keepdims 不能丢，否则广播形状会对不上，可能得到完全错乱的矩阵。

1.2 直接写一组手工参数

为了避免调试时被随机参数干扰，后面所有手工验证我都用下面这组参数：

pi = np.array([0.6, 0.4]) A = np.array([[0.7, 0.3], [0.4, 0.6]]) B = np.array([[0.5, 0.2, 0.3], [0.1, 0.6, 0.3]])

含义是：晴天比雨天稍常见；晴天之后 70% 还是晴天，雨天之后 60% 还是雨天；晴天最可能散步，雨天最可能宅家。这里方向非常关键：A 的行是当前状态，列是下一状态，也就是A[i][j] = P(q_{t+1}=s_j | q_t=s_i)，和很多教科书的花括号公式一致。代码里写矩阵乘法时，行方向搞反会导致整个模型完全错乱，建议动手前先在注释里标清楚方向。

1.3 观测序列必须整数编码

obs 数组里的每个元素必须是一个整数，并且满足 0 <= obs[t] < D。如果你的原始数据是字符串，要先建立映射：

symbol2id = {"散步": 0, "宅家": 1, "购物": 2} obs = np.array([symbol2id[x] for x in ["散步", "宅家", "购物", "散步"]])

这里踩过的坑包括：D 定义成 2，但 obs 里有 2 导致 B[:, 2] 越界；或者字符串没映射直接放进代码报 TypeError。还有一点：如果是连续观测（比如语音的浮点特征），就不能用离散发射矩阵，得换成高斯发射或混合高斯发射，那是另一个话题，本文先把离散版本讲透。

2. 前向算法代码实现：从公式到矩阵乘法的每一步

2.1 从枚举到动态规划

前向算法解决"观测序列概率"问题：给定模型 λ=(pi,A,B)，计算 P(O|λ)。朴素做法是枚举所有可能的状态序列，再把每条路径的生成概率加起来。状态数 K、序列长度 T 时，有 K^T 条路径，直接枚举显然不现实。前向算法用动态规划把复杂度降到 O(T*K^2)。在 HMM 类里加一个 forward 方法：

def forward(self, obs): T = len(obs) alpha = np.zeros((T, self.K)) alpha[0] = self.pi * self.B[:, obs[0]] for t in range(1, T): alpha[t] = self.B[:, obs[t]] * (alpha[t - 1] @ self.A) return alpha, alpha[-1].sum()

逐行拆开看。第一行self.pi * self.B[:, obs[0]]是按元素相乘：pi 是初始状态概率，B[:, obs[0]] 是每个隐藏状态发射第一个观测的概率，乘出来的每个元素就是"序列一开始进入状态 i，并且观测到 obs[0]"的联合概率。第二行到循环体是递推，核心就一句alpha[t - 1] @ self.A，这个矩阵乘法的结果是一个长度 K 的向量，第 j 个元素等于 Σ_i alpha[t-1][i] * A[i][j]，可以理解为"把所有前一时刻状态 i 的累计概率，按照转移概率搬到新状态 j 上"。最后再乘self.B[:, obs[t]]，把 t 时刻看到 obs[t] 这件事加进去。

2.2 用"账户余额"理解 alpha

我说个更好记的类比：alpha[t] 就像每个隐藏状态在 t 时刻的账户余额。前一天晚上每个状态各有一笔余额；过了一夜，余额按转移矩阵 A 重新分配给今天可能出现的状态；等今天观测结果公布，再根据发射概率把余额进一步调整。到了序列末尾，把所有状态的余额相加，就是整个观测序列出现的总概率。这个方法的价值是它避免了枚举路径，所有中间结果都是软概率，也就是同时考虑了所有可能的隐藏状态路径。

有一点必须提前说：每一步都在乘小于 1 的小数，序列稍微一长，浮点数就可能下溢到 0。0.7 的 200 次方大约就成了 1e-31，而真实模型的概率可能比这还小，所以 alpha 矩阵在长序列下经常全零。这个问题后面单独讲。

2.3 工程上必须处理的缩放

解决下溢的经典办法是每一步把 alpha 归一化，让和变成 1，同时单独记录每一步的缩放系数。观测序列的总概率等于所有缩放系数的乘积，取对数后就是 log 似然。代码可以这样写：

def forward_scaled(self, obs): T = len(obs) alpha = np.zeros((T, self.K)) scales = np.zeros(T) alpha[0] = self.pi * self.B[:, obs[0]] scales[0] = alpha[0].sum() alpha[0] = alpha[0] / scales[0] for t in range(1, T): alpha[t] = self.B[:, obs[t]] * (alpha[t - 1] @ self.A) scales[t] = alpha[t].sum() alpha[t] = alpha[t] / scales[t] return alpha, scales, np.log(scales).sum() def log_likelihood(self, obs): return self.forward_scaled(obs)[2]

log 似然的好处是数值稳定，而且可以用它做收敛判断和模型比较。注意返回的 alpha 是缩放后的版本，已经不是原始概率，但它保留了相对大小，后续 Baum-Welch 训练和 Viterbi 都在这种缩放版本上继续推。

3. 后向算法和维特比解码：从总概率到最优路径

3.1 后向算法的方向最容易错

后向算法的定义是：beta[t][i] 表示已知 t 时刻处于状态 i，后续观测 obs[t+1:] 出现的概率。t=T 时没有后续观测，所以最后一行全是 1。

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

最容易写错的就是最后一行的矩阵方向。把公式写开就是 beta[t][j] = Σ_i A[j][i] * B[i][obs[t+1]] * beta[t+1][i]。这里 A 的第 j 行与向量相乘，正好得到从状态 j 出发的所有路径。如果你随手转置成self.A.T @ ...，形状不会报错，但结果完全不对，而且这种 bug 很难靠看打印发现。我的检查习惯是每一步都盯形状：beta[t] 必须是长度 K 的向量，中间的self.B[:, obs[t+1]] * beta[t+1]也是长度 K 的向量，所以self.A必须放左边。

3.2 维特比：把求和换成取最大值

前向和后向算的是所有路径的总概率，维特比则要回答：给定观测序列，最可能的隐藏状态路径是哪一条。思路同样是动态规划，区别是把"求和"变成"取最大值"，并且额外记录每一步是从哪个状态跳过来的。代码：

def viterbi(self, obs): T = len(obs) delta = np.zeros((T, self.K)) psi = np.zeros((T, self.K), dtype=int) delta[0] = self.pi * self.B[:, obs[0]] psi[0] = -1 for t in range(1, T): trans = delta[t - 1][:, None] * self.A delta[t] = self.B[:, obs[t]] * trans.max(axis=0) psi[t] = trans.argmax(axis=0) 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[-1].max()

trans的形状是 K×K，trans[i][j] = delta[t-1][i] * A[i][j]，表示"前一步最优走到状态 i 的得分，再转移到状态 j 的得分"。trans.max(axis=0)是对每个目标状态 j 取最大值，也就是从任意 i 跳到 j 的最佳得分；trans.argmax(axis=0)记录这个最佳得分来自哪个 i，塞进 psi 矩阵。

3.3 回溯路径为什么会写错

psi[t][j] 的含义是：在 t 时刻到达状态 j 的最优路径中，t-1 时刻的状态是谁。所以最终回溯要从后往前：先找到整个序列最后一个时刻得分最高的状态，然后循环往回找前驱。循环里的写法是states[t] = psi[t + 1, states[t + 1]]，不是psi[t][states[t - 1]]。我把这个错踩过一次，当时路径总是反着来或者整体移位，定位了很久才发现回溯方向写反。

未做缩放的维特比在长序列下也会下溢。delta 全部变成 0 时，argmax等于随机挑了一个状态，回溯出来的路径毫无意义。工程实现里一般直接用 log 域版本。

3.4 log 域维特比是工程标配

def viterbi_log(self, obs): T = len(obs) log_pi = np.log(self.pi) log_A = np.log(self.A) log_B = np.log(self.B) delta = np.zeros((T, self.K)) psi = np.zeros((T, self.K), dtype=int) delta[0] = log_pi + log_B[:, obs[0]] for t in range(1, T): tmp = delta[t - 1][:, None] + log_A delta[t] = log_B[:, obs[t]] + tmp.max(axis=0) psi[t] = tmp.argmax(axis=0) 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[-1].max()

log 域把乘法变成加法，数值上稳很多。注意log(0)会直接报错，如果发射概率矩阵里有严格 0，建议在训练或初始化时加一个很小的平滑项，比如B[B == 0] = 1e-10再归一化。不要为了省事写成np.log(self.B + 1e-10)，虽然也能用，但会让原本精确的模型参数被无差别平移，后面 EM 表现会受影响。

4. 参数估计：Baum-Welch 算法的完整实现与训练工程

4.1 EM 思想在 HMM 里长什么样

没有隐藏状态标签时，怎么从无标注观测序列里学出 pi、A、B？最经典的就是 Baum-Welch，本质上就是 EM 算法。思路是先假定当前参数是对的，用它去"猜测"每个时刻最可能处于哪个状态、每个相邻时刻最可能在哪两个状态之间转移；但这些猜测不是硬划分，而是带着后验概率的软计数。拿到这些软计数后，重新估计参数，再重复。每一步迭代都能保证观测序列的 log 似然单调上升，最终收敛到局部最优。这也是为什么 4.4 节要强调多随机初始化。

4.2 gamma 和 xi 的计算

定义两个核心变量：gamma_t(i) 表示已知整个观测序列时，t 时刻处于状态 i 的后验概率；xi_t(i,j) 表示 t 时刻处于 i、t+1 时刻处于 j 的后验概率。计算时用到前向和后向结果：

alpha, _ = self.forward(obs) beta = self.backward(obs) ll = alpha[-1].sum() gamma = alpha * beta / ll xi = np.zeros((T - 1, self.K, self.K)) for t in range(T - 1): xi[t] = alpha[t][:, None] * self.A * self.B[:, obs[t + 1]] * beta[t + 1] / ll

alpha[t][i] 乘 beta[t][i] 再除以 ll，就是"t 时刻处于 i 且生成完整观测序列"的后验概率。xi 的公式更直接：alpha[t][i] 管住前缀，A[i][j] 管住转移，B[j][obs[t+1]] 管住下一时刻发射，beta[t+1][j] 管住后缀，最后除以 ll 做归一。每个变量的名字跟着经典记号走，因为大批论文和工具包都用这个记号，熟了对读资料很有帮助。

4.3 参数更新的代码实现

得到 gamma 和 xi 后，M 步的更新十分直观：

def baum_welch(self, obs, epochs=100, tol=1e-6): T = len(obs) last_ll = -np.inf for epoch in range(epochs): alpha, _ = self.forward(obs) beta = self.backward(obs) ll = alpha[-1].sum() gamma = alpha * beta / ll xi = np.zeros((T - 1, self.K, self.K)) for t in range(T - 1): xi[t] = alpha[t][:, None] * self.A * self.B[:, obs[t + 1]] * beta[t + 1] / ll new_pi = gamma[0].copy() new_pi /= new_pi.sum() denom_A = gamma[:-1].sum(axis=0) new_A = xi.sum(axis=0) / denom_A[:, None] new_B = np.zeros_like(self.B) for k in range(self.D): mask = (np.array(obs) == k) new_B[:, k] = gamma[mask].sum(axis=0) / gamma.sum(axis=0) self.pi = new_pi self.A = new_A self.B = new_B if abs(ll - last_ll) < tol: break last_ll = ll return ll

new_pi 是第一个时刻的 gamma，因为初始状态只能出现在 t=0。new_A 的分子是所有相邻时刻从 i 到 j 的期望转移次数，分母是除了最后时刻以外处于状态 i 的期望总次数。new_B 稍微麻烦一点：对于每一个观测符号 k，把所有观测到 k 的时刻对应的 gamma 行汇总，再除以状态 j 在所有时刻的期望总次数。分母不能错误写成gamma[mask].sum(axis=0)，那会丢掉"状态 j 在非 k 时刻的占比"，导致发射概率矩阵每一行不归一。

4.4 训练时的工程细节和避坑思路

这个实现是教学版，没有做缩放，序列短没问题，序列一长就可能在除法里冒出 NaN。工程上要做的第一件事就是加缩放：用缩放后的 alpha 和对应的缩放系数去计算 gamma 与 xi，或者干脆整个写到 log 空间里重算。比较稳妥的做法是每次迭代后检查np.isfinite(ll)，不是有限值就立刻停，别让 NaN 一路传染下去。

另一个工程习惯是多起点随机初始化。EM 只保证收敛到局部最优，一次随机初始化可能得到一个很差的模型。我通常写一个外层循环：

def train_with_restarts(obs, K, D, n_restarts=5, epochs=100): best_ll, best_model = -np.inf, None for seed in range(n_restarts): model = HMM(K, D, random_seed=seed) ll = model.baum_welch(obs, epochs=epochs) if ll > best_ll: best_ll, best_model = ll, model return best_model

还有一个隐蔽问题：EM 更新里隐含"E 步必须用旧参数算 gamma/xi，M 步结束后才更新"。如果你在计算 xi 的循环里顺手改了self.A，或者提前修改了self.B，迭代的单调性就没了。我遇到"似然不升反降"时，第一反应就是看变量是否被复用和翻新，而不是怀疑公式。

5. 代码调试与常见问题速查

5.1 用暴力枚举验证前向算法

HMM 代码看起来简单，但方向、维度、下标一错就是静默失败。最可靠的验证办法是拿一个小模型暴力枚举所有状态路径，和动态规划结果对比。状态只有 2、序列长度只有 4 时，路径只有 16 条，手算或者暴力枚举都毫无压力。

from itertools import product def brute_force_hmm(pi, A, B, obs): K = len(pi) total = 0.0 for states in product(range(K), repeat=len(obs)): p = pi[states[0]] for t, o in enumerate(obs): p *= B[states[t], o] if t > 0: p *= A[states[t - 1], states[t]] total += p return total

测试时把手动参数塞进 model：

obs_test = [0, 1, 2, 0] model = HMM(2, 3) model.pi = np.array([0.6, 0.4]) model.A = np.array([[0.7, 0.3], [0.4, 0.6]]) model.B = np.array([[0.5, 0.2, 0.3], [0.1, 0.6, 0.3]]) print(model.forward(obs_test)[1]) print(brute_force_hmm(model.pi, model.A, model.B, obs_test))

输出的两个值应该在浮点误差范围内相等。如果不等，优先看 B 的维度、A 的方向、obs 的编码；这个检查能过滤掉八成新手 bug。

5.2 数据格式速查表

很多运行时报错其实和数据格式有关。我整理一张速查表：

5.3 下溢问题的现场排查

你可以通过几个信号识别下溢：model.forward(obs)[1]返回 0，可你手工估算概率不至于那么小；viterbi的分数在不同运行里没有任何变化但不是你预期的路径；log 似然从某个 epoch 开始变成-inf。碰到这些不要急着调参，先把概率计算从普通域切到 log 域再跑一遍。同样一组参数，log 域版本一般能给出稳定的结果。

5.4 EM 似然不行为时的检查路径

如果训练时 log 似然不升反降，第一个要查的是 E 步和 M 步的变量隔离问题：alpha、beta、gamma、xi 必须全部用旧参数计算完毕，再更新 self.pi、self.A、self.B。我在代码里专门用了 new_pi、new_A、new_B 三个新变量，就是防止把 M 步更新提前引入了 E 步。

第二个要查的是计算似然时有没有复用被缩放的 alpha。如果你把forward_scaled返回的 alpha 当成原始概率，直接拿去算 ll，ll 就会变成接近 1 的值而不是真实似然。教学版forward返回的才是原始概率，工程版要自己维护好缩放系数并取对数似然。

第三个要查的是 B 的更新是否让每一行重新归一化。发射概率矩阵每一行必须仍是概率分布，更新完最好立刻打印model.B.sum(axis=1)，如果某一行不是 1，说明分母或广播写错了。

5.5 完整可运行 Demo 与最后的建议

把前面所有方法放进同一个 HMM 类后，可以用一段小数据跑完整流程：

obs_seq = np.array([0, 1, 2, 0, 1, 0, 2, 1, 0, 1, 1, 2, 0]) model = HMM(2, 3, random_seed=42) print("训练前 log P(O):", model.log_likelihood(obs_seq)) best_states, score = model.viterbi_log(obs_seq) print("Viterbi 路径:", best_states, "log score:", score) model.baum_welch(obs_seq, epochs=50) print("训练后 log P(O):", model.log_likelihood(obs_seq)) print("训练后的 A:") print(model.A)

你可能会发现训练后的 log P(O) 比训练前高不少，这是正常现象，说明 EM 在持续改善模型。但不代表这是全局最优，用不同随机种子多跑几次，挑 log 似然最高的模型才算完整流程。观测序列太短或状态数设置不合理时，模型还可能收敛到重复状态（两个状态学得一模一样），这时候要么减 K，要么增大数据量。

这套代码我建议你先直接用玩具数据跑通，再动手改参数。改完任何一部分，都用 5.1 节那种暴力枚举去校验一下，否则你很难分清是模型参数的问题还是代码实现的问题。把 HMM 的三大计算问题和这套工程习惯练熟之后，再去读更大规模序列模型的源码，会轻松不少。