【机器学习算法推导】隐马尔科夫模型HMM及相关算法

1.隐马尔可夫模型HMM

1.1 简介

  隐马尔科夫模型是与概率论密切相关的一个模型，隐马尔科夫模型（Hidden Markov Model，以下简称HMM）在语音识别，自然语言处理等方面都有所应用。它被用来解决这样的问题：有两组序列Q和V，Q无法被观察到，称为隐藏序列。但是与V有着某种概率上的关系，而V是可以直接被观察到的，称为状态序列。

1.2 定义

  令$Q$={$q_1,q_2,...,q_N$}表示所有隐藏状态的集合，$V$={$v_1,v_2,...,v_M$}表示所有可观察状态的集合，N表示隐藏状态数，M表示可观察状态数。对于一个长度为T的序列来说，S表示隐藏状态序列，O表示观察序列，即$S$={$s_1,s_2,...,s_T$},$O$={$o_1,o_2,...o_T$}，其中$s_t\in Q$,$o_t\in V$。
  一个HMM模型需要由三部分组成，分别是隐藏状态的初始概率分布$\Pi$,状态转移概率矩阵$A$，以及观测状态概率矩阵$B$。
  HMM模型假定了一个隐藏状态只跟前一个隐藏状态有关。令时刻t的状态$s_t=q_i$，时刻t+1的状态为$s_{t+1}=q_j$即从时刻t到时刻t+1的HMM状态转移概率$a_{ij}=P(s_{t+1}=q_j|s_t=q_i)$，所有的$a_{ij}$组成了状态转移概率矩阵A。
  HMM模型同时也假定了一个被观测到的状态，只跟当前的隐藏状态有关。令时刻t的隐藏状态$s_t=q_j$，对应的观察状态$o_t=v_k$即时刻t的观察状态$v_k$在隐藏状态$q_j$下的生成概率为$b_j(k)=P(o_t=v_k|s_t=q_j)$，所有的$b_j(k)组成了观测状态概率矩阵B。

1.3 实例

  用概率论中经典的“取球模型”来类比的话，假设我们有3个盒子，每个盒子里都有红色和白色两种球，这三个盒子里球的数量分别是：
$$\begin{array}{cccc}& \text{盒子1}& \text{盒子2}& \text{盒子3}\\ 白球数& 5& 4& 7\\ 红球数& 5& 6& 3\end{array}$$
  根据盒子中球的分布，可以得到观测状态概率矩阵B：
$$A=\left(\begin{array}{cc}0.5& 0.5\\ 0.4& 0.6\\ 0.7& 0.3\end{array}\right)$$
  按照下面的方法从盒子里抽球，开始的时候，从第一个盒子抽球的概率是0.2，从第二个盒子抽球的概率是0.4，从第三个盒子抽球的概率是0.4。则初始概率分布$$\Pi =(0.2,0.4,0.4{)}^T$$  以这个概率抽一次球后，将球放回。然后从当前盒子转移到下一个盒子进行抽球。规则是：如果当前抽球的盒子是第一个盒子，则以0.5的概率仍然留在第一个盒子继续抽球，以0.2的概率去第二个盒子抽球，以0.3的概率去第三个盒子抽球。如果当前抽球的盒子是第二个盒子，则以0.5的概率仍然留在第二个盒子继续抽球，以0.3的概率去第一个盒子抽球，以0.2的概率去第三个盒子抽球。如果当前抽球的盒子是第三个盒子，则以0.5的概率仍然留在第三个盒子继续抽球，以0.2的概率去第一个盒子抽球，以0.3的概率去第二个盒子抽球。说了这么多，最终得到的就是状态转移概率矩阵A：
$$A=\left(\begin{array}{ccc}0.5& 0.2& 0.3\\ 0.3& 0.5& 0.2\\ 0.2& 0.3& 0.5\end{array}\right)$$
  其中
  观察集合是: V = { 红 ， 白 } ， M = 2 V=\{红，白\}，M=2 V={红，白}，M=2
  状态集合是： Q = { 盒 子 1 ， 盒 子 2 ， 盒 子 3 } ， N = 3 Q=\{盒子1，盒子2，盒子3\}，N=3 Q={盒子1，盒子2，盒子3}，N=3

