《统计学习方法》——隐马尔可夫模型(下)

学习算法

HMM的学习，在有观测序列的情况下，根据训练数据是否包含状态序列，可以分别由监督学习算法和无监督学习算法实现。

监督学习算法

监督学习算法就比较简单，基于已有的数据利用极大似然估计法来估计隐马尔可夫模型的参数。分为以下几步：

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

设样本中时刻$t$处于状态$i$，时刻$t+1$转到到状态$j$的频数为$A_{ij}$，那么状态转移概率$a_{ij}$的估计是
$${\hat{a}}_{ij}=\frac{A_{ij}}{{\sum }_j^N{A}_{ij}},i=1,2,\cdots ,N;j=1,2,\cdots ,N\text{\hspace{1em}(10.30)}$$
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}},j=1,2,\cdots ,N;k=1,2,\cdots ,M\text{\hspace{1em}(10.31)}$$
3.初始状态概率${\pi }_i$的估计${\hat{\pi }}_i$为$S$个样本中初始状态为$q_i$的频率

一般没有这么多标注的训练数据，因此通常采用的是无监督学习方法，下面介绍一种。

Baum-Welch算法

假设给定训练数据只包含$S$个长度为$T$的观测序列 { O 1 , ⋯ , O S } \{O_1,\cdots,O_S\} {O1​,⋯,OS​}而没有对应的状态序列，目标是学习隐马尔可夫模型$\lambda =(A,B,\pi )$的参数。我们将观测序列数据看作观测数据$O$，状态序列数据看作不可观测的隐数据$I$，那么隐马尔可夫模型事实上是一个含有隐变量的概率模型
$$P(O|\lambda )=\sum _{I}P(O,I|\lambda )=\sum _{I}P(O|I,\lambda )P(I|\lambda )\text{\hspace{1em}(10.32)}$$
它的参数学习可以由EM算法实现。

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

所有观测数据写成$O=(o_1,\cdots ,o_T)$，所有隐数据写成$I=(i_1,\cdots ,i_T)$，完全数据是$(O,I)=(o_1,o_2,\cdots ,o_T,i_1,i_2,\cdots ,i_T)$。完全数据的对数似然函数是$\log P(O,I|\lambda )$。

2.EM算法的E步：求$Q$函数$Q(\lambda ,\bar{\lambda })$

按照$Q$函数的定义，
$$\begin{array}{rlr}Q(\lambda ,\bar{\lambda })& =E_I[\log P(O,I|\lambda )|O,\bar{\lambda }]& \\ & =\sum _I\log P(O,I|\lambda )P(I|O,\bar{\lambda })& \\ & =\sum _I\log P(O,I|\lambda )\frac{P(O,I|\bar{\lambda })}{P(O|\bar{\lambda })}& \\ & =\frac{1}{P(O|\bar{\lambda })}\sum _I\log P(O,I|\lambda )P(O,I|\bar{\lambda })& 与I无关可以提出去\end{array}$$
其中$\bar{\lambda }$是隐马尔可夫模型参数的当前估计值，$\lambda$是要极大化的隐马尔可夫模型参数。

$\bar{\lambda }$是一个常量，上式中$P(O|\bar{\lambda })$对于$\lambda$而言是一个常数项，省去该项就得到了：
$$Q(\lambda ,\bar{\lambda })=\sum _I\log P(O,I|\lambda )P(O,I|\bar{\lambda })\text{\hspace{1em}(10.33)}$$
而根据式$(10.12)$：
$$P(O,I|\lambda )={\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1,i_2}{b}_{i_2}(o_2)\cdots a_{i_{T-1},i_T}{b}_{i_T}(o_T)$$
我们要求的就是$\pi ,A,B$，因此对上式按这三个参数分开。

于是函数$Q(\lambda ,\bar{\lambda })$可以写成：

