HMM和Viterbi算法

背景

试图实现一个基于隐马尔科夫模型的拼音输入法，根据用户输入的拼音序列转换为对应的汉字，完成中文的输入。

输入法基本需求

• 通常一个拼音会对应多个汉字，如何确定在当前状态下，输入的拼音对应哪个汉字；
• 拼音序列通常可以对应多种可能的汉字组合，如何确定哪一种组合是相对最好的一种组合，如对于拼音序列yi chang，有异常、一场、宜昌等多种汉字组合，如何进行选择；

HMM简介

模型定义

HMM（Hidden Markov Models，隐马尔科夫模型）是一种基本的统计模型，可以应用在语音识别、自然语言处理、模式识别等很多领域。

HMM的一个前提是马尔科夫假设，即假设模型的当前状态仅依赖于前面的几个状态。马尔科夫假设极大的简化了问题，但也可能丢失一些重要信息。一个马尔科夫过程是指状态间的转移仅依赖于前n个状态的过程，该过程即为n阶马尔科夫模型，其中n是影响下一个状态选择的前n个状态。（通常所说的HMM，指一阶HMM）

HMM是在一个标准的马尔科夫过程中引入一组隐藏状态，以及观察状态与隐藏状态之间的一些概率关系，描述了一个含有隐藏状态的马尔科夫过程。

使用HMM模型时，通常问题有以下两个主要特征：

1. 问题是基于序列的，如时间序列、状态序列；
2. 问题中有两类数据，一类数据序列是可以观测到的，即观察序列，另一类数据序列是不能观测到的，即隐藏状态序列，简称状态序列；

例如在打算要做拼音输入法任务中，可以将拼音看做是隐藏状态，而拼音转化的文字结果作为一系列的观察状态，可以看到，引入隐藏状态的同时，观察序列与隐藏过程也具有一定的概率关系。

可以参考上图，其中$Z_i​$为隐藏状态序列，$X_i​$为隐藏状态生成的观察状态序列。隐藏状态序列$Z_i​$满足马尔科夫过程的要求，且观察状态序列$X_i​$与$Z_i​$之间存在概率关系，即模型中的$X_i​$与$X_i+1​$是存在关联的。

模型参数

首先定义一些基本符号：

• $Q=\lbrace q_1,q_2,…,q_N \rbrace$，Q是所有可能的隐藏状态集合，其中N为可能的隐藏状态数，对应所有可能的汉字的状态数量；
• $V=\lbrace v_1,v_2,…,v_M \rbrace$，V是所有可能的观察状态集合，其中M为可能的观察状态数，对应所有可能的拼音的状态数量；
• $I=\lbrace i_1,i_2,…,i_T \rbrace​$，I 是长度为T的隐藏状态序列，对应用户输入的拼音序列所对应的汉字序列；
• $O=\lbrace o_1,o_2,…o_T \rbrace​$，O 是对应的观察序列，对应用户输入的拼音序列；

HMM除上述状态集合外，还包括三组概率集合，用一个三元组$\lambda​$来刻画HMM，可写作$\lambda = (A,B,\pi)​$。

• A 是隐藏状态转移概率分布，通常用矩阵表示，称为状态转移矩阵：

  $$A=[a_{ij}]_{N\times N} \quad where\quad a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$$

  即，$a_{ij}​$ 是在时刻 t 处于状态 $q_i​$ 的条件下，在时刻 $t+1​$ 转移到状态 $q_j​$ 的概率，对应汉字到汉字之间的转移概率；

• B 是观察状态发射概率分布，通常用矩阵表示，称为混淆矩阵或发射矩阵：

  $$B=[b_{ik}]_{N\times M} \quad where\quad b_{ik}=P(o_t=v_k|i_t=q_i)$$

  即，$b_{ik}​$ 是在时刻 t 下的隐藏状态 $q_i​$ 到观察状态 $v_k​$ 的发射概率，对应汉字到拼音之间的发射概率；

• π 是初始状态概率，用向量表示：

  $$\pi = (\pi_i)_N \quad where \quad \pi=P(i_1=q_i)$$

  即，$\pi_i​$ 是在时刻 t=1 时处于隐藏状态 $q_i​$ 的概率；

可以看到，HMM（一阶）有两个基本假设：

1. 齐次假设：

   $$P(i_t|i_{t-1},o_{t-1},i_{t-2},o_{t-2},...,i_1,o_1)=P(i_t|i_{t-1})$$

   即，任意时刻的隐藏状态只依赖于它前一个时刻的隐藏状态；

