HMM模型及相关算法

HMM模型及相关算法

写在前面：本文主要参考了刘建平Pinard老师的博客，做了一定程度的归纳，其中也有一些自己的理解，包括图和公式，希望对大家学习有所帮助。

一、HMM定义和前置知识

1、条件独立的判定

贝叶斯网络的head-to-tail结构形式如下图：

还是分c未知跟c已知这两种情况：
c未知时，有：$P(a,b,c)=P(a)\ast P(c|a)\ast P(b|c)$，但无法推出$P(a,b)=P(a)P(b)$，即c未知时，a、b不独立。
c已知时，有：$P(a,b|c)=P(a,b,c)/P(c)$，且根据$P(a,c)=P(a)\ast P(c|a)=P(c)\ast P(a|c)$，可化简得到：
$$\begin{gathered} P(a,b|c)=P(a,b,c)/P(c) \\ =P(a)\ast P(c|a)\ast P(b|c)/P(c) \\ =P(a,c)\ast P(b|c)/P(c) \\ =P(a|c)\ast P(b|c) \end{gathered}$$
在c给定的条件下，a，b被阻断(blocked)，是独立的，称之为head-to-tail条件独立。

2、动态规划问题

引入递归问题：已知$a_n-a_{n-1}=2,a_1=1$ 求解a的通项公式

求解使用错位相加：
$$\begin{gathered} a_n-a_{n-1}=2 \\ {a}_{n-1}-a_{n-2}=2 \\ ... \\ {a}_2-a_1=2 \\ -------------------------------------- \\ {a}_n-a_{n-1}+a_{n-1}-a_{n-2}...+a_2-a_1=2\ast (n-1) \\ {a}_n=2n-1 \\ {S}_n=a_1+a_2...+a_n=(2n-1+1)\ast n/2=n^2 \end{gathered}$$
其实所谓递归问题通常就是这样一个求通项或者求通项和的问题。

例子：动态规划在分词上的应用

例句：“经常有意见分歧”
词典： [“经常”, “经”， “有”, “有意见”，“意见”，“分歧”，“见”， “意”， “见分歧”， “分”]
概率：[ 0.1, 0.05, 0.1, 0.1, 0.2, 0.2, 0.05, 0.05, 0.05, 0.1]
-log(x):[ 2.3, 3, 2.3, 2.3, 1.6, 1.6, 3, 3, 3, 2.3]

动态规划求解过程如下：
$$\begin{gathered} f(8)=f(7)+20=f(6)+1.6=f(5)+3 \\ f(7)=f(6)+2.3 \\ f(6)=f(5)+3=f(4)+1.6=f(3)＋2.3 \\ f(5)=f(4)+3 \\ f(4)=f(3)+2.3 \\ f(3)=f(2)+20=f(1)+3 \\ f(2)=3 \\ f(1)=0 \end{gathered}$$
有上式：

可得最佳路径

1	2	3	4	5	6	7	8
0	3	2.3	4.6	7.6	4.6	6.9	6.2

最佳路径：8<–6<–3<–1:[经常，有意见，分歧]

动态规划的特点，后一个时间节点的路径可以用前面路径节点递推得到。

3、隐马尔科夫（HMM）模型

隐马尔科夫模型（Hidden Markov Model，以下简称HMM）是比较经典的机器学习模型了，它在语言识别，自然语言处理，模式识别等领域得到广泛的应用。当然，随着目前深度学习的崛起，尤其是RNN，LSTM等神经网络序列模型的火热，HMM的地位有所下降。但是作为一个经典的模型，学习HMM的模型和对应算法，对我们解决问题建模的能力提高以及算法思路的拓展还是很好的。

1. 什么样的问题需要HMM模型

首先我们来看看什么样的问题解决可以用HMM模型。使用HMM模型时我们的问题一般有这两个特征：

１）我们的问题是基于序列的，比如时间序列，或者状态序列。

