隐马尔可夫模型(二)

马尔可夫模型的学习，根据训练数据是包括观测序列和对应的状态序列还是只有观测序列，可以分别由监督学习和非监督学习实现。本部分首先介绍监督学习，后介绍非监督学习Baum-Welch算法

监督学习方法

假设已给训练数据包含$s$个长度相同的观测序列和对应的状态序列 { ( O 1 , I 1 ) , ( O 2 , I 2 ) , ⋯   , ( O s , I s ) } \left\{\left(O_{1}, I_{1}\right),\left(O_{2}, I_{2}\right), \cdots,\left(O_{s}, I_{s}\right)\right\} {(O1​,I1​),(O2​,I2​),⋯,(Os​,Is​)}那么可以用极大似然估计法来估计隐马尔可夫模型的参数。具体方法如下：

1. 转移概率$a_{ij}$的估计
   设样本中时刻$t$处于状态$i$转移到时刻$t+1$状态的频数为$A_{ij}$,那么状态转移概率$a_{ij}$的估计是$${\hat{a}}_{ij}=\frac{A_{ij}}{{\sum }_{j=1}^N{A}_{ij}},i=1,2,\cdots ,N;j=1,2,\cdots ,N$$
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(1)$$
3. 初始状态概率${\pi }_i$的估计${\hat{\pi }}_i$为$S$个样本中初始状态为$q_i$的频率
   由于监督学习需要使用训练数据，而人工标注训练数据的代价往往会很高，有时会利用非监督学习方法。

aum-Welch算法
假设给定训练数据只包含$S$个长度为$T$的观测序列 { O 1 , O 2 , ⋯   , O s } \left\{O_{1}, O_{2}, \cdots, O_{s}\right\} {O1​,O2​,⋯,Os​}而没有对应的状态序列，目标是学习隐马尔可夫模型$\lambda =(A,B,\pi )$的参数，我们将观测序列数据看作观测书记$O$，状态序列数据看作不可观测的隐数据$I$，那么隐马尔可夫模型事实上是一个含有隐变量的概率模型$$P(O;\lambda )=\sum _{I}P(O|I;\lambda )P(I;\lambda )(2)$$
其参数学习可有由$\text{EM}$算法实现

1. 确定完全数据的对数似然函数
   所有观测数据写成$O=\left(o_1,o_2,\cdots ,o_T\right)$，所有隐数据写成$I=\left(i_1,i_2,\cdots ,i_T\right)$，完全数据是$(O,I)=\left(o_1,o_2,\cdots ,o_T,i_1,i_2,\cdots ,i_T\right)$，完全数据的对数似然函数是$\log P(O,I;\lambda )$
2. $\text{EM}$算法中的$\text{E}$步：求$Q$函数$Q(\lambda ,\overline{\lambda })$
   按照$Q$函数的定义 $Q(\lambda ,\overline{\lambda })=E_t[\log P(O,I;\lambda )|O;\overline{\lambda }]$
   略去对$\lambda$而言的常数因子$1/P(O;\overline{\lambda })$,得$$Q(\lambda ,\overline{\lambda })=\sum _I\log P(O,I;\lambda )P(O,I;\overline{\lambda })(3)$$其中$\overline{\lambda }$是隐马尔可夫模型参数的当前估计值,$\lambda$是要极大化的隐马尔可夫模型参数$$\begin{gathered} P(O,I;\lambda )=P(O|I;\lambda )P(I;\lambda ) \\ ={\pi }_{i_1}{b}_{i_1}\left(o_1\right)a_{i_1{i}_2}{b}_{i_2}\left(o_2\right)\cdots a_{i_{T-1}{i}_T}{b}_{i_T}\left(o_T\right) \end{gathered}$$于是函数$Q(\lambda ,\overline{\lambda })$可以写成：$$\begin{gathered} Q(\lambda ,\overline{\lambda })=\sum _I\log {\pi }_{i_1}P(O,I;\overline{\lambda }) \\ +\sum _I\left(\sum _{t=1}^{T-1}\log a_{i_t{i}_{t+1}}\right)P(O,I;\overline{\lambda })+\sum _I\left(\sum _{t=1}^T\log b_{i_t}\left(o_t\right)\right)P(O,I;\overline{\lambda })(4) \end{gathered}$$
   式中求和都是对所有训练数据的序列总长度$T$进行的
