隐马尔科夫模型

HMM 算法小结

HMM 算法背景&概念简介

基本概念

隐马尔科夫模型是关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测从而产生观测随机序列的过程，隐藏的马尔科夫链随机生成的状态的序列，称为状态序列；每个状态生成一个规则，而由此产生的观测的随机序列称为观测序列。序列的每一个位置又可以看作是一个时刻。

HMM 由初始概率分布、状态转移概率分布和观测概率分布确定。

HMM 的形式定义如下：

设 $Q$ 是所有可能的状态的集合，$V$ 是所有可能的观测的集合：

$$Q = \{q_1,q_2,\cdots,q_N \}, \ V = \{v_1,v_2,\cdots, v_M \}$$

其中，$N$ 是可能的状态数，$M$ 是可能的观测数。

$I$ 是长度为 $T$ 的状态序列， $O$ 是对应的观测序列：

$$I = (i_1,i_2,\cdots,i_T), \ O = (o_1,o_2,\cdots, o_T)$$

$A$ 是状态转移概率矩阵：

$$A = [a_{ij}]_{N \times N}$$

其中，

$$a_{ij} = P(i_{t+1} = q_j | i_t = q_i), \ i = 1,2,\cdots,N; \ j = 1,2,\cdots,N$$

是在时刻 $t$ 处于状态 $q_i$ 的条件下在时刻 $t+1$ 转移到状态 $q_j$ 的概率。

$B$ 是观测概率矩阵：

$$B = [b_j(k)]_{N \times M}$$

其中， $$b_j(k) = P(o_t = v_k | i_t = q_j), \ k = 1,2,\cdots,M; \ j = 1,2,\cdots, N$$

是在时刻 $t$ 处于状态 $q_j$ 的条件下生成观测 $v_k$ 的概率。

$\pi$ 是初始态概率向量： $$\pi = (\pi_i)$$ 其中， $$\pi_i = P(i_1 = q_i), \ i = 1,2,\cdots,N$$ 是时刻 $t=1$ 处于状态 $q_i$ 的概率。

隐马尔科夫模型由初始状态概率向量 $\pi$、状态转移矩阵 $A$ 和观测概率矩阵 $B$ 决定。$\pi$ 和 $A$ 决定状态序列，$B$ 决定观测序列。因此，隐马尔科夫模型 $\lambda$ 可以用三元符号表示，即

$$\lambda = (A,B,\pi)$$

$A,B,\pi$ 称为隐马尔科夫模型的三要素。

状态转移概率矩阵 $A$ 与初始状态概率向量 $\pi$ 确定了隐藏的马尔科夫链，生成不可观测的状态序列。观测概率矩阵 $B$ 确定了如何从状态生成观测，与状态序列综合确定了如何产生观测序列。

基本假设

  1. 齐次马尔科夫性假设，即假设隐藏的马尔科夫链在任意时刻 $t$ 的状态只依赖于前一时刻的状态，与其他时刻的状态及观测无关，业余时刻 $t$ 无关：
$$P(i_t | i_{t-1}, o_{t-1}, \cdots, i_1, o_1) = P(i_t|i_{t-1}), \ t = 1,2,\cdots, T$$
  1. 观测独立性假设，即假设任意时刻的观测只依赖于该时刻的马尔科夫链的状态，与其他观测以状态无关。

实例

假设有 4 个盒子，每个盒子里装有红、白两种颜色的球，盒子里的红、白球数由下表列出

按照下面方法抽球，产生一个球的颜色的观测序列：

  1. 开始，从 4 个盒子里以等概率随机选取 1 个盒子，规则是：

    • 如果当前盒子是盒子 1，那么下一个盒子一定是盒子 2；
    • 如果当前是盒子 2 或 3，那么分别以概率 0.4 和 0.6 转移到左边或右边的盒子；
    • 如果当前是盒子 4，那么各以 0.5 的概率停留在 4 或转移到盒子 3；

  2. 确定转移的盒子后，再从这个盒子里随机抽一个球，记录其颜色，放回；

  3. 如此下去，重复进行 5 次，得到一个球的颜色的观测序列：
$$O = (红，红，白，白，红)$$

在这个过程中，观察者只能观测到球的颜色的序列，观测不到球是从哪个盒子取出的，也就是观测不到盒子的序列。

在这个例子中有两个随机序列，一个是盒子的序列（状态序列），一个是球的颜色的观测序列（观测序列）。前者是隐藏的，只有后者是可观测的。这是一个隐马尔科夫模型的例子。根据所给条件，可以明确状态集合、观测集合、序列长度以及模型的三要素。

