15、非参数隐马尔可夫模型：原理、算法与应用

非参数隐马尔可夫模型：原理、算法与应用

1. 引言

在许多统计建模和机器学习问题中，序列数据是核心。例如文本由单词序列组成，金融数据常是价格序列，语音信号可表示为短期功率谱系数序列等。虽然可以直接对时间序列元素间的关系进行建模，但很多情况下数据存在潜在的隐藏结构。

隐马尔可夫模型（HMM）是一种用于处理序列数据的重要模型，它能捕捉这种隐藏结构。HMM 假设存在一个马尔可夫链 (s_1, s_2, \cdots, s_T)，每个 (s_t) 处于 (K) 个可能状态之一，状态在时间 (t) 的分布仅通过时间 (t - 1) 的状态依赖于之前的状态，由一个 (K \times K) 的随机转移矩阵 (\pi) 描述，其中 (\pi_{ij} = P(s_t = j|s_{t - 1} = i))。通常，我们并不直接观察马尔可夫链，而是观察到依赖于状态相关参数 (\theta_{s_t}) 参数化的观测模型 (F) 的观测值 (y_t)。HMM 诱导的概率分布可以表示为：
[p(y_{1:T}, s_{1:T}|K, \pi, \theta) = \prod_{t = 1}^{T} p(s_t|s_{t - 1})p(y_t|s_t) = \prod_{t = 1}^{T} \pi_{s_{t - 1}, s_t}F(y_t; \theta_{s_t})]

设计 HMM 时，一个关键问题是如何选择状态数 (K)。如果我们对生成观测序列的底层物理过程有先验知识，并且知道该过程中的状态数，那么可以将 (K) 设置为该值。然而，在许多应用中，底层状态数是未知的，需要从数据中推断。

传统的贝叶斯方法在处理 HMM 时存在统计和计算上的局限性。统计上，假设有限数量的状态能充分建模序列往往是不合理的；计算上，比较不同有限状态数的 HMM 需要近似难以处理的边际似然。

为了克服这些问题，引入了非参数贝叶斯方法，特别是无限隐马尔可夫模型（iHMM）。iHMM 定义了一个具有可数无限个隐藏状态的马尔可夫链，对于任何有限观测序列，只会访问有限个状态。随着序列长度的增加和新情况的出现，可以从无限状态池中招募新状态。

2. 从 HMM 到贝叶斯 HMM

在实际应用中，我们通常已知序列 (y_{1:T})，并希望了解隐藏表示 (s_{1:T}) 以及参数 (\pi)、(\theta) 和 (K)。根据可用的领域知识，我们可以区分以下三个计算问题：
- (\pi)、(\theta)、(K) 已知 ：在这种情况下，我们只需要根据观测值 (y_{1:T}) 推断 (s_{1:T})。可以应用贝叶斯规则找到 (s_{1:T}) 的后验分布：
[p(s_{1:T}|K, \pi, \theta, y_{1:T}) = \frac{p(y_{1:T}, s_{1:T}|K, \pi, \theta)}{p(y_{1:T}|K, \pi, \theta)} \propto \prod_{t = 1}^{T} p(s_t|s_{t - 1})p(y_t|s_t)]
计算这个分布可以使用动态规划的前向 - 后向算法。
- (K) 已知，(\pi)、(\theta) 待学习 ：当只有隐藏状态数 (K) 和观测值 (y_{1:T}) 已知时，我们通常希望学习最佳参数 (\theta) 和 (\pi) 以及隐藏表示 (s_{1:T})。常见的标准是最大似然和最大后验目标。期望最大化（EM）等算法可以搜索最大似然和最大后验解，但通常只能找到局部最优估计。
- (\pi)、(\theta)、(K) 待学习 ：在这种情况下，除了隐藏表示 (s_{1:T}) 和其他参数 (\pi)、(\theta) 外，还需要找到 (K) 的统计上有意义的值。使用最大似然准则不是一个好主意，因为更多的状态总是会导致更好的数据拟合，但可能会出现过拟合问题。赤池信息准则（AIC）和贝叶斯信息准则（BIC）可以用于调整最大似然估计，通过惩罚参数数量来避免过拟合。

