企业官网建设流程全解析

如果你已经把前向算法和后向算法啃下来了，那HMM参数估计这块，也就是Baum-Welch算法，其实就剩一层窗户纸了。

标题写的是第2讲，默认你已经知道隐马尔可夫模型是什么：一条状态马尔可夫链，每个状态再按一定概率发射观测值。今天要面对的问题是：给定一堆观测序列，怎么把一个完整的模型参数λ=(A, B, π)给估出来。

这问题严格来说属于EM算法在离散状态空间上的一个特例。很多人学到这里被满屏的希腊字母劝退，其实核心思想在生活里很常见——你只能看到结果，过程是隐藏的，那就只能用结果去反推过程的统计规律。比如你看不见天气，但能从一个人每天是否带伞推断出天气变化的规律，这就是HMM参数估计在做的事。

这篇内容会让你搞清楚两件事：第一，为什么参数估计不能像普通马尔可夫链那样直接计数；第二，Baum-Welch具体怎么迭代、代码怎么落地、实际运行时有哪些看不见的坑。适合正在学HMM、准备做序列建模，或者想把EM思想落到真实代码里的朋友参考。

1. 参数估计问题到底难在哪

1.1 状态可见与不可见的本质差别

如果HMM里的状态序列Z=(z₁, z₂, ..., z_T)是直接可见的，那参数估计就退化成一个简单的统计问题。

比如我有N=3个状态，M=2个观测符号，跑了T=1000步，每一步我都知道自己在哪个状态、发射了什么观测。那么转移矩阵A和发射矩阵B的估计就是纯粹的频率计算：

• π_i的估计 = 第一步处于状态i的次数 / 总序列数
• a_ij的估计 = 从状态i转移到状态j的次数 / 从状态i出发的总次数
• b_j(k)的估计 = 在状态j下观测到k的次数 / 状态j出现的总次数

这就是最大似然估计，没有任何技巧难度，统计完事。

问题在于，真实场景下状态序列是不可见的。我能看到的只有观测序列O=(o₁, o₂, ..., o_T)，比如一段语音的声学特征、一连串股票涨跌符号、一列基因碱基序列。而导致这些观测的隐状态路径，可能有N^T种之多。

这直接导致似然函数变成了一个没法化简的求和：

P(O|λ) = Σ_Z P(O|Z,λ) · P(Z|λ)

里面套着一个对所有长度为T的隐状态路径求和的运算。

1.2 对数似然的log-sum困境

做最大似然估计时，我们一般喜欢对似然函数取对数，因为连乘变成连加，求导方便。但这里有个尴尬的地方：虽然P(O|Z,λ)和P(Z|λ)本身都是连乘形式，可以取对数化简，可外层还有个对Z的求和。

取对数之后，log里面依然是一个巨大的求和式。log(a+b+c)没法拆成log(a)+log(b)+log(c)，这个数学结构决定了直接对λ求偏导、令导数为零，得不到闭式解。

这是HMM参数估计和普通马尔可夫链的本质区别。它不是一个简单的计数问题，而是一个在隐变量存在条件下的最优化问题。

1.3 EM算法搭的桥

解决这类含隐变量问题的标准工具就是EM算法，全称Expectation-Maximization算法。HMM里专门有个历史悠久的名字叫Baum-Welch算法，其实就是EM算法在HMM参数估计上的具体实现。

EM的核心思路是分两步走：

• E步（期望步）：利用当前参数λ的旧值，计算隐状态的后验分布，也就是在已知观测序列的前提下，推算每个时刻最可能处于哪些状态。
• M步（最大化步）：把这些后验概率当成“软计数”来用，重新估计一组让期望似然最大化的新参数λ。

这实际上是在玩一个“先猜后修正”的游戏。初始参数可以不那么准，但只要迭代足够多次，似然值就会单调上升，最终收敛到某个局部最优解。注意我说的是局部最优，这点后面还会专门聊。

2. E步的核心：前向后向概率与两个关键辅助量

2.1 从前向到后向的完整拼图