盒子对应状态，状态的集合是：

$$Q = \{盒1，盒2，盒3，盒4 \}, \ N = 4$$

球的颜色对应观测。观测的集合是：

$$V = {红，白}, \ M = 2$$

状态序列和观测序列长度 $T = 5$。

初始概率分布为

$$\pi = (0.25,0.25,0.25,0.25)^T$$

状态转移概率分布（transition matrix）为：

观测概率分布为

观测序列的生成过程

根据 HMM 定义，可以将一个长度为 $T$ 的观测序列 $O = (o_1,o_2,\cdots, o_T)$ 的生成过程描述如下。

算法：

输入：$\lambda = (A,B,\pi)$，$T$ 输出：$O = (o_1,o_2,\cdots, o_T)$

  1. 按照初始状态分布 $\pi$ 产生状态 $i_1$;
  2. 令 $t=1$;
  3. 按照状态 $i_t$ 的观测概率分布 $b_{i_t}(k)$ 生成 $o_t$；
  4. 按照状态 $i_t$ 的状态转移概率分布 $\{a_{{i_t}i_{t+1}} \}$ 产生状态 $i_{t+1}$， $i_{t+1} = 1,2, \cdots, N$；
  5. 令 $t=t+1$；如果 $t<T$，转步骤 3；否则停止。

3个基本问题

  1. 概率计算问题。 在给定模型 $\lambda = (A,B,\pi)$ 和观测序列 $O$，计算在模型 $\lambda$ 下观测序列 $O$ 出现的概率 $P(O|\lambda)$.
  2. 学习问题。 已知观测序列 $O = (o_1,o_2, \cdots, o_T)$，估计模型 $\lambda$ 的参数，使得在该模型下观测序列概率 $P(O|\lambda)$ 最大。也就是用极大似然估计来估计参数。
  3. 预测问题，也称为解码（decoding）问题。 已知模型 $\lambda$ 和观测序列 $O$，求对给定观测序列条件概率 $P(I|O)$ 最大的状态序列 $I = (i_1,i_2,\cdots, i_T)$。即给定观测序列，求最有可能的对应的状态序列（隐藏）。

HMM 算法详解

此处介绍之前说的第一个问题：观测序列概率 $P(O|\lambda)$ 的计算方法。前向和后向算法。

首先直接计算法，也就是直接用概率公式计算。通过列举所以可能的长度为 $T$ 的状态序列 $I$，求各个状态序列 $I$ 与观测序列 $O$ 的联合概率 $P(O,I|\lambda)$，然后对所有可能的状态序列求和，得到 $P(O|\lambda)$。这样计算量极大，复杂度达到 $O(TN^T)$ 阶的，所以不可行。

1. 前向算法

定义前向概率。

给定隐马尔科夫模型 $\lambda$，定义到时刻 $t$ 部分观测序列为 $o_1, o_2, \cdots, o_t$，且状态为 $q_i$ 的概率为前向概率，记作

$$\alpha_t(i) = P(o_1,o_2,\cdots,o_t,i_t = q_i|\lambda)$$

可以递推地求得前向概率 $\alpha_t(i)$ 以及观测序列概率 $P(O|\lambda)$

算法：

输入：$\lambda$，$O$ 输出：观测序列概率 $P(O|\lambda)$

  1. 初值
$$\alpha_1(i) = \pi_ib_i(o_1), \ i = 1,2,\cdots,N$$

  1. 递推

    对 $t=1,2,\cdots,T-1$，

$$\alpha_{t+1}(i) = \Bigg[\sum_{j=1}^{N}\alpha_t(j)a_{ji} \Bigg]b_i(o_{t+1}), \ i=1,2,\cdots,N$$
  1. 终止
$$P(O|\lambda) = \sum_{i=1}^{N}\alpha_T(i)$$

解释：

步骤 1，初始化前向概率，也就是初始时刻（t=1）的状态 $i_1 = q_i$ 和观测 $o_1$ 的联合概率。

