机器学习（10）：隐马尔可夫模型（HMM）

最近自己会把自己个人博客中的文章陆陆续续的复制到CSDN上来，欢迎大家关注我的 个人博客，以及我的github。

本文将讲解有关隐马尔可夫模型的相关知识，主要包括隐马尔可夫模型的概念、两个基本假设、三个基本问题，以及解决这三个基本问题的前向-后向算法、Baum-Welch 算法和 Viterbi 算法等内容。

本文主要是依据李航老师的《统计学习方法》和邹博老师的机器学习教学视频总结编写的。文中所用到的有关机器学习的基本概念和方法可以参考本人博客中该系列之前的文章，或者直接上网搜索相关的内容。以下文章所列出的公式以及公式的推导读者们最好是在草稿本上自己推导一遍。由于本人水平所限，文章中难免有错误和不当之处，欢迎大家多多批评指正！

一、隐马尔可夫模型

如上图所示，隐马尔可夫模型（HMM）一个关于时序的概率模型，红色圆圈是隐藏的不可观测的序列，被称为状态序列，其中 $i_t$ 的值只取决于它的前一个状态 $i_{t-1}$；绿色的矩形是可观测的序列，被称为观测序列。序列的每一个位置可以看作是一个时刻。

举个例子，在天气预报中，每天的天气情况就是 HMM 中的观测序列，而每一天对应的气压情况就是 HMM 中的状态序列，我们假设它是不可观测的，并且假设天气与气压有关。

---

为了后面描述的方便，下面先来引入几个定义。

设 Q 是所有可能的状态（红色圆圈）的结合，V 是所有可能的观测（绿色矩形）的集合：
$$Q=(q_1,q_2,...,q_N)$$

$$V=(v_1,v_2,...,v_M)$$

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

设 I 是长度为 T 的状态序列，O 是对应的观测序列：
$$I=(i_1,i_2,...,i_T)$$

$$O=(o_1,o_2,...,o_T)$$

HMM可用符号 $\lambda$ 表示，该模型可由初始概率分布 $\pi$、状态转移概率分布 A 和观测概率分布 B 确定，即 $\lambda =(A,B,\pi )$。

**(1) 初始概率分布 $\pi$ **就是图中第一个红色圆圈 $i_1$ 取每个值的概率：
$$\pi =({\pi }_1,{\pi }_2,...,{\pi }_N)$$

$${\pi }_i=P(i_1=q_i),i=1,2,...N$$

其中 $q_i$ 就是上面提到过的可能的状态。

**(2) 状态转移概率分布 A **就是在时刻 t 的红色圆圈 $i_t$ 的状态为 $q_i$ 的条件下，下一个时刻的红色圆圈 $i_{t+1}$ 的状态变为 $q_j$ 的概率：
$$A=[a_{ij}{]}_{N\times N}$$

$$a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$$

**(3) 观测概率分布 B **是指在时刻 t 红色圆圈 $i_t$ 的状态是 $q_i$ 的条件下，该时刻的观测值为 $v_k$ 的概率：
$$B=[b_j(k){]}_{N\times M}$$

$$b_j(k)=P(o_t=v_k|i_t=q_j)$$

沿用上面天气预报的例子，在该例子中，每天的气压情况只与前一天的气压情况有关，则如果我们得知了第一天的气压情况（初始概率分布 $\pi$），又得知了每相邻两天气压的变化规律（状态转移概率分布 A）以及每天气压和天气之间的情况（观测概率分布 B），那么我们就可以预测第 T 天的天气情况了。

二、HMM 的两个基本假设和三个基本问题

1. 两个基本假设

(1) 齐次马尔科夫性假设：

假设隐藏的马尔科夫链（即图中红色圆圈形成的链条）在任意时刻 t 的状态只依赖于其前一个时刻。

(2) 观测独立假设：

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

2. 三个基本问题

(1) 概率计算问题：

给定模型 $\lambda =(A,B,\pi )$ 和观测序列 $O=(o_1,o_2,...,o_T)$ 计算在模型 $\lambda$ 下观测序列 O 出现的概率 $P(O|\lambda )$。也就是计算模型和观测序列之间的匹配程度。

(2) 学习问题：