２）我们的问题中有两类数据，一类序列数据是可以观测到的，即观测序列；而另一类数据是不能观察到的，即隐藏状态序列，简称状态序列。

下面我们用精确的数学符号来表述我们的HMM模型。

2. HMM模型的定义

对于HMM模型，首先我们假设Q 是所有可能的隐藏状态的集合，V 是所有可能的观测状态的集合，即：
$$\begin{gathered} Q=q_1,q_2,...,q_N, \\ V=v_1,v_2,...v_M \end{gathered}$$
其中，N 是可能的隐藏状态数，M是所有的可能的观察状态数。

对于一个长度为T 的序列，I 对应的状态序列, O 是对应的观察序列，即：
$$\begin{gathered} I=i_1,i_2,...,i_T, \\ O=o_1,o_2,...o_T \end{gathered}$$
其中，任意一个隐藏状态$i_T\in Q$ ,任意一个观察状态$o_T\in V$

HMM模型做了两个很重要的假设如下：

（1）齐次马尔科夫链假设。
即任意时刻的隐藏状态只依赖于它前一个隐藏状态。当然这样假设有点极端，因为很多时候我们的某一个隐藏状态不仅仅只依赖于前一个隐藏状态，可能是前两个或者是前三个。但是这样假设的好处就是模型简单，便于求解。如果在时刻t 的隐藏状态是$i_t=q_i$ ,在时刻t+1的隐藏状态是 $i_{t+1}=q_j$ , 则从时刻t 到时刻t+1的HMM状态转移概率 $a_{ij}$ 可以表示为：
$$a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$$
这样$a_{ij}$ 可以组成马尔科夫链的状态转移矩阵A :
$$A=[a_{ij}{]}_{N\times N}$$
（2）观测独立性假设。
即任意时刻的观察状态只仅仅依赖于当前时刻的隐藏状态，这也是一个为了简化模型的假设。如果在时刻t的隐藏状态是$i_t=q_j$ , 而对应的观察状态为$o_t=v_k$ , 则该时刻观察状态$v_k$ 在隐藏状态$q_j$下生成的概率为$b_j(k)$,满足：
$$b_j(k)=P(o_t=v_k|i_t=q_j)$$
这样$b_j(k)$ 可以组成观测状态生成的概率矩阵B :
$$B=[b_j(k){]}_{N\times M}$$
除此之外，我们需要一组在时刻t=1的隐藏状态概率分布$\Pi$:
$$\Pi =[\pi (i){]}_N,其中\pi (i)=P(i_1=q_i)$$
一个HMM模型，可以由隐藏状态初始概率分布$\Pi$, 状态转移概率矩阵A 和观测状态概率矩阵B决定。$\Pi ,A$ 决定状态序列，B 决定观测序列。因此，HMM模型参数可以由一个三元组λ 表示如下：
$$\lambda =(A,B,\Pi )$$

二、HMM中的三个问题

（1） 评估观察序列概率。

即给定模型$\lambda =(A,B,\Pi )$ 和观测序列$O=o_1,o_2,...o_T$，计算在模型λ下观测序列O出现的概率$P(O|\lambda )$。这个问题的求解需要用到前向后向算法。这个问题是HMM模型三个问题中最简单的。

首先我们回顾下HMM模型的问题。这个问题是这样的。我们已知HMM模型的参数$\lambda =(A,B,\Pi )$ 。其中A 是隐藏状态转移概率的矩阵，B 是观测状态生成概率的矩阵， Π是隐藏状态的初始概率分布。同时我们也已经得到了观测序列$O=o_1,o_2,...o_T$ ,现在我们要求观测序列O 在模型λ 下出现的条件概率$P(O|\lambda )$ 。

我们知道所有的隐藏状态之间的转移概率和所有从隐藏状态到观测状态生成概率，那么我们是可以暴力求解的。

我们可以列举出所有可能出现的长度为T 的隐藏序列$I=i_1,i_2,...,i_T$ ,分布求出这些隐藏序列与观测序列$O=o_1,o_2,...o_T$ 的联合概率分布$P(O,I|\lambda )$ ，这样我们就可以很容易的求出边缘分布$P(O|\lambda )$ 了。

