隐马尔科夫模型

隐马尔科夫模型

1 隐马尔科夫模型概述

1.1 隐马尔科夫模型的定义

隐马尔科夫模型是关于时序的概率模型，描述由隐藏的马尔科夫链生成不可观测的状态序列，再由状态序列生成可观测的观测序列的过程。

设隐马尔科夫模型的状态集合$\mathcal{Y}$为 { s 1 , s 2 , . . . , s N } \{s_1,s_2,...,s_N\} {s1​,s2​,...,sN​}，观测集合$\mathcal{X}$为 { o 1 , o 2 , . . . , o M } \{o_1,o_2,...,o_M\} {o1​,o2​,...,oM​}。一个隐马尔科夫模型的变量可以分为两组，一组是状态变量$Y=(y_1,y_2,...,y_T)$，其中$y_i$表示系统在第$i$时刻的状态，状态是不可观测的；一组是观测变量$X=(x_1,x_2,...,x_T)$，其中$x_i$表示系统在第$i$时刻的观测值。

为确定一个隐马尔科夫模型，需要以下三组参数

• 状态转移概率：模型在状态间转换的概率，表示为矩阵$A=[a_{ij}{]}_{N\times N}$，其中
  $$a_{ij}=P(y_{t+1}=s_j\mid y_t=s_i)$$

• 输出观测概率：模型在当前状态获得各个观测值的概率，表示为矩阵$B=[b_j(k){]}_{N\times M}$，其中
  $$b_j(k)=P(x_t=o_k\mid y_t=s_j)$$

• 初始状态概率：模型在初识时刻各状态出现的概率，通常表示为$\pi =({\pi }_1,{\pi }_2,...,{\pi }_N)$，其中
  $${\pi }_i=P(y_1=s_i)$$

因此，隐马尔科夫模型可用三元组$\lambda =(A,B,\pi )$表示。

1.2 观测序列的生成过程

给定隐马尔科夫模型$\lambda$，观测序列的生成过程描述如下：

$(1)$ 按照初始状态概率选择$y_1$，令$t=1$；

$(2)$ 根据输出观测概率矩阵$B$和$y_t$选择观测值$x_t$；

$(3)$ 根据状态转移矩阵$B$和$y_t$选择状态$y_{t+1}$；

$(4)$ 如果$t<T$，设置$t=t+1$，转到第(2)步，否则停止。

1.3 隐马尔科夫模型的三个基本问题

隐马尔科夫模型有$3$个基本问题。

$(1)$ 概率计算问题。给定模型$\lambda =(A,B,\pi )$和观测序列$X=(x_1,x_2,...,x_T)$，计算在模型$\lambda$下观测$X$出现的概率$P(X\mid \lambda )$

$(2)$ 参数估计问题。给定观测序列$X=(x_1,x_2,...,x_T)$，估计模型$\lambda =(A,B,\pi )$的参数，使得在该模型下观测概率$P(X\mid \lambda )$最大。

$(3)$ 状态预测问题。给定模型$\lambda =(A,B,\pi )$和观测序列$X=(x_1,x_2,...,x_T)$，求最有可能的状态序列，即使得$P(Y\mid X)$最大的状态序列$Y=(y_1,y_2,...,y_T)$。

下面将逐一介绍这三种问题的解法。

2 概率计算

2.1 直接计算法

给定模型$\lambda =(A,B,\pi )$和观测序列$X=(x_1,x_2,...,x_T)$，计算在模型$\lambda$下观测$X$出现的概率$P(X\mid \lambda )$，最直接的方法就是枚举所有可能的状态序列$Y=(y_1,y_2,...,y_T)$，计算联合概率$P(X,Y\mid \lambda )$，然后对所有可能的状态序列求和，计算$P(X\mid \lambda )$。