3. $\text{EM}$算法中的$\text{M}$步:极大化$Q$函数$Q(\lambda ,\overline{\lambda })$求模型参数$A,B,\pi$
   由于要极大化的参数在上式$(4)$中单独地出现在$3$个项中，所以只需对各项分别极大化$(1)$上式$(4)$第一项可以写成$$\sum _I\log {\pi }_{i_0}P(O,I;\overline{\lambda })=\sum _{i=1}^N\log {\pi }_{i}P\left(O,i_1=i;\overline{\lambda }\right)$$注意到${\pi }_i$满足约束条件${\sum }_{i=1}^N{\pi }_i=1$，利用拉格朗日乘子法，写出拉格朗日函数：$$\sum _{i=1}^N\log {\pi }_{i}P\left(O,i_1=i;\overline{\lambda }\right)+\gamma \left(\sum _{i=1}^N{\pi }_i-1\right)$$对其求偏导数并令结果为0 $$\frac{\partial }{\partial {\pi }_i}\left[\sum _{i=1}^N\log {\pi }_{i}P\left(O,i_1=i;\overline{\lambda }\right)+\gamma \left(\sum _{i=1}^N{\pi }_i-1\right)\right]=0(5)$$得 $$P\left(O,i_1=i;\overline{\lambda }\right)+\gamma {\pi }_i=0$$对$i$求和得到$\gamma$ $$\gamma =-P(O;\overline{\lambda })$$ 代入上式得$${\pi }_i=\frac{P\left(O,i_1=i;\overline{\lambda }\right)}{P(O;\overline{\lambda })}(6)$$
   $(2)$ 式$(4)$第二项可以写成$$\sum _I\left(\sum _{t=1}^{T-1}\log a_{i_t,i_{t+1}}\right)P(O,I;\overline{\lambda })=\sum _{i=1}^N\sum _{j=1}^N\sum _{i=1}^{N-1}\log a_{ij}P\left(O,i_t=i,i_{t+1}=j;\overline{\lambda }\right)4$$ 类似第一项，应用具体约束条件${\sum }_{j=1}^N{a}_{ij}=1$的拉格朗日乘子法可以求出$$a_{ij}=\frac{{\sum }_{i=1}^{T-1}P\left(O,i_t=i,i_{t+1}=j;\overline{\lambda }\right)}{{\sum }_{i=1}^{T-1}P\left(O,i_t=i;\overline{\lambda }\right)}$$ $(3)$ 式$(4)$第三项为$$\sum _I\left(\sum _{t=1}^T\log b_{i_t}\left(o_t\right)\right)P(O,I;\overline{\lambda })=\sum _{j=1}^N\sum _{t=1}^T\log b_j\left(o_t\right)P\left(O,i_t=j;\overline{\lambda }\right)(7)$$同样用拉格朗日乘子法，约束条件是${\sum }_{k=1}^M{b}_j(k)=1$，注意只有在$o_t=v_k$时$b_j\left(o_t\right)$对$b_j(k)$的偏导数才不为$0$，以$I\left(o_t=v_k\right)$表示，求得$$b_j(k)=\frac{{\sum }_{i=1}^{T}P(O,i_t=j;\overline{\lambda })I\left(o_t=v_k\right)}{{\sum }_{i=1}^{T}P\left(O,i_t=j;\overline{\lambda }\right)}(8)$$
   Baum-Welch模型参数估计公式
   将式$(6)-(8)$中的各概率分别用${\gamma }_t(i),{\xi }_i(i,j)$，则可将相应的公式写成：$$a_{ij}=\frac{{\sum }_{i=1}^{T-1}{\xi }_i(i,j)}{{\sum }_{i=1}^{T-1}{\gamma }_t(i)}(9)$$ $$b_j(k)=\frac{{\sum }_{t=1,o_k=y_k}^T{\gamma }_t(j)}{{\sum }_{i=1}^T{\gamma }_t(j)}(10)$$ $${\pi }_i={\gamma }_1(i)(11)$$式$(9)-(11)$就是Baum-Welch算法，它是$\text{EM}$算法在隐马尔可夫模型学习中的具体实现。
   预测算法