步骤 2，是前向概率的递推公式，计算到时刻 $t+1$ 部分观测序列为 $o_1,o_2,\cdots, o_t,o_{t+1}$，且在时刻 $t+1$ 处于状态 $q_i$ 的前向概率。

  1. $\alpha_t(j)$ 指的是到 $t$ 时刻观测到 $o_1,o_2,\cdots,o_t$ 并在时刻 $t$ 处于状态 $q_j$ 的前向概率

  2. $\alpha_t(j)a_{ji}$ 就是到 $t$ 时刻观测到 $o_1,o_2,\cdots,o_t$ 并在时刻 $t$ 处于状态 $q_j$ 而在时刻 $t+1$ 到达状态 $q_i$ 的联合概率。

  3. $\sum_{j=1}^{N}\alpha_t(j)a_{ji}$ 就等于对在时刻 $t$ 的所有 $N$ 个状态下求和，最终得到时刻 $t$ 观测为 $o_1,o_2,\cdots,o_t$ 并在时刻 $t+1$ 处于状态 $q_i$ 的联合概率。

  4. 最后再与观测概率 $b_i(o_{t+1})$ 的乘积恰好是到时刻 $t+1$ 观测到 $o_1,o_2,\cdots, o_t,o_{t+1}$ 并在时刻 $t+1$ 处于状态 $q_i$ 的前向概率 $\alpha_{t+1}(i)$

步骤（3）给出 $P(O|\lambda)$ 的计算公式，因为

$$ \begin{aligned} \alpha_T(i) &= P(o_1,o_2,\cdots,o_T,i_T = q_i|\lambda) \end{aligned} $$

所以

$$P(O|\lambda) = \sum_{i=1}^{N}\alpha_T(i)$$

如上图所示，实际上前向算法是根据状态序列的路径结构递推计算 $P(O|\lambda)$ 的算法。其实通过转移矩阵 $A$ 和观测概率分布 $B$ 来分别对从时刻 $t$ 到 $t+1$ 的状态和观测值的递推，从而局部计算前向概率，然后再递推到全局，得到 $P(O|\lambda)$。

2. 后向算法

给定隐马尔科夫模型 $\lambda$，定义在时刻 $t$ 状态为 $q_i$ 的条件下，从 $t+1$ 到 $T$ 的部分观测序列为 $o_{t+1}, o_{t+2}, \cdots, o_T$ 的概率为后向概率，记作

$$\beta_t(i) = P(o_{t+1}, o_{t+2}, \cdots, o_T|i_t = q_i,\lambda)$$

可以用递推的方法求得后向概率 $\beta_t(i)$ 及观测序列概率 $P(O|\lambda)$。

算法：

输入：$\lambda$，$O$ 输出：$P(O|\lambda)$

  1. 初始化： $\beta_T(i) =1, \ i=1,2,\cdots,N$

  2. 对 $t = T-1,T-2,\cdots,1$ $$\beta_t(i) = \sum_{j=1}^{N}a_{ij}b_j(o_{t+1})\beta_{t+1}(j), \ i=1,2,\cdots,N$$

  3. 计算 $P(O|\lambda)$
    $$P(O|\lambda) = \sum_{i=1}^N \pi_ib_i(o_1)\beta_1(i)$$

解释：

步骤 1， 初始化后向概率，对最终时刻的所有状态 $q_i$ 规定 $\beta_T(i) = 1$。

步骤 2， 后向概率的递推公式，为了计算在时刻 $t$ 状态为 $q_i$ 的条件下，时刻 $t+1$ 之后的观测序列为 $o_{t+1}, o_{t+2}, \cdots, o_T$ 的后向概率 $\beta_t(i)$，只需考虑在时刻 $t+1$ 所有可能的 $N$ 个状态 $q_j$ 的转移概率 （$a_{ij}$），以及在此状态下的观测 $o_{t+1}$ 的观测概率（$b_j(o_{t+1})$），然后考虑状态 $q_j$ 之后的观测序列的后向概率 （$\beta_{t+1}(j)$）。

步骤 3 求 $P(P|\lambda)$ 的思路与步骤 2 一致，只是初始概率 $\pi_i$ 代替转移概率。

利用前向概率和后向概率的定义可以将观测序列概率 $P(O|\lambda)$ 统一写成

$$P(O|\lambda) = \sum_{i=1}^{N}\sum_{j=1}^{N} \alpha_t(i)\alpha_{ij}b_j(o_{t+1})\beta_{t+1}(j), \ t=1,2,\cdots,T-1$$

3. 概率与期望值的计算

