一、隐马尔科夫模型HMM

隐马尔科夫模型HMM（一）HMM模型基础

隐马尔科夫模型（Hidden Markov Model，以下简称HMM）是比较经典的机器学习模型了，它在语言识别，自然语言处理，模式识别等领域得到广泛的应用。当然，随着目前深度学习的崛起，尤其是RNN，LSTM等神经网络序列模型的火热，HMM的地位有所下降。但是作为一个经典的模型，学习HMM的模型和对应算法，对我们解决问题建模的能力提高以及算法思路的拓展还是很好的。本文是HMM系列的第一篇，关注于HMM模型的基础。

1. 什么样的问题需要HMM模型

　　　　首先我们来看看什么样的问题解决可以用HMM模型。使用HMM模型时我们的问题一般有这两个特征：１）我们的问题是基于序列的，比如时间序列，或者状态序列。２）我们的问题中有两类数据，一类序列数据是可以观测到的，即观测序列；而另一类数据是不能观察到的，即隐藏状态序列，简称状态序列。

　　　　有了这两个特征，那么这个问题一般可以用HMM模型来尝试解决。这样的问题在实际生活中是很多的。比如：我现在在打字写博客，我在键盘上敲出来的一系列字符就是观测序列，而我实际想写的一段话就是隐藏序列，输入法的任务就是从敲入的一系列字符尽可能的猜测我要写的一段话，并把最可能的词语放在最前面让我选择，这就可以看做一个HMM模型了。再举一个，我在和你说话，我发出的一串连续的声音就是观测序列，而我实际要表达的一段话就是状态序列，你大脑的任务，就是从这一串连续的声音中判断出我最可能要表达的话的内容。

　　　　从这些例子中，我们可以发现，HMM模型可以无处不在。但是上面的描述还不精确，下面我们用精确的数学符号来表述我们的HMM模型。

2. HMM模型的定义

　　　　对于HMM模型，首先我们假设QQ是所有可能的隐藏状态的集合，VV是所有可能的观测状态的集合，即：

Q={q1,q2,...,qN},V={v1,v2,...vM}Q={q1,q2,...,qN},V={v1,v2,...vM}

　　　　其中，NN是可能的隐藏状态数，MM是所有的可能的观察状态数。

　　　　对于一个长度为TT的序列，II对应的状态序列, OO是对应的观察序列，即：

I={i1,i2,...,iT},O={o1,o2,...oT}I={i1,i2,...,iT},O={o1,o2,...oT}

　　　　其中，任意一个隐藏状态it∈Qit∈Q,任意一个观察状态ot∈Vot∈V

　　　　HMM模型做了两个很重要的假设如下：

　　　　1） 齐次马尔科夫链假设。即任意时刻的隐藏状态只依赖于它前一个隐藏状态，这个我们在MCMC(二)马尔科夫链中有详细讲述。当然这样假设有点极端，因为很多时候我们的某一个隐藏状态不仅仅只依赖于前一个隐藏状态，可能是前两个或者是前三个。但是这样假设的好处就是模型简单，便于求解。如果在时刻tt的隐藏状态是it=qiit=qi,在时刻t+1t+1的隐藏状态是it+1=qjit+1=qj, 则从时刻tt到时刻t+1t+1的HMM状态转移概率aijaij可以表示为：

aij=P(it+1=qj|it=qi)aij=P(it+1=qj|it=qi)

 　　　　这样aijaij可以组成马尔科夫链的状态转移矩阵AA:

A=[aij]N×NA=[aij]N×N

　　　　2） 观测独立性假设。即任意时刻的观察状态只仅仅依赖于当前时刻的隐藏状态，这也是一个为了简化模型的假设。如果在时刻tt的隐藏状态是it=qjit=qj, 而对应的观察状态为ot=vkot=vk, 则该时刻观察状态vkvk在隐藏状态qjqj下生成的概率为bj(k)bj(k),满足：

bj(k)=P(ot=vk|it=qj)bj(k)=P(ot=vk|it=qj)

　　　　这样bj(k)bj(k)可以组成观测状态生成的概率矩阵BB:

B=[bj(k)]N×MB=[bj(k)]N×M

　　　　除此之外，我们需要一组在时刻t=1t=1的隐藏状态概率分布ΠΠ:

Π=[π(i)]N其中π(i)=P(i1=qi)Π=[π(i)]N其中π(i)=P(i1=qi)

　　　　一个HMM模型，可以由隐藏状态初始概率分布ΠΠ, 状态转移概率矩阵AA和观测状态概率矩阵BB决定。Π,AΠ,A决定状态序列，BB决定观测序列。因此，HMM模型可以由一个三元组λλ表示如下：

λ=(A,B,Π)λ=(A,B,Π)

3.一个HMM模型实例

　　　　下面我们用一个简单的实例来描述上面抽象出的HMM模型。这是一个盒子与球的模型，例子来源于李航的《统计学习方法》。

　　　　假设我们有3个盒子，每个盒子里都有红色和白色两种球，这三个盒子里球的数量分别是：

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

　　　　按照下面的方法从盒子里抽球，开始的时候，从第一个盒子抽球的概率是0.2，从第二个盒子抽球的概率是0.4，从第三个盒子抽球的概率是0.4。以这个概率抽一次球后，将球放回。然后从当前盒子转移到下一个盒子进行抽球。规则是：如果当前抽球的盒子是第一个盒子，则以0.5的概率仍然留在第一个盒子继续抽球，以0.2的概率去第二个盒子抽球，以0.3的概率去第三个盒子抽球。如果当前抽球的盒子是第二个盒子，则以0.5的概率仍然留在第二个盒子继续抽球，以0.3的概率去第一个盒子抽球，以0.2的概率去第三个盒子抽球。如果当前抽球的盒子是第三个盒子，则以0.5的概率仍然留在第三个盒子继续抽球，以0.2的概率去第一个盒子抽球，以0.3的概率去第二个盒子抽球。如此下去，直到重复三次，得到一个球的颜色的观测序列:

O={红，白，红}O={红，白，红}