$$\begin{array}{rlr}Q(\lambda ,\bar{\lambda })& =\sum _I\log P(O,I|\lambda )P(O,I|\bar{\lambda })& \\ & =\sum _I\log [{\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1,i_2}{b}_{i_2}(o_2)\cdots a_{i_{T-1},i_T}{b}_{i_T}(o_T)]P(O,I|\bar{\lambda })& \\ & =\sum _I[\log ({\pi }_{i_1})+\log (a_{i_1,i_2}{a}_{i_2,i_3}{a}_{i_{T-1},i_T})+\log (b_{i_1}(o_1)b_{i_2}(o_2)\cdots b_{i_T}(o_T)]P(O,I|\bar{\lambda })& 根据三个参数分开\\ & =\sum _I\log ({\pi }_{i_1})P(O,I|\bar{\lambda })+\sum _I\left(\sum _{t=1}^{T-1}\log a_{i_t,i_{t+1}}\right)P(O,I|\bar{\lambda })+\sum _I\left(\sum _{t=1}^T{b}_{i_t}(o_t)\right)P(O,I|\bar{\lambda })& \end{array}\text{\hspace{1em}(10.34)}$$
3.EM算法的M步： 极大化$Q$函数$Q(\lambda ,\bar{\lambda })$ 求模型参数$A,B,\pi$

我们上一步已经把要极大化的参数单独地分开成三个项，所以只需要对各项分别极大化。

(1)式$(10.34)$第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 \left(\sum _{i=1}^N{\pi }_i-1\right)$$
上式对${\pi }_i$求偏导并令结果为0：
$$\frac{\partial }{\partial {\pi }_i}\left[\sum _{i=1}^N\log ({\pi }_i)P(O,i_1=i|\bar{\lambda })+\gamma \left(\sum _{i=1}^N{\pi }_i-1\right)\right]=0\text{\hspace{1em}(10.35)}$$
得
$$\frac{1}{{\pi }_i}P(O,i_1=i|\bar{\lambda })+\gamma =0\Rightarrow P(O,i_1=i|\bar{\lambda })+{\pi }_i\gamma =0$$
对$i$求和得到$\gamma$：
$$\begin{array}{rl}{\pi }_i\gamma & =-P(O,i_1=i|\bar{\lambda })\\ \sum _i{\pi }_i\gamma & =\sum _i-P(O,i_1=i|\bar{\lambda })\\ \gamma & =-P(O|\bar{\lambda })\end{array}$$
代入式$(10.35)$得
$$\begin{array}{rl}P& (O,i_1=i|\bar{\lambda })+{\pi }_i\gamma =0\\ {\pi }_i& =-\frac{P(O,i_1=i|\bar{\lambda })}{\gamma }\\ {\pi }_i& =\frac{P(O,i_1=i|\bar{\lambda })}{P(O|\bar{\lambda })}\end{array}\text{\hspace{1em}(10.36)}$$
(2) 式$(10.34)$的第2项可以写成
$$\sum _I\left(\sum _{t=1}^{T-1}\log a_{i_t,i_{t+1}}\right)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 })$$
$a_{ij}$满足约束条件${\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 \left(\sum _{i=1}^N{a}_{ij}-1\right)$$
上式对$a_{ij}$求偏导并令结果为0得：
$$\begin{gathered} \frac{1}{a_{ij}}\sum _{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\bar{\lambda })+\gamma =0 \\ {a}_{ij}\gamma =-\sum _{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\bar{\lambda }) \end{gathered}$$
两边对$j$求和：
$$\begin{gathered} \sum _j{a}_{ij}\gamma =-\sum _j\sum _{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\bar{\lambda }) \\ \gamma =-\sum _{t=1}^{T-1}P(O,i_t=i|\bar{\lambda }) \end{gathered}$$
代入得
$$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 })}\text{\hspace{1em}(10.37)}$$
(3) 式$(10.34)$的第3项为
$$\sum _I\left(\sum _{t=1}^T{b}_{i_t}(o_t)\right)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$，回顾一下，有$M$个观测变量；$b_j(k)$表示状态为$q_j$的情况下观测为$v_k$的概率。

注意，只有在$o_t=v_k$时$b_j(o_t)$对$b_j(k)$的偏导数才不为0，以$I(o_t=v_k)$表示。