更原则性的方法是对模型进行全贝叶斯分析。贝叶斯分析将参数 (\pi)、(\theta) 视为未知量，并将它们作为随机变量引入模型。这需要添加先验分布，例如 (p(\theta|H)) 和 (p(\pi|\alpha))，并扩展全联合分布：
[p(y_{1:T}, s_{1:T}, \pi, \theta|K) = p(\pi|\alpha)p(\theta|H) \left(\prod_{t = 1}^{T} p(s_t|s_{t - 1})p(y_t|s_t) \right) = p(\pi|\alpha)p(\theta|H) \left(\prod_{t = 1}^{T} \pi_{s_{t - 1}, s_t}F(y_t; \theta_{s_t}) \right)]

在贝叶斯分析中，我们可以计算边际似然 (p(y_{1:T}|K))，用于比较、选择或平均不同的 (K) 值。然而，解析计算 HMM 的边际似然是难以处理的。以下是三种处理这种难处理性的方法：
- 马尔可夫链蒙特卡罗（MCMC）技术 ：可以分为显式估计边际似然的方法和在不同 (K) 之间切换的方法。例如退火重要性采样和桥采样，但这些方法计算准确的边际似然估计可能很昂贵。可逆跳 MCMC 方法可以在不同大小的模型之间“跳跃”。
- 状态合并算法 ：如果隐藏状态 (s_{1:T}) 被观察到，参数 (\pi) 和 (\theta) 变得独立，并且在假设先验和似然是共轭的情况下，可以解析计算边际似然。通过选择一个好的状态序列，可以积分出其他参数并计算边际似然的近似值。
- 变分贝叶斯（VB）推理 ：计算边际似然的下界。VB 方法不使用参数 (\pi)、(\theta) 的点估计，而是使用这些参数的近似后验分布，并且不使用状态序列 (s_{1:T}) 的点估计，而是使用该量的完整分布。

3. 无限隐马尔可夫模型

传统的贝叶斯处理 HMM 需要为每个 (K) 值定义大量的 HMM。而无限隐马尔可夫模型（iHMM）引入了一个概率模型，其中 (K) 是无限的，任何有限序列 (s_{1:T}) 只访问有限个状态。iHMM 隐式地定义了 HMM 中访问状态数的分布，而不是定义 (K) 个不同的 HMM。

3.1 分层波利亚瓮方案

波利亚瓮方案是一类通过向瓮中填充彩色球的隐喻来定义离散概率分布的算法。分层波利亚瓮模型将波利亚瓮方案扩展到时间序列设置。

考虑一个可数无限个参数为 (\alpha) 的波利亚瓮，每个对应 iHMM 可能处于的一个状态，称为转移瓮。除了彩色球，瓮本身也被着色，球和瓮的颜色都对应于状态。

在每个时间步 (t)，记录前一个球的颜色 (s_{t - 1})，然后根据波利亚瓮方案从颜色为 (s_{t - 1}) 的瓮中抽取球。将 (s_t) 设置为添加到瓮 (s_{t - 1}) 中的额外球的颜色。转移概率可以表示为：
[p(s_t = j|s_{t - 1} = i, \alpha) = \frac{n_{ij}}{\sum_{j’} n_{ij’} + \alpha}]

如果需要一个新颜色的球，我们引入一个参数为 (\gamma) 的额外神谕波利亚瓮。从神谕瓮中抽取球，记录其颜色 (s_t)，并将两个相同颜色的球放回神谕瓮，一个相同颜色的球放回转移瓮 (s_{t - 1})。

超参数 (\gamma) 和 (\alpha) 在确定 iHMM 中的状态数方面起着重要作用。(\gamma) 控制添加全新颜色球的频率，即访问新状态的概率；(\alpha) 控制查询神谕的频率。

为了将分层波利亚瓮方案完成到一个完整的隐马尔可夫模型，我们需要引入一个机制来根据状态序列 (s_{1:T}) 生成观测值 (y_{1:T})。为此，我们引入：
- 一个基础分布 (H)，从中为每个状态 (k \in {1, \cdots, \infty}) 抽取输出似然参数 (\theta_k)。
- 一个似然模型 (F)，它接受状态 (s_t) 及其参数 (\theta_{s_t}) 并生成数据点 (y_t)。

