隐含马尔科夫(HMM)模型算法推导

隐含马尔科夫(HMM)模型

1.简介

​ 隐含马尔科夫模型是一种关于时序的有向图概率模型，可用于最高气温预测、分词等时序数据建模问题。首先介绍下这一类问题的特点，用$o_t$表示t时刻观测值，记观测值序列为O={o1,o2,…,on}O=\{o_1,o_2,\dots,o_n\}O={o1​,o2​,…,on​}，某一个时刻t的观测值与前n个时刻观测值有一定关系；根据常识来看，某天的最高气温跟前两天的最高气温有一定关系，这就是一个时序数据建模问题。解决这类问题最简单的方法就是用多项式来拟合，以前n个时刻观测值作为输入，以t时刻观测值输出。虽然模型很简单，但其效果并不是很好，因为真实情况比这个复杂得多。隐含马尔科夫模型能够更好地解决这类问题。

​ 隐含马尔科夫模型并不直接使用观测值作为输入，而是假设观测值序列O是由对应的状态序列S产生的，模型最后预测的结果也是状态序列，这样就能够求解出隐变量。状态S就是隐变量，满足马尔科夫假设，即满足任意t时刻的状态$s_t$与其前n个状态相关；并且观测值序列O由状态序列S产生。为了简便，特作出两个假设：

• 齐次马尔科夫假设：任意t时刻状态$s_t$只与其前一个状态$s_{t-1}$有关。
• 观测独立性假设：任意t时刻观测值$o_t$只与该时刻状态$s_t$有关，观测值之间相互独立。

特别说明，这里的状态S一定是离散的，观测值O既可以是离散的也可以是连续的，为了便于说明，采用离散的形式来讨论。

​ 假设Q={q1,q2,…,qM}Q=\{q_1,q_2,\dots,q_M\}Q={q1​,q2​,…,qM​}表示所有可能状态的集合；V={v1,v2,…,vN}V=\{v_1,v_2,\dots,v_N\}V={v1​,v2​,…,vN​}表示所有可能观测值；状态转移概率矩阵A，其中$a_{ij}$表示状态$q_i$转移到状态$q_j$的概率，即$a_{ij}=P(s_t=q_j|s_{t-1}=q_i)$
$$A={\left[\begin{array}{cccc}{a}_{11}& a_{12}& \dots & a_{1M}\\ a_{21}& a_{22}& \dots & a_{2M}\\ ⋮& ⋮& \dots & ⋮\\ a_{M1}& a_{M2}& \dots & a_{MM}\end{array}\right]}_{M\times M}$$
观测概率矩阵B，其中$b_{ij}$表示状态$q_i$产生观测值$v_j$的概率，即$b_{ij}=P(o_t=v_j|s_t=q_i)$
$$B={\left[\begin{array}{cccc}{b}_{11}& b_{12}& \dots & b_{1N}\\ b_{21}& b_{22}& \dots & b_{2N}\\ ⋮& ⋮& \dots & ⋮\\ b_{M1}& b_{M2}& \dots & b_{MN}\end{array}\right]}_{M\times N}$$
时刻1的状态概率向量$\Pi$，其中${\pi }_i$表示时刻1状态$q_i$出现的概率，即${\pi }_i=P(s_1=q_i)$
$$\Pi =[{\pi }_1,{\pi }_2,\dots ,{\pi }_M{]}_{1\times M}$$
于是，隐含马尔科夫模型可记为$\lambda =(A,B,\Pi )$，其中$A,B,\Pi$称为隐含马尔科夫模型三要素。

​ 训练隐含马尔科夫模型需要解决3个基本问题:

• 概率计算问题：给出观测序列O，如何计算似然函数概率$P(O|\lambda )$
• 学习问题：给定观测序列O，如何最大化似然函数求解参数$\lambda$
• 预测问题：求解得到模型后，如何预测得出最优状态序列

2.概率计算问题