利用拉格朗日乘子法，写出拉格朗日函数：
$$\sum _{j=1}^N\sum _{t=1}^T\log b_j(o_t)P(O,i_t=j|\bar{\lambda })+\gamma \left(\sum _{k=1}^M{b}_j(k)-1\right)$$
上式对$b_j(k)$求偏导并令结果为0得：
$$\begin{gathered} \frac{1}{b_j(k)}\sum _{t=1}^{T}P(O,i_t=j|\bar{\lambda })I(o_t=v_k)+\gamma =0 \\ {b}_j(k)\gamma =-\sum _{t=1}^{T}P(O,i_t=j|\bar{\lambda })I(o_t=v_k) \end{gathered}$$
两边对$k$求和，注意右边本来通过指示函数$I(o_t=v_k)$限制$o_t=v_k$，对$k$求和的话，相当于整个指示函数没了。
$$\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 })}\text{\hspace{1em}(10.38)}$$

Baum-Welch模型参数估计公式

将式$(10.36)$~式$(10.38)$中的各概率分别用式$(10.23),(10.25)$中的${\gamma }_t(i),{\xi }_t(i,j)$表示，则可以写成：
$$\begin{array}{rl}{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 })}\\ & =\frac{{\sum }_{t=1}^{T-1}P(i_t=i,i_{t+1}=j|O,\bar{\lambda })P(O|\bar{\lambda })}{{\sum }_{t=1}^{T-1}P(O,i_t=i|\bar{\lambda })}\\ & =\frac{{\sum }_{t=1}^{T-1}P(i_t=i,i_{t+1}=j|O,\bar{\lambda })}{{\sum }_{t=1}^{T-1}P(O,i_t=i|\bar{\lambda })/P(O|\bar{\lambda })}\\ & =\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)}\end{array}\text{\hspace{1em}(10.39)}$$

$$\begin{array}{rl}{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 })}\\ & =\frac{{\sum }_{t=1}^{T}P(i_t=j|O,\bar{\lambda })I(o_t=v_k)}{{\sum }_{t=1}^{T}P(i_t=j|O,\bar{\lambda })}\\ & =\frac{{\sum }_{t=1,o_t=v_k}^T{\gamma }_t(j)}{{\sum }_{t=1}^T{\gamma }_t(j)}\end{array}\text{\hspace{1em}(10.40)}$$

$${\pi }_i={\gamma }_1(i)\text{\hspace{1em}(10.41)}$$

然后就可以基于这种形式得到Baum-Welch算法的步骤。

预测算法

预算算法是用来解决预测问题的，即给定模型参数和观测序列，求出最有可能的状态序列。

近似算法

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

(从式$(10.24)$)给定隐马尔可夫模型$\lambda$和观测序列$O$，在时刻$t$处于状态$q_i$的概率${\gamma }_t(i)$是
$${\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(i){\beta }_t(i)}\text{\hspace{1em}(10.42)}$$
在每一时刻$t$最有可能的状态$i_t^{\ast }$是
$$i_t^{\ast }=\arg \underset{1\le i\le N}{\max }[{\gamma }_t(i)],t=1,2,\cdots ,T\text{\hspace{1em}(10.43)}$$
从而得到状态序列$I^{\ast }=(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast })$。

近似算法的优点是计算简单，缺点是不能保证预测的状态序列整体是最有可能的状态序列。

维特比算法

维特比算法实际是用动态规划解决隐马尔可夫模型预测问题，即用动态规划求概率最大路径(最优路径)。这时一条路径对应着一个状态序列。

依据动态规划原理，我们从时刻$t=1$开始，递推地计算在时刻$t$状态为$i$的各条部分路径的最大概率，直到得到时刻$t=T$状态为$i$的各条路径的最大概率。然后，在时刻$t=T$的最大概率即为最优路径的概率$P^{\ast }$，最优路径的终结点$i_T^{\ast }$也同时得到。然后从终结点开始，由后向前逐步求得结点$i_{T-1}^{\ast },\cdots ,i_1^{\ast }$，得到最优路径。