3.2 分层狄利克雷过程

分层狄利克雷过程为 iHMM 提供了另一种解释。从狄利克雷过程 (G \sim DP(\alpha, H)) 中抽取的是一个离散分布，它将概率质量分配到 (H) 域中的可数个原子上。狄利克雷过程可以表示为：
[G = \sum_{i = 1}^{\infty} \beta_i \delta_{\theta_i}]
其中 (\beta_i = \hat{\beta} i \prod {l = 1}^{i - 1} (1 - \hat{\beta} l))，(\hat{\beta}_l \sim Beta(1, \alpha))，(\theta_i \sim H)，(\delta {\theta_i}) 是在 (\theta_i) 处的点质量。

分层狄利克雷过程的思想是，一个狄利克雷过程的基础分布本身是从另一个狄利克雷过程中抽取的：
[G_0|\gamma, H \sim DP(\gamma, H)]
[G_j|\alpha, G_0 \sim DP(\alpha, G_0)]

所有 (G_j) 共享相同的原子，并且 (G_j) 的概率质量根据棍子重新断裂过程分布。我们可以使用分层狄利克雷过程的构造来构建 iHMM，将无限维转移矩阵的行设置为可数无限个 (G_j) 的棍子。

虽然不明显，但分层波利亚瓮模型和分层狄利克雷过程定义了马尔可夫链 (s_{1:T}) 上的相同分布。

以下是分层波利亚瓮模型和分层狄利克雷过程的比较表格：
| 模型 | 特点 | 状态生成方式 | 超参数作用 |
| — | — | — | — |
| 分层波利亚瓮模型 | 通过向瓮中添加彩色球来定义状态转移 | 根据波利亚瓮方案从转移瓮和神谕瓮中抽取球 | (\gamma) 控制添加新颜色球的频率，(\alpha) 控制查询神谕的频率 |
| 分层狄利克雷过程 | 基于狄利克雷过程的分层结构 | 将无限维转移矩阵的行设置为 (G_j) 的棍子 | (\gamma) 影响基础分布的离散程度，(\alpha) 影响子狄利克雷过程的分布 |

4. 推理

iHMM 引入时提出了两种不同的推理算法：一种是每迭代时间复杂度为 (O(KT^2)) 的吉布斯采样器，另一种是时间复杂度为 (O(KT)) 的近似吉布斯采样器。分层狄利克雷过程框架对近似吉布斯采样器进行了修改，使其成为精确算法，同时保留了 (O(KT)) 的复杂度。

4.1 塌缩吉布斯采样器

塌缩吉布斯采样器的一次迭代包括两个步骤：
- 采样 (s_{1:T}) ：如果观察到 (\beta) 和 (s_{1:T})，转移矩阵的行变得条件独立。我们可以将每个转移视为从波利亚瓮方案中抽取，并使用波利亚瓮的知识重新采样状态序列。对于每个 (s_t)，其重新采样的概率为：
[p(s_t = k|s_{-t}, \alpha, \beta, \gamma) \propto \begin{cases}
\frac{(n_{s_{t - 1}, k}^{-t} + \alpha \beta_k)(n_{k, s_{t + 1}}^{-t} + \alpha \beta_{s_{t + 1}})}{n_{k \cdot}^{-t} + \alpha} & \text{if } k \leq K, k \neq s_{t - 1} \
\frac{(n_{s_{t - 1}, k}^{-t} + \alpha \beta_k)(n_{k, s_{t + 1}}^{-t} + 1 + \alpha \beta_{s_{t + 1}})}{n_{k \cdot}^{-t} + 1 + \alpha} & \text{if } k = s_{t - 1} = s_{t + 1} \
\frac{(n_{s_{t - 1}, k}^{-t} + \alpha \beta_k)(n_{k, s_{t + 1}}^{-t} + \alpha \beta_{s_{t + 1}})}{n_{k \cdot}^{-t} + 1 + \alpha} & \text{if } k = s_{t - 1} \neq s_{t + 1} \
\alpha \beta_k \beta_{s_{t + 1}} & \text{if } k = K + 1
\end{cases}]
- 采样 (\beta) ：(\beta) 是 HDP 基础分布原子上的分布。引入一个辅助变量 (m_{ij})，表示从状态 (i) 查询神谕时返回标签为 (j) 的球的神谕调用次数。通过采样 (m_{ij})，可以隐式采样 (c_j = \sum_{i} m_{ij})，从而重新采样 (\beta)。