　　　　 注意在这个过程中，观察者只能看到球的颜色序列，却不能看到球是从哪个盒子里取出的。

　　　　那么按照我们上一节HMM模型的定义，我们的观察集合是:

V={红，白}，M=2V={红，白}，M=2

　　　　我们的状态集合是：

Q={盒子1，盒子2，盒子3}，N=3Q={盒子1，盒子2，盒子3}，N=3

　　　　而观察序列和状态序列的长度为3.

　　　　初始状态分布为：

Π=(0.2,0.4,0.4)TΠ=(0.2,0.4,0.4)T

　　　　状态转移概率分布矩阵为：

 

A=⎛⎝⎜0.50.30.20.20.50.30.30.20.5⎞⎠⎟A=(0.50.20.30.30.50.20.20.30.5)

 　　　　观测状态概率矩阵为：

 

B=⎛⎝⎜0.50.40.70.50.60.3⎞⎠⎟B=(0.50.50.40.60.70.3)

4. HMM观测序列的生成

　　　　从上一节的例子，我们也可以抽象出HMM观测序列生成的过程。

　　　　输入的是HMM的模型λ=(A,B,Π)λ=(A,B,Π),观测序列的长度TT

　　　　输出是观测序列O={o1,o2,...oT}O={o1,o2,...oT}

　　　　生成的过程如下：

　　　　1）根据初始状态概率分布ΠΠ生成隐藏状态i1i1

　　　　2) for t from 1 to T

　　　　　　a. 按照隐藏状态itit的观测状态分布bit(k)bit(k)生成观察状态otot

　　　　　　b. 按照隐藏状态itit的状态转移概率分布aitit+1aitit+1产生隐藏状态it+1it+1

　　　　所有的otot一起形成观测序列O={o1,o2,...oT}O={o1,o2,...oT}

5. HMM模型的三个基本问题

 　　　　HMM模型一共有三个经典的问题需要解决：

　　　　1） 评估观察序列概率。即给定模型λ=(A,B,Π)λ=(A,B,Π)和观测序列O={o1,o2,...oT}O={o1,o2,...oT}，计算在模型λλ下观测序列OO出现的概率P(O|λ)P(O|λ)。这个问题的求解需要用到前向后向算法，我们在这个系列的第二篇会详细讲解。这个问题是HMM模型三个问题中最简单的。

　　　　2）模型参数学习问题。即给定观测序列O={o1,o2,...oT}O={o1,o2,...oT}，估计模型λ=(A,B,Π)λ=(A,B,Π)的参数，使该模型下观测序列的条件概率P(O|λ)P(O|λ)最大。这个问题的求解需要用到基于EM算法的鲍姆-韦尔奇算法， 我们在这个系列的第三篇会详细讲解。这个问题是HMM模型三个问题中最复杂的。

　　　　3）预测问题，也称为解码问题。即给定模型λ=(A,B,Π)λ=(A,B,Π)和观测序列O={o1,o2,...oT}O={o1,o2,...oT}，求给定观测序列条件下，最可能出现的对应的状态序列，这个问题的求解需要用到基于动态规划的维特比算法，我们在这个系列的第四篇会详细讲解。这个问题是HMM模型三个问题中复杂度居中的算法。

隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率

在隐马尔科夫模型HMM（一）HMM模型中，我们讲到了HMM模型的基础知识和HMM的三个基本问题，本篇我们就关注于HMM第一个基本问题的解决方法，即已知模型和观测序列，求观测序列出现的概率。

1. 回顾HMM问题一：求观测序列的概率

　　　　首先我们回顾下HMM模型的问题一。这个问题是这样的。我们已知HMM模型的参数λ=(A,B,Π)λ=(A,B,Π)。其中AA是隐藏状态转移概率的矩阵，BB是观测状态生成概率的矩阵， ΠΠ是隐藏状态的初始概率分布。同时我们也已经得到了观测序列O={o1,o2,...oT}O={o1,o2,...oT},现在我们要求观测序列OO在模型λλ下出现的条件概率P(O|λ)P(O|λ)。

　　　　乍一看，这个问题很简单。因为我们知道所有的隐藏状态之间的转移概率和所有从隐藏状态到观测状态生成概率，那么我们是可以暴力求解的。

　　　　我们可以列举出所有可能出现的长度为TT的隐藏序列I={i1,i2,...,iT}I={i1,i2,...,iT},分布求出这些隐藏序列与观测序列O={o1,o2,...oT}O={o1,o2,...oT}的联合概率分布P(O,I|λ)P(O,I|λ)，这样我们就可以很容易的求出边缘分布P(O|λ)P(O|λ)了。

　　　　具体暴力求解的方法是这样的：首先，任意一个隐藏序列I={i1,i2,...,iT}I={i1,i2,...,iT}出现的概率是：

P(I|λ)=πi1ai1i2ai2i3...aiT−1iTP(I|λ)=πi1ai1i2ai2i3...aiT−1iT

　　　　对于固定的状态序列I={i1,i2,...,iT}I={i1,i2,...,iT}，我们要求的观察序列O={o1,o2,...oT}O={o1,o2,...oT}出现的概率是：

P(O|I,λ)=bi1(o1)bi2(o2)...biT(oT)P(O|I,λ)=bi1(o1)bi2(o2)...biT(oT)

　　　　则OO和II联合出现的概率是：

P(O,I|λ)=P(I|λ)P(O|I,λ)=πi1bi1(o1)ai1i2bi2(o2)...aiT−1iTbiT(oT)P(O,I|λ)=P(I|λ)P(O|I,λ)=πi1bi1(o1)ai1i2bi2(o2)...aiT−1iTbiT(oT)

　　　　然后求边缘概率分布，即可得到观测序列OO在模型λλ下出现的条件概率P(O|λ)P(O|λ)：

P(O|λ)=∑IP(O,I|λ)=∑i1,i2,...iTπi1bi1(o1)ai1i2bi2(o2)...aiT−1iTbiT(oT)P(O|λ)=∑IP(O,I|λ)=∑i1,i2,...iTπi1bi1(o1)ai1i2bi2(o2)...aiT−1iTbiT(oT)