既然E步要求“在已知观测序列的条件下计算隐状态的后验分布”，那就绕不开前向算法和后向算法这两个基础工具。

前向概率α_t(i)定义为：t时刻处于状态i，且观测到o₁到o_t这一段的联合概率。它可以从t=1开始逐步递推：

α₁(i) = π_i · b_i(o₁)

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

后向概率β_t(i)定义为：在t时刻处于状态i的前提下，从t+1时刻开始一直观测到o_T的概率。它从序列末尾往前递推：

β_T(i) = 1

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

前向和后向配合起来，就能得到P(O|λ)的两种等价表达：

P(O|λ) = Σ_i α_T(i) = Σ_i α_t(i) · β_t(i)（任意时刻t）

这个等式很关键。它意味着前向和后向在任意时间点上都能相遇并拼出整个观测序列的似然值。

2.2 γ_t(i)：单点状态占据概率

定义γ_t(i)为t时刻处于状态i的后验概率。根据贝叶斯公式：

γ_t(i) = P(z_t=i|O,λ) = [α_t(i) · β_t(i)] / P(O|λ)

直观理解就是：t时刻“经前向传播过来的信息”和“经后向传播过来的信息”在状态i上碰头，归一化后就是该时刻处于状态i的可能性。

这个量在M步里会被用来重新估计π和B。因为π是初始状态概率，只需要看t=1时的γ₁(i)；而π_i的估计就是γ₁(i)本身，这就是为什么E步算完，π的重估几乎顺带就完成了。

2.3 ξ_t(i,j)：相邻时刻的状态迁移概率

比γ更精细一点，我们还需要知道t时刻处于状态i、t+1时刻处于状态j的后验概率，记作ξ_t(i,j)。

当初我在学的时候，总是记不住ξ的公式，后来发现拆开看就特别清楚：

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

分子里四个因子各司其职：

• α_t(i)告诉你t时刻之前的信息。
• a_ij告诉你从i跳到j的先验概率。
• b_j(o_{t+1})告诉你跳到j之后能解释o_{t+1}的能力。
• β_{t+1}(j)告诉你j之后的信息。

这四个因子一乘，再除以观测似然，就把“从i跳到j”这件事的后验概率算出来了。

γ_t(i)和ξ_t(i,j)之间还有一个恒等关系：

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

这个关系在M步化简公式时会用上。

2.4 E步的计算量分析

如果你需要自己实现，要意识到E步并不是只跑一次前向和后向就完事。对于一个长度为T、状态数为N的观测序列，每一轮迭代都要：

• 跑一遍前向递推：时间复杂度O(T·N²)
• 跑一遍后向递推：同样是O(T·N²)
• 再算所有时刻的所有γ和ξ：也是O(T·N²)

空间复杂度方面，α、β、γ、ξ都需要存下来，每个都是T×N或T×N×N的量级。N=10、T=10000时，仅ξ矩阵就是10⁶量级的浮点数，内存压力不大；但如果N到几百，T到几十万，直接存储会变得很吃力。实际工程里可能需要对序列做分块处理，这个问题等实操部分再细说。

3. M步的三个重估公式

3.1 从软计数到硬公式

有了γ和ξ，M步就变成了一件很机械的事情：用期望的“软计数”替代真实的“硬计数”，重新估计π、A、B。

先看初始状态概率π：

π_i_new = γ₁(i)

这个公式的含义是：在所有可能的状态序列中，第一步落在状态i的期望次数占总序列数的比例。如果有K条独立观测序列，那就对所有序列的γ₁(i)取平均。

再看转移概率a_ij：

a_ij_new = Σ_t ξ_t(i,j) / Σ_t γ_t(i)

分子是所有时刻从i转移到j的期望次数，分母是从i转移出去的期望总次数。两者一比，就是转移概率的最大似然估计。这里t从1到T-1求和，因为最后一次转移发生在T-1到T之间。

最后看发射概率b_j(k)：

b_j(k)_new = Σ_t γ_t(j)·1{o_t=k} / Σ_t γ_t(j)