1.4 生成观测序列

  假设我们存在一个HMM模型$\lambda =(A,B,\Pi )$，且已知观测序列的长度T，需要生成一个观测序列 O = { o 1 , o 2 , . . . o T } O=\{o_1,o_2,...o_T\} O={o1​,o2​,...oT​}，其生成过程如下：

• 根据初始概率分布$\Pi$生成隐藏状态$s_1$
• 根据当前隐藏状态$s_t$以及观察状态概率$b_{s_t}(k)$生成观察状态$o_t$。
• 根据当前隐藏状态$s_t$以及状态转移概率$a_{s_t{s}_{t+1}}$生成隐藏状态$s_{t+1}$
• 重复T次步骤2和3，得到所有的$o_t$组成观察序列 O = { o 1 , o 2 , . . . o T } O=\{o_1,o_2,...o_T\} O={o1​,o2​,...oT​}

1.5 HMM的三个基本问题

• 评估观察序列概率。即给定模型λ=(A,B,Π)和观测序列O={o1,o2,…oT}，计算在模型λ下观测序列O出现的概率P(O|λ)。这个问题的求解需要用到基于动态规划的前向后向算法。
• 模型参数学习问题。即给定观测序列O={o1,o2,…oT}，估计模型λ=(A,B,Π)的参数，使该模型下观测序列的条件概率P(O|λ)最大。这个问题的求解需要用到基于EM算法的鲍姆-韦尔奇算法。
• 预测问题，也称为解码问题。即给定模型λ=(A,B,Π)和观测序列O={o1,o2,…oT}，求给定观测序列条件下，最可能出现的对应的状态序列，这个问题的求解需要用到基于动态规划的维特比算法。

2.前向后向算法

  问题描述：给定模型λ=(A,B,Π)和观测序列O={o1,o2,…oT}，计算在模型λ下观测序列O出现的概率P(O|λ)。

2.1 暴力求法

  任意一个隐藏序列 S = { s 1 , s 2 , . . . , s T } S=\{s_1,s_2,...,s_T\} S={s1​,s2​,...,sT​}出现的概率是：
$$P(S|\lambda )={\Pi }_{s_1}{a}_{s_1{s}_2}{a}_{s_2{s}_3}...a_{s_{T-1}{s}_T}$$
  在已知S的情况下，我们要求的观察序列 O = { o 1 , o 2 , . . . o T } O=\{o_1,o_2,...o_T\} O={o1​,o2​,...oT​}出现的概率为
$$P(O|S,\lambda )=b_{s_1}(o_1)b_{s_2}(o_2)...b_{s_T}(o_T)$$
  则S和O联合出现的概率为
$$P(O,S|\lambda )=P(S|\lambda )P(O|S,\lambda )={\Pi }_{s_1}{a}_{s_1{s}_2}{a}_{s_2{s}_3}...a_{s_{T-1}{s}_T}{b}_{s_1}(o_1)b_{s_2}(o_2)...b_{s_T}(o_T)$$
  然后求边缘概率分布，即可得到观测序列O在模型λ下出现的条件概率P(O|λ)：
$$P(O|\lambda )=\sum _{S}P(O,S|\lambda )$$
  然而暴力算法存在一个大问题，当隐藏状态较多时，预测状态有$N^T$种组合，上述的算法太过耗时。我们一般使用前向算法或者后向算法来求解这个问题。

