隐马尔科夫模型的概率计算算法——解决评估问题

1. 问题背景与定位

在初识隐马尔可夫模型中，我们已经知道一个完整的隐马尔可夫模型 (HMM) 由三要素 $\lambda =(A,B,\pi )$ 定义。现在，我们面临 HMM 的第一个核心问题：评估 (Evaluation)。

评估问题描述：给定一个 HMM 模型 $\lambda =(A,B,\pi )$ 和一个观测序列 $O=(o_1,o_2,\dots ,o_T)$，我们如何计算这个模型产生出该特定观测序列的概率 $P(O|\lambda )$？

例如，在语音识别中，我们有两个 HMM 模型：${\lambda }_{\text{你好}}$ 和 ${\lambda }_{\text{再见}}$。当系统接收到一段语音（观测序列 $O$）时，我们可以分别计算 $P(O|{\lambda }_{\text{你好}})$ 和 $P(O|{\lambda }_{\text{再见}})$。哪个概率更高，就说明这段语音更可能对应哪个词语。因此，高效计算 $P(O|\lambda )$ 是许多 HMM 应用的基础。

接下来，我们将探索解决这个问题的三种方法：从最直观但最笨拙的“直接计算法”，到作为标准解决方案的“前向算法”与“后向算法”。

2. 直接计算法

2.1 直观理解

最符合直觉的方法是什么？就是把所有可能发生的情况都列出来，计算每一种情况的概率，然后把它们全部加起来。

具体来说，一个观测序列 $O=(o_1,\dots ,o_T)$ 是由一个我们看不见的隐藏状态序列 $Q=(q_1,\dots ,q_T)$ 生成的。我们不知道真实的隐藏状态序列是什么，但我们可以枚举所有可能的隐藏状态序列。

对于任何一条特定的隐藏状态序列，例如 $Q=(s_{i_1},s_{i_2},\dots ,s_{i_T})$，计算它和观测序列 $O$ 同时发生的联合概率 $P(O,Q|\lambda )$ 是很直接的：

$$P(O,Q|\lambda )=P(Q|\lambda )\times P(O|Q,\lambda )$$

让我们把这两项拆开：

1. $P(Q|\lambda )$: 隐藏状态序列自身的概率。根据马尔可夫链的性质，它等于：

   $$P(Q|\lambda )=P(q_1=s_{i_1})\times P(q_2=s_{i_2}|q_1=s_{i_1})\times \cdots \times P(q_T=s_{i_T}|q_{T-1}=s_{i_{T-1}})$$

   用我们的模型参数表示就是：

   $$P(Q|\lambda )={\pi }_{i_1}\cdot a_{i_1{i}_2}\cdot a_{i_2{i}_3}\cdots a_{i_{T-1}{i}_T}$$

2. $P(O|Q,\lambda )$: 在这个特定的隐藏状态序列下，生成观测序列 $O$ 的概率。根据观测独立性假设，它等于：

   $$P(O|Q,\lambda )=P(o_1|q_1=s_{i_1})\times P(o_2|q_2=s_{i_2})\times \cdots \times P(o_T|q_T=s_{i_T})$$

   用我们的模型参数表示就是：

   $$P(O|Q,\lambda )=b_{i_1}(o_1)\cdot b_{i_2}(o_2)\cdots b_{i_T}(o_T)$$

将两者相乘，就得到了一条特定路径的概率。

2.2 最终计算与问题

为了得到最终的 $P(O|\lambda )$，我们需要将所有可能的隐藏状态序列 $Q$ 的概率加起来：

$$P(O|\lambda )=\sum _{\text{所有可能的}Q}P(O,Q|\lambda )=\sum _{i_1,i_2,\dots ,i_T}{\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1{i}_2}{b}_{i_2}(o_2)\cdots a_{i_{T-1}{i}_T}{b}_{i_T}(o_T)$$

问题出在哪里？ 这个方法的计算量是灾难性的。

• 隐藏状态共有 $N$ 种。
• 序列长度为 $T$。
• 那么，所有可能的隐藏状态序列 $Q$ 的总数是 $N^T$。
• 计算每一条序列的概率，需要进行 $O(T)$ 次乘法。
• 总的计算复杂度大约是 $O(T\cdot N^T)$。

