拟合隐马尔可夫模型

使用 Baum-Welch 期望最大化算法将隐马尔可夫模型（HMM）拟合到序列观测数据，通过 Viterbi 解码最可能的隐状态序列，并通过信息准则选择最优隐状态数目。

适用场景

• 观测到一系列发射值，但底层生成状态不可直接观测
• 怀疑数据由在有限数量状态之间切换的系统生成
• 需要将时间序列分割为隐含阶段（如市场状态、语音音素、生物序列注释）
• 想要计算观测序列在生成模型下的概率
• 需要给定观测的最可能隐状态序列（解码）
• 比较不同隐状态数目的模型以获得最佳复杂度-拟合权衡

输入

必需

输入	类型	描述
observations	sequence/matrix	观测数据序列（单变量或多变量）
n_hidden_states	integer	要拟合的隐状态数（或用于模型选择的范围）
emission_type	string	发射分布族："gaussian"、"discrete"、"poisson"、"multinomial"

可选

输入	类型	默认值	描述
initial_params	dict	random/heuristic	初始转移矩阵、发射参数和初始概率
n_restarts	integer	10	随机重启次数以缓解局部最优
max_iterations	integer	500	每次重启的最大 EM 迭代次数
convergence_tol	float	1e-6	EM 的对数似然收敛阈值
state_range	list of ints	[n_hidden_states]	用于模型选择的状态数范围
covariance_type	string	"full"	高斯发射："full"、"diagonal"、"spherical"
regularization	float	1e-6	添加到协方差矩阵对角线的小常数以防止奇异

步骤

第 1 步：定义隐状态和观测模型

1.1. 指定隐状态数 K（或在第 5 步中用于模型选择的候选范围）。

1.2. 根据数据类型选择发射分布族：

• 连续数据：高斯（单变量或多变量）
• 计数数据：泊松或负二项
• 类别数据：离散/多项式

1.3. 定义模型组件：

• 转移矩阵 A，大小 K x K：A[i,j] = P(z_t = j | z_{t-1} = i)
• 发射参数 theta_k，每个状态 k：分布特定（如高斯的均值和协方差）
• 初始状态分布 pi：pi[k] = P(z_1 = k)

1.4. 验证观测数据格式正确：序列中无缺失值，维度一致，且长度相对于参数数量足够。

预期结果： 清晰指定的 HMM 架构，具有 K 个状态、已选择的发射族和长度 T >> K^2 的干净观测数据。

失败处理： 如果数据包含缺失值，插补或删除受影响的段。如果 T 相对于 K 太小，减少 K 或获取更多数据。

第 2 步：初始化参数

2.1. 为每次 n_restarts 重启生成初始参数：

• 转移矩阵：随机行随机矩阵（每行从 Dirichlet 分布抽取）或略微扰动的均匀矩阵。
• 发射参数：使用 K-means 聚类观测值来初始化均值；为高斯发射计算聚类方差。
• 初始分布：均匀或与 K-means 的聚类大小成比例。

2.2. 对于第一次重启，使用 K-means 知情的初始化（通常是最强的起始）。对于后续重启，使用随机扰动。

2.3. 验证所有初始参数有效：

• 转移矩阵行和为 1，所有条目为正。
• 发射参数在有效域中（如协方差矩阵正定）。
• 初始分布和为 1。

预期结果： n_restarts 组有效的初始参数，至少有一个数据驱动的初始化。

失败处理： 如果 K-means 无法收敛，使用更多重启的纯随机初始化。如果协方差矩阵奇异，将正则化常数添加到对角线。

第 3 步：运行 Baum-Welch EM 进行参数估计

3.1. E 步（前向-后向算法）：

• 使用递推计算前向概率 alpha[t,k] = P(o_1,...,o_t, z_t=k | model)：
  • alpha[1,k] = pi[k] * b_k(o_1)
  • alpha[t,k] = sum_j(alpha[t-1,j] * A[j,k]) * b_k(o_t)
