隐马尔可夫模型hidden Markov model

本文是《统计学习方法》李航著学习笔记。
为了叙述方便，将hidden Markov model简称HMM。HMM是一种用于标注问题的生成模型，模型工作过程：“隐藏的马尔科夫链”随机生成“不可观测的状态序列”，“每个状态”生成“一个观测”，从而得“观测序列”。在标注问题中，给定“最终的观测序列”，预测其对应的“状态序列”，也称为“标记序列”。可以把整个过程想象成一个网络层，各个状态点和观测点表示网络中的结点。

下面的叙述中主要涉及：
（a）介绍HMM的构成，包括初始状态概率、状态转移概率、观测概率；
（b）在给定模型和观测序列情况下，计算该观测序列出现的概率；
（c）在给定观测序列的情况下，以该观测序列出现概率最大为目标，利用极大似然或者EM算法计算模型参数，就是模型学习过程；
（d）在给定模型和观测序列的情况下，以状态序列的条件概率最大为目标，预测该观测序列对应的状态序列。
其中，（c）在状态序列未知情况下，需要利用EM算法迭代求解HMM参数时，就是以（b）为理论基础实施计算的。

HMM的构成

记HMM的数学表示为$\lambda$，它由初始概率分布$\pi$，状态转移概率分布$A$和观测概率分布$B$确定，即$\lambda=(A,B,\pi)$。

为了给出$\lambda$的具体数学表示，设所有可能的状态的集合为

$$Q=\{q_{1},q_{2},\cdots,q_{N}\}$$
所有的观测的集合为
$$V=\{v_{1},v_{2},\cdots,v_{M}\}$$
考虑长度为 $T$的状态序列$I=(i_{1},i_{2},\cdots,i_{T})$，其对应的观测序列 $O=(o_{1},o_{2},\cdots,o_{T})$，其中对于 $k=1,2,\cdots,T$， $i_{k}$取自 $Q$，$o_{k}$取自 $V$，则
$$初始状态概率向量\pi=(\pi_{i})，\pi_{i}=P(i_{1}=q_{i})\\ 状态转移概率矩阵A=[a_{ij}]_{N\times N}，a_{ij}=P(i_{t+1}=q_{j}|i_{t}=q_{i})\\ 观测概率矩阵B=[b_{j}(k)]_{N\times M}，b_{j}(k)=P(o_{t}=v_{k}|i_{t}=q_{j})$$

观测序列的生成：
输入：HMM$\lambda=(A,B,\pi)$，观测序列长度$T$
输出：观测序列$O=(o_{1},o_{2},\cdots,o_{T})$
（1）t=1，由初始状态概率向量$\pi$产生状态$i_{1}$
（2）按照状态$i_{t}$的观测概率矩阵$B$中的$b_{i_{t}}(k)$生成$o_{t}$
（3）按照状态$i_{t}$的状态转移概率矩阵$A$中的$a_{i_{t}i_{t+1}}$产生状态$i_{t+1}$
（4）置$t=t+1$；$if：t<T，turn：(2)；else：stop$
即$i_{1}\rightarrow o_{1}\rightarrow i_{2}\rightarrow o_{2}\cdots\rightarrow i_{T}\rightarrow o_{T}$；这里需要注意，状态序列与观测序列的生成过程满足齐次马尔科夫性和观测独立性假设。

观测序列概率计算

给定HMM参数，可以通过直接法、前向算法和后向算法计算观测序列的概率。
（一）直接法：利用联合分布和边缘分布的关系，对状态序列进行遍历求和，只是理论上可行，计算量在$O(TN^{T})$阶。
（二）前向算法：根据前向概率向后递推求观测序列概率
前向概率：给定HMM参数$\lambda=(A,B,\pi)$，定义到时刻$t$部分观测序列为$o_{1},o_{2},\cdots,o_{t}$，且此时状态$i_{t}=q_{i}$的概率为前向概率：

$$\alpha_{t}(i)=P(o_{1},o_{2},\cdots,o_{t},i_{t}=q_{i}|\lambda)$$

观测序列概率的前向算法：
输入：HMM$\lambda=(A,B,\pi)$，观测序列$O$
输出：观测序列概率$P(O|\lambda)$
（1）初始前向概率

