HMM参数估计详解：从Baum-Welch原理到Python实现-尧图手机网站定制

做隐马尔可夫模型（HMM）系列文章，第一篇我们搞清楚了马尔可夫链和HMM三大基本问题的大框架，也把符号约定好了。今天这篇是第二篇，专门讲参数估计，也就是给定一串观测序列，怎么反推出模型参数 λ=(π, A, B)。这一块是把HMM从“玩具模型”推向真实应用的关键一步，无论是语音识别、股票状态划分还是用户行为建模，都要用到这个环节。适合已经了解HMM基本概念、准备上手写代码的读者，我会把推导逻辑和代码实现一起讲明白，而不是只丢公式。

1. 参数估计到底在求什么：已知观测序列，反推 λ=(π, A, B)

1.1 模型里藏着的三组“看不见的尺子”

HMM 之所以叫“Hidden”，就是因为我们只能看到观测序列 O，看不到真正的状态序列 I。模型里需要估计的，正是控制这个“看不见过程”的三组参数：

初始概率 π：t=1 时刻处于每个隐藏状态的概率。比如第一天的天气，可能是晴天也可能是雨天。
状态转移矩阵 A：从状态 i 转移到状态 j 的概率。今天晴明天雨、今天雨明天依旧下，都由它刻画。
发射概率矩阵 B：在某状态下产生某个观测值的概率。就像晴天更容易让人去散步，雨天更容易让人待在家里打扫。

这三组参数合起来就是 λ=(π, A, B)。参数估计的目标只有一个：找到一组 λ，使得当前观测序列 O 出现的概率 P(O|λ) 最大。换句话说，去猜那根控制黑箱的旋钮应该拧到什么位置。

1.2 一个最小例子把问题钉住

为了后面不迷路，我们先立一个贯穿全文的例子。

假设有两个隐藏状态：晴天(0)、雨天(1)。观测值有三种：散步(0)、购物(1)、打扫(2)。我们只能看到连续几天的行为记录，比如 O=[0,1,1,0,2,1]，但看不到到底哪天下雨哪天晴。我们的任务，就是靠这串行为记录，反推出晴雨之间怎么切换，以及晴天下各种行为出现的概率、雨天下各种行为出现的概率。

这个例子很小，但麻雀虽小五脏俱全。它把参数估计的核心困境暴露得很清楚：

观测序列是“果”，状态序列是“因”。我们只看得到果，却要同时推断因和因的规则。

这就是 HMM 参数估计比普通概率估计难的地方：普通最大似然只猜一组参数，而HMM要同时猜“隐藏过程”和“隐藏过程产生观测”的两层结构。

2. 为什么直接最大似然不行：隐变量的锅

2.1 观测序列的联合概率公式

如果直接做最大似然，我们写出的观测概率是：

P(O|λ) = Σ_I P(I|λ) · P(O|I, λ)

展开就是：

P(O|λ) = Σ_{i1} Σ_{i2} ... Σ_{iT} [ π_{i1} b_{i1}(o1) a_{i1,i2} b_{i2}(o2) ... a_{i_{T-1},i_T} b_{i_T}(o_T) ]

这个式子本身是对的。问题在于，状态序列 I=(i1,...,iT) 的组合数有 N^T 种。就算只有两个状态、一百个时间步，也是 2^100 种路径，直接求和连计算机都扛不住，更不用说求导优化了。

这就像一个人想从北京徒步到广州，理论上走哪条路都能到，但你不可能先枚举完所有路线再挑最短的那条。直接最大似然看起来很美，实际上一根毛都抓不到。

2.2 EM算法切入：把“缺失”的状态序列当作补全数据

既然全量状态序列太可怕，那就换个思路：先猜一组初始参数，然后根据这组参数去猜状态的后验分布，再用“猜出来的状态分布”去更新参数，如此反复迭代。这就是 EM 算法的思想，而 Baum-Welch 算法就是 EM 在 HMM 下的特例。

EM 的核心逻辑是交替执行两步：

E步（期望步）：在固定当前参数 λ 的前提下，计算隐藏变量（状态序列）的后验概率，进而得到一组“期望统计量”，比如某个时刻处于状态 i 的概率、某两个状态相邻转移的概率。
M步（极大步）：把这组期望统计量当作充分统计量，重新估计 π、A、B，使似然函数在这组统计量下尽量大。

直觉上可以这样理解：我们没有完整的路径数据，那就让当前参数帮我们“脑补”每条路径的贡献权重，然后把这些权重带进统计计数里。反复做，参数就会越猜越准。