根据贝叶斯公式
$$P(X,Y\mid \lambda )=P(X\mid Y,\lambda )P(Y\mid \lambda )$$
其中
$$\begin{array}{rl}P(Y\mid \lambda )& ={\pi }_{y_1}{a}_{y_1{y}_2}{a}_{y_2{y}_3},...,a_{y_{T-1}{y}_T}\\ P(X\mid Y,\lambda )& =b_{y_1}(x_1)b_{y_2}(x_2),...,b_{y_T}(x_T)\end{array}$$
因此
$$P(X,Y\mid \lambda )={\pi }_{y_1}{b}_{y_1}(x_1)a_{y_1{y}_2}{b}_{y_2}(x_2),...,a_{y_{T-1}{y}_T}{b}_{y_T}(x_T)$$
然后，通过对所有可能的状态序列求和，就可以得到观测序列$X$的概率$P(X\mid \lambda )$，即
$$\begin{array}{rl}P(X\mid \lambda )& =\sum _{Y}P(X,Y\mid \lambda )\\ & =\sum _{y_1,y_2,...,y_T}{\pi }_{y_1}{b}_{y_1}(x_1)a_{y_1{y}_2}{b}_{y_2}(x_2),...,a_{y_{T-1}{y}_T}{b}_{y_T}(x_T)\end{array}$$
这种方法计算量很大，不可行。下面介绍更有效的算法：前向-后向算法。

2.2 前向算法

前向概率：给定隐马尔科夫模型$\lambda$，定义到时刻$t$观测序列为$x_1,x_2,...,x_t$且状态为$s_i$的概率为前向概率，记作
$${\alpha }_t(i)=P(x_1,x_2,...x_t,y_t=s_i\mid \lambda )$$

观测序列概率的前向算法：

输入：隐马尔科夫模型$\lambda$，观测序列$X$；

输出：观测序列概率$P(X\mid \lambda )$。

$(1)$ 初值
$${\alpha }_1(i)={\pi }_i{b}_i(x_1)i=1,2,...,N$$
$(2)$ 递推，对$t=1,2,...,T-1$
$${\alpha }_{t+1}(i)=(\sum _{j=1}^N{\alpha }_t(j)a_{ji})b_i(x_{t+1})i=1,2,...,N$$
$(3)$ 终止
$$P(X\mid \lambda )=\sum _{j=1}^N{\alpha }_T(j)$$

2.3 后向算法

后向概率：给定隐马尔科夫模型$\lambda$，定义时刻$t$状态为$s_i$的条件下，从$t+1$到$T$观测序列为$x_{t+1},x_{t+2},...,x_T$的概率为后向概率，记作
$${\beta }_t(i)=P(x_{t+1},x_{t+2},...,x_T\mid y_t=s_i,\lambda )$$

观测序列概率的后向算法：

输入：隐马尔科夫模型$\lambda$，观测序列$X$；

输出：观测序列概率$P(X\mid \lambda )$。

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

2.4 一些概率与期望值的计算

利用前向概率和后向概率，可以方便的对一些概率进行计算。