利用前向概率和后向概率，可以得到关于单个状态和两个状态概率的计算公式。

  1. 给定模型 $\lambda$ 和观测 $O$，在时刻 $t$ 处于状态 $q_i$ 的概率。记为

    $$\gamma_t(i) = P(i_t = q_i|O,\lambda)$$

    可以通过前后向概率计算。事实上，

    $$\gamma_t(i) = P(i_t = q_i|O,\lambda) = \frac{P(i_t = q_i,O|\lambda)}{P(O|\lambda)}$$

    前向概率：

    $$\alpha_t(i) = P(o_1,o_2,\cdots,o_t,i_t = q_i|\lambda)$$

    后向概率：

    $$\beta_t(i) = P(o_{t+1}, o_{t+2}, \cdots, o_T|i_t = q_i,\lambda)$$

    由定义可知：

    $$\alpha_t(i)\beta_t(i) = P(i_t = q_i,O|\lambda)$$

    于是得到

    $$\gamma_t(i) = \frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^{N}\alpha_t(j)\beta_t(j)}$$

  2. 给定模型 $lambda$ 和观测 $O$，在时刻 $t$ 处于状态 $q_i$ 且在时刻 $t+1$ 处于状态 $q_j$ 的概率。记为

    $$\xi_t(i,j) = P(i_t = q_i, i_{t+1 = q_j} | O,\lambda)$$

    可以通过前向后向概率计算：

    $$ \begin{aligned} \xi_t(i,j) &= \frac{P(i_t = q_i, i_{t+1 = q_j} ,O| \lambda)}{P(O|\lambda)} \\ &= \frac{P(i_t = q_i, i_{t+1 = q_j} ,O| \lambda)}{\sum_{i=1}^{N}\sum_{j=1}^{N}P(i_t = q_i, i_{t+1 = q_j} ,O| \lambda)} \end{aligned} $$

  3. 将 $\gamma_t(i)$ 和 $\xi_{t}(i,j)$ 对各个时刻 $t$ 求和，可以得到一些有用的期望值。

    1. 在观测 $O$ 下状态 $i$ 出现的期望值： $$\sum_{t=1}^{T}\gamma_t(i)$$
    2. 在观测 $O$ 下由状态 $i$ 转移的期望值： $$\sum_{t=1}^{T-1}\gamma_t(i)$$
    3. 在观测 $O$ 下由状态 $i$ 转移到状态 $j$ 的期望值： $$\sum_{t=1}^{T-1}\xi_t(i,j)$$

HMM 算法学习&实现

监督学习方法

假设已给训练数据包含 $S$ 个长度相同的观测序列和对应的状态序列 $\{ (O_1,I_1), (O_2,I_2), \cdots, (O_S,I_S) \}$，那么可以利用极大似然估计来估计隐马尔科夫模型的参数。

  1. 转移概率 $a_{ij}$ 的估计 设样本中时刻 $t$ 处于状态 $t+1$ 转移到状态 $j$ 的频数为 $A_{ij}$ ，那么状态转移概率 $a_{ij}$ 的估计是

    $$\hat a_{ij} = \frac{A_{ij}}{\sum_{j=1}^{N}A_{ij}}, \ i=1,2,\cdots,N; \ j = 1,2,\cdots,N$$

  2. 观测概率 $b_j(k)$ 的估计 设样本中状态为 $j$ 并观测为 $k$ 的频数为 $B_{jk}$，那么状态为 $j$ 观测为 $k$ 的概率 $b_j(k)$ 的估计是

    $$\hat b_j(k) = \frac{B_{jk}}{\sum_{k=1}^{M}B_{jk}}, \ j=1,2,\cdots,N; \ k = 1,2,\cdots, M$$

  3. 初始状态概率 $\pi_i$ 的估计 $\hat \pi_i$ 为 $S$ 个样本中初始状态为 $q_i$ 的频率

无监督学习（Baum-Welch 算法）

假设给定数据只包含了 $S$ 个长度为 $T$ 的观测序列 $\{ O_1, O_2, \cdots, O_S \}$ 而没有对应的状态序列，目标是学习隐马尔科夫模型 $\lambda = (A,B,\pi)$ 的参数。我们将观测是序列数据看作观测数据 $O$，状态序列数据看作不可观测的隐数据 $I$，那么 HMM 事实上是一个含隐变量的概率模型。

$$P(O|\lambda) = \sum_{I}P(O|I,\lambda)P(I|\lambda)$$