分子里那个1{o_t=k}是指示函数，只有当t时刻观测恰好是k时才计入。这个公式读起来就是：在所有处于状态j的时刻里，观测到符号k的期望比例。

3.2 重估公式的收敛保证

你可能会问：这几个公式看起来就是简单的频率计算，凭什么说它一定能让似然值增大？

关键在于EM算法的理论保证。Baum-Welch本质上是EM的特例，因此满足EM的单调性定理：每一轮迭代之后的似然P(O|λ_new)不会低于迭代前的P(O|λ_old)。

实际跑代码的时候，我建议你在每一轮迭代后都计算一次对数似然log P(O|λ)，并打印出来。正常情况下它会一直往上走，只是增幅会逐渐变小。如果在某一轮似然值突然下降了，那几乎可以肯定代码里存在bug，最常见的锅是ξ或γ的计算公式写错了。

3.3 多个观测序列的处理方式

很多教程都只讲单序列的情况，但我实测下来，真实项目里几乎都是多个观测序列一起喂给模型。

比如我要做用户行为序列建模，手里有一万个用户的会话记录，每条长度还不一样。这时候的处理方式是：

E步每条序列单独跑前向后向，得到各自的γ和ξ，但要把所有序列的γ和ξ在分子分母上分别累加起来。M步用累加后的总数做估计。

换句话说，π的重估是所有序列第一步状态的期望平均，A和B的重估是跨序列的总期望计数之比。

这种处理方式天然保证了不同长度序列的权重均衡——长序列和短序列在计算中都以“序列”为单位参与期望计数，不会因为某条序列特别长就主导整个参数估计。

4. 数值稳定性与工程实现细节

4.1 下溢问题的根源

在实际写代码之前，有一个必须先解决的问题：数值下溢。

前向算法的递推式是不断相乘的过程，每一步都乘上转移概率和发射概率，而这些值都小于1。当序列长度T达到几百甚至几千时，α_t(i)的值会小到超出双精度浮点数的表示范围，变成0。

后向算法也是一样。如果α、β算出来全是0，那γ和ξ就全是0/0，整个EM迭代直接报废。

解决这个问题的标准做法是引入缩放因子（scaling），也叫归一化前向后向算法。

4.2 缩放法的具体操作

前向算法的缩放思路很简单：每一轮递推后，把α向量除以当前时刻所有状态概率之和。

具体来说：

1. 计算α₁(i) = π_i · b_i(o₁)。
2. 求c₁ = Σ_i α₁(i)，然后让α₁(i) = α₁(i) / c₁。

这样α₁向量的和就是1，不会跑到下溢界。 3. 继续递推α_{t+1}(j)，再除以c_{t+1} = Σ_j α_{t+1}(j)归一化。

后向算法也做同样的缩放，每一轮除以对应的c_t。

这么做的妙处在于，缩放因子c_t的乘积恰好等于观测似然的对数：

log P(O|λ) = Σ_t log(c_t)

这是因为每一次缩放都相当于把前向概率归一化了一次，累积起来，缩放因子乘积就是这个序列的似然值。

这个结果不光解决了下溢问题，还顺带给EM迭代提供了一个天然的收敛监控指标。

4.3 Python实现Baum-Welch

下面给你一份可以直接跑的Python实现。我不会用一个现成的库包装给你，而是把每一步都摊开，因为参数估计的调试过程本身就是学习的一部分。

