隐马尔可夫模型(HMM explained)

隐马尔可夫模型

隐马尔可夫模型（Hidden Markov Model，HMM），是上世纪七十年代由美国数学家Leonard E. Baum等人提出来的，可以用于解决标注（tagging）问题，属于一种生成模型（generative model），即直接生成观测序列和状态序列的联合概率分布。

1 描述

1.1 术语解释：

• 状态序列（state sequence）：模型假设由一个隐藏的马尔可夫链生成一个不可观测的序列，称这个序列为状态序列，其中序列的每个时间步都代表着一个不可被观测到的状态
• 观测序列（observation sequence）：状态序列的每个状态都会对应生成一个结果，这个结果是可以观测得到的，将这个由状态序列生成的结果序列称为观测序列

1.2 形式定义

序列

• 状态序列：$I=(i_1,i_2,...,i_T)$
• 观测序列：$O=(o_1,o_2,...,o_T)$
  （其中T为序列长度，也可以理解为总的时间步）

集合

• 状态集合： Q = { q 1 , q 2 , . . . , q N } Q = \{q_1, q_2, ..., q_N\} Q={q1​,q2​,...,qN​}
• 观测结果集合：$V=(v_1,v_2,...,v_M)$
  （其中N为可能的状态数，M为可能的观测结果数目）

概率矩阵

• 状态转移概率矩阵：$A=[a_{ij}{]}_{N\times N}$
  矩阵的每个元素代表从一个状态$q_i$转移到另一个状态$q_j$的概率，即$a_{ij}=Pr(i_{t+1}=q_j|i_t=q_i)$
• 观测概率矩阵：$B=[b_j(k){]}_{N\times M}$
  矩阵每个元素代表从一个状态$q_j$生成观测结果$v_k$的概率，即$b_j(k)=Pr(o_t=v_k|i_t=q_j)$

初始状态矩阵概率

• 初始状态概率向量： $\pi = (\pi_i) $
  向量的元素代表在第一时刻（t=1），处于状态$q_i$的概率，即${\pi }_i=Pr(i_1=q_i)$

有了以上的形式定义，我们可以使用$\lambda =(A,B,\pi )$来表示模型，这三个元素可以决定一个隐马尔可夫模型，因为由初始状态$\pi$和状态转移概率矩阵$A$可以确定隐藏的状态序列$I$，然后通过得出的这个状态序列$I$以及观测概率矩阵$B$，可以确定如何从这些状态生成观测结果，也就是确定观测序列$O$。

2 例子

2.1 问题描述

假设有三个盒子，每个盒子里有红白两种颜色的球，盒子里的红白球数目如下：

盒子标号	1	2	3
红球数目	5	3	6
白球数目	5	7	4

抽取球的规则为：首先从三个盒子里面以概率选择一个盒子，然后在这个盒子里面随机选择一个球，记录球的颜色后放回，并以当前盒子随机转移到下一个盒子。

挑选初始盒子的规则为：0.5的概率挑选盒子1，0.25的概率挑选盒子2，0.25的概率挑选盒子3

转移盒子的规则为：如果当前是盒子1，则一定转移到盒子2，如果当前是盒子2，则0.4的概率转移到1，以0.6的概率转移到3，如果当前为盒子3，则以0.5的概率转移到1，以0.5的概率转移到2

重复这个过程五次，最后得到颜色的观测序列为： O = { 红 ， 红 ， 白 ， 白 ， 红 } O= \{红， 红， 白， 白， 红\} O={红，红，白，白，红}。注意，这里我们只能观测到选出球的颜色，而盒子的序列是不可观测的，也就是说盒子的序列为隐藏的状态序列。

2.2 模型定义

• 问题的集合如下：

状态的集合为： Q = { 盒 子 1 ， 盒 子 2 ， 盒 子 3 } , N = 3 Q = \{盒子1， 盒子2， 盒子3\}, N=3 Q={盒子1，盒子2，盒子3},N=3

观测的集合为： V = { 红 ， 白 } , M = 2 V = \{红， 白\}, M=2 V={红，白},M=2

• 概率矩阵：

状态转移概率矩阵：
$$A=\left(\begin{array}{ccc}0& 1& 0\\ 0.4& 0& 0.6\\ 0.5& 0.5& 0\end{array}\right)$$