2.2 前向算法

  算法比较好的同学也许已经发现，上述的过程可以使用动态规划来完成。动态规划的本质是递推，也就是说，我们需要将问题分解为子问题，从子问题的最优解中一步一步，根据递推式，推出整个问题的最优解。
  这里给出动态规划中的子问题，也就是前向概率的定义。假设在时刻t，隐藏状态为$q_i$，观测状态为$o_1,o_2,...o_t$的概率为前向概率$f_t(i)=P(o_1,o_2,...,o_t|s_t=q_i,\lambda )$。
  考虑时刻t+1的情况，因为HMM的隐藏状态只与前一个隐藏状态有关，既然这样，我们就考虑所有前一个状态，让它们都向t+1时刻的某个状态转移，就可以得到相应的概率，这样的话，我们就得到了前向概率的递推式:
$$f_{t+1}(i)=b_i(o_{t+1})\sum _{j=1}^N{f}_t(j)a_{ji}$$
  有了递推式，我们就可以推算出最终的目标，得到观测序列O在模型λ下出现的条件概率P(O|λ)。前向算法过程如下：

1. 初始化时刻1各个状态的前向概率$f_1(i)={\Pi }_i{b}_i(o_1),i=1,2,3,...,N$
2. 使用递推式$f_{t+1}(i)=b_i(o_{t+1}){\sum }_{j=1}^N{f}_t(j)a_{ji}$递推时刻2,3,…,T的前向概率。
3. 计算观测序列O在模型λ下出现的条件概率$P(O|\lambda )=\sum _{i=1}^N{f}_T(i)$

2.3 后向算法

  后向算法和前向算法类似，不同的地方在于后向算法是从后往前递推。假设在时刻t，隐藏状态为$q_i$，观测状态为$o_{t+1},o_{t+2},...o_T$的概率为后向概率$g_t(i)=P(o_{t+1},o_{t+2},...o_T|s_t=q_i,\lambda )$。
  考虑时刻t+1的情况，假定我们已经知道当前所有隐藏状态的后向概率$g_{t+1}(j)$，在时刻t的某个状态i下，由i到达观测状态为$o_{t+1},o_{t+2},...,o_T$的概率为$a_{ij}{b}_j(o_{t+1})g_{t+1}(j)$，由于这这样的j可以有多个，我们将时刻t+1下所有的状态j考虑进来，就得到了后向概率的递推式:
$$g_t(i)=\sum _{j=1}^N{a}_{ij}{b}_j(o_{t+1})g_{t+1}(j)$$
  有了递推式，我们就可以推算出最终的目标，得到观测序列O在模型λ下出现的条件概率P(O|λ)。后向算法过程如下：

1. 初始化时刻T各个状态的前向概率$g_T(i)=1,i=1,2,3,...,N$
2. 使用递推式$g_t(i)={\sum }_{j=1}^N{a}_{ij}{b}_j(o_{t+1})g_{t+1}(j)$递推时刻T-1,T-2,…,1的后向概率。
3. 计算观测序列O在模型λ下出现的条件概率$P(O|\lambda )=\sum _{i=1}^N{\Pi }_i{b}_i(o_1)g_1(i)$

3.鲍勃-韦尔奇算法

  给定若干训练样本O={o1,o2,…oT}，估计模型λ=(A,B,Π)的参数，使该模型下观测序列的条件概率P(O|λ)最大。这个问题的求解需要用到基于EM算法的鲍姆-韦尔奇算法。
  联合分布$P(O,S|\lambda )$的表达式如下：
$$P(O,S|\lambda )=\prod _{d=1}^D{\Pi }_{i_1^{(d)}}{b}_{i_1^{(d)}}(o_1^{(d)})a_{i_1^{(d)}{i}_2^{(d)}}{b}_{i_2^{(d)}}(o_2^{(d)})...a_{i_{T-1}^{(d)}{i}_T^{(d)}}{b}_{i_T^{(d)}}(o_T^{(d)})$$

  E-Step的期望表达式为：
$$L(\lambda ,\bar{\lambda })=\sum _{S}P(S|O,\bar{\lambda })\log P(O,S|\lambda )$$
  因为$P(S|O,\bar{\lambda })=\frac{P(O,S|\bar{\lambda })}{P(O|\bar{\lambda })}$，分母是常数，所以在M-Step中我们要最大化的表达式为：
