统计学习方法之隐马尔可夫模型

隐马尔可夫模型

1. 基本概念

1.1 定义

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测从而产生观测随机序列的过程。隐马尔科夫链生成的状态序列，称为状态序列；每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列。序列的每一个位置又可看作为一个时刻。

隐马尔可夫模型的形式定义如下：

设$Q$是所有可能的状态的集合，$V$是所有可能的观测的集合：
Q = { q 1 , q 2 , . . . , q N } , V = { v 1 , v 2 , . . . , v M } Q=\{q_1,q_2,...,q_N\}, \quad V=\{v_1,v_2,...,v_M\} Q={q1​,q2​,...,qN​},V={v1​,v2​,...,vM​}
其中，$N$是可能的状态数，$M$是可能的观测数。

$I$是长度为$T$的状态序列，$O$是对应的观测序列
I = { i 1 , i 2 , . . . , i T } , O = { o 1 , o 2 , . . , o T } I=\{i_1,i_2,...,i_T\}, \quad O=\{o_1,o_2,..,o_T\} I={i1​,i2​,...,iT​},O={o1​,o2​,..,oT​}

$A$是状态转移概率矩阵
$$A=[a_{ij}{]}_{N\times N}$$
其中，
$$a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$$
是在时刻$t$处于状态$q_i$的条件下在时刻$t+1$转移到状态$q_j$的概率

$B$是观测概率矩阵
$$B=[b_j(k){]}_{N\times M}$$
其中，
$$b_j(k)=P(o_t=v_k|i_t=q_j)$$
是在时刻$t$处于状态$q_j$的条件下生成观测$v_k$的概率

$\pi$是初始状态概率向量
$$\pi =({\pi }_i)$$
其中，
$${\pi }_i=P(i_1=q_i)$$
是时刻$t=1$处于状态$q_i$的概率

隐马尔可夫模型由初始状态概率向量$\pi$、状态转移概率矩阵$A$和概念概率矩阵$B$决定。因此，隐马尔可夫模型$\lambda$可以用三元符合表示：
$$\lambda =(A,B,\pi )$$

注意：隐马尔可夫模型作了两个基本假设：

1. 齐次马尔可夫性假设： 隐马尔科夫链在任意时刻$t$的状态只依赖于前一时刻的状态，与其他时刻的状态和观测无关
   $$P(i_t|i_{t-1},o_{t-1},...,i_1,o_1)=P(i_t|i_{t-1})$$
2. 观测独立性假设： 任意时刻的观测只依赖于该时刻的马尔科夫链的状态，与其他观测和状态无关
   $$P(o_t|i_T,o_T,...,i_{t+1},o_{t+1},i_t,i_{t-1},o_{t-1},...,i_1,o_1)=P(o_t|i_t)$$

1.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 )$最大。即用EM算法估计参数
3. 预测问题： 已知模型$\lambda =(A,B,\pi )$和观测序列$O=(o_1,o_2,...,o_T)$，求对给定观测序列条件概率$P(I|O)$最大的状态序列$I=(i_1,i_2,...,i_T)$。即给定观测序列，求最有可能的状态序列

2. 概率计算算法

2.1 直接计算法

给定模型$\lambda =(A,B,\pi )$和观测序列$O=(o_1,o_2,...,o_T)$，计算观测序列$O$出现的概率$P(O|\lambda )$。

最直接的方法：通过列举所有可能长度为$T$的状态序列$I=(i_1,i_2,...,i_T)$，求各个状态序列$I$与观测序列$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 }_{i_1}{a}_{i_1{i}_2}{a}_{i_2{i}_3}...a_{i_{T-1}{i}_T}$$
2. 对固定的状态序列$I=(i_1,i_2,...,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)$$
3. $O和I$同时出现的联合概率为：
   $$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)$$
4. 对所有可能的状态$I$求和，得到观测序列$O$的概率$P(O|\lambda )$：
   $$P(O|\lambda )=\sum _{I}P(O|I,\lambda )P(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)$$

由于计算量很大，是$O(T{N}^T)$阶，故不可行

2.2 前向算法

前向概率： 给定隐马尔可夫模型$\lambda$，定义导时刻$t$部分观测序列为$o_1,o_2,...,o_t$且状态为$q_i$的概率为前向概率，记作
$${\alpha }_t(i)=P(o_1,o_2,...,o_t,i_t=q_i|\lambda )$$