　　　　虽然上述方法有效，但是如果我们的隐藏状态数NN非常多的那就麻烦了，此时我们预测状态有NTNT种组合，算法的时间复杂度是O(TNT)O(TNT)阶的。因此对于一些隐藏状态数极少的模型，我们可以用暴力求解法来得到观测序列出现的概率，但是如果隐藏状态多，则上述算法太耗时，我们需要寻找其他简洁的算法。

　　　　前向后向算法就是来帮助我们在较低的时间复杂度情况下求解这个问题的。

2. 用前向算法求HMM观测序列的概率

　　　　前向后向算法是前向算法和后向算法的统称，这两个算法都可以用来求HMM观测序列的概率。我们先来看看前向算法是如何求解这个问题的。

　　　　前向算法本质上属于动态规划的算法，也就是我们要通过找到局部状态递推的公式，这样一步步的从子问题的最优解拓展到整个问题的最优解。

　　　　在前向算法中，通过定义“前向概率”来定义动态规划的这个局部状态。什么是前向概率呢, 其实定义很简单：定义时刻tt时隐藏状态为qiqi, 观测状态的序列为o1,o2,...oto1,o2,...ot的概率为前向概率。记为：

αt(i)=P(o1,o2,...ot,it=qi|λ)αt(i)=P(o1,o2,...ot,it=qi|λ)

　　　　既然是动态规划，我们就要递推了，现在我们假设我们已经找到了在时刻tt时各个隐藏状态的前向概率，现在我们需要递推出时刻t+1t+1时各个隐藏状态的前向概率。

　　　　从下图可以看出，我们可以基于时刻tt时各个隐藏状态的前向概率，再乘以对应的状态转移概率，即αt(j)ajiαt(j)aji就是在时刻tt观测到o1,o2,...oto1,o2,...ot，并且时刻tt隐藏状态qjqj, 时刻t+1t+1隐藏状态qiqi的概率。如果将想下面所有的线对应的概率求和，即∑j=1Nαt(j)aji∑j=1Nαt(j)aji就是在时刻tt观测到o1,o2,...oto1,o2,...ot，并且时刻t+1t+1隐藏状态qiqi的概率。继续一步，由于观测状态ot+1ot+1只依赖于t+1t+1时刻隐藏状态qiqi, 这样[∑j=1Nαt(j)aji]bi(ot+1)[∑j=1Nαt(j)aji]bi(ot+1)就是在在时刻t+1t+1观测到o1,o2,...ot，ot+1o1,o2,...ot，ot+1，并且时刻t+1t+1隐藏状态qiqi的概率。而这个概率，恰恰就是时刻t+1t+1对应的隐藏状态ii的前向概率，这样我们得到了前向概率的递推关系式如下：

αt+1(i)=[∑j=1Nαt(j)aji]bi(ot+1)αt+1(i)=[∑j=1Nαt(j)aji]bi(ot+1)

　　　　我们的动态规划从时刻1开始，到时刻TT结束，由于αT(i)αT(i)表示在时刻TT观测序列为o1,o2,...oTo1,o2,...oT，并且时刻TT隐藏状态qiqi的概率，我们只要将所有隐藏状态对应的概率相加，即∑i=1NαT(i)∑i=1NαT(i)就得到了在时刻TT观测序列为o1,o2,...oTo1,o2,...oT的概率。

　　　　下面总结下前向算法。

　　　　输入：HMM模型λ=(A,B,Π)λ=(A,B,Π)，观测序列O=(o1,o2,...oT)O=(o1,o2,...oT)

　　　　输出：观测序列概率P(O|λ)P(O|λ)

　　　　1) 计算时刻1的各个隐藏状态前向概率：

α1(i)=πibi(o1),i=1,2,...Nα1(i)=πibi(o1),i=1,2,...N

　　　　2) 递推时刻2,3,...T2,3,...T时刻的前向概率：

αt+1(i)=[∑j=1Nαt(j)aji]bi(ot+1),i=1,2,...Nαt+1(i)=[∑j=1Nαt(j)aji]bi(ot+1),i=1,2,...N

　　　　3) 计算最终结果：

P(O|λ)=∑i=1NαT(i)P(O|λ)=∑i=1NαT(i)

　　　　从递推公式可以看出，我们的算法时间复杂度是O(TN2)O(TN2)，比暴力解法的时间复杂度O(TNT)O(TNT)少了几个数量级。

3. HMM前向算法求解实例

　　　　这里我们用隐马尔科夫模型HMM（一）HMM模型中盒子与球的例子来显示前向概率的计算。

　　　　我们的观察集合是:

V={红，白}，M=2V={红，白}，M=2

　　　　我们的状态集合是：

Q={盒子1，盒子2，盒子3}，N=3Q={盒子1，盒子2，盒子3}，N=3

　　　　而观察序列和状态序列的长度为3.

　　　　初始状态分布为：

Π=(0.2,0.4,0.4)TΠ=(0.2,0.4,0.4)T

　　　　状态转移概率分布矩阵为：

 

A=⎛⎝⎜0.50.30.20.20.50.30.30.20.5⎞⎠⎟A=(0.50.20.30.30.50.20.20.30.5)

 　　　　观测状态概率矩阵为：

 

B=⎛⎝⎜0.50.40.70.50.60.3⎞⎠⎟B=(0.50.50.40.60.70.3)

　　　　球的颜色的观测序列:

O={红，白，红}O={红，白，红}

　　　　按照我们上一节的前向算法。首先计算时刻1三个状态的前向概率：

　　　　时刻1是红色球，隐藏状态是盒子1的概率为：

α1(1)=π1b1(o1)=0.2×0.5=0.1α1(1)=π1b1(o1)=0.2×0.5=0.1

　　　　隐藏状态是盒子2的概率为：

α1(2)=π2b2(o1)=0.4×0.4=0.16α1(2)=π2b2(o1)=0.4×0.4=0.16

　　　　隐藏状态是盒子3的概率为：