1)暴力求解

首先，任意一个隐藏序列$I=i_1,i_2,...,i_T$ 出现的概率是：
$$P(I|\lambda )={\pi }_{i1}{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)$$
则O和I联合出现的概率是：
$$P(O,I|\lambda )=P(I|\lambda )P(O|I,\lambda )={\pi }_{i1}{b}_{i1}(o_1)a_{i_1{i}_2}{b}_{i_2}(o_2)...a_{i_{T-1}{i}_T}{b}_{i_T}(o_T)$$
然后求边缘概率分布，即可得到观测序列O 在模型λ下出现的条件概率$P(O|\lambda )$ ：
$$P(O|\lambda )=\sum _{I}P(O,I|\lambda )=\sum _{i_1,i_2,...i_T}{\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)$$
虽然上述方法有效，但是如果我们的隐藏状态数N非常多的那就麻烦了，此时我们预测状态有$N^T$种组合，算法的时间复杂度是$O(T{N}^T)$ 阶的。因此对于一些隐藏状态数极少的模型，我们可以用暴力求解法来得到观测序列出现的概率，但是如果隐藏状态多，则上述算法太耗时，我们需要寻找其他简洁的算法。

前向后向算法就是来帮助我们在较低的时间复杂度情况下求解这个问题的。

2)用前向算法求HMM观测序列的概率

前向算法本质上属于动态规划的算法，也就是我们要通过找到局部状态递推的公式，这样一步步的从子问题的最优解拓展到整个问题的最优解。

在前向算法中，通过定义“前向概率”来定义动态规划的这个局部状态。什么是前向概率呢, 其实定义很简单：定义时刻t 时隐藏状态为$q_i$ , 观测状态的序列为$o_1,o_2,...o_t$ 的概率为前向概率。记为：
$${\alpha }_t(i)=P(o_1,o_2,...o_t,i_t=q_i|\lambda )$$
既然是动态规划，就要找递推了，现在假设我们已经找到了在时刻t 时各个隐藏状态的前向概率，现在我们需要递推出时刻t+1时各个隐藏状态的前向概率。

从下图可以看出，我们可以基于时刻t时各个隐藏状态的前向概率，再乘以对应的状态转移概率，即${\alpha }_t(j)a_{ji}$ 就是在时刻t 观测到$o_1,o_2,...o_t$，并且时刻t 隐藏状态$q_j$ , 时刻t+1 隐藏状态$q_i$ 的概率。如果将想下面所有的线对应的概率求和，即 ${\sum }_{j=1}^N{\alpha }_t(j)a_{ji}$ 就是在时刻t观测到$o_1,o_2,...o_t$ ，并且时刻t+1隐藏状态$q_i$ 的概率。继续一步，由于观测状态$o_{t+1}$只依赖于t+1时刻隐藏状态$q_i$ , 这样$[{\sum }_{j=1}^N{\alpha }_t(j)a_{ji}]b_i(o_{t+1})$ 就是在在时刻t+1观测到$o_1,o_2,...o_t，o_{t+1}$，并且时刻t+1隐藏状态$q_i$的概率。而这个概率，恰恰就是时刻t+1对应的隐藏状态i 的前向概率，这样我们得到了前向概率的递推关系式如下：
$${\alpha }_{t+1}(i)=[\sum _{j=1}^N{\alpha }_t(j)a_{ji}]b_i(o_{t+1})$$

我们的动态规划从时刻1开始，到时刻T 结束，由于${\alpha }_T(i)$ 表示在时刻T 观测序列为$o_1,o_2,...o_T$，并且时刻T 隐藏状态qi的概率，我们只要将所有隐藏状态对应的概率相加，即${\sum }_{i=1}^N{\alpha }_T(i)$ 就得到了在时刻T观测序列为$o_1,o_2,...o_T$的概率。

下面总结下前向算法：