首先导入两个变量$\delta$和$\Psi$，定义在时刻$t$状态为$i$的所有单个路径$(i_1,i_2,\cdots ,i_t)$中概率最大值为
$${\delta }_t(i)=\underset{i_1,\cdots ,i_{t-1}}{\max }P(i_t=i,i_{t-1},\cdots ,i_1,o_t,\cdots ,o_1|\lambda ),i=1,2,\cdots ,N\text{\hspace{1em}(10.44)}$$
由定义可得变量$\delta$的递推公式：
$$\begin{array}{rlr}{\delta }_{t+1}(i)& =\underset{i_1,\cdots ,i_t}{\max }P(i_{t+1}=i,i_t,\cdots ,i_1,o_{t+1},\cdots ,o_1|\lambda )& \\ & =\underset{i_1,\cdots ,i_{t-1},i_t}{\max }P(i_{t+1}=i,i_t,\cdots ,i_1,o_{t+1},\cdots ,o_1|\lambda )& \\ & =\underset{1\le j\le N}{\max }\underset{i_1,\cdots ,i_{t-1}}{\max }P(i_{t+1}=i,i_t=j,\cdots ,i_1,o_{t+1},\cdots ,o_1|\lambda )& \\ & =\underset{1\le j\le N}{\max }\underset{i_1,\cdots ,i_{t-1}}{\max }P(o_{t+1},i_{t+1}=i|i_t=j,\cdots ,i_1,o_t,\cdots ,o_1,\lambda )P(i_t=j,\cdots ,i_1,o_t,\cdots ,o_1|\lambda )& \\ & =\underset{1\le j\le N}{\max }{\delta }_t(j)P(o_{t+1},i_{t+1}=i|i_t=j,\cdots ,i_1,o_t,\cdots ,o_1,\lambda )& \\ & =\underset{1\le j\le N}{\max }{\delta }_t(j)P(o_{t+1},i_{t+1}=i|i_t=j,\lambda )& D-划分\\ & =\underset{1\le j\le N}{\max }{\delta }_t(j)P(o_{t+1}|i_{t+1}=i,i_t=j,\lambda )P(i_{t+1}=i|i_t=j,\lambda )& \\ & =\underset{1\le j\le N}{\max }{\delta }_t(j)P(o_{t+1}|i_{t+1}=i,\lambda )P(i_{t+1}=i|i_t=j,\lambda )& \\ & =\underset{1\le j\le N}{\max }{\delta }_t(j)b_i(o_{t+1})a_{ji}& \\ & =\underset{1\le j\le N}{\max }[{\delta }_t(j)a_{ji}]b_i(o_{t+1}),i=1,2,\cdots ,N;t=1,2,\cdots ,T-1& \end{array}\text{\hspace{1em}(10.45)}$$
$P(o_{t+1},i_{t+1}=i|i_t=j,\cdots ,i_1,o_t,\cdots ,o_1,\lambda )=P(o_{t+1},i_{t+1}=i|i_t=j,\lambda )$利用了D-划分：

维特比算法很像前向算法除了前者取上个时刻路径概率的最大值，而后者取的是求和。

该递推公式除了通过推导，还可以通过画图来理解。如上图所示，已知$t$时刻各个状态下的${\delta }_t(j),1\le j\le N$。那么$${\delta }_{t+1}(i)$$即在$t$时刻状态为$j$的${\delta }_t(j)$乘以由状态$j$转移到$t+1$时刻状态$i$的概率的最大者，再乘以由状态$i$观测到输出$o_{t+1}$的概率。

定义在时刻$t$状态为$i$的所有单个路径$(i_1,i_2,\cdots ,i_t)$中概率最大的路径的第$t-1$个结点为
$${\Psi }_t(i)=\arg \underset{1\le j\le N}{\max }[{\delta }_t(j)a_{ji}],i=1,2,\cdots ,N\text{\hspace{1em}(10.46)}$$
算法 10.5 (维特比算法)

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

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

(1) 初始化
$$\begin{gathered} {\delta }_1(i)={\pi }_i{b}_i(o_1),i=1,2,\cdots ,N \\ {\Psi }_1(i)=0,i=1,2,\cdots ,N \end{gathered}$$
在时刻$t=1$状态为$i$的概率最大值即为初始概率乘以$b_i(o_1)$，${\Psi }_1(i)=0$表示未知，或者说后面回溯时的终止条件。