其参数学习可以由 EM 算法实现

  1. 确定完全数据的对数似然函数

    所有观测数据写成 $O = (o_1,o_2,\cdots,o_T)$，所有隐数据写成 $I = (i_1,i_2,\cdots,i_T)$，完全数据是 $(O,I) = (o_1,o_2,\cdots,o_T,i_1,i_2,\cdots,i_T)$。完全数据的对数似然函数是 $\log P(O,I|\lambda)$。

  2. EM 算法的 E 步：求 $Q$ 函数

    $$Q(\lambda,\overline \lambda) = \sum_{I} \log P(O,I|\lambda)P(O,I|\overline \lambda)$$ 其中， $\overline \lambda$ 是隐马尔科夫模型参数的当前估计值，$\lambda$ 是要极大化的隐马尔科夫模型参数。

    $$P(O,I|\lambda) = \pi_{i_1}b_{i_1}(o_1)a_{i_1i_2}b_{i_2}(o_2) \cdots a_{i_{T-1}i_T}b_{i_T}(o_T)$$

    于是函数 $Q(\lambda,\overline \lambda)$ 可以写成：

    $$ \begin{aligned} Q(\lambda,\overline \lambda) &= \sum_{I} \log \pi_{i_1} P(O,I|\overline \lambda) + \sum_{I} \Bigg(\sum_{t=1}^{T-1} \log a_{i_ti_{t+1}} \Bigg) P(O,I|\overline \lambda)+ \\ & \sum_{I} \Bigg(\sum_{t=1}^{T} \log b_{i_t}(o_t) \Bigg) P(O,I|\overline \lambda) \end{aligned} $$

    式中求和都是对所有数据的序列总长度 $T$ 进行的。

  3. EM 算法的 M 步：极大化 $Q$ 函数 $Q(\lambda,\overline \lambda)$ 求模型参数 $A,B,\pi$

    由于要极大化的参数在上面的式子中单独的出现在 3 个项中，所以只需要对各项分别极大化即可。

    • 第 1 项可以写成：

      $$\sum_{I} \log \pi_{i_1} P(O,I|\overline \lambda) = \sum_{i=1}^{N} \log \pi_i P(O,i_1 = i|\overline \lambda)$$

      注意到 $\pi_i$ 满足约束条件 $\sum_{i=1}^{N} \pi_i = 1$，利用拉格朗日乘子法，写出拉格朗日函数：

      $$\sum_{i=1}^{N} \log \pi_iP(O,i_1=i|\overline \lambda) + \gamma \Bigg(\sum_{i=1}^{N}\pi_i-1 \Bigg)$$

      对其求偏导并令其结果等于 0

      $$\frac{\partial}{\partial\pi_i} \Bigg [\sum_{i=1}^{N} \log \pi_i P(O,i_1=i|\overline \lambda)+ \gamma \Bigg (\sum_{i=1}^{N} \pi_i -1\Bigg)\Bigg] = 0$$

      得到

      $$P(O,i_1 = i|\overline \lambda) + \gamma \pi_i = 0$$

      对 $i$ 求和得到

      $$ \begin{aligned} \sum_{i=1}^{N} P(O,i_1 &= i|\overline \lambda) = - \sum_{i=1}^{N} \gamma \pi_i \\ \gamma &= -\sum_{i=1}^{N} P(O,i_1 = i|\overline \lambda) \\ &= -P(O|\overline \lambda) \end{aligned} $$

      带入上面的式子得到

      $$\pi_i = \frac{P(O,i_1=i|\overline \lambda)}{P(O|\overline \lambda)}$$

    • 第 2 项可以写成

      $$\sum_{I} \Bigg(\sum_{t=1}^{T-1} \log a_{i_ti_{t+1}} \Bigg) P(O,I|\overline \lambda) = \sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T-1} \log a_{ij} P(O,i_t=i,i_{t+1} = j |\overline \lambda)$$

      类似第 1 项，应用具有条件 $\sum_{j=1}^{N}a_{ij} = 1$ 的拉格朗日乘子法可以求出

      拉格朗日函数

      $$\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T-1} \log a_{ij} P(O,i_t=i,i_{t+1} = j |\overline \lambda)+\gamma(\sum_{j=1}^{N}a_{ij} -1) $$ 求偏导等于 0

      $$P(O,i_t=i,i_{t+1} = j |\overline \lambda) + \gamma a_{ij} = 0$$

      求出 $a_{ij}$

      $$a_{ij} = \frac{ \sum_{t=1}^{T-1}P(O,i_t=i,i_{t+1} = j |\overline \lambda) }{\sum_{t=1}^{T-1}P(O,i_t=i|\overline \lambda)}$$

    • 第 3 项为

      $$\sum_{I} \Bigg(\sum_{t=1}^{T} \log b_{i_t}(o_t) \Bigg) P(O,I|\overline \lambda) = \sum_{j=1}^{N}\sum_{t=1}^{T} \log b_j(o_t)P(O,i_t=j|\overline \lambda)$$

      同样使用拉格朗日乘子法，约束条件为 $\sum_{k=1}^{M}b_j(k)=1$。注意，只有在 $o_t = v_k$ 时 $b_j(o_t)$ 对 $b_j(k)$ 的偏导数才不为 0，以 $I(o_t = v_k)$ 表示。求得

      $$b_j(k) = \frac{\sum_{t=1}^{T}P(O,i_t=j|\overline\lambda)I(o_t=v_k)}{\sum_{t=1}^{T}P(O,i_t=j|\overline \lambda)}$$