​ 给定观测序列O={o1,o2,…,oT}O=\{o_1,o_2,\dots,o_T\}O={o1​,o2​,…,oT​}与参数模型$\lambda =(A,B,\Pi )$有
$$\begin{array}{rl}P(O|\lambda )& =\sum _{S\in Q^T}P(O,S|\lambda )\\ & =\sum _{s_1,s_2,\dots ,s_T}P(o_1,o_2,\dots ,o_T,s_1,s_2,\dots ,s_T|\lambda )\\ & =\sum _{s_1,s_2,\dots ,s_T}{\pi }_{s_1}{b}_{s_1{o}_1}{a}_{s_1{s}_2}{b}_{s_2{o}_2}{o}_2\dots a_{s_{T-1}{s}_T}{b}_{s_T{o}_T}\end{array}$$
​ 其中$s_1,s_2,\dots ,s_T$为任意长度为T的状态序列，并且观测值之间相互独立。显然如果直接计算，那么计算复杂度会非常高，接下来介绍前向后向算法来计算上述概率表达式。

2.1 前向算法

​ 前向算法做出如下定义：
$${\alpha }_t(i)=P(o_1,o_2,\dots ,o_t,s_t=q_i|\lambda ),t=1,2,\dots ,T$$
​ 那么，初值为${\alpha }_1(i)=P(o_1,s_1=q_i|\lambda )={\pi }_i{b}_{q_1{o}_1}$，递推公式为
$$\begin{array}{rl}{\alpha }_{t+1}(i)& =P(o_1,o_2,\dots ,o_t,o_{t+1},s_{t+1}=q_i|\lambda )\\ & =b_{i{o}_{t+1}}\cdot \sum _{j=1}^M{\alpha }_t(j)\cdot a_{ji}\end{array}$$
​ 于是通过迭代的方式可以求得${\alpha }_T(i),i=1,2,\dots ,M$，即求得
$$P(O|\lambda )=\sum _{i=1}^{M}P(o_1,o_2,\dots ,o_T,s_T=q_i|\lambda )=\sum _{i=1}^M{\alpha }_T(i)$$

2.2 后向算法

​ 后向算法做出如下定义：
$${\beta }_t(i)=P(o_{t+1},\dots ,o_T|s_t=q_i,\lambda ),t=T,T-1,\dots ,1$$
​ 那么，初值为${\beta }_T(i)=1$，递推公式为
$$\begin{array}{rl}{\beta }_{t-1}(i)& =P(o_t,o_{t+1},\dots ,o_T|s_{t-1}=q_i,\lambda )\\ & =\sum _{j=1}^M{\alpha }_{ij}\cdot b_{j{o}_t}\cdot {\beta }_t(j)\end{array}$$
​ 于是通过迭代的方式可以求得${\beta }_1(i),i=1,2,\dots ,M$，即求得
$$P(O|\lambda )=\sum _{i=1}^{M}P(o_1,o_2,\dots ,o_T|s_1=q_i,\lambda )=\sum _{i=1}^M{\pi }_i\cdot b_{i{o}_1}\cdot {\beta }_1(i)$$

3.学习问题

​ 隐含马尔科夫模型只能够含有隐变量S，可使用EM算法进行求解，求解后得到的算法称鲍姆-韦尔奇(Baum-Welch)算法。算法推导过程如下：