$$\alpha_{1}(i)=P(o_{1},i_{1}=q_{i}|\lambda)=\pi_{i}b_{i}(o_{1})，i=1,2,\cdots,N$$
（2）向后递推
$$\alpha_{t+1}(i)=P(o_{1},o_{2},\cdots,o_{t+1},i_{t+1}=q_{i}|\lambda)\\ =\sum_{j=1}^{N}\Big[P(o_{1},o_{2},\cdots,o_{t},i_{t}=q_{j}|\lambda)a_{ji}\Big]b_{i}(o_{t+1})\\ =\sum_{j=1}^{N}\Big[\alpha_{t}(j)a_{ji}\Big]b_{i}(o_{t+1})$$
（3）整体观测序列
$$P(O|\lambda)=\sum_{i=1}^{N}P(o_{1},o_{2},\cdots,o_{T},i_{T}=q_{i}|\lambda)=\sum_{i=1}^{N}\alpha_{T}(i)$$
即， 在计算$\alpha_{t+1}(i)$时，不关注$t$时刻的$i_{t}$处于哪个状态，而$\alpha_{t}(i)$中包含$i_{t}=q_{j}$的状态约束，要想通过$\alpha_{t}(i)$递推得$\alpha_{t+1}(i)$，必须遍历所以状态，去掉某个状态$q_{j}$的影响。可类似看作网络层中，后一层每个节点的计算，需要用到前一层的每个节点。该算法的计算量是$O(N^2 T)$。

（三）后向算法：根据后向概率向前递推求观测序列概率
后向概率：给定HMM参数$\lambda=(A,B,\pi)$，定义在时刻$t$状态为$i_{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)$$

观测序列概率的后向算法：
输入：HMM$\lambda=(A,B,\pi)$，观测序列$O$
输出：观测序列概率$P(O|\lambda)$
（1）规定初始后向概率

$$\beta_{T}(i)=1$$
（2）向前递推
$$\beta_{t}(i)=P(o_{t+1},o_{t+2},\cdots,o_{T}|i_{t}=q_{i},\lambda)\\ =\sum_{j=1}^{N}a_{ij}b_{j}(o_{t+1})P(o_{t+2},o_{t+3},\cdots,o_{T}|i_{t+1}=q_{j},\lambda)\\ =\sum_{j=1}^{N}a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j)$$
（3）整体观测序列
$$P(O|\lambda)=P(o_{1},o_{2},\cdots,o_{T}|\lambda)=P(o_{1})P(o_{2},o_{3},\cdots,o_{T}|\lambda)\\ =\sum_{i=1}^{N}P(i_{1}=q_{i})P(o_{1}|i_{1}=q_{i})P(o_{2},o_{3},\cdots,o_{T}|i_{1}=q_{i},\lambda)=\sum_{i=1}^{N}\pi_{i}b_{i}(o_{1})\beta_{1}(i)$$

最后，在给定HMM$\lambda=(A,B,\pi)$和观测序列$O$的条件下，可利用前向概率和后验概率计算单个状态和两个相邻状态：
（1）时刻$t$处于状态$q_{i}$：

$$\gamma_{i}(t)=P(i_{t}=q_{i}|O,\lambda)=\frac{P(i_{t}=q_{i},O|\lambda)}{P(O|\lambda)}=\frac{P(i_{t}=q_{i},O|\lambda)}{\sum\limits_{j=1}^{N}P(i_{t}=q_{j},O|\lambda)}\\ =\frac{\alpha_{t}(i)\beta_{t}(i)}{\sum\limits_{j=1}^{N}\alpha_{t}(j)\beta_{t}(j)}=\frac{\sum\limits_{j=1}^{N}\alpha_{t}(i)a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j)}{\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}\alpha_{t}(i)a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j)}$$
（2）时刻 $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)}{P(O|\lambda)}\\ =\frac{P(i_{t}=q_{i},i_{t+1}=q_{j},O|\lambda)}{\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}P(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\limits_{i=1}^{N}\sum\limits_{j=1}^{N}\alpha_{t}(i)a_{ij}b_{j}(o_{t+1})\beta_{t+1}(j)}$$

HMM生成，参数学习算法

（1）给定观测序列的同时，也给定了状态序列，用极大似然估计得$\lambda=(A,B,\pi)$。有关MLE的内容参考http://blog.csdn.net/cymy001/article/details/78016109

