隐马尔科夫模型

一、序言

重新复习隐马尔科夫模型，重点是HMM模型的三个问题及前向、后向和维特比算法。

二、基本概念

2.1 定义

隐马尔可夫模型由初始概率分布、状态转移概率分布以及观测概率分布确定。隐马尔可夫模型的形式定义如下:
设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$$
$A$是状态转移概率矩阵:
$$A=[a_{ij}{]}_{n\times n}$$
其中,
$$a_{ij}=P(i_t=q_j|i_{t-1}=q_i),i=1,...n;j=1,...,n$$
是在时刻$t-1$处于状态$q_i$的条件下在时刻$t$转移到状态$q_j$的概率。
$B$是观测概率矩阵:
$$B=[b_{jk}{]}_{n\times m}$$
其中，
$$b_{kj}=P(o_t=v_k|i_t=q_j),k=1,...m,j=1,...,n$$
是在时刻$t$处于状态$q_j$的条件下生成观测$v_k$的概率。
$\pi$是初始状态概率向量:
$$\pi ={\pi }_i$$
其中，
$${\pi }_i=P(i_1=q_i),i=1,...n$$
是初始时刻$t=1$处于状态$q_i$的概率。

2.2 例子

举个例子，当观察到屋外艳阳高照，那么肯定是晴天；若是半乌云密布，则是阴天；若是电闪雷鸣，则是雨天。艳阳高照，乌云密布，电闪雷鸣是我们能直接观察到的，对应着上面定义的观测序列。
而它们对应的天气状态分别是晴天、阴天和雨天，则是状态序列，因为我们先观察到外边的环境是艳阳高照，乌云密布，电闪雷鸣，然后再推测出是晴天、阴天还是雨天。
如下图所示，上面的是一条隐马尔科夫链，下面对应着其随机生成的状态序列。

如下图所示，是一个完整的 HMM 模型。

状态集合 $Q=q_1,q_2,q_3$，其中$q_1=艳阳高照$，$q_2=乌云密布$，$q_3=电闪雷鸣$。
观测集合 $V=v_1,v_2,v_3$，其中$v_1=晴天$，$v_2=阴天$，$v_3=雨天$。
状态转移概率矩阵 $A$：

观测概率矩阵 $B$：

初始状态概率 $\pi$：

以上数据是随便写的。

2.3 基本假设

三、三个问题

只看这个可能有点晦涩，下面就例子说的通俗一下：

3.1 概率计算问题

评估问题，即概率计算问题，是三个问题中最简单的。给定 HMM 模型 $\lambda$，也就是已经知道状态转移概率矩阵 $A$、观测概率矩阵 $B$ 和初始状态概率 $\pi$，同时给出观测序列 $O=o_1,o_2,...,o_t$，求在该 HMM 模型下这个观测序列生成的概率。例如求接下来三天的观测天气是（阴天，雨天，晴天）的概率。解决算法：前向-后向算法。

3.2 学习问题

学习问题是三个问题中最复杂的一个。这个问题中只给出观测序列 $O=o_1,o_2,...,o_t$，让求 HMM 模型 $\lambda$ 的三个参数： $A$、 $B$ 和 $\pi$。例如，给出观测天气是（阴天，雨天，晴天），根据观测序列求一个 HMM 模型。解决算法：$Baum-Welch$ 算法（EM算法）。

3.3 预测问题

预测问题，也称为解码问题。给定 HMM 模型 $\lambda$ 和观测序列 $O=o_1,o_2,...,o_t$，求在该 HMM 模型下最有可能生成这个观测序列的隐状态序列。例如，观测天气是（阴天，雨天，晴天），求最有可能对应该观测序列的状态序列是（艳阳高照，乌云密布，电闪雷鸣），还是（乌云密布，电闪雷鸣，艳阳高照），或者是其他的某个状态序列。解决算法：$Viterbi$ 算法（一种动态规划）。

四、三个问题解决算法

4.1 概率计算算法

