统计学习方法第十章隐马尔科夫模型

隐马尔科夫模型（hidden Markov model，HMM）是可用于标注问题的统计学模型，描述由隐藏的马尔科夫链随机生成观测序列的过程，属于生成模型。隐马尔科夫模型在语音识别、自然语言处理、生物信息、模式识别等领域有着广泛的应用。

10.1 隐马尔科夫模型的基本概念

10.1.1 隐马尔科夫模型的定义

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

隐马尔科夫模型由初始概率分布、状态转移概率分布以及观测概率分布确定。隐马尔科夫模型的形式定义如下：

设 $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=\left[a_{ij}\right]_{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=\left[b_j(k)\right]_{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. 齐次马尔科夫性假设

$$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$$

b. 观测独立性假设

$$P(o_t|i_T,o_T,i_{T-1},o_{T-1},\cdots,i_{t+1},o_{t+1},i_t,i_{t-1},o_{t-1},\cdots,i_1,o_1)=P(o_t|i_t)$$

隐马尔科夫模型可以用于标注问题，这时状态对应着标记。标注问题是给定观测的序列预测其对应的标记序列。可以假设标注问题的数据是由隐马尔科夫模型生成的。这样我们可以利用隐马尔科夫模型的学习与预测算法进行标注。

10.1.2 观测序列的生成过程

算法 10.1（观测序列的生成)
输入：隐马尔科夫模型 $\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_ti_{t+1}}\}$ 产生状态 $i_{t+1}$，$i_{t+1}=1,2,\cdots,,N$
(5) 令 $t=t+1$；如果 $t<T$，转 (3)；否则，终止

10.1.3 隐马尔科夫模型的 3 个基本问题

$\blacksquare$ 概率计算问题

给定模型 $\lambda=(A,B,\pi)$ 和观测序列 $O=(o_1,o_2,\cdots,o_T)$，计算在模型 $\lambda$ 下观测序列 $O$ 出现的概率 $P(O|\lambda)$ 。

$\blacksquare$ 学习问题

已知观测序列 $O=(o_1,o_2,\cdots,o_T)$，估计模型参数 $\lambda=(A,B,\pi)$，使得在该模型下观测序列 $P(O|\lambda)$ 最大。即用极大似然估计的方法估计参数。

$\blacksquare$ 预测问题

也称为解码（decoding）问题。已知模型 $\lambda=(A,B,\pi)$ 和观测序列 $O=(o_1,o_2,\cdots,o_T)$，求对给定观测序列条件概率 $P(I|O)$ 最大的状态序列 $I=(i_1,i_2,\cdots,i_T)$。即给定观测序列，求最可能的对应的状态序列。

10.2 概率计算算法

10.2.1 前向算法

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

$$\alpha_i(t)=P(o_1,o_2,\cdots,o_t,t_i=q_i|\lambda)$$

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

算法 10.2（观测序列概率的前向算法）
输入：隐马尔科夫模型 $\lambda$，观测序列 $O$
输出：观测序列概率 $P(O|\lambda)$
(1) 初值

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

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

利用前向概率计算 $P(O|\lambda)$ 的计算量是 $O(N^2T)$ 阶的，而不是直接计算的 $O(TN^T)$ 阶。

10.2.2 后向算法

定义 10.3（后向概率）给定隐马尔科夫模型 $\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)$$

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

算法 10.3（观测序列概率的后向算法）
输入：隐马尔科夫模型 $\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}^Na_{ij}b_j(o_{t+1})\beta_{i+1}(j),\ i=1,2,\cdots,N$$
(3)
$$P(O|\lambda)=\sum_{i=1}^N\pi_ib_i(o_1)\beta_1(i)$$

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

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

当 $t=1$ 或 $t=T-1$ 时，分别为后向概率和前向概率。

根据定义，我们可得如下等式

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

10.2.3 一些概率与期望值的计算

$\blacksquare$ 给定模型 $\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)}=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N\alpha_t(j)\beta_t(j)}$$

$\blacksquare$ 给定模型 $\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{eqnarray} \xi_t(i,j) &=& P(i_t=q_i,i_{t+1}=q_j|O,\lambda)=\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}^NP(i_t=q_i,i_{t+1}=q_j,O|\lambda)} \\ &=& \frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)} \end{eqnarray}$$

$\blacksquare$ 将 $\gamma_t(i)$ 和 $\xi_t(i,j)$ 对各个时刻 $t$ 求和，可以得到一些有用的期望值

