机器学习笔记之隐马尔科夫模型

一、序言

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

二、基本概念

2.1 定义

definition
隐马尔可夫模型由初始概率分布、状态转移概率分布以及观测概率分布确定。隐马尔可夫模型的形式定义如下:
设Q是所有可能的状态的集合,V是所有可能的观测的集合。
$Q={q_{1},q_{2},...,q_{n}}, \quad V={v_{1},v_{2},...v_{m}}$
其中, $n$ 是可能的状态数, $m$ 是可能的观测数。
设 $I$ 是长度为 $t$ 的状态序列, $O$ 是对应的观测序列:
$I={i_{1},i_{2},...,i_{t}}, \quad O={o_{1},o_{2},...,o_{t}}$
$A$ 是状态转移概率矩阵:
$A=[a_{ij}]_{n×n}$
其中,
$a_{ij}=P(i_{t}=q_{j}|i_{t-1}=q_{i}), \quad i=1,...n;j=1,...,n$
是在时刻 $t-1$ 处于状态 $q_{i}$ 的条件下在时刻 $t$ 转移到状态 $q_{j}$ 的概率。
$B$ 是观测概率矩阵:
$B=[b_{jk}]_{n×m}$
其中，
$b_{kj}=P(o_{t}=v_{k}|i_{t}=q_{j}), \quad k=1,...m,j=1,...,n$
是在时刻 $t$ 处于状态 $q_{j}$ 的条件下生成观测 $v_{k}$ 的概率。
$\pi$ 是初始状态概率向量:
$\pi=\pi_{i}$
其中，
$\pi_{i}=P(i_{1}=q_{i}), \quad i=1,...n$
是初始时刻 $t=1$ 处于状态 $q_{i}$ 的概率。
definition

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(TN^{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×0.5=0.10$
$t=1$ 时刻取出红球，隐状态是盒子2的概率：
$\alpha_{1}(2) = \pi_{2}b_{2o_{1}} = 0.4×0.4=0.16$
$t=1$ 时刻取出红球，隐状态是盒子3的概率：
$\alpha_{1}(3) = \pi_{3}b_{3o_{1}} = 0.4×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×0.5+0.16×0.3+0.28*0.2)×0.5=0.154×0.5=0.077$
$t=2$ 时刻取出白球，隐状态是盒子2的概率：
$\alpha_{2}(2) = [\sum_{i=1}^3 \alpha_{1}(i)a_{i2}]b_{2o_{2}} = (0.10×0.2+0.16×0.5+0.28*0.3)×0.6=0.184×0.6=0.1104$
$t=2$ 时刻取出白球，隐状态是盒子3的概率：
$\alpha_{2}(3) = [\sum_{i=1}^3 \alpha_{1}(i)a_{i3}]b_{3o_{2}} = (0.10×0.3+0.16×0.2+0.28*0.5)×0.3=0.202×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_{jo_{t+1}}\beta_{t+1}(i)$
如下图所示：

最后，观测概率 $P(O|\lambda) = \sum_{i=1}^n \pi_{i}b_{io_{1}}\beta_{1}(i)$ 。
其实，观测概率 $P(O|\lambda)$ 还可以这么写：
$P(O|\lambda) = \sum_{i=1}^n \sum_{j=1}^n \alpha_{t}(i)a_{ij}b_{jo_{t+1}}\beta_{t+1}(j)$
是不是其实很好理解。

4.1.3 一些概率与期望的计算

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

给定模型 $\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)}$
给定模型 $\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_{jo_{t+1}}\beta_{t+1}(j)$
所以
$\xi_{t}(i,j) = q_{j},O|\lambda) = \frac{\alpha_{t}(i)a_{ij}b_{jo_{t+1}}\beta_{t+1}(j)}{\sum_{i=1}^N \sum_{j=1}^N \alpha_{t}(i)a_{ij}b_{jo_{t+1}}\beta_{t+1}(j)}$

4.2 学习算法

目的：

给定观测序列 $O={o_{1},o_{2},...,o_{T}}$ 和状态序列 $I={i_{1},i_{2},...,i_{T}}$ ，求HMM 模型 $\lambda = (A,B,\pi)$ 的三个参数。
给定观测序列 $O={o_{1},o_{2},...,o_{T}}$ ，求HMM 模型 $\lambda = (A,B,\pi)$ 的三个参数。

解决方法：

监督算法
$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}}$
所以
$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}}))] \\ = \underbrace{ \sum_{I} ln \pi_{i_{1}} P(I,O|\overline \lambda)}_{式1} + \underbrace{ \sum_{I} (\sum_{t=1}^{T-1} ln a_{i_{t}i_{t+1}}) P(I,O|\overline \lambda)}_{式2} + \underbrace{ \sum_{I} (\sum_{t=1}^{T} ln b_{i_{t}o_{t}}) P(I,O|\overline \lambda)}_{式3}$

（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}^Na_{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} lnb_{jo_{t}}P(O,i_{t}=q_{j}|\overline \lambda)$
　同样有约束条件 $\sum_{k=1}^M b_{jk}=1$ ，要注意的是只有在 $o_{t}=v_{k}$ 时 $b_{jo_{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^*={i_{1}^*,i_{2}^*,...,i_{T}^*}$ ，也就是最大化 $P(I^*|O)$ 。
解决：Viterbi 算法。
其实维特比算法就用动态规划的方法求概率最大路径，计算过程中的每条路径都对应着一个状态序列。计算过程中将最优路径经过的点都保存下来。得到最优路径后，由后向前逐步求得最优结点，这就是维特比算法。
过程：
因为计算过程很简单，就直接给出书中的截图了。
首先导入两个变量 $\delta$ 和 $\psi$ 。定义在时刻 $t$ 状态为 $i$ 的所有单个路径 $(i_{1} ,i_{2},...,i_{t} )$ 中概率最大值为：

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

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