$1.$ 给定参数$\lambda$和观测$X$，在时刻$t$处于状态$s_i$的概率，记为
$${\gamma }_t(i)=P(y_t=s_i\mid X,\lambda )$$
根据贝叶斯公式
$${\gamma }_t(i)=\frac{P(y_t=s_i,X\mid \lambda )}{P(X\mid \lambda )}=\frac{P(y_t=s_i,X\mid \lambda )}{{\sum }_{i=1}^{N}P(y_t=s_i,X\mid \lambda )}$$
根据前向概率与后向概率的定义有
$$P(y_t=s_i\mid X,\lambda )={\alpha }_t(i){\beta }_t(i)$$
因此得到
$${\gamma }_t(i)=\frac{{\alpha }_t(i){\beta }_t(i)}{{\sum }_{j=1}^T{\alpha }_t(j){\beta }_t(j)}$$
$2.$ 给定参数$\lambda$和观测$X$，在时刻$t$处于状态$s_i$且在时刻$t+1$处于状态$s_j$的概率，记为
$${\xi }_t(i,j)=P(y_t=s_i,y_{t+1}=s_j\mid X,\lambda )$$
根据贝叶斯公式
$${\xi }_t(i,j)=\frac{P(y_t=s_i,y_{t+1}=s_j,X\mid \lambda )}{P(X,\mid \lambda )}=\frac{P(y_t=s_i,y_{t+1}=s_j,X\mid \lambda )}{{\sum }_{i=1}^N{\sum }_{j=1}^{N}P(y_t=s_i,y_{t+1}=s_j,X\mid \lambda )}$$
而
$$P(y_t=s_i,y_{t+1}=s_j,X\mid \lambda )={\alpha }_t(i)a_{ij}{b}_j(x_{t+1}){\beta }_{t+1}(j)$$
因此得到
$${\xi }_t(i,j)=\frac{{\alpha }_t(i)a_{ij}{b}_j(x_{t+1}){\beta }_{t+1}(j)}{{\sum }_{i=1}^N{\sum }_{j=1}{\alpha }_t(i)a_{ij}{b}_j(x_{t+1}){\beta }_{t+1}(j)}$$
$3.$ 一些有用的期望：

$(1)$ 在观测$X$下状态$s_i$出现的期望值
$$\sum _{t=1}^T{\gamma }_t(i)$$
$(2)$ 在观测$X$下由状态$s_i$转移的期望值
$$\sum _{t=1}^{T-1}{\gamma }_t(i)$$
$(3)$ 在观测$X$下由状态$s_i$转移到状态$s_j$的期望值
$$\sum _{t=1}^{T-1}\xi (i,j)$$

3 参数估计

通过观测序列或者观测序列和状态序列来估计隐马尔科夫模型参数的算法，称为隐马尔科夫模型的学习算法

3.1 监督学习

假设给定的训练数据包含若干个观测序列和对应的状态序列，那么可以利用极大似然估计来估计参数

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

设样本中时刻$t$处于状态$i$时刻$t+1$处于状态$j$的频数是$A_{ij}$，那么$a_{ij}$的估计是
$$a_{ij}=\frac{A_{ij}}{{\sum }_{j=1}^N{A}_{ij}}$$
$2.$ 观测概率$b_j(k)$的计算

设样本中状态$j$且观测为$k$的频数是$B_{jk}$，那么$b_j(k)$的估计是
$$b_j(k)=\frac{B_{jk}}{{\sum }_{k=1}^M{B}_{jk}}$$
$3.$ 初始状态${\pi }_i$的估计为所有样本中初始状态为$s_i$的概率

3.2 非监督学习

对于非监督学习，此时的训练数据将只包含观测序列，状态序列是隐变量，这时应该使用$\mathrm{B}\mathrm{a}\mathrm{u}\mathrm{m}-\mathrm{W}\mathrm{e}\mathrm{l}\mathrm{c}\mathrm{h}$算法，也即$\mathrm{E}\mathrm{M}$算法。

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

所有观测数据为$X=(x_1,x_2,...,x_T)$，所有隐数据为$Y=(y_1,y_2,...,y_T)$，完全数据是$(X,Y)=(x_1,x_2,...,x_T,y_1,y_2,...,y_T)$，完全数据的对数似然函数是$\log P(X,Y\mid \lambda )$

$2.$ $\mathrm{E}\mathrm{M}$算法的$Q$函数
$$\begin{array}{rl}Q(\lambda ,\bar{\lambda })& =\sum _{Y}P(Y\mid X,\bar{\lambda })P(X,Y\mid \lambda )\\ & =\sum _Y\frac{P(X,Y\mid \bar{\lambda })}{P(X\mid \bar{\lambda })}P(X,Y\mid \lambda )\end{array}$$
其中$\bar{\lambda }$表示当前估计值，由于$P(X\mid \bar{\lambda })$是常数项，略去后可写成
$$Q(\lambda ,\bar{\lambda })=\sum _{Y}P(X,Y\mid \bar{\lambda })P(X,Y\mid \lambda )$$
由于
$$P(X,Y,\mid \lambda )={\pi }_{y_1}{b}_{y_1}(x_1)a_{y_1{y}_2}{b}_{y_2}(x_2)...a_{y_{T-1}{y}_T}{b}_{y_T}(x_T)$$