α1(3)=π3b3(o1)=0.4×0.7=0.28α1(3)=π3b3(o1)=0.4×0.7=0.28

　　　　现在我们可以开始递推了，首先递推时刻2三个状态的前向概率：

　　　　时刻2是白色球，隐藏状态是盒子1的概率为：

α2(1)=[∑i=13α1(i)ai1]b1(o2)=[0.1∗0.5+0.16∗0.3+0.28∗0.2]×0.5=0.077α2(1)=[∑i=13α1(i)ai1]b1(o2)=[0.1∗0.5+0.16∗0.3+0.28∗0.2]×0.5=0.077

　　　　隐藏状态是盒子2的概率为：

α2(2)=[∑i=13α1(i)ai2]b2(o2)=[0.1∗0.2+0.16∗0.5+0.28∗0.3]×0.6=0.1104α2(2)=[∑i=13α1(i)ai2]b2(o2)=[0.1∗0.2+0.16∗0.5+0.28∗0.3]×0.6=0.1104

　　　　隐藏状态是盒子3的概率为：

α2(3)=[∑i=13α1(i)ai3]b3(o2)=[0.1∗0.3+0.16∗0.2+0.28∗0.5]×0.3=0.0606α2(3)=[∑i=13α1(i)ai3]b3(o2)=[0.1∗0.3+0.16∗0.2+0.28∗0.5]×0.3=0.0606

　　　　继续递推，现在我们递推时刻3三个状态的前向概率：

　　　　时刻3是红色球，隐藏状态是盒子1的概率为：

α3(1)=[∑i=13α2(i)ai1]b1(o3)=[0.077∗0.5+0.1104∗0.3+0.0606∗0.2]×0.5=0.04187α3(1)=[∑i=13α2(i)ai1]b1(o3)=[0.077∗0.5+0.1104∗0.3+0.0606∗0.2]×0.5=0.04187

　　　　隐藏状态是盒子2的概率为：

α3(2)=[∑i=13α2(i)ai2]b2(o3)=[0.077∗0.2+0.1104∗0.5+0.0606∗0.3]×0.4=0.03551α3(2)=[∑i=13α2(i)ai2]b2(o3)=[0.077∗0.2+0.1104∗0.5+0.0606∗0.3]×0.4=0.03551

　　　　隐藏状态是盒子3的概率为：

α3(3)=[∑i=13α3(i)ai3]b3(o3)=[0.077∗0.3+0.1104∗0.2+0.0606∗0.5]×0.7=0.05284α3(3)=[∑i=13α3(i)ai3]b3(o3)=[0.077∗0.3+0.1104∗0.2+0.0606∗0.5]×0.7=0.05284

　　　　最终我们求出观测序列:O={红，白，红}O={红，白，红}的概率为：

P(O|λ)=∑i=13α3(i)=0.13022P(O|λ)=∑i=13α3(i)=0.13022

4. 用后向算法求HMM观测序列的概率

　　　　熟悉了用前向算法求HMM观测序列的概率，现在我们再来看看怎么用后向算法求HMM观测序列的概率。

　　　　后向算法和前向算法非常类似，都是用的动态规划，唯一的区别是选择的局部状态不同，后向算法用的是“后向概率”，那么后向概率是如何定义的呢？

　　　　定义时刻tt时隐藏状态为qiqi, 从时刻t+1t+1到最后时刻TT的观测状态的序列为ot+1,ot+2,...oTot+1,ot+2,...oT的概率为后向概率。记为：

βt(i)=P(ot+1,ot+2,...oT|it=qi,λ)βt(i)=P(ot+1,ot+2,...oT|it=qi,λ)

　　　　后向概率的动态规划递推公式和前向概率是相反的。现在我们假设我们已经找到了在时刻t+1t+1时各个隐藏状态的后向概率βt+1(j)βt+1(j)，现在我们需要递推出时刻tt时各个隐藏状态的后向概率。如下图，我们可以计算出观测状态的序列为ot+2,ot+3,...oTot+2,ot+3,...oT， tt时隐藏状态为qiqi, 时刻t+1t+1隐藏状态为qjqj的概率为aijβt+1(j)aijβt+1(j), 接着可以得到观测状态的序列为ot+1,ot+2,...oTot+1,ot+2,...oT， tt时隐藏状态为qiqi, 时刻t+1t+1隐藏状态为qjqj的概率为aijbj(ot+1)βt+1(j)aijbj(ot+1)βt+1(j), 则把下面所有线对应的概率加起来，我们可以得到观测状态的序列为ot+1,ot+2,...oTot+1,ot+2,...oT， tt时隐藏状态为qiqi的概率为∑j=1Naijbj(ot+1)βt+1(j)∑j=1Naijbj(ot+1)βt+1(j)，这个概率即为时刻tt的后向概率。

　　　　这样我们得到了后向概率的递推关系式如下：

βt(i)=∑j=1Naijbj(ot+1)βt+1(j)βt(i)=∑j=1Naijbj(ot+1)βt+1(j)

　　　　现在我们总结下后向算法的流程,注意下和前向算法的相同点和不同点：

　　　　输入：HMM模型λ=(A,B,Π)λ=(A,B,Π)，观测序列O=(o1,o2,...oT)O=(o1,o2,...oT)

　　　　输出：观测序列概率P(O|λ)P(O|λ)

　　　　1) 初始化时刻TT的各个隐藏状态后向概率：

βT(i)=1,i=1,2,...NβT(i)=1,i=1,2,...N

　　　　2) 递推时刻T−1,T−2,...1T−1,T−2,...1时刻的后向概率：

βt(i)=∑j=1Naijbj(ot+1)βt+1(j),i=1,2,...Nβt(i)=∑j=1Naijbj(ot+1)βt+1(j),i=1,2,...N

　　　　3) 计算最终结果：

P(O|λ)=∑i=1Nπibi(o1)β1(i)P(O|λ)=∑i=1Nπibi(o1)β1(i)

　　　　此时我们的算法时间复杂度仍然是O(TN2)O(TN2)。

5. HMM常用概率的计算