$$\bar{\lambda }=\arg \underset{\lambda }{\max }\sum _{d=1}^D\sum _{S}P(O,S|\bar{\lambda })(\log {\pi }_{s_1}+\sum _{t=1}^{T-1}\log a_{s_t,s_{t+1}}+\sum _{t=1}^T{b}_{s_t}(o_t))............(1)$$
  式子1分别对A，B，$\Pi$求导，得到相应的迭代公式。在开始推导前，我们先定义一些变量。
  由前向概率和后向概率定义可知：
$$\begin{gathered} P(s_t=q_i,O|\lambda )=f_t(i)g_t(i) \\ P(s_t=q_i,s_{t+1}=q_j,O|\lambda )=f_t(i)a_{ij}{b}_j(o_{t+1})g_{t+1}(j) \end{gathered}$$
  给定模型λ和观测序列O,在时刻t处于状态$q_i$的概率记为:
$${\gamma }_t(i)=P(s_t=q_i|O,\lambda )=\frac{P(s_t=q_i,O|\lambda )}{P(O|\lambda )}=\frac{f_t(i)g_t(i)}{\sum _{j=1}^N{f}_t(j)g_t(j)}$$
  给定模型λ和观测序列O,在时刻t处于状态$q_i$，且时刻t+1处于状态$q_j$的概率记为:
$${\xi }_t(i,j)=P(s_t=q_i,s_{t+1}=q_j|O,\lambda )=\frac{P(s_t=q_i,s_{t+1}=q_j,O|\lambda )}{P(O|\lambda )}=\frac{f_t(i)a_{ij}{b}_j(o_{t+1})g_{t+1}(j)}{\sum _{r=1}^N\sum _{s=1}^N{f}_t(r)a_{rs}{b}_s(o_{t+1})g_{t+1}(s)}$$

3.1 对${\pi }_i$求导

  式子1对${\Pi }_i$求偏导，式子中后面部分的$\sum _{t=1}^{T-1}\log a_{s_t,s_{t+1}}+\sum _{t=1}^T{b}_{s_t}(o_t)$为常数，直接无视，求导的式子变成
$$\bar{\lambda }=\arg \underset{\lambda }{\max }\sum _{d=1}^D\sum _{S}P(O,S|\bar{\lambda })\log {\pi }_{s_1}............(2)$$
  由于${\Pi }_i$是初始概率分布，所以还满足$\sum _{i=1}^N{\pi }_i=1$。因此需要用到拉格朗日乘子法来求。式子2的拉格朗日函数为:
$$L({\pi }_i)=\sum _{d=1}^D\sum _{S}P(O,S|\bar{\lambda })\log {\pi }_{s_1}+\gamma (\sum _{i=1}^N{\pi }_i-1)............(3)$$
  对${\pi }_i$进行求导，并令偏导数为0，得到
$$\sum _{d=1}^{D}P(O,s_1^{(d)}=i|\bar{\lambda })+\gamma {\pi }_i=0.........(4)$$
  由于i可以取值范围为1到N，因此我们得到N个方程，将所有方程相加后得到
$$\sum _{d=1}^{D}P(O|\bar{\lambda })+{\pi }_i=0..........(5)$$
  联立4和5得到
$$\begin{array}{rl}{\pi }_i& =\frac{\sum _{d=1}^{D}P(O,s_1^{(d)}=i|\bar{\lambda })}{\sum _{d=1}^{D}P(O|\bar{\lambda })}\\ & =\frac{\sum _{d=1}^{D}P(O,s_1^{(d)}=i|\bar{\lambda })}{DP(O|\bar{\lambda })}\\ & =\frac{\sum _{d=1}^{D}P(s_1^{(d)}=i|O^{(d)},\bar{\lambda })}{D}\\ & =\frac{\sum _{d=1}^D{\gamma }_i^{(d)}}{D}..........(6)\end{array}$$
  至此就得到了$\pi$的迭代公式。