目的：给定 HMM 模型 $\lambda =(A,B,\pi )$ 和观测序列 $O=o_1,o_2,...,o_T$，求在该HMM 模型下生成该观测序列的概率$P(O|\lambda )$。

4.1.1 直接计算法（暴力法）

首先，在给定 HMM 模型下，生成 一个 隐状态序列 $I=i_1,i_2,...,i_T$ 的概率为：
$$P(I|\lambda )={\pi }_{i_1}{a}_{i_1{i}_2}{a}_{i_2{i}_3}...a_{i_{T-1}{i}_T}$$
然后，在该状态序列，生成对应的观测序列 $O=o_1,o_2,...,o_T$ 的改立为：
$$P(O|I,\lambda )=b_{i_1{o}_1}{b}_{i_2{o}_2}...b_{i_T{o}_T}$$
最后，在给定 HMM 模型下，生成状态序列 $I$ 和观测序列 $O$ 的联合概率为：
$$P(O,I|\lambda )=P(O|I,\lambda )P(I|\lambda )={\pi }_{i_1}{b}_{i_1{o}_1}{a}_{i_1{i}_2}{b}_{i_2{o}_2}...a_{i_{T-1}{i}_T}{b}_{i_T{o}_T}$$
综上是 HMM 模型生成一个状态序列，再生成观测序列的概率。只要对所有不同的状态序列 $I$ 求和，就是要求的给定观测序列的概率$P(O|\lambda )$：
$$P(O|\lambda )=\sum _{I}P(O|I,\lambda )P(I|\lambda )=\sum _{I_1,I_2,...}{\pi }_{i_1}{b}_{i_1{o}_1}{a}_{i_1{i}_2}{b}_{i_2{o}_2}...a_{i_{T-1}{i}_T}{b}_{i_T{o}_T}$$
使用该算法原理简单，但是计算量巨大。时间复杂度：$O(T{N}^T)$。

4.1.2 前向算法

4.1.2.1 详解

前向概率 ${\alpha }_t(i)$ 如下图所示：

其实，前向算法可以看做是动态规划。
注意看呦，${\alpha }_t(i)=P(o_1,o_2,...,o_t,i_t=q_t|\lambda )$ 这不就是 暴力法中第三步 求的状态序列 $I$ 和观测序列 $O$ 的联合概率吗？
我们只要把 $t$ 时刻中所有状态序列 $q_i\in (q_1,q_2,...,q_n)$ 做累加，然后乘上$t+1$ 时刻 $q_i$ 对应的观测概率，即 $[{\sum }_{j=1}^n{\alpha }_t(j)a_{ji}]b_i{o}_{t+1}$，就得到 $t+1$ 时刻的状态序列 $I$ 和观测序列 $O$ 的联合概率，即前向概率 ${\alpha }_{t+1}(i)$。
如下图所示：

所以，只要计算出 $t=1$ 时刻的前向概率 ${\alpha }_1(i)$，往后依次递推就可以了。例如 ${\alpha }_1(i)={\pi }_{i_1}{b}_{i_1{o}_1}$，${\alpha }_2(i)={\alpha }_1(1)b_1{o}_i+{\alpha }_1(2)b_2{o}_i+...+{\alpha }_1(n)b_n{o}_i$。

综上：

4.1.2.2 例子

观测集合为：
$$V=(红，白)$$
状态集合为：
$$Q=(盒子1,盒子2,盒子3)$$
观测序列为：
$$O=(红，白，红)$$
状态转移概率矩阵为 $A$ :
第 $i$ 行表示选择第 $i$ 个盒子，第 $j$ 列表示转移到第 $j$ 个盒子， 比如：$A_{23}$ 表示上一次选择第二个盒子，这次选择第三个盒子的概率为 0.2。
观测概率矩阵 $B$：
第 $i$ 行表示选择的是第 $i$ 个盒子，第 $j$ 列表示从该盒子取到 $j$ 号球， 比如：$B_{31}$ 表示从第二个盒子取出球的概率为 0.7。