已知观测序列 $O=(o_1,o_2,...,o_T)$，估计模型 $\lambda =(A,B,\pi )$ 的参数，使得在该模型下观测序列出现的概率 $P(O|\lambda )$ 最大。也就是如何构造一个模型去最好的拟合观测数据。

(3) 预测问题：

又称解码问题。已知模型 $\lambda =(A,B,\pi )$ 和观测序列 $O=(o_1,o_2,...,o_T)$，求对给定观测序列条件下，概率 $p(I|O)$ 最大的状态序列 $I=(i_1,i_2,...,i_T)$。也就是给定观测序列，求最有可能的对应的状态序列。

针对以上三个基本问题，对每个问题的求解都提出了一种比较高效的解决方法，下面就是来分别介绍这些算法。

三、概率计算问题

再放一遍最开始的图，以便和公式对照着看。

1. 直接计算法

要解决的问题是计算模型和观测序列之间的匹配程度 $P(O|\lambda )$。可以先求出在给定模型 $\lambda$ 的前提下，生成状态序列 I 的概率 $P(I|\lambda )$，然后再求出模型和状态序列已知的情况下，生成观测序列 O 的概率 $P(O|I,\lambda )$。这样就可以利用公式 $P(O,I|\lambda )=P(O|I,\lambda )\cdot P(I|\lambda )$ 来计算了。其实也就是算出了 $P(O|\lambda )$，只不过还多求出了状态序列。

但是由于对于一个已知的观测序列 $O=(o_1,o_2,...,o_T)$ 来说，它可能是由多个不同的状态序列来生成的。先要想明白这句话是为什么。因为——每个状态 $i_t$ 都以一定的概率生成对应的观测 $o_t$，而每个状态 $i_t$ 又以一定的概率生成下一个状态 $o_{t+1}$。

所以可以通过枚举所有可能的长度为 T 的状态序列 $I=(i_1,i_2,...,i_T)$，然后对概率求和就可以得到最终的结果了。其公式如下：
$$P(O|\lambda )=\sum _{I}P(O|I,\lambda )\cdot P(I|\lambda )$$
但是这种暴力枚举的方式时间复杂度太高了，所以人们又提出了“前向-后向算法”来解决这个问题。

2. 前向算法

前向-后向算法利用了类似于动态规划的思路来求解 $P(O|\lambda )$，前向算法从前往后计算概率，后向算法从后往前计算概率。

首先定义前向概率：给定隐马尔可夫模型 $\lambda$，定义到时刻 t 为止的部分观测序列为 $O=(o_1,o_2,...,o_t)$，并且第 t 个时刻的状态为 $q_i$ 的概率为前向概率。记作：
$${\alpha }_t(i)=P(o_1,o_2,...,o_T,i_t=q_i|\lambda )$$
用前向概率计算 $P(O|\lambda )$ 过程如下：

(1) 初值
$${\alpha }_1(i)={\pi }_i{b}_i(o_1),i=1,2,...,N$$
${\alpha }_1(i)$ 表示第 1 个状态为 $q_i$，到此为止观测序列为 $o_1$ 的概率。就需要先以概率 ${\pi }_i$ 生成状态 $q_i$。再由该状态以概率 $b_{i,o_1}$ 生成观测 $o_1$ 。

(2) 递推

对 $t=1,2,...,T-1$
$${\alpha }_{t+1}(i)=[\sum _{j=1}^N{\alpha }_t(j)a_{ji}]b_i(o_{t+1}),i=1,2,...,N$$
在已知 ${\alpha }_t(j)$，即已知到前 t 个时刻为止的观测序列为 $o_1,o_2,...,o_t$，且第 t 时刻的状态为 $q_j$ 时，求 ${\alpha }_{t+1}(i)$ 的概率，即求到前 t+1 个时刻为止的观测序列为 $o_1,o_2,...,o_t,o_{t+1}$，且第 t+1 时刻的状态为 $q_i$ 的概率。需要先以上的基础上，以概率 $a_{ji}$ 由第 j 个状态 $q_j$ 生成下一个状态 $q_i$，然后再由以概率 $b_{i,o_{t+1}}$ 由下一个状态 $q_i$ 生成对应的观测 $o_{t+1}$。由于 t 时刻的状态 $q_j$ 可能有多种情况，所以还需要对 j 进行求和。