可以递推求得前向概率${\alpha }_t(i)$及观测序列概率$P(O|\lambda )$

前向算法：

$输入：隐马尔可夫模型\lambda ，观测序列O$

$输出：观测序列概率P(O|\lambda )$

1. 初值：
   $${\alpha }_1(i)=P(o_1,i_1=q_i|\lambda )=P(o_1|i_1=q_i,\lambda )P(i_1=q_i|\lambda )={\pi }_i{b}_i(o_1)$$
2. 递推，对$t=1,2,...,T-1$
   $$\begin{gathered} {\alpha }_{t+1}(i)=P(o_1,o_2,...,o_t,o_{t+1},i_{t+1}=q_i|\lambda )=\sum _{j=1}^{N}P(o_1,o_2,...,o_t,o_{t+1},i_t=q_j,i_{t+1}=q_i|\lambda ) \\ =\sum _{j=1}^{N}P(o_1,o_2,...,o_t,i_t=q_j|\lambda )P(o_{t+1},i_{t+1}=q_i|\lambda ,o_1,o_2,...,o_t,i_t=q_j) \\ =\sum _{j=1}^{N}P(o_1,o_2,...,o_t,i_t=q_j|\lambda )P(o_{t+1}|i_{t+1}=q_i,\lambda )P(i_{t+1}=q_i|\lambda ,i_t=q_j) \\ =\sum _{j=1}^N{\alpha }_t(j)a_{ji}{b}_i(o_{t+1}) \end{gathered}$$
3. 终止
   $$P(O|\lambda )=P(o_1,o_2,...,o_T|\lambda )=\sum _{i=1}^{N}P(o_1,o_2,...,o_T,i_t=q_i|\lambda )=\sum _{i=1}^N{\alpha }_T(i)$$

2.3 后向算法

后向概率： 给定隐马尔可夫模型$\lambda$，定义在时刻$t$状态为$q_i$的条件下，从$t+1$到$T$的部分观测序列$o_{t+1},o_{t+2},...,o_T$的概率为后向概率，记作：
$${\beta }_t(i)=P(o_{t+1},o_{t+2},...,o_T|i_t=q_i,\lambda )$$

可以用递推的方法求后向概率${\beta }_t(i)$及观测序列概率$P(O|\lambda )$

后向算法：

$输入：隐马尔可夫模型\lambda ，观测序列O$

$输出：观测序列概率P(O|\lambda )$

1. 初值
   $${\beta }_T(i)=1$$
2. 递推，对$t=T-1,T-2,...,1$
   $$\beta =\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)$$
3. 终止
   $$P(O|\lambda )=\sum _{i=1}^N{\pi }_i{b}_i(o_1){\beta }_1(i)$$

3. 学习算法

3.1 监督学习方法

假设已给训练数据包含$S$个长度相同的观测序列和对应的状态序列 { ( O 1 , I 1 ) , ( O 2 , I 2 ) , . . . , ( O S , I S ) } \{(O_1,I_1),(O_2,I_2),...,(O_S,I_S)\} {(O1​,I1​),(O2​,I2​),...,(OS​,IS​)}，那么可以利用极大似然估计法来估计隐马尔可夫模型的参数。

1. 转移概率$a_{ij}$的估计

   设样本中时刻$t$处于状态$i$时刻$t+1$转移到状态$j$的频数为$A_{ij}$，那么状态转移概率$a_{ij}$的估计是
   $${\hat{a}}_{ij}=\frac{A_{ij}}{{\sum }_{j=1}^N{A}_{ij}}$$

2. 观测概率$b_j(k)$的估计
   设样本中状态为$j$并观测为$k$的频数为$B_{jk}$，那么状态为$j$观测为$k$的概率$b_j(k)$的估计是：
   $${\hat{b}}_j(k)=\frac{B_{jk}}{{\sum }_{k=1}^M{B}_{jk}}$$

3. 初始状态概率${\pi }_i$的估计${\hat{\pi }}_i$为$S$个样本中初始状态为$q_i$的频率

3.2 Baum-Welch算法

假设给定训练数据只包含$S$个长度为$T$的观测序列 { O 1 , O 2 , . . . , O S } \{O_1,O_2,...,O_S\} {O1​,O2​,...,OS​}而没有对应的状态序列，目标是学习隐马尔可夫模型$\lambda =(A,B,\pi )$的参数。