（1) 计算初值 $t=1$
$t=1$ 时刻取出红球，隐状态是盒子1的概率：
$${\alpha }_1(1)={\pi }_1{b}_{1o_1}=0.2\times 0.5=0.10$$
$t=1$ 时刻取出红球，隐状态是盒子2的概率：
$${\alpha }_1(2)={\pi }_2{b}_{2o_1}=0.4\times 0.4=0.16$$
$t=1$ 时刻取出红球，隐状态是盒子3的概率：
$${\alpha }_1(3)={\pi }_3{b}_{3o_1}=0.4\times 0.7=0.28$$
（2) 递推计算 $t=2$
$t=2$ 时刻取出白球，隐状态是盒子1的概率：
$${\alpha }_2(1)=[\sum _{i=1}^3{\alpha }_1(i)a_{i1}]b_{1o_2}=(0.10\times 0.5+0.16\times 0.3+0.28\ast 0.2)\times 0.5=0.154\times 0.5=0.077$$
$t=2$ 时刻取出白球，隐状态是盒子2的概率：
$${\alpha }_2(2)=[\sum _{i=1}^3{\alpha }_1(i)a_{i2}]b_{2o_2}=(0.10\times 0.2+0.16\times 0.5+0.28\ast 0.3)\times 0.6=0.184\times 0.6=0.1104$$
$t=2$ 时刻取出白球，隐状态是盒子3的概率：
$${\alpha }_2(3)=[\sum _{i=1}^3{\alpha }_1(i)a_{i3}]b_{3o_2}=(0.10\times 0.3+0.16\times 0.2+0.28\ast 0.5)\times 0.3=0.202\times 0.3=0.0606$$
（3) 递推计算 $t=3$
$t=3$ 时刻取出红球，隐状态是盒子1的概率：
$${\alpha }_3(1)=[\sum _{i=1}^3{\alpha }_2(i)a_{i1}]b_{1o_2}=0.04187$$
$t=3$ 时刻取出红球，隐状态是盒子2的概率：
$${\alpha }_3(2)=[\sum _{i=1}^3{\alpha }_2(i)a_{i2}]b_{2o_2}=0.03551$$
$t=3$ 时刻取出红球，隐状态是盒子2的概率：
$${\alpha }_3(3)=[\sum _{i=1}^3{\alpha }_2(i)a_{i3}]b_{3o_2}=0.05284$$
（4) 终止
$$P(O|\lambda )=\sum _{i=1}^3{\alpha }_3(i)=0.13022$$

4.1.3 后向算法

其实后向算法和前向算法类似，只不过是从后往前递推。

后向概率 $\beta t(i)$ 如下图所示：

首先，定义最后时刻的 ${\beta }_T(i)=1$。
然后 ，对于 $t=T-1,T-2,...,1$，后向概率 ${\beta }_t(i)$ 就等于 $t$ 时刻的状态 $i_t=q_i$ 转移到时刻 $t+1$ 的状态 $i_{t+1}=q_j$ 的概率 × $t+1$ 时刻状态 $i_{t+1}$ 对应的观测状态 $o_{t+1}$ 的概率 × $t+1$ 时刻的后向概率 ${\beta }_{t+1}(i)$。即：
$${\beta }_t(i)=\sum _{j=1}^n{a}_{ij}{b}_{j{o}_{t+1}}{\beta }_{t+1}(i)$$
如下图所示：

最后，观测概率 $P(O|\lambda )={\sum }_{i=1}^n{\pi }_i{b}_{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)$$
是不是其实很好理解。

4.1.3 一些概率与期望的计算

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

1. 给定模型 $\lambda$ 和观测 $O$，在时刻 $t$ 处于状态 $q_i$ 的概率。记
   $${\gamma }_t(i)=P(i_t=q_i|O,\lambda )=\frac{P(i_t=q_i,O|\lambda )}{P(O|\lambda )}$$
   由前向概率 ${\alpha }_t(i)$ 和后向概率 ${\beta }_t(i)$ 定义可知:
   $${\alpha }_t(i){\beta }_t(i)=P(i_t=q_t|O,\lambda )$$
   于是得到:
   $${\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)}$$
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 )=\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 )}$$
   而
   $$P(i_t=q_i,i_{t+1}=q_j,O|\lambda )={\alpha }_t(i)a_{ij}{b}_{j{o}_{t+1}}{\beta }_{t+1}(j)$$
   所以
   $${\xi }_t(i,j)=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)}$$
   