输入：HMM模型$\lambda =(A,B,\Pi )$，观测序列$O=(o_1,o_2,...o_T)$

输出：观测序列概率$P(O|\lambda )$
1、 计算时刻1的各个隐藏状态前向概率：
$${\alpha }_1(i)={\pi }_i{b}_i(o_1),i=1,2,...N$$
2、 递推时刻2,3,…T 时刻的前向概率：
$${\alpha }_{t+1}(i)=[\sum _{j=1}^N{\alpha }_t(j)a_{ji}]b_i(o_{t+1}),i=1,2,...N$$
3、计算最终结果：
$$P(O|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$$
从递推公式可以看出，我们的算法时间复杂度是$O(T{N}^2)$ ，比暴力解法的时间复杂度$O(T{N}^T)$ 少了几个数量级。

与之相对的还有后向算法也是类似的思想，这里就不详细说了。
$${\beta }_t(i)=P(o_{t+1},o_{t+2},...o_T|i_t=q_i,\lambda )$$
输入：HMM模型$\lambda =(A,B,\Pi )$ ，观测序列$O=(o_1,o_2,...o_T)$
输出：观测序列概率$P(O|\lambda )$

1、初始化时刻T 的各个隐藏状态后向概率：
$${\beta }_T(i)=1,i=1,2,...N$$
2、递推时刻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$$
3、计算最终结果：
$$P(O|\lambda )=\sum _{i=1}^N{\pi }_i{b}_i(o_1){\beta }_1(i)$$
利用前向概率和后向概率，我们可以计算出HMM中单个状态概率公式。

给定模型λ和观测序列O ,在时刻t 处于状态$q_i$ 的概率记为:
$${\gamma }_t(i)=P(i_t=q_i|O,\lambda )=P(i_t=q_i,O|\lambda )P(O|\lambda )$$
利用前向概率和后向概率的定义可知：
$$P(i_t=q_i,O|\lambda )={\alpha }_t(i){\beta }_t(i)$$

（2）模型参数学习问题

即给定观测序列$O=o_1,o_2,...o_T$ ，估计模型$\lambda =(A,B,\Pi )$ 的参数，使该模型下观测序列的条件概率$P(O|\lambda )$ 最大。这个问题的求解需要用到基于EM算法。

对于m个样本观察数据$x=(x^{(1)},x^{(2)},...x^{(m)})$ 中，找出样本的模型参数λ, 极大化模型分布的对数似然函数如下：
$$\lambda :=\underset{\lambda }{\mathrm{argmax}}\sum _{i=1}^{m}logP(x^{(i)}|\lambda )$$
如果我们得到的观察数据有未观察到的隐含数据$z=(z^{(1)},z^{(2)},...z^{(m)})$ ，此时我们的极大化模型分布的对数似然函数如下：
$$\lambda :=\underset{\lambda }{\mathrm{argmax}}\sum _{i=1}^{m}log\sum _{z(i)}P(x^{(i)},z^{(i)}|\lambda )$$
上面这个式子是没有 办法直接求出λ 的。因此需要一些特殊的技巧，我们首先对这个式子进行缩放如下：
$$\begin{gathered} \sum _{i=1}^{m}log\sum _{z(i)}P(x^{(i)},z^{(i)}|\lambda )=\sum _{i=1}^{m}log\sum _{z(i)}{Q}_i(z^{(i)})\frac{P(x^{(i)},z^{(i)}|\lambda )}{Q_i(z^{(i)})}(1) \\ \ge \sum _{i=1}^m\sum _{z(i)}{Q}_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)}|\lambda )}{Q_i(z^{(i)})}(2) \end{gathered}$$
上面第(1)式引入了一个未知的新的分布$Q_i(z^{(i)})$，第(2)式用到了Jensen不等式：
$$log\sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_{j}log{y}_j,{\lambda }_j\ge 0,\sum _j{\lambda }_j=1$$
或者说由于对数函数是凹函数，所以有:
$$f(E(x))\ge E(f(x)),如果f(x)是凹函数$$
此时，如果要满足Jensen不等式的等号，则有：
$$\frac{P(x^{(i)},z^{(i)}|\lambda )}{Q_i(z^{(i)})}=c,c为常数$$
由于$Q_i(z^{(i)})$是一个分布，所以满足：
$$\sum _z{Q}_i(z^{(i)})=1$$
从上面两式，我们可以得到：
$$Q_i(z^{(i)})=\frac{P(x^{(i)},z^{(i)}|\lambda )}{{\sum }_{z}P(x^{(i)},z^{(i)}|\lambda )}=\frac{P(x^{(i)},z^{(i)}|\lambda )}{P(x^{(i)}|\lambda )}=P(z^{(i)}|x^{(i)},\lambda ))$$
如果$Q_i(z^{(i)})=P(z^{(i)}|x^{(i)},\lambda ))$ , 则第(2)式是我们的包含隐藏数据的对数似然的一个下界。如果我们能极大化这个下界，则也在尝试极大化我们的对数似然。即我们需要最大化下式：
$$\underset{\lambda }{\mathrm{argmax}}\sum _{i=1}^m\sum _{z(i)}{Q}_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)}|\lambda )}{Q_i(z^{(i)})}(3)$$
(3)式等价为：
$$\underset{\lambda }{\mathrm{argmax}}\sum _{i=1}^m(\sum _{z(i)}{Q}_i(z^{(i)})logP(x^{(i)},z^{(i)}|\lambda )-\sum _{z(i)}{Q}_i(z^{(i)})log{Q}_i(z^{(i)}))(4)$$
去掉上式中为常数的部分(后面减掉的${\sum }_{z(i)}{Q}_i(z^{(i)})log{Q}_i(z^{(i)})$)，这部分求$\lambda$过程中是不用考虑的。