塌缩吉布斯采样器的计算复杂度为 (O(TK + K^2))。然而，由于状态序列变量在马尔可夫链中通常是强相关的，吉布斯采样在存在强相关性时往往混合缓慢。

4.2 束采样器

束采样器是一种辅助变量 MCMC 算法，用于解决 iHMM 中状态序列强相关的问题。其核心思想是引入辅助变量 (u_{1:T})，使得在给定 (u_{1:T}) 的条件下，具有正概率的轨迹数量是有限的。

• 采样 (\pi) 和 (\theta) ：在采样 (u_{1:T}) 之前，采样转移和似然模型参数。给定状态序列 (s_{1:T})，((\pi_{k1} \cdots \pi_{kK}, \sum_{k’ = K + 1}^{\infty} \pi_{kk’})) 的条件分布为狄利克雷分布。每个 (\theta_k) 在给定 (s_{1:T})、(y_{1:T}) 和其先验分布 (H) 的条件下是独立的。
• 采样 (u_{1:T}) ：对于每个时间步 (t)，引入一个辅助变量 (u_t)，其条件分布为 (u_t \sim Uniform(0, \pi_{s_{t - 1}s_t}))。
• 采样 (s_{1:T}) ：状态序列 (s_{1:T}) 在给定辅助变量 (u_{1:T}) 的条件下，只有当 (\pi_{s_{t - 1}s_t} \geq u_t) 时才有非零概率。通过动态规划可以计算所有这样轨迹的条件分布。
• 采样 (\beta) ：丢弃转移矩阵 (\pi)，使用与塌缩吉布斯采样器相同的算法采样 (\beta)。

束采样器的最坏情况复杂度为 (O(TK^2))，但在实际应用中，由于只需要考虑转移概率大于 (u_t) 的状态，通常可以显著减少计算量。

以下是束采样器的流程 mermaid 图：

graph TD;
    A[初始化] --> B[采样 \(\pi\) 和 \(\theta\)];
    B --> C[采样 \(u_{1:T}\)];
    C --> D[采样 \(s_{1:T}\)];
    D --> E[采样 \(\beta\)];
    E --> F{是否达到迭代次数};
    F -- 是 --> G[结束];
    F -- 否 --> B;

5. 示例：无监督词性标注

词性标注是为句子中的单词标注其适当词性标签的任务。HMM 常用于词性标注，将单词作为观测值，隐藏状态序列对应于未知的词性标签。

在这个示例中，我们使用 iHMM 自动从数据中学习词性标签的数量。设置 iHMM 的参数为 (\gamma = 5) 和 (\alpha = 0.8)，使用简单的多项分布作为似然，对称狄利克雷分布作为基础分布 (H)。

在华尔街日报（WSJ）部分的宾州树库的第 0 节上进行训练，该数据集包含 1917 个句子，共 50282 个单词标记和 7904 个单词类型。初始化采样器为 50 个状态，并将每个 (s_t) 均匀随机分配到 50 个状态之一，运行 50000 次迭代。

通过分析前 10000 次迭代中状态数的演变，可以看到 iHMM 采样器在不同大小的模型之间无缝切换。

以下是前五个最常见状态的前五个单词的表格：
| 状态 ID | 频率 | 前五个单词及频率 |
| — | — | — |
| 9 | 8699 | of (1152), to (1000), in (829), and (439), for (432) |
| 12 | 4592 | the (1943), a (944), its (173), an (157), this (77) |
| 8 | 4216 |, (2498), and (260), ’s (146), said (136), – (95) |
| 18 | 3502 | company (93), year (85), market (52), U.S. (46), time (32) |
| 17 | 3470 | it (204), that (117), he (112), ” (110), they (101) |

从表格中可以看出，状态 9 分配了很多权重给介词，与真实标注方案中的介词类（IN）对应良好。这表明 iHMM 能够自动学习到有意义的词性标签。