4.2 学习算法

目的：

1. 给定观测序列 $O=o_1,o_2,...,o_T$ 和状态序列 $I=i_1,i_2,...,i_T$，求HMM 模型 $\lambda =(A,B,\pi )$ 的三个参数。
2. 给定观测序列 $O=o_1,o_2,...,o_T$，求HMM 模型 $\lambda =(A,B,\pi )$ 的三个参数。

解决方法：

1. 监督算法
2. $Baum-Welch$ 算法

4.2.1 监督算法

第二步求观测概率应该是 $b_{jk}$，因为懒，就直接截图了。

4.2.2 Baum-Welch 算法

现在已经知道的是观测数据 $O=o_1,o_2,...,o_T$， 设隐状态数据为 $I=i_1,i_2,...,i_T$，那么完全数据是 $(O,I)=(o_1,o_2,...,o_T,i_1,i_2,...,i_T)$。完全数据的对数似然函数是 $logP(O,I|\lambda )$。

既然 $Baum-Welch$ 算法使用的就是 $EM$ 算法，那么就要走两个步骤：
（1） $E$ 步
求出联合分布 $P(O,I|\lambda )$ 基于条件概率 $(I|O,\overline{\lambda })$ 的期望，其中 $\overline{\lambda }$ 为 HMM 模型参数的当前估计值，$\lambda$ 为极大化的 HMM 模型参数。
（2） $M$ 步
最大化这个期望，得到更新的模型参数λ。接着不停的进行EM迭代，直到模型参数的值收敛为止。

公式推导：
（1） $E$ 步：求$Q$ 函数
根据 $EM$ 的 $Q$ 函数定义，即这里要求的联合分布的期望为：
$$Q(\lambda ,\overline{\lambda })=\sum _{I}P(I|O,\overline{\lambda })lnP(O,I|\lambda )=\sum _{I}lnP(O,I|\lambda )\frac{P(I,O|\overline{\lambda })}{P(O,\overline{\lambda })}$$
$P(O,\overline{\lambda })$ 表示上次求出的参数与观测数据的联合概率，没有什么影响，所以：
$$Q(\lambda ,\overline{\lambda })=\sum _{I}lnP(O,I|\lambda )P(I,O|\overline{\lambda })$$
而
$$P(O,I|\lambda )={\pi }_{i_1}{b}_{i_1{o}_1}{a}_{i_1{i}_2}{b}_{i_2{o}_2}...a_{i_{T-1}{i}_T}{b}_{i_T{o}_T}$$
所以
$$\begin{gathered} Q(\lambda ,\overline{\lambda })=\sum _{I}P(I,O|\overline{\lambda })[ln{\pi }_{i_1}+ln(a_{i_1{i}_2}+...+a_{i_{T-1}{i}_T}+ln(b_{i_1{o}_1}+...+b_{i_T{o}_T}))] \\ =\underset{式1}{\underbrace{\sum _{I}ln{\pi }_{i_1}P(I,O|\overline{\lambda })}}+\underset{式2}{\underbrace{\sum _I(\sum _{t=1}^{T-1}ln{a}_{i_t{i}_{t+1}})P(I,O|\overline{\lambda })}}+\underset{式3}{\underbrace{\sum _I(\sum _{t=1}^{T}ln{b}_{i_t{o}_t})P(I,O|\overline{\lambda })}} \end{gathered}$$