Baum-Welch 模型参数估计公式

将上面几个式子的各概率分别用 $\gamma_t(i), \xi_t(i,j)$ 表示，则可将相应的公式写成：

$$a_{ij} = \frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)}$$$$b_j(k) = \frac{\sum_{t=1,o_t = v_k}^{T}\gamma_t(j)}{\sum_{t=1}^{T}\gamma_t(j)}$$$$\pi_i = \gamma_1(i)$$

其中， $\gamma_t(i), \xi_t(i,j)$ 分别由前面推导的公式给出。

算法：

HMM 预测算法

近似算法

在每个时刻 $t$ 选择在该时刻最有可能出现的状态 $i_t^*$，从而得到一个状态序列 $I^* = (i_1^*,i_2^*, \cdots, i_T^*)$，将它作为预测的结果。

给定隐马尔科夫模型 $\lambda$ 和观测序列 $O$，在时刻 $t$ 处于状态 $q_i$ 的概率 $\gamma_t(i)$ 是 $$\gamma_t(i) = \frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^{N}\alpha_t(j)\beta_t(j)}$$

在每一时刻 $t$ 最有可能的状态 $i_t^*$ 是 $$i_t^* = \arg \max_{1 \leq i \leq N} [\gamma_t(i)], \ t = 1,2,\cdots,T$$

从而得到状态序列 $I^* = (i_1^*,i_2^*, \cdots, i_T^*)$。

优点： 计算简单 缺点： 不能保证预测的状态序列整体是最有可能的状态序列，有可能对于某些 $i,j,a_{ij} = 0$。

维特比算法

维特比算法（Viterbi）实际是用动态规划解隐马尔科夫模型预测问题，即用动态规划求概率最大路径（最优）。这时一条路径对应着一个状态序列。

根据动态规划原理，最优路径具有这样的特性：如果最优路径在时刻 $t$ 通过结点 $i_t^*$，那么这一路径从结点 $i_t^*$ 到终点 $i_T^*$ 的部分路径来说，必须是最优的。

那么我们只需要从时刻 $t=1$ 开始，递推地计算咋时刻 $t$ 状态为 $i$ 的各条部分路径的最大概率，直至得到时刻 $t=T$ 状态为 $i$ 的各条路径的最大概率。时刻 $t=T$ 的最大概率即为最优路径的概率 $P^*$，最优路径的终结点 $i_T^*$ 也同时得到。之后，为了找出最优路径的各个结点，从终结点 $i_T^*$ 开始，由后向前逐步求得结点 $i_{T-1}^*, \cdots, i_1^*$，得到最优路径 $I^* = (i_1^*,i_2^*, \cdots, i_T^*)$。这就是维特比算法。

先导入两个变量 $\delta$ 和 $\Psi$。定义在时刻 $t$ 状态为 $i$ 的所有单个路径 $(i_1,i_2,\cdots,i_t)$ 中概率最大值为

$$\delta_t(i) = \max_{i_1,i_2,\cdots,i_{t-1} } P(i_t = i,i_{t-1}, \cdots, i_1, o_t, \cdots, o_1 |\lambda), \ i=1,2,\cdots,N$$

由定义可得变量 $\delta$ 的递推公式： $$ \begin{aligned} \delta_{t+1}(i) &= \max_{i_1,i_2,\cdots,i_t } P(i_{t+1} = i,i_t, \cdots, i_1, o_{t+1}, \cdots, o_1 |\lambda) \\ &= \max_{1 \leq j \leq N} [\delta(j) a_{ji}] b_i(o_{t+1}), \ i=1,2,\cdots,N; \ t=1,2,\cdots,T-1 \end{aligned} $$