于是$Q$函数可以写成
$$\begin{array}{rl}Q(\lambda ,\bar{\lambda })& =\sum _Y\log {\pi }_{y_1}P(X,Y\mid \bar{\lambda })+\sum _Y(\sum _{t=1}^{T-1}\log a_{y_t{y}_{t+1}})P(X,I,\mid \bar{\lambda })\\ & +\sum _Y(\sum _{t=1}^T\log b_{y_t}(x_t))P(X,I,\mid \bar{\lambda })\end{array}$$

$3.$ 极大化$Q$函数

由于$Q$函数中的$3$个参数都是单独出现的，因此可以对每一项分别极大化。

$(1)$ 初始状态概率的极大化
$$\sum _Y\log {\pi }_{y_1}P(X,Y\mid \bar{\lambda })=\sum _{i=1}^N\log {\pi }_{i}P(X,y_1=i\mid \bar{\lambda })$$
这里可以理解为合并同类项，由于存在约束条件${\sum }_{i=1}^N{\pi }_i=1$，利用拉格朗日乘子法，拉格朗日函数为
$$L(\pi ,\gamma )=\sum _{i=1}^N\log {\pi }_{i}P(X,y_1=i\mid \bar{\lambda })+\gamma (\sum _{i=1}^N{\pi }_i-1)$$
对${\pi }_i$求偏导
$$\frac{\partial L(\pi ,\lambda )}{\partial {\pi }_i}=\frac{P(X,y_1=i\mid \bar{\lambda })}{{\pi }_i}+\gamma$$
令偏导数等于$0$
$$P(X,y_1=i\mid \bar{\lambda })+\gamma {\pi }_i=0$$
对$i$求和得
$$\gamma =-P(X\mid \bar{\lambda })$$
因此
$${\pi }_i=\frac{P(X,y_1=i\mid \bar{\lambda })}{P(X\mid \bar{\lambda })}$$

$(2)$ 状态转移概率的极大化

$$\sum _Y(\sum _{t=1}^{T-1}\log a_{y_t{y}_{t+1}})P(X,I,\mid \bar{\lambda })=\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}\log a_{ij}P(X,y_t=i,y_{t+1}=j\mid \bar{\lambda })$$
由于存在约束条件${\sum }_{j=1}^N{a}_{ij}=1$，拉格朗日函数为
$$L(a,\gamma )=\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}\log a_{ij}P(X,y_t=i,y_{t+1}=j\mid \bar{\lambda })+\sum _{i=1}^N{\gamma }_i(\sum _{j=1}^N{a}_{ij}-1)$$
对$a_{ij}$求偏导
$$\frac{\partial L(a,\lambda )}{\partial a_{ij}}=\frac{{\sum }_{t=1}^{T-1}P(X,y_t=i,y_{t+1}=j\mid \bar{\lambda })}{a_{ij}}+{\gamma }_i$$
令偏导数等于$0$
$$\sum _{t=1}^{T-1}P(X,y_t=i,y_{t+1}=j\mid \bar{\lambda })+a_{ij}{\gamma }_i=0$$
对所有的$j$求和的
$${\gamma }_i=-\sum _{t=1}^{T-1}P(X,y_t=i\mid \bar{\gamma })$$
因此
$$a_{ij}=\frac{{\sum }_{t=1}^{T-1}P(X,y_t=i,y_{t+1}=j\mid \bar{\lambda })}{{\sum }_{t=1}^{T-1}P(X,y_t=i\mid \bar{\lambda })}$$
$(3)$ 输出观测概率的极大化