a. 在观测 $O$ 下状态 $i$ 出现的期望值

$$\sum_{t=1}^T\gamma_t(i)$$

b. 在观测 $O$ 下由状态 $i$ 转移的期望值

$$\sum_{t=1}^{T-1}\gamma_t(i)$$

c. 在观察 $O$ 下由状态 $i$ 转移到状态 $j$ 的期望值

$$\sum_{t=1}^{T-1}\xi_t(i,j)$$

10.3 学习算法

隐马尔科夫模型的学习，根据训练数据是否包括观测序列和对应的状态序列还是只有观测序列，可以分别由监督学习与非监督学习实现。

10.3.1 监督学习方法

假设已给训练数据包含 $S$ 个长度相同的观测序列和对应的状态序列

$$\{(O_1,I_1),(O_2,I_2),\cdots,(O_S,I_S)\}$$

那么可以用极大似然估计法来估计隐马尔科夫模型的参数。

$\blacksquare$ 转移概率 $a_{ij}$ 的估计

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

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

$\blacksquare$ 观测概率 $b_j(k)$ 的估计

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

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

$\blacksquare$ 初始状态概率 $\pi_i$ 的估计 $\hat{\pi}_i$ 为 $S$ 个样本中初始状态为 $q_i$ 的频率

由于监督学习需要使用训练数据，而人工标注训练数据往往代价很高，有时就会利用非监督学习的方法。

10.3.2 Baum-Welch 算法

假设给定训练数据只包含 $S$ 个长度为 $T$ 的观测序列 $\{O_1,O_2,\cdots,O_S\}$ 而没有对应的状态序列，目标是学习隐马尔科夫模型的参数。

我们将观测序列数据看做观测数据 $O$，状态序列数据看做不可观测的隐数据 $I$，那么隐马尔科夫模型实际上是一个含有隐变量的概率模型

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

它的参数学习可以由 EM 算法实现。

$\blacksquare$ 确定完全数据的对数似然函数

所有观测数据写成 $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)$

$\blacksquare$ EM 算法的 E 步，求 $Q$ 函数 $Q(\lambda,\bar{\lambda})$

$$Q(\lambda,\bar{\lambda})=\sum_I\log P(O,I|\lambda)P(O,I|\bar{\lambda})$$

其中，$\bar{\lambda}$ 是隐马尔科夫模型的当前估计值，$\lambda$ 是要极大化的隐马尔科夫模型参数，$Q$ 函数其中略去了对 $\lambda$ 而言的常数因子 $1/P(O|\bar{\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)$$

于是

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

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

$\blacksquare$ EM 算法的 M 步：极大化 $Q(\lambda,\bar{\lambda})$ 求模型参数 $\lambda=(A,B,\pi)$

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

a. 第一项可写成

$$\sum_I\log\pi_{i_1}P(O,I|\bar{\lambda})=\sum_{i=1}^N\log\pi_iP(O,i_1=q_i|\bar{\lambda})$$

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

$$\sum_{i=1}^N\log\pi_iP(O,i_1=q_i|\bar{\lambda})+\gamma\left(\sum_{i=1}^N\pi_i-1\right)$$

求偏导并令结果为 0，得

$$P(O,i_1=q_i|\bar{\lambda})+\gamma\pi_i=0$$

对 $i$ 求和得到 $\gamma$

$$\gamma=-P(O|\bar{\lambda})$$

带入得

$$\pi_i=\frac{P(O,i_1=q_i|\bar{\lambda})}{P(O|\bar{\lambda})}$$

b. 第二项可以写成

$$\sum_I\left(\sum_{t=1}^{T-1}\log a_{i_ti_{t+1}}\right)P(O,I|\bar{\lambda}) \\ =\sum_{i=1}^N\sum_{j=1}^N\sum_{t=1}^{T-1}\log a_{ij}P(O,i_t=q_i,i_{t+1}=q_j|\bar{\lambda})$$

注意约束条件 $\sum_{j=1}^Na_{ij}=1$，利用拉格朗日乘子法求出

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

c. 第三项可以写成

$$\sum_I\left(\sum_{t=1}^T\log b_{i_t}(o_t)\right)P(O,I|\bar{\lambda}) \\ =\sum_{j=1}^N\sum_{t=1}^T\log b_j(o_t)P(O,i_t=q_j|\bar{\lambda})$$

注意约束条件 $\sum_{k=1}^Mb_{j}(k)=1$，注意，只有在 $o_t=v_k$ 时，$b_j(o_t)$ 对 $b_j(k)$ 的偏导才不是 0，利用拉格朗日乘子法求出

$$b_j(k)=\frac{\sum_{t=1}^TP(O,i_t=q_j|\bar{\lambda})\mathbb{I}(o_t=v_k)}{\sum_{t=1}^TP(O,i_t=q_j|\bar{\lambda})}$$

10.3.3 Baum-Welch 模型参数估计公式

结合 10.2 节，得到

$$ 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) $$