(2) 递推。对$t=2,3,\cdots ,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,\cdots ,N \\ {\Psi }_t(i)=\arg \underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}],i=1,2,\cdots ,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,T-2,\cdots ,1$

上一步得到的$T$时刻的终结点，然后从$T-1$时刻利用$i_{T-1}^{\ast }={\Psi }_T(i_T^{\ast })$求得$T$时刻状态为$i_T^{\ast }$概率最大路径的前一个时刻最优结点$i_{T-1}^{\ast }$，以此类推：
$$i_t^{\ast }={\Psi }_{t+1}(i_{t+1}^{\ast })$$
求得最优路径$I^{\ast }=(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast })$。

下面通过书上的例子，并画图来理解一下维特比算法。

例 10.3 已知模型参数$\lambda =(A,B,\pi )$和观测序列$O=(红，白，红)$，试求最优状态序列，即最优路径。

按照上面介绍的维特比算法一步一步来求。

对于$t=1$，即初始化时，对所有的状态$i(i=1,2,3)$求状态为$i$观测$o_1$为红球的概率${\delta }_1(i)$：
$$\begin{gathered} {\delta }_1(1)={\pi }_1{b}_1(o_1)=0.2\times 0.5=0.10 \\ {\delta }_1(2)={\pi }_2{b}_2(o_1)=0.4\times 0.4=0.16 \\ {\delta }_1(3)={\pi }_3{b}_3(o_1)=0.4\times 0.7=0.28 \end{gathered}$$
且${\Psi }_1(i)=0,i=1,2,3$。

此时，如下图所示：

对于$t=2$，对所有的状态$i(i=1,2,3)$求状态为$i$观测$o_2$为白球的概率${\delta }_2(i)$。

对于$i=1$有：
$${\delta }_2(1)=\underset{1\le j\le 3}{\max }[{\delta }_1(j)a_{j1}]b_1(o_2)=\max [0.1\times 0.5,0.16\times 0.3,0.28\times 0.2]\times 0.5=0.056\times 0.5=0.028$$

类似地，对于$i=2$有：
$${\delta }_2(2)=\underset{1\le j\le 3}{\max }[{\delta }_1(j)a_{j2}]b_2(o_2)=\max [0.1\times 0.2,0.16\times 0.5,0.28\times 0.3]\times 0.6=0.084\times 0.6=0.0504$$

对于$i=3$有：
$${\delta }_2(3)=\underset{1\le j\le 3}{\max }[{\delta }_1(j)a_{j3}]b_3(o_2)=\max [0.1\times 0.3,0.16\times 0.2,0.28\times 0.5]\times 0.3=0.14\times 0.3=0.042$$

这样我们得到了在时刻$t=2$时，转移到各个状态的最优路径，这里恰巧它们都是从$t=1$时状态$q_3$出发的。

同理，对于$t=3$，对所有的状态$i(i=1,2,3)$求状态为$i$观测$o_3$为红球的概率${\delta }_3(i)$。
$$\begin{gathered} {\delta }_3(1)=\underset{1\le j\le 3}{\max }[{\delta }_2(j)a_{j1}]b_1(o_3)=\max [0.028\times 0.5,\overline{0.0504}\times 0.3,0.042\times 0.2]\times 0.5=0.01512\times 0.5=0.000756 \\ {\delta }_3(2)=\underset{1\le j\le 3}{\max }[{\delta }_2(j)a_{j2}]b_2(o_3)=\max [0.028\times 0.2,\overline{0.0504}\times 0.5,0.042\times 0.3]\times 0.4=0.02520\times 0.4=0.001008 \\ {\delta }_3(3)=\underset{1\le j\le 3}{\max }[{\delta }_2(j)a_{j3}]b_3(o_3)=\max [0.028\times 0.3,0.0504\times 0.2,\overline{0.042}\times 0.5]\times 0.7=0.02100\times 0.7=0.001470 \end{gathered}$$
最终我们得到的结果如下：

回溯得到最优路径的状态序列$(q_3,q_3,q_3)$。

参考

1. 统计学习方法 第二版
2. 概率导论
3. 维基百科
4. PRML
5. MLAPP