（2） $M$ 步：极大化 $Q$，求模型参数 $A,B,\pi$
　1）求 ${\pi }_i$：
　既然是求极值，肯定是要求导了。对于 ${\pi }_i$ 来说，满足约束条件 ${\sum }_{t=1}^N{\pi }_i=1$。现在就变成了带约束条件的求极值，直接上拉格朗日乘子法。
　式 1 可以写成：
　$$\sum _{I}ln{\pi }_{i_1}P(I,O|\overline{\lambda })=\sum _{i=1}^{N}ln{\pi }_{i}P(O,i_1=q_i|\overline{\lambda })$$
　拉格朗日函数：
　$$L=\sum _{i=1}^{N}ln{\pi }_{i}P(O,i_1=q_i|\overline{\lambda })+\gamma (\sum _{i=1}^N{\pi }_i-1)$$
　首先把求和 $\sum$ 去掉，只对单个的 ${\pi }_i$ 求偏导并等于 0：
　$$\frac{\partial L}{\partial {\pi }_i}=\frac{P(O,i_1=q_i|\overline{\lambda })}{{\pi }_i}+\gamma =0$$
　等价于：
　$$\frac{\partial L}{\partial {\pi }_i}=P(O,i_1=q_i|\overline{\lambda })+\gamma {\pi }_i=0$$
　然后再添上对 $i$ 的求和 $\sum$，可得到：
　$$\gamma =-P(O|\overline{\lambda })$$
　带入到第三项公式，可得：
　$${\pi }_i=\frac{P(O,i_1=q_i|\overline{\lambda })}{P(O|\overline{\lambda })}$$

2）求 $a_{ij}$：
　式 2 可以写成：
　$$\sum _I(\sum _{t=1}^{T-1}ln{a}_{i_t{i}_{t+1}})P(O,I|\overline{\lambda })=\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}ln{a}_{ij}P(O,i_t=q_i,i_{t+1}=q_j|\overline{\lambda })$$
　同样有约束条件 ${\sum }_{j=1}^N{a}_{ij}=1$，最后可以得到：
　$$a_{ij}=\frac{{\sum }_{t=1}^{T-1}P(O,i_1=q_i,i_{t+1}=q_j|\overline{\lambda })}{P(O,i_t=q_i|\overline{\lambda })}$$
　3）求 $b_{ij}$：
　式 3 可以写成：
　　$$\sum _I(\sum _{t=1}^{T}ln{b}_{i_t{o}_t})P(I,O|\overline{\lambda })=\sum _{j=1}^N\sum _{t=1}^{T-1}ln{b}_{j{o}_t}P(O,i_t=q_j|\overline{\lambda })$$
　同样有约束条件 ${\sum }_{k=1}^M{b}_{jk}=1$，要注意的是只有在 $o_t=v_k$ 时 $b_{j{o}_t}$ 对 $b_{jk}$ 的偏导数才不为 0，以 $I(o_t=v_k)$ 表示，最后可以得到：
　$$b_{jk}=\frac{{\sum }_{t=1}^{T}P(O,i_t=q_j|\overline{\lambda })I(o_t=v_k)}{{\sum }_{t=1}^{T}P(O,i_t=q_j|\overline{\lambda })}$$

参数估计公式：
得到参数后，可以用 4.1.3 节的 ${\gamma }_t(i),{\xi }_t(i,j)$ 表示：

算法总结：

4.3 预测算法

目的：给定 HMM 模型 $\lambda =(A,B,\pi )$ 和观测序列 $O=o_1,o_2,...,o_T$，求在该观测序列下，最可能对应的状态序列 $I^{\ast }=i_1^{\ast },i_2^{\ast },...,i_T^{\ast }$，也就是最大化 $P(I^{\ast }|O)$。
解决：Viterbi 算法。
其实维特比算法就用动态规划的方法求概率最大路径，计算过程中的每条路径都对应着一个状态序列。计算过程中将最优路径经过的点都保存下来。得到最优路径后，由后向前逐步求得最优结点，这就是维特比算法。
过程：
因为计算过程很简单，就直接给出书中的截图了。
首先导入两个变量 $\delta$和 $\psi$。定义在时刻 $t$ 状态为 $i$ 的所有单个路径 $(i_1,i_2,...,i_t)$ 中概率最大值为：


过程是不是很好理解？如果还不好理解，就继续看个例子。
例子：

整个计算过程如下图所示，

Reference

统计学习方法 李航