　　　　利用前向概率和后向概率，我们可以计算出HMM中单个状态和两个状态的概率公式。

　　　　1）给定模型λλ和观测序列OO,在时刻tt处于状态qiqi的概率记为:

γt(i)=P(it=qi|O,λ)=P(it=qi,O|λ)P(O|λ)γt(i)=P(it=qi|O,λ)=P(it=qi,O|λ)P(O|λ)

　　　　利用前向概率和后向概率的定义可知：

P(it=qi,O|λ)=αt(i)βt(i)P(it=qi,O|λ)=αt(i)βt(i)

　　　　于是我们得到：

γt(i)=αt(i)βt(i)∑j=1Nαt(j)βt(j)γt(i)=αt(i)βt(i)∑j=1Nαt(j)βt(j)

　　　　2）给定模型λλ和观测序列OO,在时刻tt处于状态qiqi，且时刻t+1t+1处于状态qjqj的概率记为:

ξt(i,j)=P(it=qi,it+1=qj|O,λ)=P(it=qi,it+1=qj,O|λ)P(O|λ)ξt(i,j)=P(it=qi,it+1=qj|O,λ)=P(it=qi,it+1=qj,O|λ)P(O|λ)

　　　　而P(it=qi,it+1=qj,O|λ)P(it=qi,it+1=qj,O|λ)可以由前向后向概率来表示为:

P(it=qi,it+1=qj,O|λ)=αt(i)aijbj(ot+1)βt+1(j)P(it=qi,it+1=qj,O|λ)=αt(i)aijbj(ot+1)βt+1(j)

　　　　从而最终我们得到ξt(i,j)ξt(i,j)的表达式如下：

ξt(i,j)=αt(i)aijbj(ot+1)βt+1(j)∑r=1N∑s=1Nαt(r)arsbs(ot+1)βt+1(s)ξt(i,j)=αt(i)aijbj(ot+1)βt+1(j)∑r=1N∑s=1Nαt(r)arsbs(ot+1)βt+1(s)

 　　　　3) 将γt(i)γt(i)和ξt(i,j)ξt(i,j)在各个时刻tt求和，可以得到：

　　　　在观测序列OO下状态ii出现的期望值∑t=1Tγt(i)∑t=1Tγt(i)

　　　　在观测序列OO下由状态ii转移的期望值∑t=1T−1γt(i)∑t=1T−1γt(i)

　　　　在观测序列OO下由状态ii转移到状态jj的期望值∑t=1T−1ξt(i,j)∑t=1T−1ξt(i,j)

　　　　上面这些常用的概率值在求解HMM问题二，即求解HMM模型参数的时候需要用到。我们在这个系列的第三篇来讨论求解HMM参数的问题和解法。

隐马尔科夫模型HMM（三）鲍姆-韦尔奇算法求解HMM参数

在本篇我们会讨论HMM模型参数求解的问题，这个问题在HMM三个问题里算是最复杂的。在研究这个问题之前，建议先阅读这个系列的前两篇以熟悉HMM模型和HMM的前向后向算法，以及EM算法原理总结，这些在本篇里会用到。在李航的《统计学习方法》中，这个算法的讲解只考虑了单个观测序列的求解，因此无法用于实际多样本观测序列的模型求解，本文关注于如何使用多个观测序列来求解HMM模型参数。

1. HMM模型参数求解概述

　　　　HMM模型参数求解根据已知的条件可以分为两种情况。

　　　　第一种情况较为简单，就是我们已知D个长度为T的观测序列和对应的隐藏状态序列，即{(O1,I1),(O2,I2),...(OD,ID)}是已知的，此时我们可以很容易的用最大似然来求解模型参数。

　　　　假设样本从隐藏状态qi转移到qj的频率计数是Aij,那么状态转移矩阵求得为：A=[aij],其中aij=AijN∑s=1Ais

　　　　假设样本隐藏状态为qj且观测状态为vk的频率计数是Bjk,那么观测状态概率矩阵为：B=[bj(k)],其中bj(k)=BjkM∑s=1Bjs

　　　　假设所有样本中初始隐藏状态为qi的频率计数为C(i),那么初始概率分布为：Π=π(i)=C(i)N∑s=1C(s)

　　　　可见第一种情况下求解模型还是很简单的。但是在很多时候，我们无法得到HMM样本观察序列对应的隐藏序列，只有D个长度为T的观测序列，即{(O1),(O2),...(OD)}是已知的，此时我们能不能求出合适的HMM模型参数呢？这就是我们的第二种情况，也是我们本文要讨论的重点。它的解法最常用的是鲍姆-韦尔奇算法，其实就是基于EM算法的求解，只不过鲍姆-韦尔奇算法出现的时代，EM算法还没有被抽象出来，所以我们本文还是说鲍姆-韦尔奇算法。

2. 鲍姆-韦尔奇算法原理

　　　　鲍姆-韦尔奇算法原理既然使用的就是EM算法的原理，那么我们需要在E步求出联合分布P(O,I|λ)基于条件概率P(I|O,¯λ)的期望，其中¯λ为当前的模型参数，然后再M步最大化这个期望，得到更新的模型参数λ。接着不停的进行EM迭代，直到模型参数的值收敛为止。

　　　　首先来看看E步，当前模型参数为¯λ, 联合分布P(O,I|λ)基于条件概率P(I|O,¯λ)的期望表达式为：L(λ,¯λ)=∑IP(I|O,¯λ)logP(O,I|λ)

　　　　在M步，我们极大化上式，然后得到更新后的模型参数如下：　¯λ=argmaxλ∑IP(I|O,¯λ)logP(O,I|λ)

　　　　通过不断的E步和M步的迭代，直到¯λ收敛。下面我们来看看鲍姆-韦尔奇算法的推导过程。

3. 鲍姆-韦尔奇算法的推导

　　　　我们的训练数据为{(O1,I1),(O2,I2),...(OD,ID)}，其中任意一个观测序列Od={o(d)1,o(d)2,...o(d)T},其对应的未知的隐藏状态序列表示为：Id={i(d)1,i(d)2,...i(d)T}