观测概率矩阵：
$$B=\left(\begin{array}{cc}0.5& 0.5\\ 0.3& 0.7\\ 0.6& 0.4\end{array}\right)$$

• 初始概率向量：

$$\pi =(0.5,0.25,0.25)$$

根据以上定义，可以确定隐马尔可夫模型$\lambda =(A,B,\pi )$

2.3 三个基本问题

Rabiner(1989)¹提出隐马尔可夫模型有下面三个基本问题：

1. Likelihood：即概率计算问题，给定模型$\lambda =(A,B,\pi )$和观测序列$O$，计算在模型下观测序列出现的概率$Pr(O|\lambda )$
2. Learning：即学习问题，给定观测序列$O$，估计模型参数$A$ $B$ 和$\pi$
3. Decoding：即预测问题，给定模型$\lambda =(A,B,\pi )$和观测序列$O$，求最可能的状态序列$I$

问题一：Likelihood Computation

1 直接计算的方法

概率计算问题定义如下：

给定模型$\lambda =(A,B,\pi )$和观测序列$O=(o_1,o_2,...,o_T)$，计算在模型下观测序列出现的概率$Pr(O|\lambda )$

一个最简单的方法是按照概率公式直接计算，方法是列举所有可能的状态序列$I=(i_1,i_2,...,i_T)$，然后求每个状态序列产生观测序列的概率，最后对这些可能的状态序列求和，就得出观测序列的概率了。

• 具体过程如下：

首先计算状态序列的概率：
$$Pr(I|\lambda )={\pi }_{i_1}{a}_{i_1{i}_2}{a}_{i_2{i}_3}...a_{i_{T-1}{i}_T}$$

然后计算给定状态序列下，观测序列的概率：

$$Pr(O|I,\lambda )=b_{i_1}(o_1)b_{i_2}(o_2)...b_{i_T}(o_T)$$

所以可以求出状态序列和观测序列的联合概率（所以说HMM是生成模型）：

$$\begin{gathered} Pr(O,I|\lambda )=Pr(I|\lambda )Pr(O|I,\lambda ) \\ ={\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1{i}_2}{b}_{i_{\mathrm{２}}}(o_{\mathrm{２}})...a_{i_{T-1}{i}_{T}2}{b}_{i_T}(o_T) \end{gathered}$$

最后将可能的状态序列求和就可以得出最终的观测序列概率：
$$\begin{gathered} Pr(O|\lambda )=\sum _{I}Pr(I|\lambda )Pr(O|I,\lambda ) \\ =\sum _{i_1,i_2,...,i_T}{\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1{i}_2}{b}_{i_{\mathrm{２}}}(o_{\mathrm{２}})...a_{i_{T-1}{i}_{T}2}{b}_{i_T}(o_T) \end{gathered}$$
显然这种方法需要遍历T时间步的可能状态序列，而每个时间步都可能是$N$个状态之一，所以总的计算量是$O(T{N}^T)$的，显然是不可行的。

这里引入一种动态规划（dynamic programming，DP）的算法来解决计算量问题，称之为前向后向算法（forward-backward algorithm）。动态规划算法使用表来存储中间值以减少重复计算，以空间换时间，这里的programming本来指的是表格的意思，DP可直译为动态的表格。

2 前向后向算法

1) 前向算法

定义在时刻t，给定隐马尔可夫模型$\lambda$，观测序列为$o_1,o_2,...,o_t$且状态为$q_i$的概率为前向概率，即：
$${\alpha }_t(i)=Pr(o_1,o_2,...,o_t,q_i|\lambda )$$
显然前向概率可以递推的求解，从而可以递推的求出观测序列概率$P(O|\lambda )$，算法具体过程如下：

1. 初始化： ${\alpha }_1(j)={\pi }_j{b}_j(o_1)$
2. 递推：${\alpha }_{t+1}(i)=\left[{\sum }_{j=1}^N{\alpha }_t(j)\right]b_i(o_{t+1})$
3. 终止：$Pr(O|\lambda )={\sum }_{i=1}^N{\alpha }_T(i)$

2) 后向算法

定义在时刻t，给定隐马尔可夫模型$\lambda$，状态为$q_i$的条件下，从t+1到T的部分观测序列为$o_{t+1},o_{t+2},...,o_T$的概率为后向概率，即：
$${\beta }_t(i)=Pr(o_{t+1},o_{t+2},...,o_T|i_t=q_i,\lambda )$$
可以递推的方式得到后向概率，从而获得观测序列概率$P(O|\lambda )$，算法过程如下：