$$\sum _Y(\sum _{t=1}^T\log b_{y_t}(x_t))P(X,I,\mid \bar{\lambda })=\sum _{j=1}^N\sum _{t=1}^T\log b_j(x_t)P(X,y_t=j\mid \bar{\lambda })$$
存在约束条件${\sum }_{k=1}^M{b}_j(k)=1$，拉格朗日函数为
$$L(b,\gamma )=\sum _{j=1}^N\sum _{t=1}^T\log b_j(x_t)P(X,y_t=j\mid \bar{\lambda })+\sum _{j=1}^N{\gamma }_j(\sum _{k=1}^M{b}_j(k)-1)$$
对$b_j(k)$求偏导时，只有在$x_t=o_k$时$b_j(x_t)$对$b_{jk}$的偏导数才不为$0$，用$I(x_t=o_k)$表示
$$\frac{\partial L(b,\gamma )}{\partial b_j(k)}=\frac{{\sum }_{t=1}^{T}I(x_t=o_k)P(X,y_t=j\mid \bar{\lambda })}{b_j(k)}+{\gamma }_j$$
令偏导数等于$0$
$$\sum _{t=1}^{T}I(x_t=o_k)P(X,y_t=j\mid \bar{\lambda })+b_j(k){\gamma }_j=0$$
对$k$求和得
$${\gamma }_j=-\sum _{t=1}^{T}P(X,y_t=j\mid \bar{\lambda })$$
因此
$$b_j(k)=\frac{{\sum }_{t=1}^{T}I(x_t=o_k)P(X,y_t=j\mid \bar{\lambda })}{{\sum }_{t=1}^{T}P(X,y_t=j\mid \bar{\lambda })}$$

到此就推导完了$\mathrm{E}\mathrm{M}$算法中的迭代更新公式，我们可以使用在$2.3$中定义的概率值来重写更新公式
$$\begin{array}{rl}{\pi }_i& ={\gamma }_1(i)\\ a_{ij}& =\frac{{\sum }_{t=1}^{T-1}\xi (i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)}\\ b_j(k)& =\frac{{\sum }_{t=1}^T{\gamma }_t(i)I(x_t=o_k)}{{\sum }_{t=1}^T{\gamma }_t(i)}\end{array}$$
$\mathrm{B}\mathrm{a}\mathrm{u}\mathrm{m}-\mathrm{W}\mathrm{e}\mathrm{l}\mathrm{c}\mathrm{h}$算法：

输入：观测数据$X=(x_1,x_2,...,x_T)$；

输出：隐马尔科夫模型参数。

$(1)$ 初始化，对$n=0$，选取$a_{ij}^{(0)},b_j(k{)}^{(0)},{\pi }_i^{(0)}$，得到模型${\lambda }^{(0)}=(A^{(0)},B^{(0)},{\pi }^{(0)})$。

$(2)$ 递推。对$n=1,2,...,$
$$\begin{array}{rl}{a}_{ij}^{(n+1)}& =\frac{{\sum }_{t=1}^{T-1}\xi (i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)}\\ b_j(k{)}^{(n+1)}& =\frac{{\sum }_{t=1}^T{\gamma }_t(i)I(x_t=o_k)}{{\sum }_{t=1}^T{\gamma }_t(i)}\\ {\pi }_i^{(n+1)}& ={\gamma }_1(i)\end{array}$$
右端各值按观测$X=(x_1,x_2,...,x_T)$和模型${\lambda }^{(n)}=(A^{(n)},B^{(n)},{\pi }^{(n)})$计算。

$(3)$ 终止，得到模型参数${\lambda }^{(n+1)}=(A^{(n+1)},B^{(n+1)},{\pi }^{(n+1)})$。

4 状态预测

4.1 近似算法