　　　　首先看鲍姆-韦尔奇算法的E步，我们需要先计算联合分布P(O,I|λ)的表达式如下：P(O,I|λ)=D∏d=1πi(d)1bi(d)1(o(d)1)ai(d)1i(d)2bi(d)2(o(d)2)...ai(d)T−1i(d)Tbi(d)T(o(d)T)

　　　　我们的E步得到的期望表达式为：L(λ,¯λ)=∑IP(I|O,¯λ)logP(O,I|λ)

　　　　在M步我们要极大化上式。由于P(I|O,¯λ)=P(I,O|¯λ)/P(O|¯λ),而P(O|¯λ)是常数，因此我们要极大化的式子等价于：¯λ=argmaxλ∑IP(O,I|¯λ)logP(O,I|λ)

　　　　我们将上面P(O,I|λ)的表达式带入我们的极大化式子，得到的表达式如下：¯λ=argmaxλD∑d=1∑IP(O,I|¯λ)(logπi1+T−1∑t=1logaitait+1+T∑t=1bit(ot))

　　　　我们的隐藏模型参数λ=(A,B,Π),因此下面我们只需要对上式分别对A,B,Π求导即可得到我们更新的模型参数¯λ　

 

　　　　首先我们看看对模型参数Π的求导。由于Π只在上式中括号里的第一部分出现，因此我们对于Π的极大化式子为：¯πi=argmaxπi1D∑d=1∑IP(O,I|¯λ)logπi1=argmaxπiD∑d=1N∑i=1P(O,i(d)1=i|¯λ)logπi

　　　　由于πi还满足N∑i=1πi=1，因此根据拉格朗日子乘法，我们得到πi要极大化的拉格朗日函数为：argmaxπiD∑d=1N∑i=1P(O,i(d)1=i|¯λ)logπi+γ(N∑i=1πi−1)

　　　　其中，γ为拉格朗日系数。上式对πi求偏导数并令结果为0， 我们得到：D∑d=1P(O,i(d)1=i|¯λ)+γπi=0

　　　　令i分别等于从1到N，从上式可以得到N个式子，对这N个式子求和可得：D∑d=1P(O|¯λ)+γ=0

　　　　从上两式消去γ,得到πi的表达式为：πi=D∑d=1P(O,i(d)1=i|¯λ)D∑d=1P(O|¯λ)=D∑d=1P(O,i(d)1=i|¯λ)DP(O|¯λ)=D∑d=1P(i(d)1=i|O,¯λ)D=D∑d=1P(i(d)1=i|O(d),¯λ)D

　　　　利用我们在隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率里第二节中前向概率的定义可得：P(i(d)1=i|O(d),¯λ)=γ(d)1(i)

　　　　因此最终我们在M步πi的迭代公式为：πi=D∑d=1γ(d)1(i)D

 

　　　　现在我们来看看A的迭代公式求法。方法和Π的类似。由于A只在最大化函数式中括号里的第二部分出现，而这部分式子可以整理为：D∑d=1∑IT−1∑t=1P(O,I|¯λ)logaitait+1=D∑d=1N∑i=1N∑j=1T−1∑t=1P(O,i(d)t=i,i(d)t+1=j|¯λ)logaij

　　　　由于aij还满足N∑j=1aij=1。和求解πi类似，我们可以用拉格朗日子乘法并对aij求导，并令结果为0，可以得到aij的迭代表达式为：aij=D∑d=1T−1∑t=1P(O(d),i(d)t=i,i(d)t+1=j|¯λ)D∑d=1T−1∑t=1P(O(d),i(d)t=i|¯λ)

　　　　利用隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率里第二节中前向概率的定义和第五节ξt(i,j)的定义可得们在M步aij的迭代公式为：aij=D∑d=1T−1∑t=1ξ(d)t(i,j)D∑d=1T−1∑t=1γ(d)t(i)

 

　　　　现在我们来看看B的迭代公式求法。方法和Π的类似。由于B只在最大化函数式中括号里的第三部分出现，而这部分式子可以整理为：D∑d=1∑IT∑t=1P(O,I|¯λ)logbit(ot)=D∑d=1N∑j=1T∑t=1P(O,i(d)t=j|¯λ)logbj(ot)

　　　　由于bj(ot)还满足M∑k=1bj(ot=vk)=1。和求解πi类似，我们可以用拉格朗日子乘法并对bj(k)求导，并令结果为0，得到bj(k)的迭代表达式为：bj(k)=D∑d=1T∑t=1P(O,i(d)t=j|¯λ)I(o(d)t=vk)D∑d=1T∑t=1P(O,i(d)t=j|¯λ)

　　　　其中I(o(d)t=vk)当且仅当o(d)t=vk时为1，否则为0. 利用隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率里第二节中前向概率的定义可得bj(ot)的最终表达式为：bj(k)=D∑d=1T∑t=1,o(d)t=vkγ(d)t(i)D∑d=1T∑t=1γ(d)t(i)

　　　　有了πi,aij,bj(k)的迭代公式，我们就可以迭代求解HMM模型参数了。

4. 鲍姆-韦尔奇算法流程总结

　　　　这里我们概括总结下鲍姆-韦尔奇算法的流程。

　　　　输入： D个观测序列样本{(O1),(O2),...(OD)}

　　　　输出：HMM模型参数

　　　　1)随机初始化所有的πi,aij,bj(k)

　　　　2) 对于每个样本d=1,2,...D，用前向后向算法计算γ(d)t(i)，ξ(d)t(i,j),t=1,2...T

　　　　3)  更新模型参数：

πi=D∑d=1γ(d)1(i)D

aij=D∑d=1T−1∑t=1ξ(d)t(i,j)D∑d=1T−1∑t=1γ(d)t(i)

bj(k)=D∑d=1T∑t=1,o(d)t=vkγ(d)t(i)D∑d=1T∑t=1γ(d)t(i)

　　　　4) 如果πi,aij,bj(k)的值已经收敛，则算法结束，否则回到第2）步继续迭代。

　　　　以上就是鲍姆-韦尔奇算法的整个过程。