这事实上是一个含有隐变量的概率模型
$$P(O|\lambda )=\sum _{I}P(O|I,\lambda )P(I|\lambda )$$
它的参数学习可以由$EM$算法实现

1. 确定完全数据的对数似然函数

   所有观测数据写成$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 )$

2. $EM$算法的$E$步：求$Q$函数$Q(\lambda ,\bar{\lambda })$
   $$\begin{gathered} Q(\lambda ,\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 })} \\ =\sum _{I}logP(O,I|\lambda )P(O,I|\bar{\lambda }) \end{gathered}$$
   其中，$\bar{\lambda }$是隐马尔可夫模型参数的当前估计值，$\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$函数可写成
   $$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 })$$

3. $EM$算法的$M$步：极大化$Q$函数，求模型参数$A,B,\pi$
   由于要极大化的参数在$Q$函数中单独地出现在3个项中，所以只需对各项分别极大化。

   • $Q$函数第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 })$$
     注意到${\pi }_i$满足约束条件${\sum }_{i=1}^N{\pi }_i=1$，利用拉格朗日乘子法，写出拉格朗日函数：
     $$\sum _{i=1}^{N}log{\pi }_{i}P(O,i_1=i|\bar{\lambda })+\gamma (\sum _{i+1}^N{\pi }_i-1)$$
     对其求偏导并令结果为0
     $$\frac{\partial }{\partial {\pi }_i}[\sum _{i=1}^{N}log{\pi }_{i}P(O,i_1=i|\bar{\lambda })+\gamma (\sum _{i+1}^N{\pi }_i-1)]=0$$
     得：
     $$P(O,i_1=i|\bar{\lambda })+\gamma {\pi }_i=0$$
     对$i$求和得到$\gamma$
     $$\gamma =-P(O|\bar{\lambda })$$
     所以：
     $${\pi }_i=\frac{P(O,i_1=i|\bar{\lambda })}{P(O|\bar{\lambda })}$$
   • $Q$函数第2项可写成：
     $$\sum _I(\sum _{t=1}^{T-1}log{a}_{i_t{i}_{t+1}})P(O,I|\bar{\lambda })=\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}log{a}_{ij}P(O,i_t=i,i_{t+1}=j|\bar{\lambda })$$
     注意到约束条件${\sum }_{j=1}^N{a}_{ij}=1$，利用拉格朗日乘子法，写出拉格朗日函数：
     $$\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}log{a}_{ij}P(O,i_t=i,i_{t+1}=j|\bar{\lambda })+\gamma (\sum _{j=1}^N{a}_{ij}-1)$$
     对其求偏导并令结果为0
     $$\frac{\partial }{\partial a_{ij}}[\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}log{a}_{ij}P(O,i_t=i,i_{t+1}=j|\bar{\lambda })+\gamma (\sum _{j=1}^N{a}_{ij}-1)]=0$$
     得：
     $$\sum _{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\bar{\lambda })+\gamma a_{ij}=0$$
     对$j$求和得到$\gamma$
     $$\gamma =-\sum _{t=1}^{T-1}P(O,i_t=i|\bar{\lambda })$$
     所以：
     $$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 })}$$
   • $Q$函数第3项可写成：
     $$\sum _I(\sum _{t=1}^{T}log{b}_{i_t}(o_t))P(O,I|\bar{\lambda })=\sum _{j=1}^N\sum _{t=1}^{T}log{b}_j(o_t)P(O,i_t=j|\bar{\lambda })$$
     同样使用拉格朗日乘子法，约束条件为${\sum }_{k=1}^M{b}_j(k)=1$，拉格朗日函数如下：（注意，只有在$o_t=v_k$时，$b_j(o_t)$对$b_j(k)$得偏导才不为0）
     $$\sum _{j=1}^N\sum _{t=1}^{T}log{b}_j(o_t)P(O,i_t=j|\bar{\lambda })+\gamma (\sum _{k=1}^M{b}_j(k)-1)$$
     对其求偏导并令结果为0
     $$\frac{\partial }{\partial b_j(k)}[\sum _{j=1}^N\sum _{t=1}^{T}log{b}_j(o_t)P(O,i_t=j|\bar{\lambda })+\gamma (\sum _{k=1}^M{b}_j(k)-1)]=0$$
     得：
     $$\sum _{t=1}^{T}P(O,i_t=j|\bar{\lambda })I(o_t=v_k)+\gamma b_j(k)=0$$
     对$k$求和得到$\gamma$
     $$\gamma =-\sum _{t=1}^{T}P(O,i_t=j|\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. 预测算法

4.1 近似算法

思想：在每个时刻$t$选择在该时刻最有可能出现得状态$i_t^{\ast }$，从而得到一个状态序列$I^{\ast }=(i_1^{\ast },i_2^{\ast },...,i_T^{\ast })$，将它作为预测结果。

给定隐马尔可夫模型$\lambda$和观测序列$O$，在时刻$t$处于状态$q_i$得概率${\gamma }_t(i)$为：
$$\begin{gathered} {\gamma }_t(i)=P(i_t=q_i|\lambda ,O)=\frac{P(i_t=q_i,O|\lambda )}{P(O|\lambda )}=\frac{P(o_1,o_2,...,o_t,i_t=q_i,o_{t+1},...,o_T|\lambda )}{P(O|\lambda )} \\ =\frac{P(o_1,o_2,...,o_t,i_t=q_i|\lambda )P(o_{t+1},...,o_T|\lambda ,o_1,o_2,...,o_t,i_t=q_i)}{P(O|\lambda )} \\ =\frac{P(o_1,o_2,...,o_t,i_t=q_i|\lambda )P(o_{t+1},...,o_T|\lambda ,i_t=q_i)}{P(O|\lambda )} \end{gathered}$$
由前向概率和后向概率得定义可知
$${\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)}$$
在每一时刻$t$最有可能得状态$i_t^{\ast }$为：
$$i_t^{\ast }=arg\underset{1\le i\le N}{\max }[{\gamma }_t(i)]$$
从而得到状态序列$I^{\ast }=(i_1^{\ast },i_2^{\ast },...,i_T^{\ast })$