3.2 对$a_{ij}$求导

  式子1对$a_{ij}$求偏导，直接无视常数部分，求导的式子变成
$$\sum _{d=1}^D\sum _S\sum _{t=1}^{T-1}P(O,S|\bar{\lambda })\log a_{s_t,s_{t+1}}=\sum _{d=1}^D\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}P(O,s_t^{(d)}=i,s_{t+1}^{(d)}=j|\bar{\lambda })\log a_{ij}............(7)$$
  同样存在$\sum _{i=1}^N{a}_{ij}=1$的限制，使用拉格朗日乘子法，最终得到$a_{ij}$的表达式为
$$a_{ij}=\frac{{\sum }_{d=1}^D{\sum }_{t=1}^{T-1}P(O^{(d)},s_t^{(d)}=i,s_{t+1}^{(d)}=j|\bar{\lambda })}{{\sum }_{d=1}^D{\sum }_{t=1}^{T-1}P(O^{(d)},s_t^{(d)}=i|\bar{\lambda })}=\frac{{\sum }_{d=1}^D{\sum }_{t=1}^{T-1}{\xi }_t^{(d)}(i,j)}{{\sum }_{d=1}^D{\sum }_{t=1}^{T-1}{\gamma }_t^{(d)}(i)}$$

3.3 对$b_j$求导

  同上述的步骤类似，$b_j(k)$的迭代公式为：
$$\begin{array}{rl}{b}_j(k)& =\frac{\sum _{d=1}^D\sum _{t=1,o_t^{(d)}=v_k}^{N}P(O,s_t^{(d)}=j|\bar{\lambda })}{\sum _{d=1}^D\sum _{t=1}^{N}P(O,s_t^{(d)}=j|\bar{\lambda })}\\ & =\frac{\sum _{d=1}^D\sum _{t=1,o_t^{(d)}=v_k}^N{\gamma }_t^{(d)}(j)}{\sum _{d=1}^D\sum _{t=1}^N{\gamma }_t^{(d)}(j)}\end{array}$$

3.4 算法总结

1. 随机初始化HMM的参数$\lambda$
2. 对所有的观察序列,计算${\gamma }_t^{(d)}(i),{\xi }_t^{(d)}(i,j)$,t=1,2,…,T
3. 根据迭代公式更新模型参数
4. 如果参数收敛，则算法结束，否则重复2和3。

4. 维特比算法

  问题描述：在已知模型参数$\lambda$和观察序列O的情况下，要找出最有可能的隐藏状态。可以使用维特比算法来解决，它是一种基于动态规划的算法。在维特比算法中，定义了两个局部状态。
  定义${\delta }_t(i)$表示在时刻t，状态为i的情况下，到目前为止所有可能的转移路径中概率最大的情况。递推式如下：
$${\delta }_{t+1}(i)=\underset{1\le j\le N}{\max }{\delta }_t(j)a_{ji}{b}_i(o_{t+1})$$
  ${\psi }_t(i)$记录在时刻t，状态为i的情况下，概率最大的转移路径的第t-1个节点的状态。递推式如下：
$${\psi }_t(i)=\underset{1\le j\le N}{\max }{\delta }_{t-1}{a}_{ji}$$
  维特比算法总结如下:

1. 初始化${\delta }_1(i)={\pi }_i{b}_i(o_1),{\psi }_1(i)=0,i=1,2,...,N$
2. 利用递推式${\delta }_{t+1}(i)={\max }_{1\le j\le N}{\delta }_t(j)a_{ji}{b}_i(o_{t+1})$和${\psi }_t(i)={\max }_{1\le j\le N}{\delta }_{t-1}{a}_{ji}$求解时刻2,3,…,T
3. 选择最大的${\delta }_T(i)$，利用${\psi }_t(i)$回溯得到前T-1个状态，最终得到最有可能的序列，并且最大的${\delta }_T(i)$为这个序列出现的概率。