import numpy as np class HMM: def __init__(self, N, M): """ N: 隐状态数 M: 观测符号数（离散观测） """ self.N = N self.M = M # 随机初始化，但要保证每行和为1 self.A = np.random.rand(N, N) self.A /= self.A.sum(axis=1, keepdims=True) self.B = np.random.rand(N, M) self.B /= self.B.sum(axis=1, keepdims=True) self.pi = np.random.rand(N) self.pi /= self.pi.sum() def forward_scaled(self, obs): T = len(obs) N = self.N alpha = np.zeros((T, N)) c = np.zeros(T) # t=0 alpha[0] = self.pi * self.B[:, obs[0]] c[0] = alpha[0].sum() alpha[0] /= c[0] # t=1..T-1 for t in range(1, T): alpha[t] = (alpha[t-1] @ self.A) * self.B[:, obs[t]] c[t] = alpha[t].sum() alpha[t] /= c[t] log_likelihood = np.sum(np.log(c)) return alpha, c, log_likelihood def backward_scaled(self, obs, c): T = len(obs) N = self.N beta = np.zeros((T, N)) beta[T-1] = 1.0 for t in range(T-2, -1, -1): beta[t] = self.A @ (self.B[:, obs[t+1]] * beta[t+1]) beta[t] /= c[t+1] return beta def baum_welch(self, observations, max_iter=100, tol=1e-6): """ observations: list of list，每条是观测序列，观测值为0..M-1的整数 """ old_loglik = None for it in range(max_iter): # E步积累量 pi_num = np.zeros(self.N) A_num = np.zeros((self.N, self.N)) A_den = np.zeros(self.N) B_num = np.zeros((self.N, self.M)) B_den = np.zeros(self.N) total_loglik = 0.0 for obs in observations: T = len(obs) alpha, c, loglik = self.forward_scaled(obs) beta = self.backward_scaled(obs, c) total_loglik += loglik # gamma 和 xi gamma = alpha * beta gamma /= gamma.sum(axis=1, keepdims=True) xi = np.zeros((T-1, self.N, self.N)) for t in range(T-1): xi[t] = (alpha[t][:, None] * self.A * self.B[:, obs[t+1]][None, :] * beta[t+1][None, :]) xi[t] /= xi[t].sum() # 累积统计量 pi_num += gamma[0] A_num += xi.sum(axis=0) A_den += gamma[:-1].sum(axis=0) for k in range(self.M): mask = (np.array(obs) == k) B_num[:, k] += gamma[mask].sum(axis=0) B_den += gamma.sum(axis=0) # M步更新 self.pi = pi_num / pi_num.sum() self.A = A_num / A_den[:, None] self.B = B_num / B_den[:, None] # 收敛判断 if old_loglik is not None and abs(total_loglik - old_loglik) < tol: print(f"迭代{it+1}轮后收敛，log-likelihood = {total_loglik:.4f}") break old_loglik = total_loglik if (it + 1) % 10 == 0: print(f"第{it+1}轮，log-likelihood = {total_loglik:.4f}") return self.A, self.B, self.pi

代码里有几个细节值得单独拎出来说。

xi的计算我用了向量化。虽然内部循环不可避免，但矩阵操作可以让N比较小时跑得飞快。注意xi的分子分母最后要归一化，保证每一时刻的xi之和为1。

B的累积里用了一个布尔索引的技巧：先用观测序列生成一个mask，再把这个mask广播到gamma上做行筛选。这个写法比循环遍历每个时刻要清晰得多。

还有一点是A_den的求和范围：gamma[:-1]表示去掉最后一个时刻，因为最后时刻不需要发生转移。如果你把gamma的全部时刻都加进去了，A的每一行之和可能略大于1，这个误差虽小，但累积起来会影响收敛精度。

4.4 观测序列长度对迭代速度的实际影响

按上面的实现，跑一个N=5、M=10、每条序列长度T=200、共100条序列的模拟数据，在我的机器上每一轮大约要0.1秒左右，50轮收敛基本稳定。

真正要注意的瓶颈其实不是单轮速度，而是收敛速度。如果你发现似然值在30轮之后还在缓慢爬升，可以考虑提前设一个相对容差，或者用早期停止的验证集方案，避免在训练集上把参数拟合到过深。

5. Baum-Welch实战中的五大经典翻车现场

5.1 似然值下降

理论上EM保证单调不减，但如果你在迭代中看到log-likelihood突然下降，不要第一时间怀疑数学。

最常见的具体原因是缩放因子的使用姿势不对。后向算法里β每次除以的缩放因子，不是前向那一轮算出来的c[t]，而是c[t+1]。这是由递推方向决定的，因为β_t的计算过程里用到了t+1时刻的观测和缩放信息。我当年把c[t+1]写成c[t]之后，前面几轮看起来正常，到第6轮就开始全面崩坏。