算法 10.4（Baum-Welch 算法）
输入：观测数据 $O=(o_1,o_2,\cdots,o_T)$
输出：隐马尔科夫模型参数
(1) 初始化
对 $n=0$，选取 $a_{ij}^{(0)}$，$b_j(k)^{(0)}$，$\pi_i^{(0)}$，得到模型 $\lambda^{(0)}=\left(A^{(0)},B^{(0)},\pi^{(0)}\right)$
(2) 递推。对 $n=1,2,\cdots,$

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

右端各式按观测 $O=(o_1,o_2,\cdots,o_T)$ 和模型 $\lambda^{(n)}=\left(A^{(n)},B^{(n)},\pi^{(n)}\right)$ 计算，式中
$$\gamma_t(i)=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N\alpha_t(j)\beta_t(j)}\\ \xi_t(i,j)=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}$$
(3) 终止。得到模型参数 $\lambda^{(n+1)}=\left(A^{(n+1)},B^{(n+1)},\pi^{(n+1)}\right)$

10.4 预测算法

10.4.1 近似算法

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

给定隐马尔科夫模型和观测序列，在时刻 $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^\star$ 是

$$i_t^\star=\arg\max_{1\leq i\leq N}[\gamma_t(i)],\ i=1,2,\cdots,T$$

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

近似算法的优点是计算简单，其缺点是不能保证预测的状态序列整体是最有可能的状态序列，因为预测的状态序列可能有实际不发生的部分。事实上，上述方法得到的状态序列中有可能存在转移概率为 0 的相邻状态。尽管如此，近似算法仍然是有用的。

10.4.2 维特比算法

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

根据动态规划原理，最优路径具有这样的特征：如果最优路径在时刻 $t$ 通过节点 $i_t^\star$，那么这一路径中“从节点 $i_1^\star$ 到终点 $i_T^\star$ 的部分路径”对于“从节点 $i_1^\star$ 到终点 $i_T^\star$ 的所有可能的部分路径”来说必须是最优的。

根据这一原理，我们只需从时刻 $t=1$ 开始，递推的计算在时刻 $t$ 状态为 $i$ 的各条部分路径的最大概率，直至得到时刻 $t=T$ 时状态为 $i$ 的各条路径的最大概率。时刻 $t=T$ 的最大概率即为最优路径的概率 $P^\star$，最优路径的终点 $i_T^\star$ 也就同时得到了。

之后，为了找出最优路径的各个节点，从终结点 $i_T^\star$ 开始，由后向前逐步求得节点 $i_{T-1}^\star,\cdots,i_1^\star$，得到最优路径 $I^\star=(i_1^\star,i_2^\star,\cdots,i_T^\star)$。这就是维特比算法。

首先定义在时刻 $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$$

得递推公式

$$\begin{eqnarray} \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_t(ja_{ji})]b_i(o_{t+1}),\ i=1,2,\cdots,N;\ t=1,2,\cdots,T-1 \end{eqnarray}$$

定义在时刻 $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$$

算法 10.5（维特比算法）
输入：模型 $\lambda=(A,B,\pi)$ 和观测 $O=(o_1,o_2,\cdots,o_T)$
输出：最优路径 $I^\star=(i_1^\star,i_2^\star,\cdots,i_T^\star)$
(1) 初始化

$$\delta_1(i)=\pi_ib_i(o_1),\ i=1,2,\cdots,N \\ \psi_1(i)=0,\ i=1,2,\cdots,N$$

(2) 递推，对 $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$$
(3) 终止
$$P^\star=\max_{1\leq i\leq N}\delta_T(i) \\ i_T^\star=\arg\max_{1\leq i\leq N}[\delta_T(i)]$$
(4) 最优路径回溯，对 $t=T-1,T-2,\cdots,1$
$$i_t^\star=\psi_{t+1}(i_{t+1}^\star)$$
求得最优路径 $I^\star=(i_1^\star,i_2^\star,\cdots,i_T^\star)$