6. 超越 iHMM 的扩展模型

iHMM 在实际应用中有很多局限性，为了更好地适应不同的场景，研究者们提出了一系列扩展模型。

6.1 输入 - 输出 iHMM（IO - iHMM）

传统的 iHMM 假设马尔可夫链自主演化，观测值在给定马尔可夫链的条件下是独立的。但在很多实际场景中，马尔可夫链会受到外部因素的影响。例如，机器人在室内环境中移动并拍照，若机器人有特定的行动策略，就可以将其行动作为输入引入 iHMM，这就是输入 - 输出 iHMM（IO - iHMM）。

IO - iHMM 可以看作是一个条件序列模型，通过隐藏表示 (s_{1:T}) 从输入 (e_{1:T}) 预测观测序列 (y_{1:T})，对结构化序列预测问题很有用。其转移概率依赖于输入 (e_t)，可看作是一个三维转移矩阵，即 (p(s_t|s_{t - 1}, e_t) = \pi_{e_t, s_{t - 1}, s_t})。向量 (\pi_{e, s, \cdot}) 的构建方式与 iHMM 中转移矩阵的行类似，通过 (\pi_{e, s, \cdot} \sim Stick(\alpha\beta)) 生成。

引入输入变量会增加参数数量。若用 iHMM 训练长度为 (T) 的序列，有效状态数为 (K)，需学习 (O(K^2)) 个转移矩阵参数；而用 IO - iHMM 训练相同序列，输入序列有 (E) 种可能输入符号，有效状态数为 (K)，则需学习 (O(EK^2)) 个转移矩阵参数。如果观测模型也依赖输入序列 (y_t \sim F(\theta_{s_t, e_t}))，观测模型的参数数量也会增加 (E) 倍。这意味着要么需要更多数据来学习与 iHMM 容量相同的模型，要么在相同数据量下学习容量更小的模型，或者通过耦合参数来降低模型复杂度。

IO - iHMM 的信念网络结构如下 mermaid 图所示：

graph LR;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    E1(输入 e1):::process --> S1(隐藏状态 s1):::process;
    E2(输入 e2):::process --> S2(隐藏状态 s2):::process;
    ET(输入 eT):::process --> ST(隐藏状态 sT):::process;
    S1 --> Y1(观测值 y1):::process;
    S2 --> Y2(观测值 y2):::process;
    ST --> YT(观测值 yT):::process;
    S1 --> S2;
    S2 --> ...(...):::process;
    ... --> ST;

6.2 粘性和块对角 iHMM

在很多应用中，明确建模状态之间的转移时间尺度很有用。对于 HMM 或 iHMM，转移矩阵的对角线权重控制着状态转移的频率，在状态 (i) 停留 (g) 个时间步的概率服从几何分布 (p(g) = \pi_{ii}^{g - 1}(1 - \pi_{ii}))。

• 粘性 iHMM ：最初的分层波利亚瓮模型引入了自转移超参数，增加了转移矩阵对角线的先验概率质量。在分层狄利克雷过程表示中，粘性 iHMM 从 (β \sim Stick) 中抽取基础度量，然后从以下分布中抽取转移矩阵的行：
  (\pi_{ij} \sim Dirichlet(\frac{\alpha}{\alpha + \kappa}\beta_1, \frac{\alpha}{\alpha + \kappa}\beta_2, \cdots, \frac{\alpha\beta_i + \kappa}{\alpha + \kappa}, \cdots))
  参数 (\kappa) 控制着过程的切换速率或长度尺度，这种改变使得所有状态的对角线元素权重更大，适用于先验未知段数量的分割问题。
• 块对角无限 HMM ：这是一个更通用的模型，允许将状态分组为块对角结构。当块大小为 1 时，就是粘性 iHMM；更大的块允许对状态进行无监督聚类。该模型用于从视频数据中无监督学习基于视图的对象模型，每个状态块或聚类对应一个对象，模型假设捕捉了相邻视频帧更可能对应同一对象不同视图而非不同对象的直觉。