则我们需要极大化的对数似然下界为：
$$\begin{gathered} \underset{\lambda }{\mathrm{argmax}}\sum _{i=1}^m\sum _{z^{(i)}}{Q}_i(z^{(i)})logP(x^{(i)},z^{(i)}|\lambda ) \\ \underset{\lambda }{\mathrm{argmax}}\sum _{i=1}^m\sum _{z^{(i)}}P(z^{(i)}|x^{(i)},\lambda ))logP(x^{(i)},z^{(i)}|\lambda ) \end{gathered}$$
上式也就是我们的EM算法的M步，那E步呢？注意到上式中$Q_i(z^{(i)})$是一个分布，因此${\sum }_{z^{(i)}}{Q}_i(z^{(i)})logP(x^{(i)},z^{(i)}|\lambda )$可以理解为$logP(x^{(i)},z^{(i)}|\lambda )$基于条件概率分布$Q_i(z^{(i)})$的期望。

至此，我们理解了EM算法中E步和M步的具体数学含义。

（3）预测问题，也称为解码问题

即给定模型$\lambda =(A,B,\Pi )$ 和观测序列$O=o_1,o_2,...o_T$ ，求给定观测序列条件下，最可能出现的对应的状态序列，这个问题的求解需要用到基于动态规划的维特比算法，这个问题是HMM模型三个问题中复杂度居中的算法。

在HMM模型的解码问题中，给定模型 $\lambda =(A,B,\Pi )$ 和观测序列$O=o_1,o_2,...o_T$ ，求给定观测序列O条件下，最可能出现的对应的状态序列 $I=i_1,i_2,...i_T$ ,即$P(I|O)$ 要最大化。

（1）近似解法

求出观测序列O在每个时刻t最可能的隐藏状态 $i_t$然后得到一个近似的隐藏状态序列 $I=i_1,i_2,...i_T$ 。要这样近似求解不难，给定模型λ 和观测序列O 时，在时刻t 处于状态$q_i$ 的概率是$ γ_t(i)​$ ，这个概率可以通过HMM的前向算法与后向算法计算。这样我们有：
$$i_t=\underset{1\le i\le N}{\mathrm{argmax}}[{\gamma }_t(i)],t=1,2,...T$$
近似算法很简单，但是却不能保证预测的状态序列是整体是最可能的状态序列，因为预测的状态序列中某些相邻的隐藏状态可能存在转移概率为0的情况。