近似算法的思想是，在每个时刻$t$，选择在该时刻最有可能出现的状态$y_t^{\ast }$作为该时刻的预测，也即
$$y_t^{\ast }=\arg \underset{1\le i\le N}{\max }[{\gamma }_t(i)]$$
从而得到状态序列$Y=(y_1^{\ast },y_2^{\ast },...,y_T^{\ast })$，近似算法是一种贪婪算法，优点是计算简单，缺点是不能保证得到的状态序列是最优状态序列。

4.2 维特比算法

维特比算法利用动态规划法求解隐马尔科夫模型的状态预测问题。维特比算法将一个状态序列当做一条路径，目的是求解概率最大路径。

根据动态规划原理，最优路径具有最优子结构的性质。即如果最优路径在时刻$t$通过节点$y_t^{\ast }$，那么这一路径中从节点$y_t^{\ast }$到节点$y_T^{\ast }$的部分路径，对于从$y_t^{\ast }$到$y_T^{\ast }$的所有可能路径来说，必然是最优的。否则就会得到一条比原来路径更优的路径，这是矛盾的。

因此，只需从时刻$t=1$开始，递推的计算在时刻$t$状态为$i$的各条部分路径的最大概率，直至得到时刻$T$状态为$i$的各条路径的最大概率。时刻$t=T$的最大概率即为最优路径的概率$P^{\ast }$，最优路径的终结点$y_T^{\ast }$也已经得到，之后从$y_T^{\ast }$开始，从后向前逐步求得$y_{T-1}^{\ast },...,y_1^{\ast }$，得到最优路径$Y^{\ast }=(y_1^{\ast },y_2^{\ast },...,y_T^{\ast })$。

定义在时刻$t$状态为$i$的所有单个路径$(y_1,y_2,...,y_t)$中概率最大值为
$${\delta }_t(i)=\underset{y_1,y_2,...,y_{t-1}}{\max }P(y_t=i,y_{t-1},...,y_1,x_t,...,x_1\mid \lambda ),i=1,2,...,N$$
变量$\delta$的递推公式为
$$\begin{array}{rl}{\delta }_{t+1}(i)& =\underset{y_1,y_2,...,y_t}{\max }P(y_{t+1}=i,y_t,...,y_1,x_{t+1},...,x_1\mid \lambda )\\ & =\underset{1\le j\le N}{\max }[{\delta }_t(j)a_{ji}]b_i(x_{t+1}),i=1,2,...,N,t=1,2,...,T-1\end{array}$$
定义在时刻$t$状态为$i$的所有单个路径$(y_1,y_2,...,y_t)$中概率最大的路径的第$t-1$个节点为
$${\psi }_t(i)=\arg \underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}],i=1,2,...,N$$
有了这两个变量，就可以介绍维特比算法。

维特比算法：

输入：隐马尔科夫模型$\lambda =(A,B,\pi )$和观测$X=(x_1,x_2,...,x_T)$；

输出：最优状态序列$Y^{\ast }=(y_1^{\ast },y_2^{\ast },...,y_T^{\ast })$。

$(1)$ 初始化
$$\begin{array}{rl}{\delta }_1(i)& ={\pi }_i{b}_i(x_1),i=1,2,...,N\\ {\psi }_1(i)& =0,i=1,2,...,N\end{array}$$
$(2)$ 递推。对$t=2,3,...,T$
$$\begin{array}{rl}{\delta }_t(i)& =\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]b_i(x_t),i=1,2,...,N\\ {\psi }_t(i)& =\arg \underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}],i=1,2,...,N\end{array}$$
$(3)$ 终止
$$\begin{array}{rl}{P}^{\ast }& =\underset{1\le i\le N}{\max }[{\delta }_T(i)]\\ y_T^{\ast }& =\arg \underset{1\le i\le N}{\max }[{\delta }_T(i)]\end{array}$$
$(4)$ 回溯。对$t=T-1,T-2,...,1$
$$y_t^{\ast }={\psi }_{t+1}(y_{t+1}^{\ast })$$

求得最优状态序列$Y^{\ast }=(y_1^{\ast },y_2^{\ast },...,y_T^{\ast })$。