这是一个指数级别的复杂度。假设一个简单的词性标注任务，有 $N=50$ 个词性，一句话长度为 $T=20$。那么 $50^{20}$ 是一个天文数字，即使用最快的计算机也无法完成计算。因此，直接计算法在实践中是不可行的，它只是一种理论思想。

3. 前向算法

为了解决直接计算法的指数爆炸问题，我们需要一种更聪明的方法。前向算法利用动态规划 (Dynamic Programming) 的思想，通过存储和复用中间计算结果，将复杂度从指数级降到了多项式级。

3.1 直观理解与核心定义

直接计算法的问题在于，它进行了大量的重复计算。例如，在计算很多不同的长路径时，它们的前缀部分（比如从时刻1到时刻3的路径）被反复计算了无数次。

前向算法的核心思想是：在时刻 $t$ 到达某个状态 $s_i$ 的概率，可以由在时刻 $t-1$ 到达所有可能状态的概率之和来递推得到。我们不需要关心 $t-1$ 之前具体的路径是什么，只需要一个“总的概率”即可。

为此，我们定义一个关键变量——前向概率 (Forward Probability) ${\alpha }_t(i)$。

定义: 前向概率 ${\alpha }_t(i)$ 是指，在 HMM 模型 $\lambda$下，观测序列为 $o_1,o_2,\dots ,o_t$ 且在时刻 $t$ 处于隐藏状态 $s_i$ 的联合概率。

$${\alpha }_t(i)=P(o_1,o_2,\dots ,o_t,q_t=s_i|\lambda )$$

请仔细理解这个定义：它包含了到时刻 $t$ 为止的两个条件：

1. 已经看到了观测子序列 $o_1,\dots ,o_t$。
2. 在时刻 $t$，系统正好处于状态 $s_i$。

这个 ${\alpha }_t(i)$ 已经将所有能够“在时刻 $t$ 到达状态 $s_i$ 并产生观测子序列 $o_1,\dots ,o_t$”的路径概率全部加在一起了。

3.2 算法推导与步骤

前向算法的推导过程分为三步：

1. 初始化 $t=1$ 我们来计算初始时刻的前向概率 ${\alpha }_1(i)$。 根据定义：

$${\alpha }_1(i)=P(o_1,q_1=s_i|\lambda )$$

根据概率的链式法则，上式可以分解为：

$$\begin{array}{rlrl}{\alpha }_1(i)& =P(q_1=s_i|\lambda )\times P(o_1|q_1=s_i,\lambda )& & \text{（应用链式法则）}\\ & ={\pi }_i\cdot b_i(o_1)& & \text{（根据 π 和 B 的定义）}\end{array}$$

含义：在时刻1，处于状态 $s_i$ 并观测到 $o_1$ 的概率，就是“初始时选择状态 $s_i$ 的概率”乘以“状态 $s_i$ 生成观测 $o_1$ 的概率”。

2. 递推 for $t=1,\dots ,T-1$ 这是算法的核心。假设我们已经计算出了时刻 $t$ 的所有前向概率 ${\alpha }_t(j)$（对于所有 $j=1,\dots ,N$），我们如何计算时刻 $t+1$ 的前向概率 ${\alpha }_{t+1}(i)$？

根据定义： ${\alpha }_{t+1}(i)=P(o_1,\dots ,o_t,o_{t+1},q_{t+1}=s_i|\lambda )$

我们可以通过对时刻 $t$ 的所有可能状态 $s_j$ 进行边缘化 (Marginalization) 来展开它：

$$\begin{array}{rlrl}{\alpha }_{t+1}(i)& =\sum _{j=1}^{N}P(o_1,\dots ,o_t,q_t=s_j,o_{t+1},q_{t+1}=s_i|\lambda )& & \text{（对 qt 的所有可能状态求和）}\\ & =\sum _{j=1}^{N}P(o_1,\dots ,o_t,q_t=s_j|\lambda )\times P(o_{t+1},q_{t+1}=s_i|o_1,\dots ,o_t,q_t=s_j,\lambda )& & \text{（应用链式法则）}\end{array}$$