(3) 终止
$$P(O|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$$
${\alpha }_T(i)$ 表示观测序列为 $O=(o_1,o_2,...,o_T)$，并且第 T 个时刻的状态为 $q_i$，对所有可能的状态求和就是最终要求的 $P(O|\lambda )$。

3. 后向算法

后向概率：给定隐马尔可夫模型 $\lambda$，定义在时刻 t 状态为 $q_i$ 的条件下，从 t+1 到 T 时刻的部分观测序列为 $o_{t+1},o_{t+1},...,o_T$ 的概率为后向概率，记作：
$${\beta }_t(i)=P(o_{t+1},o_{t+1},...,o_T|i_t=q_i,\lambda )$$
可以用与前向算法类似的方式递推得到 $P(O|\lambda )$：

(1) 初值
$${\beta }_T(i)=1,i=1,2,...,N$$
由于 T 是最终的时刻，不存在从 T+1 到 T 的观测序列，所以概率为 1。

(2) 递推

对于 $t=T-1,T-2,...,1$
$${\beta }_t(i)=\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j),i=1,2,...,N$$
先从时刻 t+1 的状态 $p_j$ 以概率 $b_j(o_{t+1})$ 的概率生成对应的观测 $o_{t+1}$。而状态 $p_j$ 是由时刻 t 的状态 $p_i$ 以概率 $a_{ij}$ 生成的。

(3) 终止
$$P(O|\lambda )=\sum _{i=1}^N{\pi }_i{b}_i(o_1){\beta }_1(i)$$
同理，倒着递推出最终的表达式。

4. 一些概率和期望值

下面给出某些概率和期望值的计算公式，在后面的算法中会用到其中的几个。

四、学习问题

学习问题是如何构造一个模型去最好的拟合观测数据。

如果训练数据包括观测序列和状态序列，则 HMM 的学习是监督学习；若只有观测序列，则 HMM的学习需要用 EM 算法解决，是非监督学习。

1. 监督学习方法

因为训练数据包括观测序列和状态序列，所以要求的模型的参数 $a_{ij},b_i(j)$ 和 ${\pi }_i$ 都可以通过统计训练数据中它们出现的频率来近似表示。

2. Baum-Welch 算法

（1）EM 算法的 E 步：

先写出 Q 函数：
$$Q(\lambda ,\bar{\lambda })=\sum _I\log P(O,I|\lambda )P(O,I|\bar{\lambda })$$
其中 $\bar{\lambda }$ 是隐马尔可夫模型参数的当前估计值，$\lambda$ 是要极大化的隐马尔可夫模型参数。

因为 $P(O,I|\lambda )={\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1}{b}_{i_2}(o_2)...a_{i_{T-1}{i}_T}{b}_{i_T}(o_T)$，所以 $Q(\lambda ,\bar{\lambda })$ 可以写为：
$$Q(\lambda ,\bar{\lambda })=\sum _I\log {\pi }_{i_1}P(O,I|\bar{\lambda })+\sum _I(\sum _{t=1}^{T-1}\log a_{i_t,i_{t+1}})P(O,I|\bar{\lambda })$$

$$+\sum _I(\sum _{t=1}^T\log b_{i_t}(o_t))P(O,I|\bar{\lambda })公式(1)$$

（2）EM 算法的 M 步：

这一步要极大化 Q 函数求模型的参数 $A,B,\pi$ 。

公式 (1) 中的第 1 项可以写为：
$$\sum _I\log {\pi }_{i_1}P(O,I|\bar{\lambda })=\sum _{i=1}^N\log {\pi }_{i}P(O,i_1=i|\bar{\lambda })$$
因为上式满足约束条件 ${\sum }_{i=1}^N{\pi }_i=1$，利用拉格朗日乘子法得拉氏函数为：
$$L=\sum _{i=1}^N\log {\pi }_{i}P(O,i_1=i|\bar{\lambda })+\gamma (\sum _{i=1}^N{\pi }_i-1)$$
拉氏函数对 ${\pi }_i$ 求偏导，并令结果为 0 得：
$$P(O,i_1=i|\bar{\lambda })+\gamma {\pi }_i=0公式(2)$$
上式在求导后的原式的基础上两边同乘了 ${\pi }_i$。上式对 i 求和得：
$$\gamma =-P(O|\bar{\lambda })$$
将其带入公式 (2) 得：
$${\pi }_i=\frac{P(O,i_1=i|\bar{\lambda })}{P(O|\bar{\lambda })}$$

