别再死记硬背公式了！用Python手撸一个HMM，带你直观理解前向后向算法

用Python手写HMM前向后向算法从数学公式到代码的思维跃迁很多机器学习初学者在第一次接触隐马尔可夫模型(HMM)时都会被那些复杂的概率公式搞得晕头转向。前向概率、后向概率、状态转移矩阵...这些概念在数学推导中看起来清晰明了但一到实际应用就让人摸不着头脑。本文将通过Python代码实现带你直观理解HMM的核心算法。1. HMM基础概念与问题定义隐马尔可夫模型是一种处理序列数据的强大工具广泛应用于语音识别、自然语言处理等领域。想象一个简单的天气预测场景我们无法直接观测到天气状态晴天或雨天但能观察到人们是否带伞。HMM就是要解决这种通过可见现象推测隐藏状态的问题。HMM由以下要素构成状态集合Q例如{晴天雨天}观测集合V例如{带伞不带伞}状态转移矩阵A描述天气状态间的转移概率观测概率矩阵B描述在特定天气下观察到某种行为的概率初始状态分布π系统初始时处于各状态的概率HMM需要解决的三大核心问题评估问题给定模型参数和观测序列计算该序列出现的概率前向后向算法解码问题给定观测序列找出最可能的状态序列Viterbi算法学习问题从观测序列中估计模型参数Baum-Welch算法# HMM模型基本参数定义示例 import numpy as np # 状态集合晴天、雨天 states [Sunny, Rainy] # 观测集合带伞、不带伞 observations [Umbrella, No_umbrella] # 状态转移矩阵AA[i][j]表示从状态i转移到状态j的概率 A np.array([[0.7, 0.3], # Sunny - Sunny, Sunny - Rainy [0.4, 0.6]]) # Rainy - Sunny, Rainy - Rainy # 观测概率矩阵BB[i][j]表示在状态i下观测到j的概率 B np.array([[0.1, 0.9], # Sunny下带伞、不带伞的概率 [0.8, 0.2]]) # Rainy下带伞、不带伞的概率 # 初始状态分布π pi np.array([0.6, 0.4]) # 初始为晴天、雨天的概率2. 前向算法从前往后的概率计算前向算法解决的是评估问题给定模型参数λ(A,B,π)和观测序列O计算P(O|λ)。直观理解就是计算在当前模型下出现这个观测序列的可能性有多大。前向算法通过动态规划的思想逐步计算前向概率αₜ(i)P(o₁,o₂,...,oₜ,qₜi|λ)即在时刻t处于状态i并观察到前t个观测的概率。算法步骤初始化计算第一个时刻的前向概率 α₁(i) πᵢ * bᵢ(o₁)递推对于t2到T计算 αₜ(j) [∑ᵢ αₜ₋₁(i) * aᵢⱼ] * bⱼ(oₜ)终止P(O|λ) ∑ᵢ α_T(i)def forward_algorithm(obs_seq, A, B, pi): 前向算法实现 :param obs_seq: 观测序列如[Umbrella, No_umbrella] :param A: 状态转移矩阵 :param B: 观测概率矩阵 :param pi: 初始状态分布 :return: 观测序列的概率P(O|λ) T len(obs_seq) # 序列长度 N A.shape[0] # 状态数 alpha np.zeros((T, N)) # 1. 初始化 obs_idx observations.index(obs_seq[0]) alpha[0, :] pi * B[:, obs_idx] # 2. 递推 for t in range(1, T): obs_idx observations.index(obs_seq[t]) for j in range(N): alpha[t, j] np.sum(alpha[t-1, :] * A[:, j]) * B[j, obs_idx] # 3. 终止 return np.sum(alpha[-1, :]) # 示例计算观测序列[Umbrella, No_umbrella]的概率 obs_seq [Umbrella, No_umbrella] prob forward_algorithm(obs_seq, A, B, pi) print(f观测序列的概率: {prob:.4f})提示前向算法的复杂度是O(TN²)远优于直接计算的O(TNᵀ)。这种效率提升来自于动态规划对中间结果的复用。3. 后向算法从后往前的概率计算与前向算法相反后向算法定义后向概率βₜ(i)P(oₜ₊₁,oₜ₊₂,...,o_T|qₜi,λ)即在时刻t处于状态i的条件下观察到后面所有观测的概率。算法步骤初始化β_T(i) 1递推对于tT-1到1计算 βₜ(i) ∑ⱼ aᵢⱼ * bⱼ(oₜ₊₁) * βₜ₊₁(j)终止P(O|λ) ∑ᵢ πᵢ * bᵢ(o₁) * β₁(i)def backward_algorithm(obs_seq, A, B, pi): 后向算法实现 :param obs_seq: 观测序列 :param A: 状态转移矩阵 :param B: 观测概率矩阵 :param pi: 初始状态分布 :return: 观测序列的概率P(O|λ) T len(obs_seq) N A.shape[0] beta np.zeros((T, N)) # 1. 初始化 beta[-1, :] 1 # 最后一个时刻的β设为1 # 2. 递推 for t in range(T-2, -1, -1): obs_idx observations.index(obs_seq[t1]) for i in range(N): beta[t, i] np.sum(A[i, :] * B[:, obs_idx] * beta[t1, :]) # 3. 终止 obs_idx observations.index(obs_seq[0]) return np.sum(pi * B[:, obs_idx] * beta[0, :]) # 验证前后向算法结果一致 prob_backward backward_algorithm(obs_seq, A, B, pi) print(f后向算法计算的概率: {prob_backward:.4f})前后向算法的关系算法计算方向定义应用场景前向从前向后αₜ(i)P(o₁..oₜ,qₜi⎮λ)实时滤波、在线学习后向从后向前βₜ(i)P(oₜ₊₁..o_T⎮qₜi,λ)序列补全、离线分析4. 前后向算法的联合应用将前后向算法结合可以解决更多实际问题平滑问题计算P(qₜi|O,λ)即在知道完整观测序列后某一时刻处于某状态的概率 γₜ(i) αₜ(i)βₜ(i) / P(O|λ)状态预测计算P(qₜ₊₁j|O₁..ₜ,λ)即预测下一时刻的状态def compute_gamma(obs_seq, A, B, pi): 计算平滑概率γ T len(obs_seq) N A.shape[0] # 计算前向和后向概率 alpha np.zeros((T, N)) beta np.zeros((T, N)) # 前向计算 obs_idx observations.index(obs_seq[0]) alpha[0, :] pi * B[:, obs_idx] for t in range(1, T): obs_idx observations.index(obs_seq[t]) for j in range(N): alpha[t, j] np.sum(alpha[t-1, :] * A[:, j]) * B[j, obs_idx] # 后向计算 beta[-1, :] 1 for t in range(T-2, -1, -1): obs_idx observations.index(obs_seq[t1]) for i in range(N): beta[t, i] np.sum(A[i, :] * B[:, obs_idx] * beta[t1, :]) # 计算γ P_O np.sum(alpha[-1, :]) gamma (alpha * beta) / P_O return gamma # 计算平滑概率 gamma compute_gamma(obs_seq, A, B, pi) print(各时刻处于各状态的概率) for t in range(len(obs_seq)): print(f时刻{t1}: 晴天{gamma[t,0]:.4f}, 雨天{gamma[t,1]:.4f})5. 算法优化与工程实践在实际应用中我们还需要考虑数值稳定性问题。由于概率乘积可能导致数值下溢通常采用对数空间计算def forward_algorithm_log(obs_seq, A, B, pi): 对数空间的前向算法 T len(obs_seq) N A.shape[0] log_alpha np.zeros((T, N)) # 初始化 obs_idx observations.index(obs_seq[0]) log_alpha[0, :] np.log(pi) np.log(B[:, obs_idx]) # 递推 for t in range(1, T): obs_idx observations.index(obs_seq[t]) for j in range(N): log_alpha[t, j] np.logaddexp.reduce(log_alpha[t-1, :] np.log(A[:, j])) np.log(B[j, obs_idx]) # 终止 return np.logaddexp.reduce(log_alpha[-1, :]) log_prob forward_algorithm_log(obs_seq, A, B, pi) print(f对数概率: {log_prob:.4f})性能优化技巧使用矩阵运算替代循环提前计算并缓存常用值并行化计算使用JIT编译如Numba# 矩阵运算优化的前向算法 def forward_algorithm_matrix(obs_seq, A, B, pi): T len(obs_seq) N A.shape[0] alpha np.zeros((T, N)) # 初始化 obs_idx observations.index(obs_seq[0]) alpha[0, :] pi * B[:, obs_idx] # 递推 for t in range(1, T): obs_idx observations.index(obs_seq[t]) alpha[t, :] (alpha[t-1, :] A) * B[:, obs_idx] return np.sum(alpha[-1, :]) # 比较三种实现的性能 %timeit forward_algorithm(obs_seq, A, B, pi) %timeit forward_algorithm_log(obs_seq, A, B, pi) %timeit forward_algorithm_matrix(obs_seq, A, B, pi)通过代码实现HMM算法那些抽象的数学公式变得直观可感。在天气预测的例子中我们看到如何用概率矩阵描述状态转移和观测生成以及如何通过动态规划高效计算序列概率。这种从理论到实践的转换能力正是机器学习工程师的核心竞争力之一。