5.2 状态编号翻转问题

另一个很容易遇到的现象是：多次用不同随机种子跑同一个数据，最后学到的参数看起来不太一样。比如第一次学到的状态1对应晴天、状态2对应雨天，第二次反过来状态1对应雨天、状态2对应晴天。

这不是bug，这是HMM天生的不可辨识性。状态只是集合，编号的顺序本身没有语义。你完全可以把所有状态重新编号，得到的模型和原来那个是等价的。

这意味着在评估模型时不能直接比较两个模型的状态编号，而要看它们对应的观测发射模式和转移模式。在工程上，如果你下游任务需要状态有稳定含义，可能需要额外加正则项或者用人工标注来锚定状态顺序。

5.3 初始值选择与局部最优

EM只能保证收敛到局部最优，不能保证找到全局最优。这在实际使用中很要命，因为HMM的似然面往往有很多局部极小点。

解决策略有两个比较有效。

第一是多次随机初始化。跑10次不同随机种子的Baum-Welch，保留似然值最高的那组参数。这几个并行实验成本可控，因为大多数数据量下每轮迭代都是秒级甚至亚秒级的。

第二是给初始参数一些先验知识。比如我做过手势识别，知道手势大致分成三到四个状态阶段，就可以用均匀转移矩阵作为A的起点，把B初始化为随机但每行方差适中的值。

还有一个我自己常用的技巧：先用KMeans之类的聚类方法把观测序列粗分类，然后再用聚类结果构造初始的B矩阵。这相当于给EM提供了一个更合理的起点，实测下来收敛速度明显更快。

5.4 零概率陷阱与平滑

发射矩阵B里如果某些状态下的某个观测符号从未出现过，重估公式会直接产生0概率。这看起来无害，但只要之后某条测试序列在这个状态下出现了该观测符号，整个似然值就直接变成0，一票否决。

解决方法是加平滑项，最常见的就是拉普拉斯平滑：

b_j(k)_new = (E[j,k] + δ) / (Σ_k E[j,k] + M·δ)

其中E[j,k]是重估公式里的分子期望计数，δ是一个小的正数，通常取0.1到1之间。

在数值上这相当于给每个观测-状态对加了一个虚拟的先验计数，防止概率被估成严格的0。这在处理文本类序列数据时尤其重要，因为自然语言的分布是长尾的，很多词的出现频率很低。

5.5 数据量不够，模型过拟合

HMM参数数量和状态数、观测数的平方成正比。N=10、M=1000的模型，光B矩阵就有10000个参数。如果你的观测序列总长度只有几百，那学出来的参数基本是噪声。

判断是否过拟合有一个简单办法：把数据划分成训练集和验证集，在训练集上跑Baum-Welch，在验证集上计算log-likelihood。训练集上的似然会一直上升，但验证集似然开始下降的时候，就是该收手的时候。

如果数据量实在不够，减少状态数N是最直接的缓解手段。也可以用先验结构来约束参数空间，比如设定某些转移概率必须为0，或者让B矩阵用参数化的概率分布形式来表示（例如高斯发射的连续HMM）。

6. 一个完整的模拟实验

理论说了一堆，不如跑个实验看结果。

我构造一个最简单的天气场景。假设隐状态有两个：0代表晴天，1代表雨天。转移概率是晴天转晴天0.7，晴天转雨天0.3，雨天转晴天0.4，雨天转雨天0.6。发射符号是带伞与否，0代表不带伞，1代表带伞。晴天的发射是0.9不带伞，雨天的是0.2不带伞。

用这个真实参数生成100条长度为100的观测序列，然后故意用一组和真实参数相差很大的随机参数做初始化，跑Baum-Welch。