1. 初始化： ${\beta }_T(i)=1,i=1,2,...,N$
2. 递推：${\beta }_t(i)={\sum }_{j=1}^N{a}_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j),i=1,2,...,N$
3. 终止：$Pr(O|\lambda )={\sum }_{i=1}^N{\pi }_i{b}_i(o_1){\beta }_1(i)$

问题二：Training

模型训练的问题定义如下：

给定观测序列$O=(o_1,o_2,...,o_T)$，确定模型的参数$\lambda =(A,B,\pi )$

隐马尔可夫模型的训练算法称为Baum-Welch算法，是EM算法的一个特例。EM算法是一种迭代的算法，算法分为两步，第一步求期望Expectation，第二步最大化期望Maximization。简单来讲就是初始化概率的估计值，然后使用这些估计值来计算一个更好的估计，不断的迭代提升。

MSE vs EM（https://machinelearningmastery.com/expectation-maximization-em-algorithm/）

Maximum likelihood estimation is an approach to density estimation for a dataset by searching across probability distributions and their parameters.

It is a general and effective approach that underlies many machine learning algorithms, although it requires that the training dataset is complete, e.g. all relevant interacting random variables are present. Maximum likelihood becomes intractable if there are variables that interact with those in the dataset but were hidden or not observed, so-called latent variables.

The expectation-maximization algorithm is an approach for performing maximum likelihood estimation in the presence of latent variables. It does this by first estimating the values for the latent variables, then optimizing the model, then repeating these two steps until convergence. It is an effective and general approach and is most commonly used for density estimation with missing data, such as clustering algorithms like the Gaussian Mixture Model.

EM算法详细步骤如下：

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

完全数据定义为：$(Q,I)=(o_1,o_2,...,o_T,i_1,i_2,...,i_T)$，它的对数似然函数为：$\log Pr(O,I|\lambda )$

2. E步：求$Q$函数

$Q$函数为：$Q(\lambda ,\hat{\lambda })={\sum }_I\log Pr(O,I|\lambda )Pr(O,I|\hat{\lambda })$

这里$ \hat{\lambda}$是模型\ast 当前的估计值\ast ，$\lambda$是通过计算而获得的更好的估计值。

这里的Q函数实际上省略了常数因子$ \frac{1}{Pr(O|\hat{\lambda})}$，按照定义，$Q(\lambda, \hat{\lambda}) =\mathbb{E}_I[\log Pr(O, I | \lambda)|O, \hat{\lambda}]$，展开得到$\sum_I \frac{\log Pr(O, I|\lambda)Pr(O,I|\hat{\lambda})}{Pr(O|\hat{\lambda})} = \frac{1}{Pr(O|\hat{\lambda})} \sum_I \log Pr(O, I|\lambda)Pr(O,I|\hat{\lambda})$。

其中$Pr(O,I|\lambda )={\pi }_{i_1}{b}_{i_1}(o_1)a_{i_1{i}_2}{b}_{i_2}{o}_2...a_{i_{T-1}{i}_T}{b}_{i_T}(o_T)$

所以$\log Pr(O,I|\lambda )=\log {\pi }_{i_1}+{\sum }_{t=1}^{T-1}{a}_{i_t{i}_{t+1}}+{\sum }_{t=1}^T{b}_{i_t}(o_t)$

带入到Q函数，$Q(\lambda ,\hat{\lambda })={\sum }_I\log {\pi }_{i_1}Pr(O,I|\hat{\lambda })+{\sum }_I\left({\sum }_{t=1}^{T-1}\log a_{i_t{i}_{t+1}}\right)Pr(O,I|\hat{\lambda })+{\sum }_I\left({\sum }_{t=1}^T\log b_{i_t}(o_t)\right)Pr(O,I|\hat{\lambda })$

3. M步：极大化Q函数

这个步骤就是求能够极大化Q函数的模型参数，即$A$、$B$和$\pi$。

从上面的Q函数表达式可以看出，要求的模型参数分别位于三个项中，所以可以逐个对各项最大化。这是一个约束优化问题，因为是概率，所以要满足归一性的约束条件，运用拉格朗日乘子法可以解决这个问题。