现在，我们利用 HMM 的两个核心假设来简化第二项：

• 齐次马尔可夫性: $q_{t+1}$ 只依赖于 $q_t$。所以 $P(q_{t+1}=s_i|\dots ,q_t=s_j)=P(q_{t+1}=s_i|q_t=s_j)=a_{ji}$。
• 观测独立性: $o_{t+1}$ 只依赖于 $q_{t+1}$。所以 $P(o_{t+1}|\dots ,q_{t+1}=s_i)=P(o_{t+1}|q_{t+1}=s_i)=b_i(o_{t+1})$。

将这些简化代入，第二项变为 $a_{ji}\cdot b_i(o_{t+1})$。 而第一项 $P(o_1,\dots ,o_t,q_t=s_j|\lambda )$ 正是我们定义的 ${\alpha }_t(j)$！

所以，我们得到了递推公式：

$${\alpha }_{t+1}(i)=\left[\sum _{j=1}^N{\alpha }_t(j)a_{ji}\right]\times b_i(o_{t+1})$$

含义: 要计算在时刻 $t+1$ 到达状态 $s_i$ 的总概率，我们首先要汇总所有从时刻 $t$ 的任意状态 $s_j$ 转移到 $s_i$ 的概率（方括号内的部分），然后再乘以 $s_i$ 生成观测 $o_{t+1}$ 的概率。

**3. 终止 当我们递推到最后一步，计算完所有 ${\alpha }_T(i)$ 后，我们如何得到最终的 $P(O|\lambda )$？ 根据定义，${\alpha }_T(i)=P(o_1,\dots ,o_T,q_T=s_i|\lambda )$。 我们要求的 $P(O|\lambda )=P(o_1,\dots ,o_T|\lambda )$，只需要对最后一个时刻的所有可能状态 $s_i$ 进行边缘化即可：