这一节的意义在于，它把“为什么会推导出前向后向算法”这个问题给回答了——Baum-Welch 需要一个计算后验概率的抓手，而抓手就是前向算法和后向算法。

3. 前向与后向算法：Baum-Welch的两块地基

3.1 前向变量 α_t(i)：从1到t“看得见”的贡献

在实现 Baum-Welch 之前，必须先掌握前向和后向两个递推算法。它们不仅用于计算 P(O|λ)，也是后面统计量 γ 和 ξ 的直接组成部分。

前向变量定义为：

α_t(i) = P(o1, o2, ..., ot, i_t = q_i | λ)

意思是：在参数 λ 下，观测到前 t 个结果并且 t 时刻落在状态 i 的联合概率。

递推关系非常优雅：

初始化：α_1(i) = π_i · b_i(o1)

递推：α_{t+1}(j) = [ Σ_i α_t(i) · a_ij ] · b_j(o_{t+1})

终止：P(O|λ) = Σ_i α_T(i)

它最大的价值是避免了枚举全部路径。α_t(i) 把“到 t 时刻为止、处于状态 i”的所有可能路径的概率合并成了一个值，每一层只需要 N^2 次乘法，复杂度从指数级降到 O(N^2 T)。

3.2 后向变量 β_t(i)：从t+1到T“未来”的贡献

后向变量是前向的镜像，定义成：

β_t(i) = P(o_{t+1}, ..., o_T | i_t = q_i, λ)

也就是已知 t 时刻在状态 i，未来从 t+1 到 T 的观测序列出现的概率。

递推从后往前：

初始化：β_T(i) = 1

递推：β_t(i) = Σ_j a_ij · b_j(o_{t+1}) · β_{t+1}(j)

前向和后向合起来，能得到一个非常实用的量：

P(i_t = q_i | O, λ) = α_t(i) · β_t(i) / P(O|λ)

这个量表示“在已知全部观测序列的前提下，某个时刻正处于状态 i 的概率”。这正是参数估计最需要的信息。

3.3 为什么薛定谔的“软归属”是关键

很多初学者会疑惑：既然我看不到状态，那 Baum-Welch 到底怎么更新矩阵？

关键就在这里：我们用“软归属”，而不是“硬归属”。

硬归属是：每个时刻强行指定一个状态，然后按离散计数更新参数。软归属则是：每个时刻按概率比例分配权重，比如 t 时刻有 0.7 的概率在晴天、0.3 的概率在雨天，那么统计晴天贡献时只记 0.7。

这种软归属正是 EM 的精髓。它让每个观测数据都对多个状态有影响，而不是武断地选一个。gamma 和 xi 就是软归属的两个具体指标，下一节详细拆解。

4. Baum-Welch重估公式：怎么把计数变成参数

4.1 γ_t(i)与ξ_t(i,j)的计算

接下来进入核心推导。为了重估参数，我们需要两个统计量。

第一个是 γ_t(i)，表示给定全部观测时，t 时刻处于状态 i 的概率：

γ_t(i) = P(i_t = q_i | O, λ) = α_t(i) · β_t(i) / P(O|λ)

第二个是 ξ_t(i,j)，表示给定全部观测时，t 时刻在状态 i、t+1 时刻转移到状态 j 的概率：

ξ_t(i,j) = P(i_t = q_i, i_{t+1} = q_j | O, λ)

它的计算稍微复杂一点，但图像很清晰。考虑 t 时刻在状态 i，t+1 时刻在状态 j，并且完整观测到 O 的联合概率，可以拆成四段：

α_t(i) · a_ij · b_j(o_{t+1}) · β_{t+1}(j)

也就是说：前 t 个观测由 α_t(i) 负责，从 i 跳到 j 用 a_ij，跳到 j 后产生观测 o_{t+1} 用 b_j(o_{t+1})，剩下从 t+2 到 T 的观测由 β_{t+1}(j) 负责。最后除以 P(O|λ) 归一化：

ξ_t(i,j) = α_t(i) · a_ij · b_j(o_{t+1}) · β_{t+1}(j) / P(O|λ)

γ 和 ξ 之间有个简单关系：

γ_t(i) = Σ_j ξ_t(i,j)

因为不管 t+1 时刻跑到哪个状态，只要 t 时刻处于状态 i，都被 γ_t(i) 算进去。

4.2 三个重估公式

有了 γ 和 ξ，参数重估就是“期望计数除以总计数”：

π_i 的新估计值：π_i = γ_1(i)。也就是第一时刻处于状态 i 的期望概率。
状态转移矩阵 A 的新估计值：