隐马尔科夫模型HMM（四）维特比算法解码隐藏状态序列

在本篇我们会讨论HMM模型最后一个问题的求解，即即给定模型和观测序列，求给定观测序列条件下，最可能出现的对应的隐藏状态序列。在阅读本篇前，建议先阅读这个系列的第一篇以熟悉HMM模型。

　　　　HMM模型的解码问题最常用的算法是维特比算法，当然也有其他的算法可以求解这个问题。同时维特比算法是一个通用的求序列最短路径的动态规划算法，也可以用于很多其他问题，比如之前讲到的文本挖掘的分词原理中我们讲到了单独用维特比算法来做分词。

　　　　本文关注于用维特比算法来解码HMM的的最可能隐藏状态序列。

1. HMM最可能隐藏状态序列求解概述

　　　　在HMM模型的解码问题中，给定模型λ=(A,B,Π)λ=(A,B,Π)和观测序列O={o1,o2,...oT}O={o1,o2,...oT}，求给定观测序列O条件下，最可能出现的对应的状态序列I∗={i∗1,i∗2,...i∗T}I∗={i1∗,i2∗,...iT∗},即P(I∗|O)P(I∗|O)要最大化。

　　　　一个可能的近似解法是求出观测序列OO在每个时刻tt最可能的隐藏状态i∗tit∗然后得到一个近似的隐藏状态序列I∗={i∗1,i∗2,...i∗T}I∗={i1∗,i2∗,...iT∗}。要这样近似求解不难，利用隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率中第五节的定义：在给定模型λλ和观测序列OO时，在时刻tt处于状态qiqi的概率是γt(i)γt(i)，这个概率可以通过HMM的前向算法与后向算法计算。这样我们有：

i∗t=argmax1≤i≤N[γt(i)],t=1,2,...Tit∗=argmax1≤i≤N[γt(i)],t=1,2,...T

　　　　近似算法很简单，但是却不能保证预测的状态序列是整体是最可能的状态序列，因为预测的状态序列中某些相邻的隐藏状态可能存在转移概率为0的情况。

　　　　而维特比算法可以将HMM的状态序列作为一个整体来考虑，避免近似算法的问题，下面我们来看看维特比算法进行HMM解码的方法。

2. 维特比算法概述

　　　　维特比算法是一个通用的解码算法，是基于动态规划的求序列最短路径的方法。在文本挖掘的分词原理中我们已经讲到了维特比算法的一些细节。

　　　　既然是动态规划算法，那么就需要找到合适的局部状态，以及局部状态的递推公式。在HMM中，维特比算法定义了两个局部状态用于递推。

　　　　第一个局部状态是在时刻tt隐藏状态为ii所有可能的状态转移路径i1,i2,...iti1,i2,...it中的概率最大值。记为δt(i)δt(i):

δt(i)=maxi1,i2,...it−1P(it=i,i1,i2,...it−1,ot,ot−1,...o1|λ),i=1,2,...Nδt(i)=maxi1,i2,...it−1P(it=i,i1,i2,...it−1,ot,ot−1,...o1|λ),i=1,2,...N

　　　　由δt(i)δt(i)的定义可以得到δδ的递推表达式：

δt+1(i)=maxi1,i2,...itP(it+1=i,i1,i2,...it,ot+1,ot,...o1|λ)=max1≤j≤N[δt(j)aji]bi(ot+1)(1)(2)(1)δt+1(i)=maxi1,i2,...itP(it+1=i,i1,i2,...it,ot+1,ot,...o1|λ)(2)=max1≤j≤N[δt(j)aji]bi(ot+1)

　　　　第二个局部状态由第一个局部状态递推得到。我们定义在时刻tt隐藏状态为ii的所有单个状态转移路径(i1,i2,...,it−1,i)(i1,i2,...,it−1,i)中概率最大的转移路径中第t−1t−1个节点的隐藏状态为Ψt(i)Ψt(i),其递推表达式可以表示为：

Ψt(i)=argmax1≤j≤N[δt−1(j)aji]Ψt(i)=argmax1≤j≤N[δt−1(j)aji]

　　　　有了这两个局部状态，我们就可以从时刻0一直递推到时刻TT，然后利用Ψt(i)Ψt(i)记录的前一个最可能的状态节点回溯，直到找到最优的隐藏状态序列。

3. 维特比算法流程总结

　　　　现在我们来总结下维特比算法的流程：

　　　　输入：HMM模型λ=(A,B,Π)λ=(A,B,Π)，观测序列O=(o1,o2,...oT)O=(o1,o2,...oT)

　　　　输出：最有可能的隐藏状态序列I∗={i∗1,i∗2,...i∗T}I∗={i1∗,i2∗,...iT∗}

　　　　1）初始化局部状态：

δ1(i)=πibi(o1),i=1,2...Nδ1(i)=πibi(o1),i=1,2...N

Ψ1(i)=0,i=1,2...NΨ1(i)=0,i=1,2...N

　　　　2) 进行动态规划递推时刻t=2,3,...Tt=2,3,...T时刻的局部状态：

δt(i)=max1≤j≤N[δt−1(j)aji]bi(0t),i=1,2...Nδt(i)=max1≤j≤N[δt−1(j)aji]bi(0t),i=1,2...N

Ψt(i)=argmax1≤j≤N[δt−1(j)aji],i=1,2...NΨt(i)=argmax1≤j≤N[δt−1(j)aji],i=1,2...N

　　　　3) 计算时刻TT最大的δT(i)δT(i),即为最可能隐藏状态序列出现的概率。计算时刻TT最大的Ψt(i)Ψt(i),即为时刻TT最可能的隐藏状态。

P∗=max1≤j≤NδT(i)P∗=max1≤j≤NδT(i)

i∗T=argmax1≤j≤N[δT(i)]iT∗=argmax1≤j≤N[δT(i)]

　　　　4) 利用局部状态Ψ(i)Ψ(i)开始回溯。对于t=T−1,T−2,...,1t=T−1,T−2,...,1：