若训练数据含$S$个等长的观测序列和对应的状态序列$\{(O_{1},I_{1}),(O_{2},I_{2}),\cdots,(O_{S},I_{S})\}$，则记$S$个序列的每个状态序列中，任何起始时刻，只要该起始时刻处于状态$i$，并且它的下一时刻处于状态$j$，则计入频数统计量$A_{ij}$，那么可得$a_{ij}$的MLE

$$a_{ij}=\frac{A_{ij}}{\sum\limits_{j=1}^{N}A_{ij}}$$
同样， $S$个序列的每个状态序列中，任何起始时刻，只要状态序列在该时刻处于状态$j$，观测序列在该时刻处于观测 $k$，则计入频数统计量$B_{jk}$，可得 $b_{j}(k)$的MLE
$$b_{j}(k)=\frac{B_{jk}}{\sum\limits_{k=1}^{M}B_{jk}}$$
对于初始状态概率，只需要统计状态序列集合 $\{O_{1},O_{2},\cdots,O_{S}\}$中每个序列的第一个元素状态频数，再除以 $S$，即得各个初始状态概率的MLE。

（2）只给观测序列，可以把状态序列看成隐变量，用EM算法得$\lambda=(A,B,\pi)$。有关EM算法的内容参考 http://blog.csdn.net/cymy001/article/details/78264648
目标是使含隐变量（这里是“状态变量”）的观测序列概率随着参数 $\lambda$的迭代计算，逐步增大，而这个问题可以通过降低下界的办法转化成EM算法的 $Q$函数极大化问题，即

$$\tilde{Q}(\lambda,\bar{\lambda})=E_{I}[logP(O,I|\lambda)|O,\bar{\lambda}] =\sum_{I}logP(O,I|\lambda)P(I|O,\bar{\lambda})\\ =\sum_{I}logP(O,I|\lambda)\frac{P(O,I|\bar{\lambda})}{P(O|\bar{\lambda})} \propto\sum_{I}logP(O,I|\lambda)P(O,I|\bar{\lambda})=Q(\lambda,\bar{\lambda})$$
给定待优化的变量记号：初始状态概率 $\pi_{i_{1}}$，状态转移概率 $a_{i_{t}i_{t+1}}$，观测概率 $b_{i_{t}}(o_{t})$，则
$$P(O,I|\lambda)=\pi_{i_{1}}b_{i_{1}}(o_{1})a_{i_{1}i_{2}}b_{i_{2}}(o_{2})\cdots a_{i_{T-1}i_{T}}b_{i_{T}}(o_{T})$$
于是
$$Q(\lambda,\bar{\lambda})=\sum_{I}logP(O,I|\lambda)P(O,I|\bar{\lambda})\\ =\sum_{I}log\pi_{i_{1}}P(O,I|\bar{\lambda})+ \sum_{I}\Big(\sum_{t=1}^{T-1}loga_{i_{t}i_{t+1}}\Big)P(O,I|\bar{\lambda})\\ +\sum_{I}\Big(\sum_{t=1}^{T}logb_{i_{t}}(o_{t})\Big)P(O,I|\bar{\lambda})\\ =\sum_{i=1}^{N}log\pi_{i}P(O,i_{1}=i|\bar{\lambda}) +\sum_{i=1}^{N}\sum_{j=1}^{N}\Big(\sum_{t=1}^{T-1}loga_{ij}\Big)P(O,i_{t}=i,i_{t+1}=j|\bar{\lambda})\\ +\sum_{j=1}^{N}\Big(\sum_{t=1}^{T}logb_{j}(o_{t})\Big)P(O,i_{t}=j|\bar{\lambda})$$
分别对上式和的三个部分利用 $Lagrange$乘子法求解，即
$$L(\pi_{i})=\sum_{i=1}^{N}log\pi_{i}P(O,i_{1}=i|\bar{\lambda})+\theta_{1}\Big(\sum_{i=1}^{N}\pi_{i}-1\Big)\\ L(a_{ij})=\sum_{i=1}^{N}\sum_{j=1}^{N}\Big(\sum_{t=1}^{T-1}loga_{ij}\Big)P(O,i_{t}=i,i_{t+1}=j|\bar{\lambda})+\theta_{2}(\sum_{j=1}^{N}a_{ij}-1)$$
$$L(b_{j}(k))=\sum_{j=1}^{N}\Big(\sum_{t=1}^{T}logb_{j}(o_{t})\Big)P(O,i_{t}=j|\bar{\lambda})+\theta_{3}(\sum_{k=1}^{M}b_{j}(k)-1)$$
分别对 $\pi_{i},a_{ij},b_{j}(k)$求导进行相关计算，即可得HMM的参数
$$\pi_{i}=\frac{P(O,i_{1}=i|\bar{\lambda})}{P(O|\bar{\lambda})} =\gamma_{1}(i)\\ a_{ij}=\frac{\sum\limits_{t=1}^{T-1}P(O,i_{t}=i,i_{t+1}=j|\bar{\lambda})}{\sum\limits_{t=1}^{T-1}P(O,i_{t}=i|\bar{\lambda})} =\frac{\sum\limits_{t=1}^{T-1}\xi_{t}(i,j)}{\sum\limits_{t=1}^{T-1}\gamma_{t}(i)}\\ b_{j}(k)=\frac{\sum\limits_{t=1}^{T}P(O,i_{t}=j|\bar{\lambda})I(o_{t}=v_{k})}{\sum\limits_{t=1}^{T}P(O,i_{t}=j|\bar{\lambda})} =\frac{\sum\limits_{t=1,o_{t}=v_{k}}^{T}\gamma_{t}(j)}{\sum\limits_{t=1}^{T}\gamma_{t}(j)}$$