• 计算后向概率 beta[t,k] = P(o_{t+1},...,o_T | z_t=k, model)：
  • beta[T,k] = 1
  • beta[t,k] = sum_j(A[k,j] * b_j(o_{t+1}) * beta[t+1,j])
• 计算状态后验 gamma[t,k] = P(z_t=k | O, model)：
  • gamma[t,k] = alpha[t,k] * beta[t,k] / P(O | model)
• 计算转移后验 xi[t,i,j] = P(z_t=i, z_{t+1}=j | O, model)。

3.2. M 步（参数重估计）：

• 更新转移矩阵：A[i,j] = sum_t(xi[t,i,j]) / sum_t(gamma[t,i])
• 使用加权充分统计量更新发射参数：
  • 高斯均值：mu_k = sum_t(gamma[t,k] * o_t) / sum_t(gamma[t,k])
  • 高斯协方差：加权散布矩阵加正则化
  • 离散：b_k(v) = sum_t(gamma[t,k] * I(o_t=v)) / sum_t(gamma[t,k])
• 更新初始分布：pi[k] = gamma[1,k]

3.3. 计算对数似然：log P(O | model) = log sum_k(alpha[T,k])。使用 log-sum-exp 技巧防止下溢。

3.4. 缩放： 使用缩放的前向-后向变量防止长序列的数值下溢。在每个时间步归一化 alpha 并累积对数缩放因子。

3.5. 重复 E 步和 M 步直到对数似然变化低于 convergence_tol 或达到 max_iterations。

3.6. 在所有重启中，保留最终对数似然最高的参数集。

预期结果： 各迭代间对数似然单调非递减，在 max_iterations 内收敛。最终参数有效（行随机矩阵、正定协方差）。

失败处理： 如果对数似然下降，E 步或 M 步中有 bug——验证公式。如果收敛非常慢，尝试更好的初始化或增加 max_iterations。如果协方差变得奇异，增加正则化。

第 4 步：应用 Viterbi 解码获取最可能状态序列

4.1. 初始化 Viterbi 变量：

• delta[1,k] = log(pi[k]) + log(b_k(o_1))
• psi[1,k] = 0（无前驱）

4.2. 对 t = 2,...,T 前向递推：

• delta[t,k] = max_j(delta[t-1,j] + log(A[j,k])) + log(b_k(o_t))
• psi[t,k] = argmax_j(delta[t-1,j] + log(A[j,k]))

4.3. 终止：

• z*_T = argmax_k(delta[T,k])
• 最佳路径对数概率：max_k(delta[T,k])

4.4. 对 t = T-1,...,1 回溯：

• z*_t = psi[t+1, z*_{t+1}]

4.5. 输出解码的状态序列 z* = (z*_1, ..., z*_T) 及其对数概率。

4.6. 将 Viterbi 路径概率与前向算法的总序列概率进行比较，评估最佳路径的主导程度。

预期结果： 长度为 T 的单个最可能状态序列，每个条目在 {1,...,K} 中。Viterbi 对数概率应小于等于总对数似然。

失败处理： 如果 Viterbi 路径的对数概率为负无穷，某些转移或发射概率在不应该为零的地方为零。添加下限值防止 log(0)。

第 5 步：执行模型选择（跨模型阶数的 BIC/AIC）

5.1. 对 state_range 中每个候选隐状态数 K，拟合完整的 HMM（第 2-4 步）。

5.2. 计算自由参数数 p：

• 转移矩阵：K * (K - 1)（每行是一个单纯形）
• 发射参数：取决于分布族（如 d 维全协方差高斯：K * (d + d*(d+1)/2)）
• 初始分布：K - 1

5.3. 计算信息准则：

• BIC = -2 * log_likelihood + p * log(T)
• AIC = -2 * log_likelihood + 2 * p
• AICc = AIC + 2*p*(p+1) / (T - p - 1)（小样本校正）