​ 如前所述，给定观测值序列O和模型参数$\lambda$，对于任意状态序列S有
$$P(O,S|\lambda )=P(o_1,o_2,\dots ,o_T,s_1,s_2,\dots ,s_T|\lambda )={\pi }_{s_1}{b}_{s_1{o}_1}{a}_{s_1{s}_2}{b}_{s_2{o}_2}\dots a_{s_{T-1}{s}_T}{b}_{s_T{o}_T}$$
​ 1.E步：求解Q函数
$$\begin{array}{rl}Q(\lambda ,{\lambda }_i)& =E_S[\log P(O,S|\lambda )|O,{\lambda }_i]\\ & =\sum _S\log [P(O,S|\lambda )]P(S|O,{\lambda }_i)\end{array}$$
​ 其中${\lambda }_i$为当前模型参数估计值，$\lambda$是要极大化似然函数的模型参数值。又因为$P(O|{\lambda }_i)$与$\lambda ,S$都无关，不影响Q函数最大化，所以Q函数也可以写成如下形式：
$$\begin{array}{rl}Q(\lambda ,{\lambda }_i)& =P(O|{\lambda }_i)\sum _S\log [P(O,S|\lambda )]P(S|O,{\lambda }_i)\\ & =\sum _S\log [P(O,S|\lambda )]P(O,S|{\lambda }_i)\end{array}$$
​ 于是将$P(O,S|\lambda )$表达式带入上式可得
$$\begin{array}{rl}Q(\lambda ,{\lambda }_i)& =\sum _S\log [{\pi }_{s_1}{b}_{s_1{o}_1}{a}_{s_1{s}_2}{b}_{s_2{o}_2}\dots a_{s_{T-1}{s}_T}{b}_{s_T{o}_T}]\cdot P(O,S|{\lambda }_i)\\ & =\sum _S\log {\pi }_{s_1}P(O,S|{\lambda }_i)+\sum _S\sum _{t=1}^{T-1}\log a_{s_t{s}_{t+1}}P(O,S|{\lambda }_i)+\sum _S\sum _{t=1}^T\log b_{s_t{o}_t}P(O,S|{\lambda }_i)\\ & =\sum _{j=1}^M\log {\pi }_{j}P(O,s_1=q_j|{\lambda }_i)+\sum _{j=1}^M\sum _{k=1}^M\sum _{t=1}^{T-1}\log a_{jk}P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)\\ & +\sum _{j=1}^M\sum _{t=1}^T\log b_{j{o}_t}P(O,s_t=q_j|{\lambda }_i)\end{array}$$
​ 2.M步：最大化Q函数

​ 接下来，对Q函数中三部分分别求解条件约束极值

​ （1）求解$\pi$

​ 因为$\pi$向量表示时刻1的所有可能状态，所以${\sum }_{j=1}^M{\pi }_j=1$。于是求${\sum }_{j=1}^M\log {\pi }_{j}P(O,s_1=q_j|{\lambda }_i)$在约束条件${\sum }_{j=1}^M{\pi }_j=1$下极值，其拉格朗日函数为
$$L_1({\mu }_1)=\sum _{j=1}^M\log {\pi }_{j}P(O,s_1=q_j|{\lambda }_i)+{\mu }_1(\sum _{j=1}^M{\pi }_j-1)$$
​ 求偏导可得
$$\frac{\partial L_1}{\partial {\pi }_j}=\frac{P(O,s_1=q_j|{\lambda }_i)}{{\pi }_j}+{\mu }_1$$

$$\frac{\partial L_1}{\partial {\mu }_1}=\sum _{j=1}^M{\pi }_j-1$$

​ 令以上两式均为0可得
$${\hat{\pi }}_j=\frac{P(O,s_1=q_j|{\lambda }_i)}{{\sum }_{k=1}^{M}P(O,s_1=q_k|{\lambda }_i)}$$
​ 又由前所述概率计算方法可得
$$P(O,s_t=q_j|{\lambda }_i)={\alpha }_t(j)\cdot {\beta }_t(j)$$
​

​ 所以
$${\hat{\pi }}_j=\frac{{\alpha }_1(j){\beta }_1(j)}{{\sum }_{k=1}^M{\alpha }_1(k){\beta }_1(k)}$$
​ （2）求解A

​ 因为A是状态转移概率矩阵，所以${\sum }_{k=1}^M{a}_{jk}=1$。于是求解${\sum }_{j=1}^M{\sum }_{k=1}^M{\sum }_{t=1}^{T-1}\log a_{jk}P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)$在${\sum }_{k=1}^M{a}_{jk}=1$下的条件极值，其拉格朗日函数为
$$L_2({\mu }_2)=\sum _{j=1}^M\sum _{k=1}^M\sum _{t=1}^{T-1}\log a_{jk}P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)+{\mu }_2(\sum _{k=1}^M{a}_{jk}-1)$$
​ 求偏导数可得
$$\frac{\partial L_2({\mu }_2)}{\partial a_{jk}}=\sum _{t=1}^{T-1}\frac{P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)}{a_{jk}}+{\mu }_2$$

$$\frac{\partial L_2({\mu }_2)}{\partial {\mu }_2}=\sum _{k=1}^M{a}_{jk}-1$$