2. 观测独立性假设：

   $$P(o_t|i_{T-1},o_{T-1},i_{T-2},o_{T-2},...i_1,o_1)=P(o_t|i_t)$$

   即，任意时刻的观察状态只依赖于当前时刻的隐藏状态；

用HMM解决的三类基本问题

一旦一个系统可以作为HMM被描述，则可以用来解决三个基本问题 ¹。其中前两个是模式识别的问题：给定HMM求一个观察序列的概率（评估），搜索最有可能生成一个观察序列的隐藏状态序列（解码）。第三个问题是给定观察序列生成一个HMM（学习）。

• 评估（概率计算问题）：前向算法——动态规划

  给定模型$\lambda = (A,B,\pi)​$和观察序列$O=\lbrace o_1,o_2,…o_T \rbrace​$，计算在模型 $ \lambda ​$ 参数已知的情况下，计算观察序列 $O​$ 出现的概率 $P(O|\lambda)​$；

• 学习（模型学习问题）：前向-后向算法（也叫Baum-Welch算法）——EM

  已知观察序列$O=\lbrace o_1,o_2,…o_T \rbrace​$，学习使得观察序列概率 $P(O|\lambda)​$最大的模型$\lambda = (A,B,\pi)​$参数；

  即根据观察序列用极大似然估计的方法估计参数生成隐马尔科夫模型；

• 解码（预测问题）：Viterbi算法——动态规划

  已知模型$\lambda = (A,B,\pi)$和观察序列$O=\lbrace o_1,o_2,…o_T \rbrace$，求给定观察序列的情况下条件概率$P(I|O,\lambda)$最大的（隐藏）状态序列$I=\lbrace i_1,i_2,…,i_T \rbrace$；

  即给定观察序列，搜索最有可能的对应的隐藏状态序列；

Viterbi算法

基于以上可知，拼音输入法的问题，属于解码问题，即已知模型和观察序列，求最有可能的对应的隐藏状态序列。

理论上，可以通过枚举所有的状态转移序列来求解解码问题，但效率非常低，暴力枚举的思路是枚举所有的长度T的状态序列，计算该状态序列与观察序列的联合概率。在状态种类为$N$的情况下，共有$N^T$种排列组合，每种组合计算联合概率的计算量为$T$，则总的计算复杂度为$O(TN^T)$ ，可见该方法并不可取。

因此常用维特比（Viterbi）算法来解决。

Viterbi算法是一个通用的求序列最短路径的动态规划算法。

如上图所示，Viterbi算法简单来说就是：从开始状态之后每走一步，就记录下到达该时刻每个状态时，对应该状态所有路径中的概率最大值，并且以这个最大值为基准继续向后前进。显然，如果这个最大值都不能使该状态成为最大似然估计路径上的节点的话，那么该节点的其他更小的概率值（包括对应的路径）就更不可能了。

Viterbi算法通过一种有效的方法来分析HMM模型的观察序列，并捕获最可能的隐藏序列，利用递归的方法减少计算量，并且对于观察序列的整个上下文都进行了很好的考虑，对包含噪音的序列也能进行良好的分析。

Viterbi算法详细过程

该部分内容主要参考52nlp上关于HMM的详细讲解，整理一遍便于自己详细理解。

如上文所说，通过暴力枚举的方法来找到最可能的序列的代价是非常昂贵的，因此Viterbi算法考虑用递归的方式来寻找最有可能的隐藏状态序列。所以，我们首先定义局部概率 $\delta​$ ，表示到达网格中的某个中间状态时的概率值，之后来详细了解如何在$t=1​$和$t=n (n>1)​$时如何计算其局部概率值。

定义局部概率和局部最佳路径

观察上图的网格，显示的是汉字对于观察序列（拼音序列）的一阶转移情况。对于网格中的每一个中间及终止状态，都有一个可以到达该状态的最可能路径。比如，在$t=3​$时刻的三个状态中，每一个都有一条到达该状态的最可能路径，可能如下图所示：

这些路径即为局部最佳路径，而每条局部最佳路径的概率值即为局部概率$\delta​$，用$ \delta(i,t)​$表示在 $t​$ 时刻到达状态 $i​$ 的所有路径概率中最大的概率值，局部最佳路径即对应该最大概率值的隐藏状态序列。基于此可知，在 $t=T​$ 时刻，每个状态都有一个局部概率及相应的局部最佳路径，因此，可以通过选择该时刻局部概率值最大的状态（及其对应的最佳路径）来确定全局最佳路径，即全局最佳隐藏状态序列。

计算$t=1$时刻的局部概率

局部概率指到达该状态时的最佳路径的概率值，当$t=1​$时，该路径是不存在的，因此，用$t=1​$时刻所处状态的初始概率值及相应观察状态的转移概率值来计算$t=1​$时的局部概率，即：