以下是粘性 iHMM 和块对角无限 HMM 的比较表格：
| 模型 | 特点 | 适用场景 | 参数作用 |
| — | — | — | — |
| 粘性 iHMM | 增加转移矩阵对角线权重 | 分割问题，如说话人分割任务 | (\kappa) 控制切换速率 |
| 块对角无限 HMM | 状态分组为块对角结构 | 视频数据中对象模型的无监督学习 | 块大小决定聚类程度 |

6.3 iHMM 与皮特曼 - 约尔基础分布

从波利亚瓮方案可知，基础过程（即神谕瓮）类似于状态上的平均分布，影响给定序列长度下访问的期望状态数。狄利克雷过程在波利亚瓮中对分布有特定要求，其频率 - 排名图的右尾急剧下降，意味着出现一两次的颜色数量很少。

皮特曼 - 约尔过程是狄利克雷过程的两参数推广，引入了更多灵活性。其频率 - 排名图遵循幂律分布，在语言建模中能编码更现实的先验。皮特曼 - 约尔过程的棍子构造为：通过抽取 (\hat{\beta} i \sim Beta(1 - d, \alpha + id)) 并计算 (\beta_i = \hat{\beta}_i \prod {l = 1}^{i - 1} (1 - \hat{\beta}_l)) 来构建棍子，当 (d = 0) 时，皮特曼 - 约尔过程就是狄利克雷过程。

在 iHMM 中，用皮特曼 - 约尔分布替换狄利克雷过程基础分布相对简单，但束采样器的实现需要特别小心。因为皮特曼 - 约尔过程的棍子构造中，(\pi_{ij}) 比狄利克雷过程下降得更快，虽然允许更多小的转移概率，但 (\sum_{l} \pi_{il}) 下降得很慢。束采样器需要考虑所有 (\pi_{ij} > u_t) 的转移，可能需要将转移矩阵扩展到 (K) 个状态，直到 (\sum_{l > K} \pi_{il} < u_t) 以确保考虑了所有可能的状态。不过，由于 (\pi) 中的大部分条目很小，在扩展转移矩阵时可以保持矩阵条目的排序列表，在正向过滤步骤中，每次只需要遍历列表，直到 (u_t) 大于列表中的元素。

6.4 自回归 iHMM 和切换线性动态系统

许多时间序列使用自回归（AR）模型或线性动态系统（LDS）进行建模，这些模型假设过程的动态是线性的。通过在一组固定的线性动态之间切换，可以获得更强大的非线性推广，这种由 HMM 驱动切换的模型通常称为切换线性动态系统（SLDS）。

一个扩展的 SLDS 用 iHMM 替换了其底层的有限马尔可夫链，允许从数据中学习任意数量的动态。该扩展有两个原型变体：
- AR - iHMM ：假设观测值遵循自回归动态。
- 无限 SLDS ：假设只有部分连续变量被观测到，未观测变量遵循线性动态。

这两个模型的推理是在非参数模型的有限截断上实现的，与贝叶斯 SLDS 类似，但非参数模型为转移矩阵的行引入了额外的“先验的先验”。

以下是自回归 iHMM 和无限 SLDS 的特点对比表格：
| 模型 | 观测值假设 | 隐藏变量假设 | 推理方式 |
| — | — | — | — |
| AR - iHMM | 遵循自回归动态 | 无特殊假设 | 可使用塌缩吉布斯采样器 |
| 无限 SLDS | 部分连续变量被观测 | 未观测变量遵循线性动态 | 需显式表示隐藏状态 |

7. 结论

我们回顾了无限隐马尔可夫模型（iHMM）的分层波利亚瓮构造和分层狄利克雷过程构造，并证明了这两种形式定义了相同的分布。在此基础上，介绍了塌缩吉布斯采样器和基于动态规划的束采样器，这两种推理方法的实现难度仅比基于 HMM 的对应方法略高，使得 iHMM 成为解决学习 HMM 状态数问题的有吸引力的解决方案。

iHMM 有许多有趣的扩展方式，如允许输入和输出、粘性自转移、块结构模型、皮特曼 - 约尔基础分布、自回归和切换线性结构等。考虑到隐马尔可夫模型在时间序列和序列建模中的重要作用以及非参数方法的灵活性，无限隐马尔可夫模型在未来有很大的应用和扩展潜力。