$$P(O|\lambda )=\sum _{i=1}^{N}P(o_1,\dots ,o_T,q_T=s_i|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$$

含义: 观测序列 $O$ 出现的总概率，等于在最后时刻，系统处于任何一个可能状态 $s_i$ 并完成观测 $O$ 的概率之和。

算法复杂度:

• 初始化需要 $O(N)$。
• 每一步递推，计算一个 ${\alpha }_{t+1}(i)$ 需要 $O(N)$ 次计算。要计算所有 $N$ 个状态，需要 $O(N^2)$。
• 总共有 $T-1$ 步递推。
• 总复杂度为 $O(N^{2}T)$，这是一个多项式级别的复杂度，相比 $O(T\cdot N^T)$ 是巨大的飞跃。

3.3 盒子和球模型

模型参数 $\lambda =(A,B,\pi )$:

• 状态 $Q=\{1,2,3\}$, $N=3$
• 观测 $V=\{\text{红, 白}\}$, $M=2$
• 初始概率 $\pi =(0.2,0.4,0.4)$
• 状态转移矩阵 $A=\left(\begin{array}{ccc}0.5& 0.2& 0.3\\ 0.3& 0.5& 0.2\\ 0.2& 0.3& 0.5\end{array}\right)$
• 观测概率矩阵 $B=\left(\begin{array}{cc}0.5& 0.5\\ 0.4& 0.6\\ 0.7& 0.3\end{array}\right)$ （第1列是红，第2列是白）

观测序列 $O=(\text{红, 白, 红})$, $T=3$。

1. 初始化 (t=1), 观测 $o_1=\text{红}$

• ${\alpha }_1(1)={\pi }_1\cdot b_1(o_1)=0.2\times 0.5=0.10$
• ${\alpha }_1(2)={\pi }_2\cdot b_2(o_1)=0.4\times 0.4=0.16$
• ${\alpha }_1(3)={\pi }_3\cdot b_3(o_1)=0.4\times 0.7=0.28$

2. 递推 (t=2), 观测 $o_2=\text{白}$

• 计算 ${\alpha }_2(1)$:

  $$\begin{array}{rl}{\alpha }_2(1)& =[{\alpha }_1(1)a_{11}+{\alpha }_1(2)a_{21}+{\alpha }_1(3)a_{31}]\times b_1(o_2)\\ & =[0.10\times 0.5+0.16\times 0.3+0.28\times 0.2]\times 0.5\\ & =[0.05+0.048+0.056]\times 0.5\\ & =0.154\times 0.5=0.077\end{array}$$

• 计算 ${\alpha }_2(2)$:

  $$\begin{array}{rl}{\alpha }_2(2)& =[{\alpha }_1(1)a_{12}+{\alpha }_1(2)a_{22}+{\alpha }_1(3)a_{32}]\times b_2(o_2)\\ & =[0.10\times 0.2+0.16\times 0.5+0.28\times 0.3]\times 0.6\\ & =[0.02+0.08+0.084]\times 0.6\\ & =0.184\times 0.6=0.1104\end{array}$$

• 计算 ${\alpha }_2(3)$:

  $$\begin{array}{rl}{\alpha }_2(3)& =[{\alpha }_1(1)a_{13}+{\alpha }_1(2)a_{23}+{\alpha }_1(3)a_{33}]\times b_3(o_2)\\ & =[0.10\times 0.3+0.16\times 0.2+0.28\times 0.5]\times 0.3\\ & =[0.03+0.032+0.14]\times 0.3\\ & =0.202\times 0.3=0.0606\end{array}$$

3. 递推 (t=3), 观测 $o_3=\text{红}$

• 计算 ${\alpha }_3(1)$:

  $$\begin{array}{rl}{\alpha }_3(1)& =[{\alpha }_2(1)a_{11}+{\alpha }_2(2)a_{21}+{\alpha }_2(3)a_{31}]\times b_1(o_3)\\ & =[0.077\times 0.5+0.1104\times 0.3+0.0606\times 0.2]\times 0.5\\ & =[0.0385+0.03312+0.01212]\times 0.5\\ & =0.08374\times 0.5=0.04187\end{array}$$

• 计算 ${\alpha }_3(2)$:

  $$\begin{array}{rl}{\alpha }_3(2)& =[{\alpha }_2(1)a_{12}+{\alpha }_2(2)a_{22}+{\alpha }_2(3)a_{32}]\times b_2(o_3)\\ & =[0.077\times 0.2+0.1104\times 0.5+0.0606\times 0.3]\times 0.4\\ & =[0.0154+0.0552+0.01818]\times 0.4\\ & =0.08878\times 0.4=0.035512\end{array}$$

• 计算 ${\alpha }_3(3)$:

  $$\begin{array}{rl}{\alpha }_3(3)& =[{\alpha }_2(1)a_{13}+{\alpha }_2(2)a_{23}+{\alpha }_2(3)a_{33}]\times b_3(o_3)\\ & =[0.077\times 0.3+0.1104\times 0.2+0.0606\times 0.5]\times 0.7\\ & =[0.0231+0.02208+0.0303]\times 0.7\\ & =0.07548\times 0.7=0.052836\end{array}$$

4. 终止

$$\begin{array}{rl}P(O|\lambda )& =\sum _{i=1}^3{\alpha }_3(i)\\ & ={\alpha }_3(1)+{\alpha }_3(2)+{\alpha }_3(3)\\ & =0.04187+0.035512+0.052836\\ & =0.130218\end{array}$$

4. 后向算法

后向算法是解决评估问题的另一种动态规划方法。它与前向算法对称，从序列的末尾开始，向前递推。

4.1 直观理解与核心定义

我们定义一个后向概率 (Backward Probability) ${\beta }_t(i)$。

定义: 后向概率 ${\beta }_t(i)$ 是指，在 HMM 模型 $\lambda$下，并且在时刻 $t$ 处于隐藏状态 $s_i$ 的条件下，观测到未来观测序列 $o_{t+1},o_{t+2},\dots ,o_T$ 的条件概率。

$${\beta }_t(i)=P(o_{t+1},o_{t+2},\dots ,o_T|q_t=s_i,\lambda )$$

请注意它与前向概率的区别：

• ${\alpha }_t(i)$ 是一个联合概率，包含了 $o_1..o_t$ 和 $q_t=s_i$。
• ${\beta }_t(i)$ 是一个条件概率，条件是 $q_t=s_i$，计算的是未来观测序列的概率。

4.2 算法推导与步骤

1. 初始化 $t=T$ 按照惯例，我们定义序列末尾的后向概率为一个基准值。

$${\beta }_T(i)=1,\text{for all }i=1,\dots ,N$$

含义: 在时刻 $T$ 已经处于状态 $s_i$ 的条件下，未来的观测序列（空序列）发生的概率自然是1。这是一个递归的起点。

2. 递推 for $t=T-1,\dots ,1$ 假设我们已经计算出了时刻 $t+1$ 的所有后向概率 ${\beta }_{t+1}(j)$，我们来计算时刻 $t$ 的后向概率 ${\beta }_t(i)$。

根据定义，${\beta }_t(i)=P(o_{t+1},\dots ,o_T|q_t=s_i,\lambda )$。 我们通过对时刻 $t+1$ 的所有状态 $s_j$ 进行边缘化来展开：

$$\begin{array}{rl}{\beta }_t(i)& =\sum _{j=1}^{N}P(o_{t+1},\dots ,o_T,q_{t+1}=s_j|q_t=s_i,\lambda )\\ & =\sum _{j=1}^{N}P(q_{t+1}=s_j|q_t=s_i,\lambda )\times P(o_{t+1}|q_{t+1}=s_j,q_t=s_i,\lambda )\times P(o_{t+2},\dots ,o_T|o_{t+1},q_{t+1}=s_j,q_t=s_i,\lambda )\end{array}$$

再次利用HMM假设简化：

• $P(q_{t+1}=s_j|q_t=s_i,\lambda )=a_{ij}$
• $P(o_{t+1}|\dots ,q_{t+1}=s_j,\dots )=P(o_{t+1}|q_{t+1}=s_j)=b_j(o_{t+1})$
• $P(o_{t+2},\dots ,o_T|\dots ,q_{t+1}=s_j,\dots )=P(o_{t+2},\dots ,o_T|q_{t+1}=s_j)={\beta }_{t+1}(j)$

将它们代入，得到递推公式：

$${\beta }_t(i)=\sum _{j=1}^N{a}_{ij}\cdot b_j(o_{t+1})\cdot {\beta }_{t+1}(j)$$

3. 终止 后向算法也可以用来计算总概率 $P(O|\lambda )$。我们可以在任意时刻 $t$ 将前向和后向概率结合起来。最简单的是在 $t=1$ 时：

$$P(O|\lambda )=\sum _{i=1}^{N}P(o_1,\dots ,o_T,q_1=s_i|\lambda )$$

我们将 $P(o_1,\dots ,o_T,q_1=s_i|\lambda )$ 分解为 $P(q_1=s_i)\cdot P(o_1|q_1=s_i)\cdot P(o_2,\dots ,o_T|q_1=s_i)$。 这三项正好是 ${\pi }_i$, $b_i(o_1)$, 和 ${\beta }_1(i)$。 所以：

$$P(O|\lambda )=\sum _{i=1}^N{\pi }_i\cdot b_i(o_1)\cdot {\beta }_1(i)$$

这个结果必须与前向算法得到的结果完全相同。

5. 补充概率概念

5.1 单个状态的概率：${\gamma }_t(i)$

定义: 给定模型 $\lambda$ 和 整个 观测序列 $O$，在时刻 $t$ 处于状态 $s_i$ 的概率。

$${\gamma }_t(i)=P(q_t=s_i|O,\lambda )$$

这个概率回答了“在看到了所有证据之后，我们回头看，在时刻 $t$ 系统处于状态 $s_i$ 的可能性有多大？”

推导: 根据条件概率的定义 $P(A|B)=P(A,B)/P(B)$：

$${\gamma }_t(i)=\frac{P(q_t=s_i,O|\lambda )}{P(O|\lambda )}$$

我们来分析分子 $P(q_t=s_i,O|\lambda )$。它可以被看作是“过去”、“现在”和“未来”的结合：

• “过去”：观测到 $o_1,\dots ,o_t$ 并且在时刻 $t$ 处于状态 $s_i$。这正是前向概率 ${\alpha }_t(i)$。
• “未来”：在时刻 $t$ 处于状态 $s_i$ 的条件下，观测到 $o_{t+1},\dots ,o_T$。这正是后向概率 ${\beta }_t(i)$。

因此， $P(q_t=s_i,O|\lambda )={\alpha }_t(i){\beta }_t(i)$。 代入可得：

$${\gamma }_t(i)=\frac{{\alpha }_t(i){\beta }_t(i)}{P(O|\lambda )}=\frac{{\alpha }_t(i){\beta }_t(i)}{{\sum }_{j=1}^N{\alpha }_t(j){\beta }_t(j)}$$

意义: ${\gamma }_t(i)$ 在模型训练中，可以被看作是在时刻 $t$ 访问状态 $i$ 的期望次数。对所有时刻求和 ${\sum }_{t=1}^T{\gamma }_t(i)$ 就是在整个序列中访问状态 $i$ 的总期望次数。

5.2 两个状态的转移概率：${\xi }_t(i,j)$

定义: 给定模型 $\lambda$ 和 整个 观测序列 $O$，在时刻 $t$ 处于状态 $s_i$ 且在时刻 $t+1$ 处于状态 $s_j$ 的概率。

$${\xi }_t(i,j)=P(q_t=s_i,q_{t+1}=s_j|O,\lambda )$$

这个概率回答了“在看到了所有证据之后，在时刻 $t$ 到 $t+1$ 之间，发生一次从 $s_i$到$s_j$转移的可能性有多大？”

推导: 同样地，${\xi }_t(i,j)=\frac{P(q_t=s_i,q_{t+1}=s_j,O|\lambda )}{P(O|\lambda )}$。 分子可以被分解为：

$$P(\dots )=\underset{{\alpha }_t(i)}{\underbrace{P(o_1,\dots ,o_t,q_t=s_i)}}\times \underset{a_{ij}}{\underbrace{P(q_{t+1}=s_j|q_t=s_i)}}\times \underset{b_j(o_{t+1})}{\underbrace{P(o_{t+1}|q_{t+1}=s_j)}}\times \underset{{\beta }_{t+1}(j)}{\underbrace{P(o_{t+2},\dots ,o_T|q_{t+1}=s_j)}}$$

将所有部分组合起来，得到：

$${\xi }_t(i,j)=\frac{{\alpha }_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)}{P(O|\lambda )}=\frac{{\alpha }_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)}{{\sum }_{k=1}^N{\sum }_{l=1}^N{\alpha }_t(k)a_{kl}{b}_l(o_{t+1}){\beta }_{t+1}(l)}$$