跟以前一样，根据 $\delta_t(i)$ 与转移概率相乘得到在时刻 $t$ 状态由 $j$ 转移到 $i$ 的概率，然后再乘观测概率得到在 $t+1$ 时刻状态为 $i$ 的所有单个路径中的概率最大值。

定义在时刻 $t$ 状态为 $i$ 的所有单个路径 $(i_1,i_2,\cdots,i_{t-1},i)$ 中概率最大的路径的第 $t-1$ 个结点为

$$\Psi_t(i) = \arg \max_{1 \leq j \leq N} [\delta_{t-1}(j)a_{ji}], \ i=1,2,\cdots,N$$

维特比算法：

输入：模型 $\lambda = (A,B,\pi)$ 和观测 $O = (o_1,o_2,\cdots,o_T)$ 输出：最优路径 $I^* = (i_1^*,i_2^*, \cdots, i_T^*)$

  1. 初始化
$$\delta_1(i) = \pi_ib_i(o_1), i = 1,2,\cdots,N$$$$\Psi_1(i) = 0, \ i=1,2,\cdots,N$$

  1. 递推

    对 $t=2,3,\cdots,T$

$$\delta_t(i) = \max_{1 \leq j \leq N} [\delta_{t-1}(j)a_{ji}]b_i(o_t), \ i=1,2,\cdots,N$$$$\Psi_t(i) = \arg \max_{1 \leq j \leq N} [\delta_{t-1}(j)a_{ji}], \ i=1,2,\cdots,N$$

  1. 终止

    $$P^* = \max_{1 \leq j \leq N}\delta_T(i)$$ $$i_T^* = \arg \max_{1 \leq j \leq N} [\delta_T(i)]$$

  2. 最优路径回溯

    对 $t = T-1, T-2, \cdots, 1$

    $$i_t^* = \Psi_{t+1}(i_{t+1}^*)$$

求得最优路径 $I^* = (i_1^*,i_2^*, \cdots, i_T^*)$。

实例

  1. 初始化

    $$\delta_1(i) = \pi_ib_i(o_1) = \pi_ib_i(红) , \ i=1,2,3$$

    $$ \begin{aligned} \delta_1(1) &= 0.5*0.2 = 0.1 \\ \delta_1(2) &= 0.4*0.4 = 0.16 \\ \delta_1(3) &= 0.7*0.4 = 0.28 \end{aligned}$$ 记 $\Psi_1(i) = 0, ~~i=1,2,3$

  2. 在 $t=2$ 时，对每个状态，求在 $t=1$ 时状态为 $j$ 观测为红并在 $t=2$ 时状态为 $i$ 观测 $o_2$ 为白的最大概率，记为 $\delta_2(i)$

    $$\delta_2(i) = \max_{1 \leq j \leq 3} [\delta_1(j) a_{ji}]b_i(o_2)$$

    同时，对每个状态 $i, i=1,2,3$，记录概率最大路径的前一个状态 $j$

    $$\Psi_2(i) = \arg \max_{1 \leq j \leq 3} [\delta_1(j)a_{ji}], ~~~ i = 1,2,3$$

    计算：

    $$ \begin{aligned} \delta_2(1) &= \max_{1 \leq j \leq 3} [\delta_1(j)a_{j1}]b_1(o_2) \\ &= \max_j \{0.1*0.5,0.16*0.3,0.28*0.2 \} * 0.5 \\ &= 0.028 \\ \Psi_2(1) &= 3 ~~~~~~因为~0.056 ~最大, 对应状态 ~3 \\ \delta_2(2) &= \max_{1 \leq j \leq 3} [\delta_1(j)a_{j2}]b_2(o_2) \\ &= \max_j \{ 0.1*0.2,0.16*0.5,0.28*0.3\} * 0.6 \\ &= 0.0504 \\ \Psi_2(2) &= 3 \\ \delta_2(3) &= \max_{1 \leq j \leq 3} [\delta_1(j)a_{j3}]b_3(o_2) \\ &= \max_j \{ 0.1*0.3,0.16*0.2,0.28*0.5\} * 0.3 \\ & = 0.042 \\ \Psi_2(3) &= 3 \end{aligned} $$

    同样，在 $t=3$ 时，

    $$ \begin{aligned} \delta_3(1) &= \max_{1 \leq j \leq 3} [\delta_2(j)a_{j1}]b_1(o_3) \\ &= \max_j \{0.028*0.5,0.0504*0.3,0.042*0.2 \} * 0.5 \\ &= 0.00756 \\ \Psi_3(1) &= 2 ~~~~~~因为~0.00756 ~最大, 对应状态 ~2 \\ \delta_3(2) &= \max_{1 \leq j \leq 3} [\delta_2(j)a_{j2}]b_2(o_3) \\ &= \max_j \{ 0.028*0.2,0.0504*0.5,0.042*0.3\} * 0.4 \\ &= 0.01008 \\ \Psi_3(2) &= 2 \\ \delta_3(3) &= \max_{1 \leq j \leq 3} [\delta_2(j)a_{j3}]b_3(o_3) \\ &= \max_j \{ 0.028*0.3,0.0504*0.2,0.042*0.5\} * 0.7 \\ & = 0.0147 \\ \Psi_3(3) &= 3 \end{aligned} $$

  3. 以 $P^*$ 表示最优路径的概率，则

    $$P^* = \max_{1 \leq i \leq 3} \delta_3(i) = 0.0147$$

    最优路径的终点是 $i_3^*：$

    $$i_3^* = \arg \max_i [\delta_3(i)] = 3$$

  4. 由最优路径的终点 $i_3^*$，逆向找到 $i_2^*, i_1^*$：

    $$ \begin{aligned} t=2, ~~ i_2^* = \Psi_3(3) = 3 \\ t=1, ~~ i_1^* = \Psi_2(3) = 3 \end{aligned}$$

    于是求得最优路径，也就是 $I^* = (3,3,3)$。