近似算法有点是计算简单，其缺点是不能保证预测的状态序列整体是最有可能的状态序列，因为预测的状态序列可能由实际不发生的部分，即上述方法得到的状态序列中可能存在转移概率为0的相邻状态

4.2 维特比算法

维特比算法实际是用动态规划解隐马尔可夫模型预测问题，即用动态规划求概率最大路径（最优路径）。

最优路径具有这样的特性：如果最优路径在时刻$t$通过节点$i_t^{\ast }$，那么这一路径从节点$i_t^{\ast }$到终点$i_T^{\ast }$的部分路径，对于从$i_t^{\ast }$到$i_T^{\ast }$的所有可能的部分路径来说，必须是最优的。因为假如不是这样，那么从$i_t^{\ast }$到$i_T^{\ast }$就有另一条更好的部分路径存在，如果把它和从$i_1^{\ast }$到$i_t^{\ast }$的部分路径连接起来，就会形成一条比原来路径更优的路径，这是矛盾的。

根据这一原理，我们只需从时刻$t=1$开始，递推地计算在时刻$t$状态为$i$地各条部分路径地最大概率，直至得到时刻 $t=T$ 状态为$i$的各条路径的最大概率。时刻 $t=T$ 的最大概率即为最优路径的概率$P^{\ast }$ ，最优路径的终结点$i_T^{\ast }$也同时得到。之后，为了找出最优路径的各个结点，从终结点$i_T^{\ast }$开始，由后向前逐步求得结点 $i_{T-1}^{\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$的递推公式
$${\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)$中概率最大的路径的第$t-1$个节点为
$${\psi }_t(i)=arg\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]$$

维特比算法：

$输入：模型\lambda =(A,B,\pi )和观测O=(o_1,o_2,...,o_T)$

$输出：最优路径I^{\ast }=(i_1^{\ast },i_2^{\ast },...,i_T^{\ast })$

1. 初始化
   $$\begin{gathered} {\delta }_1(i)={\pi }_i{b}_i(o_1),i=1,2,...N \\ {\psi }_1(i)=0,i=1,2,...,N \end{gathered}$$

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(o_t),i=1,2,...,N \\ {\psi }_t(i)=arg\underset{1\le j\le N}{\max }[{\delta }_{t-1}{a}_{ji}],i=1,2,...,N \end{gathered}$$

3. 终止
   $$\begin{gathered} 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)] \end{gathered}$$

4. 最优路径回溯。对$t=T-1,...,1$
   $$i_t^{\ast }={\psi }_{t+1}(i_{t+1}^{\ast })$$

求得最优路径$I^{\ast }=(i_1^{\ast },i_2^{\ast },...,i_T^{\ast })$