近似算法
近似算法，在每个时刻$t$选择在该时刻最有可能出现的状态$i_t^{\ast }$,从而得到一个状态序列$I^{\ast }=\left(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast }\right)$，将它作为预测的结果。
给定隐马尔可夫模型$\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 }_{i=1}^N{\alpha }_t(j){\beta }_t(j)}$$
在每一时刻$t$最有可能的状态$i_t^{\ast }$是$i_t^{\ast }=\arg {\max }_{1⩽i⩽N}\left[{\gamma }_t(i)\right],t=1,2,\cdots ,T$, 从而得到状态序列$I^{\ast }=\left(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast }\right)$
近似算法的优点是计算简单，其缺点是不能保证预测的状态序列整体是最有可能的状态序列，因为预测的状态序列可能有实际不发生的部分。

维比特算法
维比特算法实际是用动态规划解隐马尔可夫模型预测问题，即用动态规划求概率最大路径。这时一条路径对应着一个状态序列。
根据动态规划原理，最优路径具有这样的特性：如果最优路径在时刻$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 },\cdots ,\dot{,}{i}_1^{\ast }$,得到最优路径$I^{\ast }=\left(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast }\right)$。
首先导入两个变量$\delta$和$\psi$。定义在时刻$t$状态为$i$的所有单个路径$\left(i_1,i_2,\cdots ,i_t\right)$中概率最大值为$${\delta }_t(i)=\underset{i_1,i_2,\cdots ,i_{t-1}}{\max }P\left(i_t=i,i_{t-1},\cdots ,i_1,o_t,\cdots ,o_1;\lambda \right),i=1,2,\cdots ,N$$由定义可得变量$\delta$的递推公式：$$\begin{gathered} {\delta }_{t+1}(i)=\underset{i_1,i_2,\cdots ,i_t}{\max }P\left(i_{t+1}=i,i_t,\cdots ,i_1,o_{t+1},\cdots ,o_1;\lambda \right) \\ =\underset{1⩽j⩽N}{\max }\left[{\delta }_t(j)a_{ji}\right]b_i\left(o_{t+1}\right),i=1,2,\cdots ,N;t=1,2,\cdots ,T-1 \end{gathered}$$定义在时刻$t$状态为$i$的所有单个路径$\left(i_1,i_2,\cdots ,i_{t-1},i\right)$中概率最大的路径的第$t-1$个结点为$${\psi }_t(i)=\arg \underset{1⩽j⩽N}{\max }\left[{\delta }_{t-1}(j)a_{ji}\right],i=1,2,\cdots ,N$$维比特算法
输入：模型$\lambda =(A,B,\pi )$和观测$O=\left(o_1,o_2,\cdots ,o_T\right)$
输出：最优路径$I^{\ast }=\left(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast }\right)$
$(1)$初始化 $${\delta }_1(i)={\pi }_i{b}_i\left(o_1\right),i=1,2,\cdots ,N$$ $${\psi }_1(i)=0,i=1,2,\cdots ,N$$
$(2)$递推，对$t=2,3,\cdots ,T$
$${\delta }_t(i)=\underset{1⩽j⩽N}{\max }\left[{\delta }_{t-1}(j)a_{ji}\right]b_i\left(o_t\right),i=1,2,\cdots ,N$$ $${\psi }_t(i)=\arg \underset{1\in j\le N}{\max }\left[{\delta }_{t-1}(j)a_{ji}\right],i=1,2,\cdots ,N$$$(3)$终止$$\begin{gathered} P^{\ast }=\underset{1⩽i⩽N}{\max }{\delta }_T(i) \\ {i}_T^{\ast }=\arg \underset{1⩽i⩽N}{\max }\left[{\delta }_T(i)\right] \end{gathered}$$ $(4)$最优路径回溯。对$t=T-1,T-2,\cdots ,1$$$i_t^{\ast }={\psi }_{t+1}\left(i_{t+1}^{\ast }\right)$$求得最优路径$$I^{\ast }=\left(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast }\right)$$