而维特比算法可以将HMM的状态序列作为一个整体来考虑，避免近似算法的问题，下面我们来看看维特比算法进行HMM解码的方法。

（2）维特比算法

维特比算法是一个通用的解码算法，是基于动态规划的求序列最短路径的方法。在分词原理中已经讲到了维特比算法的一些细节。

既然是动态规划算法，那么就需要找到合适的局部状态，以及局部状态的递推公式。在HMM中，维特比算法定义了两个局部状态用于递推。

第一个局部状态是在时刻t 隐藏状态为i所有可能的状态转移路径 $i_1,i_2,...i_t$ 中的概率最大值。记为${\delta }_t(i)$ :
$${\delta }_t(i)=\underset{i_1,i_2,...i_{t-1}}{\max }P(i_t=i,i_1,i_2,...i_{t-1},o_t,o_{t-1},...o_1|\lambda ),i=1,2,...N$$
由$\delta t(i)$ 的定义可以得到δ的递推表达式：
$$\begin{gathered} {\delta }_{t+1}(i)=\underset{i1,i2,...it}{\max }P(i_{t+1}=i,i_1,i_2,...i_t,o_{t+1},o_t,...o_1|\lambda )(1) \\ =\underset{1\le j\le N}{\max }[{\delta }_t(j)a_{ji}]b_i(o_{t+1})(2) \end{gathered}$$
第二个局部状态由第一个局部状态递推得。我们定义在时刻t隐藏状态为i 的所有单个状态转移路径$(i_1,i_2,...,i_{t-1},i)$ 中概率最大的转移路径中第t−1个节点的隐藏状态为${\Psi }_t(i)$ ,其递推表达式可以表示为：
$${\Psi }_t(i)=\underset{1\le j\le N}{\mathrm{argmax}}[{\delta }_{t-1}(j)a_{ji}]$$
有了这两个局部状态，我们就可以从时刻0一直递推到时刻T，然后利用${\Psi }_t(i)$ 记录的前一个最可能的状态节点回溯，直到找到最优的隐藏状态序列。

现在来总结下维特比算法的流程：

输入：HMM模型$\lambda =(A,B,\Pi )$ ，观测序列$O=(o_1,o_2,...o_T)$

输出：最有可能的隐藏状态序列$I=i_1,i_2,...i_T$

1、初始化局部状态：
$${\delta }_1(i)={\pi }_i{b}_i(o_1),i=1,2...N$$

$${\Psi }_1(i)=0,i=1,2...N$$

2、 进行动态规划递推时刻t=2,3,…T 时刻的局部状态：
$$\begin{gathered} {\delta }_t(i)=\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]b_i(0_t),i=1,2...N \\ {\Psi }_t(i)=\underset{1\le j\le N}{\mathrm{argmax}}[{\delta }_{t-1}(j)a_{ji}],i=1,2...N \end{gathered}$$
3、计算时刻T 最大的${\delta }_T(i)$ ,即为最可能隐藏状态序列出现的概率。计算时刻T 最大的${\Psi }_t(i)$ ,即为时刻T 最可能的隐藏状态。
$$P=\underset{1\le j\le N}{\max }{\delta }_T(i)$$

$$i_T=\underset{1\le j\le N}{\mathrm{argmax}}[{\delta }_T(i)]$$

4、 利用局部状态$\Psi (i)$ 开始回溯。对于$t=T-1,T-2,...,1$ ：
$$i_t={\Psi }_{t+1}(i_{t+1})$$
计算时刻T 最大的${\delta }_T(i)$ ,即为最可能隐藏状态序列出现的概率。计算时刻T 最大的${\Psi }_t(i)$ ,即为时刻T 最可能的隐藏状态。