(1)最大化第一项

第一项可以写为：
$$\sum _I\log {\pi }_{i_1}Pr(O,I|\hat{\lambda })=\sum _{i=1}^N\log {\pi }_{i}Pr(O,i_1=i|\hat{\lambda })$$
而且${\pi }_i$需要满足约束条件${\sum }_{i=1}^N{\pi }_i=1$，所以这里需要借助拉格朗日乘子解决这个约束优化问题：
$$\sum _{i=1}^N\log {\pi }_{i}Pr(O,i_1=i|\hat{\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}Pr(O,i_1=i|\hat{\lambda })+\gamma \left(\sum _{i=1}^N{\pi }_i-1\right)\right]=0$$
得：
$$\begin{gathered} \frac{Pr(O,i_1=i|\hat{\lambda })}{{\pi }_i}+\gamma =0 \\ Pr(O,i_1=i|\hat{\lambda })+\gamma {\pi }_i=0 \end{gathered}$$
对i求和得（这里${\sum }_{i=1}^N{\pi }_i=1$）：
$$\begin{gathered} \sum _{i=1}^{N}Pr(O,i_1=i|\hat{\lambda })+\sum _{i=1}^N\gamma {\pi }_i=0 \\ Pr(O|\hat{\lambda })+\gamma =0 \end{gathered}$$
解得$\gamma$：
$$\gamma =-Pr(O|\hat{\lambda })$$
带入式$Pr(O,i_1=i|\hat{\lambda })+\gamma {\pi }_i=0$可以解出${\pi }_i$：
$${\pi }_i=\frac{Pr(O,i_1=i|\hat{\lambda })}{Pr(O|\hat{\lambda })}$$

(2)最大化第二项

第二项可以写成（这里前面两个求和是对于转移矩阵$A$的， 后面求和是对于长度为$T$的状态序列$I$的）：
$$\sum _I\left(\sum _{t=1}^{T-1}\log a_{i_t{i}_{t+1}}\right)Pr(O,I|\hat{\lambda })=\sum _{i=1}^N\sum _{j=1}^N\sum _{t=1}^{T-1}\log a_{ij}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })$$
同样是一个约束优化问题，约束条件为${\sum }_{j=1}^N{a}_{ij}=1$，应用拉格朗日乘子法可以解决这个问题，过程如下：

添加拉格朗日乘子：
$$\sum _{i=1}^N\left[\sum _{j=1}^N\sum _{t=1}^{T-1}\log a_{ij}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })+\gamma \left(\sum _{j=1}^N{a}_{ij}-1\right)\right]$$
计算其对于$a_{ij}$的偏导数并令其为0：
$$\begin{gathered} \frac{{\sum }_{t=1}^{T-1}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })}{a_{ij}}+\gamma =0 \\ \sum _{t=1}^{T-1}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })+\gamma a_{ij}=0 \end{gathered}$$
对j求和得：
$$\begin{gathered} \sum _{j=1}^N\sum _{t=1}^{T-1}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })+\sum _{j=1}^N\gamma a_{ij}=0 \\ \sum _{t=1}^{T-1}Pr(O,i_t=i|\hat{\lambda })+\gamma =0 \end{gathered}$$
解得$\gamma$：
$$\gamma =-\sum _{t=1}^{T-1}Pr(O,i_t=i|\hat{\lambda })$$
带入${\sum }_{t=1}^{T-1}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })+\gamma a_{ij}=0$解得$a_{ij}$：
$$a_{ij}=\frac{{\sum }_{t=1}^{T-1}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })}{{\sum }_{t=1}^{T-1}Pr(O,i_t=i|\hat{\lambda })}$$

(2)最大化第三项