意义: ${\xi }_t(i,j)$ 可以被看作是在时刻 $t$ 从状态 $i$ 转移到状态 $j$ 的期望次数。对时间求和 ${\sum }_{t=1}^{T-1}{\xi }_t(i,j)$ 是在整个序列中从 $i$ 转移到 $j$ 的总期望次数。

这两个量，${\gamma }_t(i)$ 和 ${\xi }_t(i,j)$，是鲍姆-韦尔奇 (Baum-Welch) 学习算法的核心，它们被用来在给定观测数据的情况下，迭代地更新模型参数 $A,B,\pi$。

6. 总结一下

本次我们深入探讨了 HMM 的概率评估问题。

• 直接计算法 思路简单，但因 $O(T\cdot N^T)$ 的指数复杂度而不可行。
• 前向算法 利用动态规划，定义了前向概率 ${\alpha }_t(i)$，自前向后递推，将复杂度降至 $O(N^{2}T)$，是解决评估问题的标准方法。
• 后向算法 是前向算法的镜像，定义了后向概率 ${\beta }_t(i)$，自后向前递推，同样可以计算总概率，并且是计算其他重要概率的必要工具。
• 组合概率 ${\gamma }_t(i)$ 和 ${\xi }_t(i,j)$ 利用前向和后向概率，让我们能深入洞察序列内部的动态，是连接评估问题和学习问题的桥梁。

至此，我们已经掌握了HMM三大问题中的第一个。下一个自然是解决隐马尔科夫模型的学习问题了。