如果你已经把前向算法和后向算法啃下来了那HMM参数估计这块也就是Baum-Welch算法其实就剩一层窗户纸了。标题写的是第2讲默认你已经知道隐马尔可夫模型是什么一条状态马尔可夫链每个状态再按一定概率发射观测值。今天要面对的问题是给定一堆观测序列怎么把一个完整的模型参数λ(A, B, π)给估出来。这问题严格来说属于EM算法在离散状态空间上的一个特例。很多人学到这里被满屏的希腊字母劝退其实核心思想在生活里很常见——你只能看到结果过程是隐藏的那就只能用结果去反推过程的统计规律。比如你看不见天气但能从一个人每天是否带伞推断出天气变化的规律这就是HMM参数估计在做的事。这篇内容会让你搞清楚两件事第一为什么参数估计不能像普通马尔可夫链那样直接计数第二Baum-Welch具体怎么迭代、代码怎么落地、实际运行时有哪些看不见的坑。适合正在学HMM、准备做序列建模或者想把EM思想落到真实代码里的朋友参考。1. 参数估计问题到底难在哪1.1 状态可见与不可见的本质差别如果HMM里的状态序列Z(z₁, z₂, ..., z_T)是直接可见的那参数估计就退化成一个简单的统计问题。比如我有N3个状态M2个观测符号跑了T1000步每一步我都知道自己在哪个状态、发射了什么观测。那么转移矩阵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(abc)没法拆成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这一段的联合概率。它可以从t1开始逐步递推α₁(i) π_i · b_i(o₁)α_{t1}(j) [Σ_i α_t(i) · a_ij] · b_j(o_{t1})后向概率β_t(i)定义为在t时刻处于状态i的前提下从t1时刻开始一直观测到o_T的概率。它从序列末尾往前递推β_T(i) 1β_t(i) Σ_j a_ij · b_j(o_{t1}) · β_{t1}(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_ti|O,λ) [α_t(i) · β_t(i)] / P(O|λ)直观理解就是t时刻“经前向传播过来的信息”和“经后向传播过来的信息”在状态i上碰头归一化后就是该时刻处于状态i的可能性。这个量在M步里会被用来重新估计π和B。因为π是初始状态概率只需要看t1时的γ₁(i)而π_i的估计就是γ₁(i)本身这就是为什么E步算完π的重估几乎顺带就完成了。2.3 ξ_t(i,j)相邻时刻的状态迁移概率比γ更精细一点我们还需要知道t时刻处于状态i、t1时刻处于状态j的后验概率记作ξ_t(i,j)。当初我在学的时候总是记不住ξ的公式后来发现拆开看就特别清楚ξ_t(i,j) [α_t(i) · a_ij · b_j(o_{t1}) · β_{t1}(j)] / P(O|λ)分子里四个因子各司其职α_t(i)告诉你t时刻之前的信息。a_ij告诉你从i跳到j的先验概率。b_j(o_{t1})告诉你跳到j之后能解释o_{t1}的能力。β_{t1}(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的量级。N10、T10000时仅ξ矩阵就是10⁶量级的浮点数内存压力不大但如果N到几百T到几十万直接存储会变得很吃力。实际工程里可能需要对序列做分块处理这个问题等实操部分再细说。3. M步的三个重估公式3.1 从软计数到硬公式有了γ和ξM步就变成了一件很机械的事情用期望的“软计数”替代真实的“硬计数”重新估计π、A、B。先看初始状态概率ππ_i_new γ₁(i)这个公式的含义是在所有可能的状态序列中第一步落在状态i的期望次数占总序列数的比例。如果有K条独立观测序列那就对所有序列的γ₁(i)取平均。再看转移概率a_ija_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_tk} / Σ_t γ_t(j)分子里那个1{o_tk}是指示函数只有当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 缩放法的具体操作前向算法的缩放思路很简单每一轮递推后把α向量除以当前时刻所有状态概率之和。具体来说计算α₁(i) π_i · b_i(o₁)。求c₁ Σ_i α₁(i)然后让α₁(i) α₁(i) / c₁。这样α₁向量的和就是1不会跑到下溢界。 3. 继续递推α_{t1}(j)再除以c_{t1} Σ_j α_{t1}(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(axis1, keepdimsTrue) self.B np.random.rand(N, M) self.B / self.B.sum(axis1, keepdimsTrue) 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) # t0 alpha[0] self.pi * self.B[:, obs[0]] c[0] alpha[0].sum() alpha[0] / c[0] # t1..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[t1]] * beta[t1]) beta[t] / c[t1] return beta def baum_welch(self, observations, max_iter100, tol1e-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(axis1, keepdimsTrue) 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[t1]][None, :] * beta[t1][None, :]) xi[t] / xi[t].sum() # 累积统计量 pi_num gamma[0] A_num xi.sum(axis0) A_den gamma[:-1].sum(axis0) for k in range(self.M): mask (np.array(obs) k) B_num[:, k] gamma[mask].sum(axis0) B_den gamma.sum(axis0) # 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迭代{it1}轮后收敛log-likelihood {total_loglik:.4f}) break old_loglik total_loglik if (it 1) % 10 0: print(f第{it1}轮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 观测序列长度对迭代速度的实际影响按上面的实现跑一个N5、M10、每条序列长度T200、共100条序列的模拟数据在我的机器上每一轮大约要0.1秒左右50轮收敛基本稳定。真正要注意的瓶颈其实不是单轮速度而是收敛速度。如果你发现似然值在30轮之后还在缓慢爬升可以考虑提前设一个相对容差或者用早期停止的验证集方案避免在训练集上把参数拟合到过深。5. Baum-Welch实战中的五大经典翻车现场5.1 似然值下降理论上EM保证单调不减但如果你在迭代中看到log-likelihood突然下降不要第一时间怀疑数学。最常见的具体原因是缩放因子的使用姿势不对。后向算法里β每次除以的缩放因子不是前向那一轮算出来的c[t]而是c[t1]。这是由递推方向决定的因为β_t的计算过程里用到了t1时刻的观测和缩放信息。我当年把c[t1]写成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参数数量和状态数、观测数的平方成正比。N10、M1000的模型光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, ppi) for _ in range(T): o np.random.choice(2, pB[z]) obs.append(o) z np.random.choice(2, pA[z]) return obs train_seqs [generate_obs(100, true_A, true_B, true_pi) for _ in range(100)] model HMM(N2, M2) 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_iter100, tol1e-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选个分数低的。另一条是走业务解读路线看学到的状态能不能对应到有业务含义的模式上。比如用户行为序列N4或5时状态恰好对应沉默期、探索期、活跃期、消费期那就比纯粹看数字有意义。两种路线可以结合使用。先用BIC圈定一个候选范围再在这个范围内用人眼检查状态语义。我自己的经验是N设得偏大比偏小更容易出问题。状态过多时部分状态会退化成分裂状态把原本一个连贯的行为模式拆成两个几乎一样的模式导致下游任务里状态的含义不稳定。最后说几句实操体会Baum-Welch这个算法从1950年代末提出到现在原理上基本没有变化但每次自己从头实现一遍还是会踩到新的坑。我个人建议如果你真的想把HMM用熟练不要只调现成库试着手写一轮完整的前向、后向和Baum-Welch。不需要用很多数据模拟数据就够。这个过程会让你对γ、ξ这些概念产生肌肉记忆之后再去学条件随机场、深度概率模型这些更复杂的东西会发现底子在HMM这里就打好了。另外一个小建议调试HMM代码时务必把每一轮的log-likelihood打印出来。如果迭代曲线平滑上升大概率代码是对的如果出现了任何异常波动先检查缩放因子传递、再检查下标错位这两个问题占了我过去踩坑的八成。参数估计不是什么玄学它是一个需要耐心观察迭代过程的工程问题。