第三项可以写为：
$$\sum _I\left(\sum _{t=1}^T\log b_{i_t}(o_t)\right)Pr(O,I|\hat{\lambda })=\sum _{j=1}^N\sum _{t=1}^T\log b_j(o_t)Pr(O,i_t=j|\hat{\lambda })$$
同样需要满足约束条件${\sum }_{k=1}^M{b}_j(k)=1$，过程如下：
$$\sum _{j=1}^N\sum _{t=1}^T\log b_j(o_t)Pr(O,i_t=j|\hat{\lambda })+\gamma \left(\sum _{k=1}^M{b}_j(k)-1\right)$$
计算其对于$b_j(k)$的偏导数并令其为0，这里左项只有在$o_t=v_k$时，$b_j(o_t)$对$b_j(k)$的偏导数才不为0，用$I(o_t=v_k)$表示，结果为：
$$\sum _{t=1}^{T}Pr(O,i_t=j|\hat{\lambda })I(o_t=v_k)+\gamma b_j(k)=0$$
对k求和得：
$$\gamma =-\sum _{t=1}^{T}Pr(O,i_t=j|\hat{\lambda })$$
代回得：
$$b_j(k)=\frac{{\sum }_{t=1}^{T}Pr(O,i_t=j|\hat{\lambda })I(o_t=v_k)}{{\sum }_{t=1}^{T}Pr(O,i_t=j|\hat{\lambda })}$$

总结：

经上面的推导，最后得出三个表达式，分别为：
$$\begin{gathered} {\pi }_i=\frac{Pr(O,i_1=i|\hat{\lambda })}{Pr(O|\hat{\lambda })} \\ {a}_{ij}=\frac{{\sum }_{t=1}^{T-1}Pr(O,i_t=i,i_{t+1}=j|\hat{\lambda })}{{\sum }_{t=1}^{T-1}Pr(O,i_t=i|\hat{\lambda })} \\ {b}_j(k)=\frac{{\sum }_{t=1}^{T}Pr(O,i_t=j|\hat{\lambda })I(o_t=v_k)}{{\sum }_{t=1}^{T}Pr(O,i_t=j|\hat{\lambda })} \end{gathered}$$
我们分别使用${\gamma }_t(i)$和${\xi }_t(i,j)$来表示它们：
$$\begin{gathered} {\pi }_i={\gamma }_1(i) \\ {a}_{ij}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)} \\ {b}_j(k)=\frac{{\sum }_{t=1,o_t=v_k}^T{\gamma }_t(j)}{{\sum }_{t=1}^T{\gamma }_t(j)} \end{gathered}$$

关于这两个概率${\gamma }_t(i)$和${\xi }_t(i,j)$的定义如下，可以通过前向后向算法求出相应的值：

1、给定模型$\lambda $和观测$O$，在时刻t处于状态$q_i $的概率为：
KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ \gamma_t(i) &=…
2、给定模型$\lambda $和观测$O$，在时刻t处于状态$q_i $的且时刻t+1处于状态$q_j$的概率为：
KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ \xi_t(i, j) &=…

以上就是Baum-Welch算法，实际上就是EM算法在学习HMM上的具体实现。

问题三：Decoding

维特比算法（Viterbi algorithm）使用动态规划（DP）解决HMM预测问题（decoding），即使用DP求概率最大路径，一条路径对应着一个状态序列。

引入两个变量$\delta$和$\Psi$。

定义在时刻$t$状态为$i$的所有单个路径$(i_1,i_2,...,i_t)$中概率最大值为：
$${\delta }_t(i)=\underset{i_1,i_2,...,i_{t-1}}{\max }Pr(i_t=i,i_{t-1},...,i_1,o_t,...,o_1|\lambda ),i=1,2,...,N$$
由定义可以得出它的递推公式：
KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ \delta_{t+1}(i…
定义在时刻$t$状态为$i$的所有单个路径$(i_1,i_2,...,i_{t-1},i)$中概率最大的路径的第$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 )$，观测序列$O=(o_q,o_2,...,o_T)$；

输出：最优路径（状态序列）$I^{\ast }=(i_1^{\ast },i_2^{\ast },...i_T^{\ast })$

(1) 初始化。
$$\begin{gathered} {\delta }_1(i)={\pi }_i{b}_i(o_1),i=1,2,...,N \\ {\Psi }_1(i)=0,i=1,2,...,N \end{gathered}$$

(2) 递推。对$t=2,3,...,T$
$$\begin{gathered} {\delta }_t(i)=\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]b_i(o_t) \\ {\Psi }_t(i)=\arg \underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}],i=1,2,...,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,...,1$
$$i_t\ast ={\Psi }_{t+1}(i_{t+1}^{\ast })$$
获得最优路径$I^{\ast }=(i_1^{\ast },i_2^{\ast },...i_T^{\ast })$

---

1. Rabiner, L. R. (1989). A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE, 77(2), 257–286. ↩︎