# Fit Hidden Markov Model

> 使用 Baum-Welch（EM）算法拟合隐马尔可夫模型，包含模型选择、用于状态序列的 Viterbi 解码以及前向-后向概率。适用于观测由不可观测的隐状态生成时、需要将 时间序列分割为隐含状态（市场状态、语音音素、生物序列）时、计算序列概率、 解码最可能的隐状态路径，或比较不同隐状态数目的模型。

- Skill: `pjt222/fit-hidden-markov-model-4` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pjt222/fit-hidden-markov-model-4`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pjt222/fit-hidden-markov-model-4/raw
- Safety review: pending (external: skill-scanner PASS, skillspector PASS)
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: AI & ML
- License: MIT
- Author: pjt222 (https://skillmd.com/u/pjt222)
- Updated: 2026-09-21
- Page: https://skillmd.com/skills/pjt222/fit-hidden-markov-model-4

---


# 拟合隐马尔可夫模型

使用 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](../model-markov-chain/SKILL.md) — 理解构成隐含层基础的转移结构的前提
- [Simulate Stochastic Process](../simulate-stochastic-process/SKILL.md) — 可用于生成合成 HMM 数据进行测试以及从拟合模型模拟进行后验预测检验