​ 令偏导数为0可得
$$\begin{array}{rl}{\hat{a}}_{jk}& =\frac{{\sum }_{t=1}^{T-1}P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)}{{\sum }_{k=1}^M{\sum }_{t=1}^{T-1}P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)}\\ & =\frac{{\sum }_{t=1}^{T-1}P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)}{{\sum }_{t=1}^{T-1}P(O,s_t=q_j|{\lambda }_i)}\end{array}$$
​ 又由前述概率计算方法可得
$$P(O,s_t=q_j,s_{t+1}=q_k|{\lambda }_i)={\alpha }_t(j)a_{jk}{b}_{k{o}_{t+1}}{\beta }_{t+1}(k)$$

$$P(O,s_t=q_j|{\lambda }_i)={\alpha }_t(j)\cdot {\beta }_t(j)$$

​ 于是
$${\hat{a}}_{jk}=\frac{{\sum }_{t=1}^{T-1}{\alpha }_t(j)a_{jk}{b}_{k{o}_{t+1}}{\beta }_{t+1}(k)}{{\sum }_{t=1}^{T-1}{\alpha }_t(j)\cdot {\beta }_t(j)}$$
​ （3）计算B

​ 因为B表示观测概率矩阵，所以${\sum }_{t=1}^T{b}_{jk}=1$。于是求解${\sum }_{j=1}^M{\sum }_{t=1}^T\log b_{j{o}_t}P(O,s_t=q_j|{\lambda }_i)$在${\sum }_{t=1}^T{b}_{j{o}_t}=1$约束条件下极值，其拉格朗日函数为
$$L_3({\mu }_3)=\sum _{j=1}^M\sum _{t=1}^T\log b_{j{o}_t}P(O,s_t=q_j|{\lambda }_i)+{\mu }_3(\sum _{t=1}^T{b}_{j{o}_t}-1)$$
​ 求偏导数得
$$\frac{\partial L_3({\mu }_3)}{\partial b_{j{o}_t}}=\frac{P(O,s_t=q_j|{\lambda }_i)}{b_{j{o}_t}}+{\mu }_3$$

$$\frac{\partial L_3({\mu }_3)}{\partial {\mu }_3}=\sum _{t=1}^T{b}_{jk}-1$$

​ 令偏导数为0可得
$$\begin{array}{rl}{\hat{b}}_{jk}& =\frac{{\sum }_{t=1}^{T}P(O,s_t=q_j|{\lambda }_i)I(o_t=v_k)}{{\sum }_{t=1}^{T}P(O,s_t=q_j|{\lambda }_i)}\end{array}$$
​ 其中k表示可能的观测值，取值范围从1到N；I(·)为0-1函数，用于筛选出符合观测值k的数据。

​ 由前述概率计算公式可得
$$P(O,s_t=q_j|{\lambda }_i)={\alpha }_t(j)\cdot {\beta }_t(j)$$
​ 于是
$${\hat{b}}_{jk}=\frac{{\sum }_{t=1}^T{\alpha }_t(j){\beta }_t(j)I(o_t=v_k)}{{\sum }_{t=1}^T{\alpha }_t(j){\beta }_t(j)}$$
​

​ 综上所述，学习问题可用上述求解出来的迭代公式解决，直至参数收敛。

4.预测问题

​ 隐含马尔科夫模型参数$\lambda =(\pi ,A,B)$，从状态序列变化上看，可以构成树形求解空间，目标是求解最优路径问题，考虑使用动态规划算法解决预测问题。维特比算法就是这样设计来解决预测问题，最优路径问题有这样的特性，最优路径的子路径也是子问题的最优路径。

因此，令${\delta }_t(i)={\max }_{s_1,s_2,\dots ,s_{t-1}}P(s_t=q_i,s_{t-1},\dots ,s_1,o_t,o_{t-1},\dots ,o_1|\lambda ),i=1,2,\dots ,M$表示状态序列满足时刻t状态为i的最大概率，由此可得递推公式
$${\delta }_{t+1}(i）=\underset{1\le j\le M}{\max }[{\delta }_t(j)a_{ji}]b_{i{o}_{t+1}}$$
定义时刻t状态为i的所有路径中最优路径中第t-1个结点为${\psi }_t(i)={\max }_{1\le j\le M}[{\delta }_{t-1}(j)a_{ji}],i=1,2,\dots ,M$。

综上所述，通过迭代可求得${\delta }_t(i),i=1,2,\dots ,M$，并且记录每次选择状态，即可得到状态序列

5.参考资料

• 《统计学习方法》- 李航