$$\delta _1(i)=\pi (i)b_{io_1}$$

计算$t>1$时刻的局部概率

考虑如图所示的网格：

计算 $t$ 时刻到达状态 $X$ 的最佳路径，显然，这条最佳路径一定会通过 $t-1$ 时刻的状态 $A$、$B$、$C$ 之中的某一个。也就是说，到达状态 $X$ 的最佳路径一定是路径(状态序列），...，A，X、（状态序列），...，B，X、（状态序列），...，C，X三条之中的一条。

我们知道在一阶马尔科夫假设下，状态 $X$ 在一个状态序列之后发生的概率只取决于之前的一个状态，也就是说，路径末端是 $AX$ 的最佳路径将是到达 $A$ 的最佳路径再紧跟 $X$ ，那么这条路径的概率即为：$P(到达状态A的最佳路径)\times P(X|A)\times P(观察状态|X)$。

由此可知，到达状态 $X$ 的最佳路径概率计算方式为：

$$P(X\ at\ time\ t)=\max_{i=A,B,C}P(i\ at\ time\ (t-1))\times P(X|i)\times P(obs\ at\ time\ t|X)$$

其中，$P(i at time (t-1))​$是$t-1​$时刻的局部概率$\delta​$，$P(X|i)​$是状态转移概率（隐藏状态到观察状态的发射概率），$P(obs at time t|X)​$是观察概率。

对上述公式进行泛化可知，在 $t​$ 时刻，观察状态是 $k_t​$ ，到达隐藏状态 $i​$ 最佳局部路径概率为：

$$\delta _t(i)=\max \limits_j(\delta _{t-1}(j)\times a_{ji}\times b_{io_t})$$

我们假设前一个状态的局部概率已知，同时利用状态转移概率和对应的观察概率，就可以从中选择当前状态的最大概率了（局部概率 $\delta$ ）。

使用反向指针记录最佳路径

考虑上图的网格，经过前面的过程我们已经得到了每个状态的局部概率，但我们最终的目标是希望得到网格中的最佳隐藏序列，也就是说，最终的目标是需要得到网格中的局部最佳路径。

回顾之前我们得到局部概率的过程，$t$ 时刻的局部概率是通过 $t-1$ 时刻的局部概率得到的，也就是说，在计算得到 $t$ 时刻的局部概率 $\delta _t(i)$ 之后，我们就可以知道这个局部概率 $\delta _t(i)$ 是由 $t-1$ 时刻的哪一个状态而得到的。因此，在这个过程中，我们通过给每一个状态赋予一个反向指针 $\psi$ 来记录，这个指针指向导致当前状态最优概率的前一时刻的某个状态。

反向指针在形式上的公式如下：

$$\psi _t(i)=argmax_j(\delta _{t-1}(j)a_{ji})$$

这个反向指针的表达式是通过转移概率（某时刻的隐藏状态，演变到下一时刻某个隐藏状态的概率）以及前一时刻的局部概率计算得到的，没有使用到观察概率（隐藏状态到观察状态的发射概率）。

经过上述过程，已经可以通过递归的方式得到网格中每个状态的最佳局部概率，以及相应的局部最佳路径。

在使用Viterbi算法的过程中，我们对于网格中的每一个状态都计算一个局部概率，同时包含一个反向指针来指向最可能到达该状态的路径。当完成整个计算过程之后，我们首先在终止状态找到最可能的状态，之后通过反向指针回溯到初始时刻，从而回溯路径对应的状态序列就是我们最终得到的最佳的隐藏状态序列了。

Viterbi算法形式化表达及计算流程

上文对Viterbi算法的流程做了详细的梳理，现在来总结一下Viterbi算法的形式化定义及完整流程：

• 输入：HMM模型 $\lambda =(A,B,\pi)$ ，观察序列 $O=(o_1,o_2,…o_T)$ ；

• 输出：最有可能的隐藏状态序列 $I^*=\lbrace i_1^*,i_2^*,…i_T^*\rbrace$ ；

• 算法流程：

  • 初始化局部状态（$t=1$）：

    $$\delta _1(i)=\pi _ib_{io_1}\ ,\ i=1,2,...N$$ $$\psi _1(i)=0\ ,\ i=1,2,...N$$

    其中， $N​$ 表示隐藏状态的长度，在拼音转汉字的过程中对应汉字的长度。

    这一步通过观察状态的初始概率，和对应的观察状态到隐藏状态的发射概率，得到初始时刻的局部概率。

  • 递推，通过动态规划递推时刻 $t=2,3,…T$ 时刻的局部状态（$t>1$）：

    $$\delta _t(i)=\max_{1\leq j\leq N}[\delta _{t-1}(j)a_{ji}]b_{io_t}\ ,\ i=1,2,...N$$ $$\psi _t(i)=\arg\max _{1\leq j\leq N}[\delta _{t-1}(j)a_{ji}]\ ,\ i=1,2,...N$$

    其中，$i$ 表示当前时刻隐藏序列的各个状态，$j$ 表示前一时刻隐藏序列的各个状态。

    当前状态节点的选取，是通过考虑所有的转移概率，包括当前的发射概率，以及前一时刻的局部概率，综合计算，最后记录最大值，同时也记录了最大概率对应的前一时刻的状态节点。

  • 终止：

    $$P^*=\max _{1\leq i\leq N}\delta _T(i)$$ $$i_T^*=\arg \max _{1\leq i\leq N}[\delta _T(i)]$$

    计算终止时刻 T 时刻最大的概率值$\delta _T(i)$，即为最佳隐藏状态序列出现的概率，计算终止时刻 T 时刻最大的$\psi _t(i)$，即为最佳的隐藏状态。

  • 回溯最优路径，对$t=T-1,T-2,…,1$：

    $$i_t^*=\psi _{t+1}(i_{t+1}^*)$$

    按照此过程回溯整个网格，回溯完成时，就得到了生成给定观察序列的最可能的隐藏状态序列 $I^*=(i_1^*,i_2^*,…i_T^*)$。

Viterbi算法代码实现

C++版

主要参考umdhmm代码实现。

HMM模型结构定义如下：

typedef struct {
	int N;		/* number of hidden states;  Q={1,2,...,N} */
	int M; 		/* number of observation symbols; V={1,2,...,M}*/
	double	**A;	/* A[1..N][1..N]. a[i][j] is the transition prob
			   of going from state i at time t to state j
			   at time t+1 */
	double	**B;	/* B[1..N][1..M]. b[j][k] is the probability of
			   of observing symbol k in state j */
	double	*pi;	/* pi[1..N] pi[i] is the initial state distribution. */
} HMM;

Viterbi算法主流程实现如下：

void Viterbi(HMM *phmm, int T, int *O, double **delta, int **psi, 
             int *q, double *pprob)
{
	int 	i, j;	/* state indices */
	int  	t;	/* time index */	

	int	maxvalind;
	double	maxval, val;

	/* 1. Initialization  */
	for (i = 1; i <= phmm->N; i++) {
		delta[1][i] = phmm->pi[i] * (phmm->B[i][O[1]]);
		psi[1][i] = 0;
	}	

	/* 2. Recursion */
	for (t = 2; t <= T; t++) {
		for (j = 1; j <= phmm->N; j++) {
			maxval = 0.0;
			maxvalind = 1;	
			for (i = 1; i <= phmm->N; i++) {
				val = delta[t-1][i]*(phmm->A[i][j]);
				if (val > maxval) {
					maxval = val;	
					maxvalind = i;	
				}
			}
			delta[t][j] = maxval*(phmm->B[j][O[t]]);
			psi[t][j] = maxvalind; 
		}
	}

	/* 3. Termination */
	*pprob = 0.0;
	q[T] = 1;
	for (i = 1; i <= phmm->N; i++) {
        if (delta[T][i] > *pprob) {
			*pprob = delta[T][i];	
			q[T] = i;
		}
	}

	/* 4. Path (state sequence) backtracking */
	for (t = T - 1; t >= 1; t--)
		q[t] = psi[t+1][q[t+1]];
}

基于HMM实现拼音输入法主流程

在拼音输入法的主要过程中，HMM模型可以通过以下过程得到：

• 对训练数据的所有内容按照单字分词，并统计每个词出现的频率，以此作为初始概率 $\pi​$ ;
• 将训练数据的所有汉字都转换成对应的拼音，统计每个拼音对应的汉字以及各自出现的频率，以此作为发射概率 $B​$ ;
• 统计训练数据中每个汉字后面出现的汉字的频率，以此作为隐藏状态的转移概率 $A$ ;

经过上述过程，即可得到一个隐马尔科夫模型 $\lambda(\pi,A,B)$ 。

在此基础上，使用Viterbi算法，对用户输入的拼音序列进行解码，即可得到最终的汉字序列。

在实际应用过程中，对于一些转移概率为0的情况，可以默认设置一个较小的转移概率来避免计算结果为0的情况；同时，对于未登录词，也可以通过设置较小的初始概率来解决概率为0的问题。

¹. http://www.52nlp.cn/hmm-learn-best-practices-four-hidden-markov-models ↩