---

同理公式 (1) 的第 2 项的约束条件为 ${\sum }_{j=1}^N{a}_{ij}=1$，第 3 项的约束条件为 ${\sum }_{k=1}^M{b}_j(k)=1$，可以分别得：
$$a_{ij}=\frac{{\sum }_{t=1}^{T-1}P(O,i_t=i,i_{t+1=j}|\bar{\lambda })}{{\sum }_{t=1}^{T-1}P(O,i_t=i|\bar{\lambda })}$$

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

而以上三个式子利用第三部分第4小节中的公式 (10.24) 和公式 (10.26) 就可以计算。

五、预测问题

预测问题是给定观测序列，求最有可能的对应的状态序列。解决预测问题有两个方法，一个是近似算法，另一个是维特比算法。

1. 近似算法

近似算法的思想是分别求出每个时刻最有可能的状态，把它们形成的状态序列作为预测结果。但是不能保证这种由局部最优解组合得到的状态序列就是全局最优解。

给定隐马尔可夫模型 $\lambda$ 和观测序列 O，在时刻 t 出于状态 $q_i$ 的概率 ${\gamma }_t(i)$ 可以由公式 (10.24) 计算得到。然后对每个时刻的状态求最大值，就得到了当前时刻最有可能的状态。

2. 维特比算法

维特比算法的本质是利用动态规划算法求解概率最优路径 $I^{\ast }=(i_1^{\ast },i_2^{\ast },...,i_T^{\ast })$，一条路径对应着一个状态序列。

最优路径有这样的特性：如果最优路径在时刻 t 的状态为 $i_t^{\ast }$，那么从 $i_t^{\ast }$ 到终点 $i_T^{\ast }$ （时刻 T 的状态）的这部分序列，对于任意的从时刻 t 到时刻 T 的状态序列来说都是最优的。不然就可以将其替换得到更优的序列。

所以可以先求得最优路径的终点 $i_T^{\ast }$，然后从后往前递推的得到状态节点 $i_{T-1}^{\ast },i_{T-2}^{\ast },...,i_1^{\ast }$，从而得到最优的路径 $I^{\ast }=(i_1^{\ast },i_2^{\ast },...,i_T^{\ast })$。

---

首先引入两个变量 $\delta$ 和 $\psi$。定义在时刻 t 状态为 i 的所有单个路径 $(i_1,i_2，...，i_t)$ 中概率最大值为：
$${\delta }_t(i)=\underset{i_1,i_2，...，i_{t-1}}{\max }P(i_t=i,i_{t-1},...,i_1,o_t,...,o_1|\lambda )$$
其递推式为：
$${\delta }_{t+1}(i)=\underset{i_1,i_2，...，i_t}{\max }P(i_{t+1}=i,i_t,...,i_1,o_{t+1},...,o_1|\lambda )$$

$$=\underset{1\le j\le N}{\max }[{\delta }_t(j)a_{ji}]b_i(o_{t+1})$$

定义在时刻 t 状态为 i 的所有单个路径 $(i_1,i_2，...，i_{t-1},i)$ 中概率最大的路径的第 t-1 个节点为：
$${\psi }_t(i)=\arg \underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]$$

---

维特比算法

(1) 初始化
$${\delta }_1(i)={\pi }_i{b}_i(o_1)$$

$${\psi }_1(i)=0$$

(2) 递推

对于 $t=2,3,...,T$，有
$${\delta }_t(i)==\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]b_i(o_t)$$

$${\psi }_t(i)=\arg \underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]$$

(3) 终止
$$P^{\ast }=\underset{1\le i\le N}{\max }{\delta }_T(i)$$

$$i_T^{\ast }=\arg \underset{1\le i\le N}{\max }[{\delta }_T(i)]$$

(4) 最优路径回溯

对于 $t=T-1,T-2,...,1$，有
$$i_t^{\ast }={\psi }_{t+1}(i_{t+1}^{\ast })$$
至此，得到最优路径（序列）$I^{\ast }=(i_1^{\ast },i_2^{\ast },...,i_T^{\ast })$。