a_ij = Σ_{t=1}^{T-1} ξ_t(i,j) / Σ_{t=1}^{T-1} γ_t(i)

发射概率矩阵 B 的新估计值：

b_j(k) = Σ_{t: o_t = v_k} γ_t(j) / Σ_{t=1}^{T} γ_t(j)

这三个公式是 Baum-Welch 的全部秘密。π 只看第一时刻，A_ij 看“从 i 跳到 j 的期望次数”除以“从 i 出发的总期望次数”，B_j(k) 看“状态 j 时观测到 v_k 的期望次数”除以“状态 j 出现的总期望次数”。

4.3 直觉：软计数除以软计数，为什么迭代收敛

重估公式的本质是“用期望代替真实计数”。设想我们真的知道状态序列，那 A 的估计无非就是“从 i 到 j 的转移次数”除以“从 i 出发的总次数”。现在状态序列未知，我们就把每个候选转移乘上它的概率权重再累加，得到“软转移次数”，分母同理。这就是为什么 Baum-Welch 常常被叫做“soft counting”。

从优化的视角看，每一轮 M 步都能保证：

P(O|λ_new) ≥ P(O|λ_old)

也就是说对数似然单调不减。这正是 EM 算法的性质。但要注意，它保证的是收敛到局部极大值，不保证全局最优。所以初始化很重要，后文会专门说。

5. 实战中躲不开的数值下溢：缩放与对数域处理

5.1 下溢为什么必然发生

理论公式推完了，如果你直接照着实现，跑几个长序列就会看到一串 NaN。原因其实很简单：

α_t(i) 是所有路径概率的加权和，每个概率都是些小于 1 的数连乘，T 一大，α_t(i) 就指数级趋向于 0。比如 T=100 的序列，即使每一步概率只有 0.5 左右，连乘结果也会低于浮点数的表示范围，变成 0。

这不是实现 bug，而是数学上必然发生的问题。所以正确实现必须做数值稳定处理，最常见的就是缩放法。

5.2 最实用的缩放法

缩放法的思路很简单：每一步算完 α_t(i) 后，统一除以一个缩放系数 c_t，让这一层所有 α_t(i) 的和等于 1。

具体操作：

每算完一层 α_t，令 c_t = Σ_i α_t(i)，再把每个 α_t(i) 都除以 c_t。
后向算法递推时，对 β_t(i) 也除以同样的 c_t，保持配套。

这样缩放并不会改变 γ_t(i)，因为 γ_t(i) 是 α_t(i)β_t(i)/P，分子分母同时被缩放，数值自动抵消。

整套 P(O|λ) 可以直接用缩放系数的连乘积恢复：

P(O|λ) = 1 / (c_1 · c_2 · ... · c_T)

这个技巧非常实用，代码里我基本都这么写，既简单又避免了 log-sum-exp 的额外复杂度。

5.3 收敛判定的细节与初始化

迭代收敛判定有几种常见做法：

对数似然差：计算 log P(O|λ)，相邻两轮差值小于某个阈值就停止。
参数差：看新旧 λ 的最大绝对差是否小于阈值。

我个人的习惯是把两者结合：先看对数似然差，再看转移矩阵是否稳定。单纯看参数差有时候会提前停，因为参数变化很小但似然还在上升；只靠似然差在局部平坦区又可能震荡。

初始化方面，几个经验：

不要让某个转移概率初始为 0。如果初始状态里 a_ij = 0，迭代过程中它很可能一直是 0，永远无法恢复。给一个小的非零值更安全。
发射矩阵也类似，避免硬 0，可以用均匀或小随机数初始化。
多随机起点跑几次，选对数似然最大的结果。HMM 的似然面非凸，单次初始化很容易掉进局部最优。

6. 跑一个完整的天气-心情参数估计Demo

6.1 数据、初始化和Python实现

回到晴雨的例子。观测序列我取：

O = [0, 1, 1, 0, 2, 1]

也就是：散步、购物、购物、散步、打扫、购物。一共 6 个时间步，很短，适合演示迭代过程。

初始参数取对称值，避免一开始就偏心：

π = [0.5, 0.5]
A = [[0.6, 0.4], [0.3, 0.7]]
B = [[0.5, 0.3, 0.2], [0.2, 0.5, 0.3]]

直接贴一段可运行的 Python 代码，实现带缩放的 Baum-Welch：