Baum_Welch算法：
输入：观测序列$O$
输出：HMM的参数$\lambda=(A,B,\pi)$
（1）初始模型参数$\lambda^{(0)}=(A^{(0)},B^{(0)},\pi^{(0)})$
（2）利用上述推导公式迭代计算

$$a_{ij}^{(n+1)}=\frac{\sum\limits_{t=1}^{T-1}\xi_{t}(i,j)}{\sum\limits_{t=1}^{T-1}\gamma_{t}(i)}， b_{j}(k)^{(n+1)}=\frac{\sum\limits_{t=1,o_{t}=v_{k}}^{T}\gamma_{t}(j)}{\sum\limits_{t=1}^{T}\gamma_{t}(j)}， \pi_{i}^{(n+1)}=\gamma_{1}(i)$$
（3）满足终止条件时，输出模型参数 $\lambda^{(n+1)}=(A^{(n+1)},B^{(n+1)},\pi^{(n=1)})$

HMM状态序列预测算法

给定HMM和观测序列，预测该观测序列对应的状态序列，也就是其对应的标签序列。
（1）贪婪的近似算法：即每一时间步都选取“当前时刻对应的状态概率最大的状态点”，这种办法会产生实际状态转移矩阵$A$中相邻转移概率为$0$的矛盾相邻转移状态。

（2）维特比算法：利用动态规划的思想，将寻找最优状态路径的问题递归成最优子状态路径问题，求解的过程是从前向后按层递推求解最大概率路径。
有关动态规划请参考http://blog.csdn.net/cymy001/article/details/78474263

定义在时刻$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_{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}(j)a_{ji}]b_{i}(o_{t+1})，i=1,2,\cdots,N$$

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

维特比算法：
输入：HMM参数$\lambda=(A,B,\pi)$，观测序列$O$
输出：最优状态序列$I^{*}$
（1）初始化

$$\delta_{1}(i)=\pi_{i}b_{i}(o_{1}),\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})，\psi_{t}(i)=arg\max_{1\leq j\leq N}[\delta_{t-1}(j)a_{ji}]，i=1,2,\cdots,N$$
（3）终止
$$P^{*}=\max_{1\leq i\leq N}\delta_{T}(i)，i_{T}^{*}=arg\max_{1\leq i\leq N}[\delta_{T}(i)]$$
（4）最优路径回溯
对 $t=T-1,T-2,\cdots,1$，求
$$i_{t}^{*}=\psi_{t+1}(i_{t+1}^{*})$$
即可得最优路径 $I^{*}$。
求解最优路径的过程，就是一个对状态转移中按层形成的完全有向图状态转移网络进行最大化过滤的过程。

如上图所示，这一个将状态变化过程展开的网络图，在前面前向概率和后向概率计算时，是对相邻层之间的状态点进行了遍历（即上图细线，第一层状态123分别会与第二层状态123连接映射），但是维比特算法中，只会追踪映射出来的对应本层概率最大的前一层点，比如绿色层状态1，可以由黑色层123三个点分别映射出对应路径，但是黑色层状态2映射到绿色层状态1的路径观测概率最大，所以只考虑这条路径。