i∗t=Ψt+1(i∗t+1)it∗=Ψt+1(it+1∗)

　　　　最终得到最有可能的隐藏状态序列I∗={i∗1,i∗2,...i∗T}I∗={i1∗,i2∗,...iT∗}

4. HMM维特比算法求解实例

　　　　下面我们仍然用隐马尔科夫模型HMM（一）HMM模型中盒子与球的例子来看看HMM维特比算法求解。

　　　　我们的观察集合是:

V={红，白}，M=2V={红，白}，M=2

　　　　我们的状态集合是：

Q={盒子1，盒子2，盒子3}，N=3Q={盒子1，盒子2，盒子3}，N=3

　　　　而观察序列和状态序列的长度为3.

　　　　初始状态分布为：

Π=(0.2,0.4,0.4)TΠ=(0.2,0.4,0.4)T

　　　　状态转移概率分布矩阵为：

 

A=⎛⎝⎜0.50.30.20.20.50.30.30.20.5⎞⎠⎟A=(0.50.20.30.30.50.20.20.30.5)

 　　　　观测状态概率矩阵为：

 

B=⎛⎝⎜0.50.40.70.50.60.3⎞⎠⎟B=(0.50.50.40.60.70.3)

　　　　球的颜色的观测序列:

O={红，白，红}O={红，白，红}

　　　　按照我们上一节的维特比算法，首先需要得到三个隐藏状态在时刻1时对应的各自两个局部状态，此时观测状态为1：

 

δ1(1)=π1b1(o1)=0.2×0.5=0.1δ1(1)=π1b1(o1)=0.2×0.5=0.1

 

δ1(2)=π2b2(o1)=0.4×0.4=0.16δ1(2)=π2b2(o1)=0.4×0.4=0.16

 

δ1(3)=π3b3(o1)=0.4×0.7=0.28δ1(3)=π3b3(o1)=0.4×0.7=0.28

 

Ψ1(1)=Ψ1(2)=Ψ1(3)=0Ψ1(1)=Ψ1(2)=Ψ1(3)=0

　　　　现在开始递推三个隐藏状态在时刻2时对应的各自两个局部状态，此时观测状态为2：

 

δ2(1)=max1≤j≤3[δ1(j)aj1]b1(o2)=max1≤j≤3[0.1×0.5,0.16×0.3,0.28×0.2]×0.5=0.028δ2(1)=max1≤j≤3[δ1(j)aj1]b1(o2)=max1≤j≤3[0.1×0.5,0.16×0.3,0.28×0.2]×0.5=0.028

 

Ψ2(1)=3Ψ2(1)=3

 

δ2(2)=max1≤j≤3[δ1(j)aj2]b2(o2)=max1≤j≤3[0.1×0.2,0.16×0.5,0.28×0.3]×0.6=0.0504δ2(2)=max1≤j≤3[δ1(j)aj2]b2(o2)=max1≤j≤3[0.1×0.2,0.16×0.5,0.28×0.3]×0.6=0.0504

 

Ψ2(2)=3Ψ2(2)=3

 

δ2(3)=max1≤j≤3[δ1(j)aj3]b3(o2)=max1≤j≤3[0.1×0.3,0.16×0.2,0.28×0.5]×0.3=0.042δ2(3)=max1≤j≤3[δ1(j)aj3]b3(o2)=max1≤j≤3[0.1×0.3,0.16×0.2,0.28×0.5]×0.3=0.042

 

Ψ2(3)=3Ψ2(3)=3

　　　　继续递推三个隐藏状态在时刻3时对应的各自两个局部状态，此时观测状态为1：

 

δ3(1)=max1≤j≤3[δ2(j)aj1]b1(o3)=max1≤j≤3[0.028×0.5,0.0504×0.3,0.042×0.2]×0.5=0.00756δ3(1)=max1≤j≤3[δ2(j)aj1]b1(o3)=max1≤j≤3[0.028×0.5,0.0504×0.3,0.042×0.2]×0.5=0.00756

 

Ψ3(1)=2Ψ3(1)=2

 

δ3(2)=max1≤j≤3[δ2(j)aj2]b2(o3)=max1≤j≤3[0.028×0.2,0.0504×0.5,0.042×0.3]×0.4=0.01008δ3(2)=max1≤j≤3[δ2(j)aj2]b2(o3)=max1≤j≤3[0.028×0.2,0.0504×0.5,0.042×0.3]×0.4=0.01008

 

Ψ3(2)=2Ψ3(2)=2

 

δ3(3)=max1≤j≤3[δ2(j)aj3]b3(o3)=max1≤j≤3[0.028×0.3,0.0504×0.2,0.042×0.5]×0.7=0.0147δ3(3)=max1≤j≤3[δ2(j)aj3]b3(o3)=max1≤j≤3[0.028×0.3,0.0504×0.2,0.042×0.5]×0.7=0.0147

 

Ψ3(3)=3Ψ3(3)=3

　　　　此时已经到最后的时刻，我们开始准备回溯。此时最大概率为δ3(3)δ3(3),从而得到i∗3=3i3∗=3

　　　　由于Ψ3(3)=3Ψ3(3)=3,所以i∗2=3i2∗=3, 而又由于Ψ2(3)=3Ψ2(3)=3,所以i∗1=3i1∗=3。从而得到最终的最可能的隐藏状态序列为：(3,3,3)(3,3,3)

5. HMM模型维特比算法总结

　　　　如果大家看过之前写的文本挖掘的分词原理中的维特比算法，就会发现这两篇之中的维特比算法稍有不同。主要原因是在中文分词时，我们没有观察状态和隐藏状态的区别，只有一种状态。但是维特比算法的核心是定义动态规划的局部状态与局部递推公式，这一点在中文分词维特比算法和HMM的维特比算法是相同的，也是维特比算法的精华所在。

　　　　维特比算法也是寻找序列最短路径的一个通用方法，和dijkstra算法有些类似，但是dijkstra算法并没有使用动态规划，而是贪心算法。同时维特比算法仅仅局限于求序列最短路径，而dijkstra算法是通用的求最短路径的方法。