np.random.seed(42) true_A = np.array([[0.7, 0.3], [0.4, 0.6]]) true_B = np.array([[0.9, 0.1], [0.2, 0.8]]) true_pi = np.array([0.5, 0.5]) def generate_obs(T, A, B, pi): obs = [] z = np.random.choice(2, p=pi) for _ in range(T): o = np.random.choice(2, p=B[z]) obs.append(o) z = np.random.choice(2, p=A[z]) return obs train_seqs = [generate_obs(100, true_A, true_B, true_pi) for _ in range(100)] model = HMM(N=2, M=2) model.A = np.array([[0.3, 0.7], [0.6, 0.4]]) model.B = np.array([[0.3, 0.7], [0.8, 0.2]]) model.pi = np.array([0.1, 0.9]) A_est, B_est, pi_est = model.baum_welch(train_seqs, max_iter=100, tol=1e-5)

因为存在状态编号翻转的可能，我做完实验后手动对比一下学到的参数。如果状态0对应晴天、状态1对应雨天，那估计结果应该接近真实值。

实际跑下来的结果让我印象很深：前向算法和后向算法在这里完全没遇到下溢问题，因为有缩放处理；EM在第20轮左右就达到了10⁻⁵量级的收敛容差；最终学到的A矩阵和真实值相差在5%以内，B矩阵误差更小。

这说明只要数据量足够、模型设定正确，Baum-Welch确实能够从纯观测序列里把参数挖出来。

7. 扩展思考与工程建议

7.1 连续观测怎么办

这篇里聊的是离散观测符号的情况。但现实中更多数据是连续值，比如语音特征、传感器数值、股票收益率。这时候发射概率不能用离散矩阵B，而要用高斯分布或高斯混合模型来建模。

参数估计框架不变，变的只是B的更新公式。原来计算的是“某个状态发射某个符号的期望计数”，现在变成了“某个状态发射某段连续值的高斯均值和协方差的加权估计”。权重的角色还是γ_t(j)。这类模型叫Gaussian-HMM或GMM-HMM，语音识别里曾经是绝对主流。

7.2 大规模数据的工程策略

当序列特别长、状态数特别大时，Baum-Welch的直接实现会面临内存和计算时间的双重压力。

工程上的常见做法是分块处理：把长序列切成多个子序列，在切分边缘处稍微重叠一点，避免丢掉跨块的状态转移信息。还可以用分段并行，把不同序列的E步分配到不同核心上执行，M步再汇聚。

如果连单序列都长到前向后向都跑不动，那就需要考虑结构化近似推断的方法了，这已经超出Baum-Welch的范畴。

7.3 模型选择：状态数N怎么定

绝大多数HMM使用者都会在某个时刻面对这个问题：到底设几个状态？

有两条路线。一条是走信息准则路线，分别训练N从2到15的模型，计算BIC或AIC，选个分数低的。另一条是走业务解读路线，看学到的状态能不能对应到有业务含义的模式上。比如用户行为序列，N=4或5时状态恰好对应沉默期、探索期、活跃期、消费期，那就比纯粹看数字有意义。

两种路线可以结合使用。先用BIC圈定一个候选范围，再在这个范围内用人眼检查状态语义。

我自己的经验是，N设得偏大比偏小更容易出问题。状态过多时，部分状态会退化成分裂状态，把原本一个连贯的行为模式拆成两个几乎一样的模式，导致下游任务里状态的含义不稳定。

最后说几句实操体会

Baum-Welch这个算法从1950年代末提出到现在，原理上基本没有变化，但每次自己从头实现一遍，还是会踩到新的坑。

我个人建议如果你真的想把HMM用熟练，不要只调现成库，试着手写一轮完整的前向、后向和Baum-Welch。不需要用很多数据，模拟数据就够。这个过程会让你对γ、ξ这些概念产生肌肉记忆，之后再去学条件随机场、深度概率模型这些更复杂的东西，会发现底子在HMM这里就打好了。

另外一个小建议：调试HMM代码时，务必把每一轮的log-likelihood打印出来。如果迭代曲线平滑上升，大概率代码是对的；如果出现了任何异常波动，先检查缩放因子传递、再检查下标错位，这两个问题占了我过去踩坑的八成。参数估计不是什么玄学，它是一个需要耐心观察迭代过程的工程问题。