5.4. 选择 BIC 最低（一致性优先）或 AIC 最低（预测优先）的模型。同时报告两者。

5.5. 汇总结果：对每个 K，显示对数似然、参数数、BIC、AIC 和收敛状态。

5.6. 如果最优 K 在 state_range 的边界上，扩展范围并重新拟合。

预期结果： BIC/AIC 中存在明确的最小值，标识最优隐状态数。所选模型应已收敛且状态含义可解释。

失败处理： 如果不存在明确的最小值（BIC 单调递减），模型可能设定错误——考虑不同的发射族。如果所有模型的对数似然都较差，数据可能不遵循 HMM 结构。

第 6 步：使用留出数据和后验解码进行验证

6.1. 将数据分为训练集和验证集（如 80/20 或使用多个序列）。

6.2. 在训练数据上拟合模型。使用前向算法在留出数据上计算对数似然（不重新拟合参数）。

6.3. 后验解码（Viterbi 的替代方法）：

• 对每个时间步，分配后验概率最高的状态：z^_t = argmax_k(gamma[t,k])
• 这最大化了正确解码状态的期望数量（与 Viterbi 最大化联合路径概率相对）。

6.4. 比较 Viterbi 和后验解码：

• 计算两个解码序列的一致率。
• 不一致的区域表示模糊的状态分配。

6.5. 评估状态可解释性：

• 检查每个状态的发射参数（均值、方差、离散分布）。
• 验证状态对应于领域上下文中有意义的状态。
• 检查状态停留时间（由 A 的对角线隐含）是否合理。

6.6. 计算每个观测的留出对数似然，并跨模型阶数比较以确认训练集模型选择。

预期结果： 留出对数似然与训练对数似然合理接近（无严重过拟合）。Viterbi 和后验解码在 90%+ 的时间步上一致。状态具有不同的、可解释的发射分布。

失败处理： 如果留出似然比训练差很多，模型过拟合——减少 K 或增加正则化。如果状态不可解释，尝试不同的初始化或不同的发射族。

验证清单

• 每次重启的 Baum-Welch 迭代间对数似然单调非递减
• 转移矩阵行随机（行和为 1，所有条目非负）
• 发射参数在有效域中（正定协方差、有效概率分布）
• Viterbi 路径对数概率不超过总序列对数概率
• BIC/AIC 曲线在所选模型阶数处显示明确的最小值
• 留出对数似然确认模型在训练集之外有泛化能力
• 前向和后向概率计算一致：P(O) = sum_k(alpha[T,k]) = sum_k(pi[k] * b_k(o_1) * beta[1,k])

常见问题

• EM 中的局部最优：Baum-Welch 算法收敛到局部最大值，不一定是全局最大值。始终使用多次随机重启并选择最佳
• 数值下溢：前向-后向概率随序列长度指数级缩小。使用对数空间计算或缩放变量防止下溢到零
• 状态过多导致过拟合：每增加一个隐状态会增加 O(K + d^2) 个参数。使用 BIC（不仅是似然）进行模型选择并在留出数据上验证
• 标签切换：隐状态仅在排列意义下可识别。跨重启比较模型时，按发射参数匹配状态，而非按索引
• 退化状态：一个状态可能坍缩为解释单个观测（方差接近零的高斯）。对协方差矩阵的正则化可防止这种情况
• 混淆 Viterbi 和后验解码：Viterbi 给出单个最佳联合路径；后验解码给出每个时间步的最佳边际状态。它们回答不同的问题，可能显著不一致
• 忽视状态停留时间：标准 HMM 中隐含的几何停留时间分布可能不适合具有长状态持续时间的数据。如果停留时间是非几何的，考虑隐半马尔可夫模型

相关技能

• Model Markov Chain — 理解构成隐含层基础的转移结构的前提
• Simulate Stochastic Process — 可用于生成合成 HMM 数据进行测试以及从拟合模型模拟进行后验预测检验