import numpy as np def baum_welch(O, pi, A, B, max_iter=50, tol=1e-4): O = np.array(O, dtype=int) T = len(O) N = len(pi) M = B.shape[1] for it in range(max_iter): # E step alpha = np.zeros((T, N)) beta = np.zeros((T, N)) scale = np.zeros(T) alpha[0] = pi * B[:, O[0]] scale[0] = alpha[0].sum() alpha[0] /= scale[0] for t in range(1, T): alpha[t] = (alpha[t-1] @ A) * B[:, O[t]] scale[t] = alpha[t].sum() alpha[t] /= scale[t] beta[-1] = 1.0 for t in range(T-2, -1, -1): beta[t] = A @ (B[:, O[t+1]] * beta[t+1]) beta[t] /= scale[t] gamma = alpha * beta gamma /= gamma.sum(axis=1, keepdims=True) xi = np.zeros((T-1, N, N)) for t in range(T-1): numerator = alpha[t][:, None] * A * B[:, O[t+1]][None, :] * beta[t+1][None, :] xi[t] = numerator / numerator.sum() # M step new_pi = gamma[0] new_A = xi.sum(axis=0) / gamma[:-1].sum(axis=0)[:, None] new_B = np.zeros_like(B) for k in range(M): mask = (O == k) new_B[:, k] = gamma[mask].sum(axis=0) new_B /= gamma.sum(axis=0)[None, :].T # log likelihood log_likelihood = -np.log(scale).sum() pi, A, B = new_pi, new_A, new_B print(f"iter {it+1}: log-likelihood = {log_likelihood:.4f}") return pi, A, B, log_likelihood pi = np.array([0.5, 0.5]) A = np.array([[0.6, 0.4], [0.3, 0.7]]) B = np.array([[0.5, 0.3, 0.2], [0.2, 0.5, 0.3]]) pi_final, A_final, B_final, llh = baum_welch([0, 1, 1, 0, 2, 1], pi, A, B) print("pi =", pi_final) print("A =", A_final) print("B =", B_final) print("log-likelihood =", llh)

这里我用了缩放后的 alpha 和 beta，每一步都归一化，确保数值不爆。gamma 和 xi 的概率由归一化后的量直接计算，最后的对数似然用 scale 累乘恢复。

6.2 迭代结果展示

跑上面代码，20 轮内对数似然就会趋于稳定。一次典型运行的输出大致如下（不同随机初始化可能略有差异，但整体规律一致）：

参数	初始化	迭代收敛后（示意）
π	[0.5, 0.5]	[0.44, 0.56] 左右
A	[[0.6, 0.4], [0.3, 0.7]]	对角线元素略升高，约 [[0.70, 0.30], [0.25, 0.75]]
B	[[0.5, 0.3, 0.2], [0.2, 0.5, 0.3]]	状态0更偏向观测0，状态1更偏向观测1，约 [[0.65, 0.23, 0.12], [0.10, 0.63, 0.27]]

注意：这只是示意数值，实际运行会因随机初始化而浮动，但方向是一致的——模型会倾向于用“状态自身持续性”和“观测偏好”来解释训练序列。

这个短例子最能说明一件事：哪怕观测序列只有 6 个点，Baum-Welch 也能给出看起来合理的参数，但这些参数不一定可靠。序列越短，估计方差越大，后面实战中要格外小心。

6.3 解读结果时容易犯的错

第一，别把“参数收敛”等同于“参数正确”。HMM 的似然面有多个峰，不同初始化可能收敛到不同的局部最优，这很正常。

第二，状态存在置换模糊性。也就是说，如果我把状态 0 和状态 1 的标签互换，观测概率的值不变。你跑出来的“状态0是晴天”，换一个随机种子可能就变成“状态0是雨天”。所以分析结果时要看结构，而不是死抠标签。

第三，关于收敛阈值。T 很小的时候，对数似然的提升空间有限，阈值设得太小会导致白跑很多轮；实用中我一般设 1e-4 到 1e-5，并且设置最大迭代次数兜底。

7. 我自己在工程中使用Baum-Welch的几条心得

7.1 短序列和小样本的陷阱

如果你只有一个很短的观测序列，Baum-Welch 估计出的参数方差会很大。拿上面的 6 点序列做例子，参数更新几乎完全跟着观测频率走，极端情况下发射概率会被某个偶然而不是规律主导。

我踩过的坑：某次项目里只有几十条用户行为序列，每条还特别短，直接跑 Baum-Welch 后得到的转移矩阵几乎全是噪声，后续的 Viterbi 解码结果一眼假。后来改成把多条序列拼接成一个长序列、或者并行做多序列联合估计，情况才稳定。多序列联合时，只需要把 gamma 和 xi 的累和分别跨序列相加即可，公式几乎不用改。