Python 代码实现

  • hmmlearn: https://hmmlearn.readthedocs.io/en/latest/

    • hmmlearn安装："pip install hmmlearn"
    • hmmlearn简介：hmmlearn实现了三种HMM模型类，按照观测状态是连续状态还是离散状态，可以分为两类。GaussianHMM和GMMHMM是连续观测状态的HMM模型，而MultinomialHMM是离散观测状态的模型，也是我们在HMM原理系列篇里面使用的模型。

  • HMM 前向后向算法理解和Python实现：https://www.cnblogs.com/gongyanzh/p/12880387.html

HMM 算法相关面经和其他资料

面经

  1. 介绍一下 HMM：

    HMM 是由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再有各个状态生成一个观测从而产生观测随机序列的过程。

  2. HMM 的三要素是什么？

    $A$：转移概率矩阵，$A_{ji}$ 描述的是从当前状态 $j$ 转移到状态 $i$ 的概率是多少。 $B$：观测概率矩阵，指的是最终从每个状态下观测的结果的分布。 $\pi$：初始状态概率向量

  3. 两个基本假设是什么？

    1. 假设隐藏的马尔科夫链在任意时刻 $t$ 的状态只依赖于前一时刻的状态。

    2. 假设任意时刻的观测只依赖于该时刻的马尔科夫链的状态。

      这样的假设简化了模型复杂度，但也同时忽略了一些重要信息。

  4. 前向算法和后向算法的区别

    前向算法先计算给定马尔科夫模型，定义到时刻 $t$ 部分观测序列且在某一状态的概率，接着通过 transition matrix 和观测概率矩阵的计算，计算出观测序列概率，也就是 $P(O|\lambda)$

    后向算法计算给定模型和在时刻 $t$ 的状态下，从 $t+1$ 到 $T$ 的部分观测序列的概率，然后往回递推最终得到观测序列概率。

  5. 解释一下维特比算法

    维特比算法用动态规划来解HMM预测问题，也就是求概率最大路径，这时对应的一个状态序列。通过引入新变量 $\delta$ 和 $\Psi$，分别表示在时刻 $t$ 状态为 $i$ 的所有单个路径中概率最大值，和最大路径的 $t-1$ 个节点。然后迭代找到最大概率路径并反向验证。

其他资料

小结

  1. 隐马尔科夫模型是关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态的序列，再由各个状态随机生成一个观测从而产生观测序列的过程。

    • HMM 可写成 $\lambda = (A,B,\pi)$，其参数分别是状态转移概率矩阵、观测概率矩阵和初始状态概率向量
    • HMM 是一个生成模型
    • HMM 可以用于标注，也就是给定观测序列预测其对应的标记序列

  2. 前后向算法通过递推来计算前后向概率，从而高效完成 HMM 的概率计算

  3. 用极大似然估计来估计参数。Baum-Welch 算法也就是 EM 算法可以高效地对 HMM 进行训练，其是一种无监督方法

  4. 维特比算法应用动态规划来高效求解最佳路径，也就是概率最大的状态序列