HMM(隐马尔科夫模型)知识点整理
在信息处理过程中,常常会遇到根据已有的序列信息推断源信息的问题,例如在自然语言处理的词性标注问题中,给定一句话,我们希望构建一个模型,能推测出这句话中每个词的词性。这个问题里,我们能够观测到的信息是单词序列,需要被推断的信息是词性序列。我们可以假设两个人在对话,尽管听者直接收到的信息是单词序列,但从说者头脑的思维转变成实际接收到的单词序列,发生了两个重要的过程。过程一,说者从他的语法库中按照一定概率分布生成了第一个词性;过程二,针对该词词性又按照一定概率生成了该词性下的第一个词汇。接下来重复过程一和二,也就变成了听者耳中的一句话。对上述过程建立概率图模型也就是隐马尔科夫模型。
一、一个简单的例子
上面的例子也许不够直观,这里举一个比较经典的例子(很多讲HMM的书或者博客都有提及)。研究概率问题,总可以想象成有一些人在玩摸球游戏。
例1.1 假设有3袋球,分别编号1,2,3。袋1中有1个红球,2个白球,袋2中有2个红球,1个白球,袋3中有3个白球。有一人甲先随机选择一个袋子,然后再从其中随机摸一个球,另一人乙可以看到甲所摸球的颜色,但无法知道其所选袋子的编号。
针对上面的例子,我们可以提出以下一些问题,
- 乙看到甲的摸球序列为“红,白,红”的概率为多少?
- 如果乙看到摸球序列为“红,白,红”,甲所选的编号序列最有可能是什么?
- 如果我进行大量的观测,能否探寻到甲选择袋子标号的规律(可能是心理因素,也可能毫无规律)?
我们不妨暂时先搁置这些问题,先对上面的例子做一些数学上的抽象分析,也就是数学建模。
二、一些概念回顾
2.1 隐变量和概率图
对于一些在随机过程中无法被直接观测到的变量,例如上面例子中的词性和袋子编号,我们称之为隐变量(hidden variable)。这些变量除非生成者告知我们,否则作为观测者无法得知,观测者唯一能做的就是推断。
接下来需要简单复习一下概率图模型的基本概念,首先回顾一下概率的链式法则,
p(x1,x2,…,xn)=p(x1)p(x2∣x1)…p(xn∣x1,…,xn−1)
联合概率分布表达了随机变量之间的全部信息,这些信息包括:
- 任意随机变量的自身属性;
- 任意两个随机变量之间的独立关系;
在数学上,图的概念刚好包含了上面两个信息:节点表征了变量自身的属性;边表征了变量与变量之间的关系。因此,很自然地联想到对一组随机变量引入图的概念,从而建立概率图模型。需要说明的是,这里,仅仅考虑有向无环图模型。
2.2 离散Markov模型
然后需要回顾的一个概念是离散Markov链,离散Markov链是一类简单的随机过程,这里不再给出严格的定义。简单来说,对于一个离散随机变量序列
x1,x2,…,xn,如果他们的联合概率分布满足下述条件,则称该随机变量序列为一个离散马尔科夫链:
p(x1,…,xn)=i=1∏np(xi∣xi−1),p(x1∣x0)=p(x1)
简单来说,就是第
n时刻的随机变量的概率分布仅与
n−1时刻相关。当然,这些随机变量有公共的状态集合
xi∈S={sk}∣k=1N。对离散Markov链建立概率图模型如下,
x1→x2→⋯→xn
直观理解就是一个不断延伸下去的线性链(如果不对
n做限制)。依附于离散Markov链的第一个重要概念是概率转移矩阵
P,该矩阵里的
(i,j)处元素为
aij=p(xn=sj∣xn−1=si)
故
P∈RN×N。这里的转移概率忽略了时间的影响,事实上存在一类时移Markov链,其概率转移矩阵是时间的函数。
接下来,自然要问,Markov链是怎样开始的呢?因而,我们需要定义Markov的初始分布
p(x1),
π=(π1,…,πn)。
三、HMM基本概念
3.1 概率图结构
有了上面的一些概念准备,我们就可以着手建立隐马尔科夫模型,我们可以假设例子中的过程一(其中的变量无法被观测,因而是隐变量)为一个离散Markov链,记作第
t时刻的隐变量为
qt,第
i时刻观测到的变量记作
ot,并且
ot仅与
qt存在条件相关性。于是可以得到下面的概率图
用黑色表示隐节点,空白节点表示可以被观测的节点。有了上面的概率图,模型的结构也就清晰起来。
3.2 模型参数
接下来的任务就是对概率图中的边引入参数。我们仍然需要引入一些符号,对于状态节点
qt,其状态空间为
S={sk}∣k=1N,对于观测节点,其状态空间为
V={vk}∣k=1M。
根据假设,隐节点序列是一个离散Markov链,于是我们对其引入转移矩阵
A=(aij)N×N和初始分布
π。接下来我们需要引入另外一个矩阵来表征概率
p(ot∣qt),即
B=(bij)N×M,称之为发射矩阵,其中的元素可以简单记作
bj(ot),表示概率
p(ot∣qt=sj)。
总结一下,一个HMM模型记作
μ=(A,B,π),其中,
-
A=(aij)N×N称为转移矩阵;
-
bj(ot),j=1,2,…,N称为发射概率;
-
π=(π1,…,πN)为初始分布;
通过引入概率图模型,我们就成功地将例1.1抽象化了。接下来,还有一步重要的工作是将例1.1中的问题用数学语言表达,也就是HMM模型中的三个基本问题,
- 给定模型
μ=(A,B,π)和观测序列
O=o1o2…oT,
p(O)为何值?
- 给定模型
μ=(A,B,π)和观测序列
O=o1o2…oT,最有可能的序列
Q=q1q2…qT是什么?
- 给定已知的序列
O=o1o2…oT和
Q=q1q2…qT,最有可能的模型
μ=(A,B,π)是什么?
接下来的三个算法分别解决了上述三个问题。
四、前向后向算法
4.1 野蛮推导
事实上,HMM的概率图模型表达的是联合概率
p(o1o2…oT,q1q2…qT)
问题1所要求的概率是
p(O∣μ),我们自然可以使用全概率公式和条件概率公式进行一系列的推导,最终得到如下的表达式:
p(O∣μ)=Q∈ST∑t=1∏Tbqt(ot)p(qt∣qt−1)=Q∈ST∑πq1bq1(o1)t=2∏Taqt−1qtbqt(ot)
4.2 前向算法
求解上式面临的最大问题是计算复杂度,需要进行约
O(TNT)次乘法和加法运算。这个计算量是巨大的,事实上,如果在序列里存在一定的递推关系,就可以可容易地降低计算复杂度。的确存在这样的递推关系,这也是前向后向算法(forward-backward algorithm)的基本思路。
假设某一个HMM发展到
t−1时刻,也就是说我们已经观测到了
o1,…,ot−1,那么,接下来会发生这样两个过程,
- 由
qt−1(尽管我们并不知道它是什么)转移到
qt;
- 由
qt发射出
ot;
也就是说,我们关心的是概率
p(o1o2…ot,qt=si∣μ),i=1,…,N
把上面的概率记作
αt(i),称为前向变量,根据上面的两个过程,递推关系也就十分明显了,
αt+1(j)=i=1∑Nαt(i)aijbj(ot+1),j=1,…,N
很显然,当HMM终止在
T时刻时,我们有
αT(i)=p(o1…oT,qT=si)
则,
p(O∣μ)=i=1∑NαT(i)
在完整给出算法以前,我们先将算法的变量用向量形式表达,定义向量
α(t)=(αt(1),…,αt(N))T
转移概率矩阵为
P=(aij)N×N,定义矩阵
C(ot)=diag(b1(ot),…,bN(ot))
下面给出矩阵形式的算法,
算法4.1(前向算法)
- 初始化
α(1)=C(o1)π;
-
Fort=2,…,Tdo
α(t)=C(ot)PTα(t−1);
- 算法终止
P(O∣μ)=i=1∑NαT(i);
上面的算法基于的是动态规划的思想,很容易计算,其算法的时间复杂度为
O(TN2),空间复杂度为
O(N2)。
4.3 后向算法
前向算法从序列的初始位置开始推理到终止状态
T,既然整个观测序列都已经知道了,那是不是可以由终止状态往前推呢?在动态规划的思想里,这两种递推结果一般是一致的。而得到后向算法的关键仍然在于定义一个存在合理递推关系的后向变量。
同样可以类似定义前向变量时的分析,假设我们从HMM序列的终止状态回溯至时刻
t(发射
ot以后,转移至
qt+1以前),我们已知的信息是观测序列
ot+1,…,oT,考虑从
t时刻往前回溯两步的过程,
-
qt以概率
bqt(ot)发射
ot;
-
qt−1转移至
qt;
也就是说,我们关心的是概率
p(ot…oT∣qt−1=si,μ)
记上面的概率为
βt−1(i),由上面的过程,可以很自然得到以下递推关系,
βt−1(i)=j=1∑Nβt(j)bj(ot)aij
则,
p(O∣μ)=j=1∑Nβ1(j)bj(o1)πj
同样,用矩阵表示,定义,
β(t)=(βt(1),…,βt(N))T
后向算法陈述如下
算法4.2(后向算法)
- 初始化
β(T)=(1,…,1)T;
-
Fort=T−1,…,1do
β(t−1)=PC(ot)β(t);
- 算法终止
p(O∣μ)=πTC(o1)β(1);
后向算法的时间复杂度和空间复杂度同样分别为
O(TN2)和
O(N2)。
4.4 联系
如果仔细观察,可以发现前向变量与后向变量之间事实上构成一种互补性,为了验证这一想法,我们做以下推导,
p(O∣μ)=i=1∑Np(o1…oT,qt=si∣μ)=i=1∑Np(ot+1…oT∣qt=si,o1…,ot,μ)p(o1,…,ot,qt∣μ)=i=1∑Np(ot…oT∣qt−1=si,μ)p(o1,…,ot,qt∣μ)=i=1∑Nαt(i)×βt(i)=α(t)Tβ(t)
通过引入动态规划算法,基本问题一可以被相对较高效地解决了。
五、维特比算法
接下来考虑基本问题二,即给定观测序列
o1…oT和模型
μ=(A,B,π),求最优的状态序列
q1q2…qT。首先我们需要定义“最优”,比较合理的想法是
Q∗=q1…qTargmaxp(Q∣O,μ)
直接对式
p(Q∣O,μ)进行时间上的递推是有些困难的,但是根据贝叶斯公式,有
p(Q∣O,μ)∝p(O,Q∣μ)
因而我们可以考虑对
p(O,Q∣μ)使用动态规划的算法,对问题进行分解,若对序列长度为
T时成立,则对子过程长度
t也应当满足
p(q1…qt,o1…ot∣μ)最大。基于上面的思路,定义如下的一个变量
δt(i)=q1…qt−1maxp(q1…qt−1,qt=si,o1…ot∣μ)
上述变量称为维特比(Viterbi)变量,类似前向算法和后向算法的分析,得到如下的递推关系(不加证明),
δt+1(j)=imax[δt(i)aij]bt+1(ot+1)
这里需要注意一个问题,由于存在转移概率,所以决定
t时刻的最优状态
qt的关键因素是
qt+1,因而,为了确定整个最优状态序列
q1…qT,需要从
qT开始回溯至
q1,因而算法需要记录下这些回溯路径,用符号
l(qt∗=si)表示,它定义为
qt∗=si时,
qt−1∗的状态。
这样,就能够对算法进行前向递推,再进行反向回溯,从而求解出最佳状态序列
q1…qT。把Viterbi算法总结如下,
算法5.1 Viterbi算法
- 初始化
δ1(i)=πi,i=1,…,N;
-
Fort=2,…,Tdo
Forj=1,…,Ndo
δt(j)=imax[δt−1(i)aij]bj(ot);recordl(qt∗=sj)=iargmax[[δt−1(i)aij];
- 递推终止
qT∗=jargmaxδT(j);
- 回溯
Fort=T−1,…,1do
qt∗=l(qt+1∗);
- 输出
q1∗…qT∗;
简单分析可知,Viterbi算法的时间复杂度和空间复杂度分别为
O(TN2)和
O(N2)。这样,基本问题二也被有效解决了。
六、Baum-Welch算法
6.1 先验or后验
上述两个算法都是基于的模型参数
μ=(A,B,π)已知,而事实上,实际问题中更多的情形是给定一些训练集,然后通过算法从中估计出合适的模型参数,也就是HMM的基本问题三,其数学形式为
μ∗=μargmaxp(μ∣Q,O)
这里的估计方式称为后验估计(posterior estimation),根据贝叶斯公式,有如下关系成立,
p(μ∣Q,O)∝p(Q,O∣μ)p(μ)
p(Q,O∣μ)称为似然函数,
p(μ)为参数的先验分布(prior distribution),这部分我们考虑极大似然估计的方法,即
μ∗=μargmaxlogp(Q,O∣μ)
6.2 监督学习
倘若足够幸运,我们搜集到了大量的序列对
(Q,O),那么,就可以构建监督学习的模型,也就是说,我们可以直接处理
对数似然函数
l(μ)=logp(Q,O∣μ)
带入参数,利用Lagrange乘子法求极值。通常情况下,序列
Q(例如词性序列)难以获得,或者获取成本极大,于是我们需要考虑无监督学习模型。
6.3 EM算法
我们依然引入监督学习中的似然函数,即
l(μ)=logp(Q,O∣μ)
这里的
Q为隐变量,无法获知,因而我们无法写出上式的具体形式。我们无法得知似然函数的具体取值(事实上它是
Q的随机变量函数),就很难通过解析或者迭代的方法一步一步地使它最大化,那么,退一步讲,如果
l(μ)的期望值一步一步最大化,是不是也就一定程度上意味着它本身也是在趋向最大呢?
上面的想法就是EM(Expectation Maximum)算法的基本思路,也就是对似然函数的期望求最大。在求期望的时候又会遇到一个新的问题,
Q的分布是依赖于参数
μ的,为了使得算法能够运行下去,我们将这里的参数设置为迭代前一时刻的参数(已知条件,通过初始化得来),而将当下的参数设为需要优化的变量,于是定义了
L函数(通常的文献里记作Q函数,这里为了区分序列Q),其表达式为
L(μ(k),μ(k−1))=Q∈ST∑p(Q,O∣μ(k−1))logp(Q,O∣μ(k))
上式完成了E阶段,即求期望阶段,接下来我们要做的就是找到上式的极大值,即M阶段。
μ(k)∗=μargmaxL(μ(k),μ(k−1))
接下来,推导上述表达式的具体形式,也就是Baum-Welch算法。
6.4 Baum-Welch算法
我们之前推导过
p(Q,O∣μ)的具体形式如下,
p(Q,O∣μ)=Q∈ST∑πq1bq1(o1)t=2∏Taqt−1qtbqt(ot)
由于
μ(k−1)已知,所以
p(Q,O∣μ(k−1))为一常数。于是,简要化简可得对数似然函数如下,
L(μ(k),μ(k−1))=Q∈ST∑p(Q,O∣μ(k−1))logπq1+Q∈ST∑t=1∑Tp(Q,O∣μ(k−1))logbqt(ot)+Q∈ST∑t=2∑Tp(Q,O∣μ(k−1))logaqt−1qt
接下来,首先考虑求参数
πi的极值点,
πi存在以下约束条件,
i=1∑Nπi=1
引入Lagrange乘子
λ,求导后可得,
∂πi∂L=p(q1=si,O∣μ(t−1))⋅πi1+λ=0i=1,2,…,N
解得,
πi(k)=∑i=1Np(q1=si,O∣μ(t−1))p(q1=si,O∣μ(t−1))
值得注意的是,在4.4部分,我们探讨了前向变量和后向变量的关系,于是,上式又可表示为
πi(k)=∑i=1Nπi(k−1)β1(i)πi(k−1)β1(i)
然后讨论参数
aij,它的约束条件为
j=1∑Naij=1i=1,…,N
同样,引入Lagrange乘子
λi,求导后可得
∂aij∂L=t=2∑Tp(qt−1=si,qt=sj,O∣μ(k−1))+λiaij=0
解得,
aij(k)=∑i=1N∑t=2Tp(qt−1=si,qt=sj,O∣μ(k−1))∑t=2Tp(qt−1=si,qt=sj,O∣μ(k−1))
利用前向变量和后向变量,我们可以推导以下关系,
p(qt−1=si,qt=sj∣O,μ)=p(O∣μ)p(qt−1=si,qt=sj,O∣μ)=p(ot+1…oT∣,qt=sj,o1…ot,qt−1=si)×p(ot,qt=sj∣qt−1=si)×p(o1…ot−1,qt−1=si)/p(O∣μ)=∑i=1N∑j=1Nαt−1(i)aijbt(ot)βt(j)αt−1(i)aijbt(ot)βt(j)
为了书写简便,我们引入两个新的符号,记
ξt(i,j)=p(qt=si,qt+1=sj∣O,μ)
γt(i)=p(qt=si∣O,μ)=j=1∑Nξt(i,j)
因此,
aij的新的参数更新式如下,
aij(k)=∑t=1T−1γt(j)∑t=1T−1ξt(i,j)
最后,我们考虑参数
bj(n),它表示概率
p(ot=vn∣qt=sj)。它满足的约束条件为
n=1∑Mbj(n)=1,j=1,…,N
利用Lagrange乘子法,得到如下方程,
∂bj(n)∂L=t=1∑Tp(qt=sj,O∣μ(k−1))1(ot=vn)+λjbj(n)=0
解得,
bj(n)=∑t=1Tp(qt=sj,O∣μ(k−1))∑t=1Tp(qt=sj,O∣μ(k−1))1(ot=vn)
带入符号,得
bj(n)=∑t=1Tγt(j)∑t=1Tγt(j)1(ot=vn)
至此,所有的参数更新式都已推导出来,下面将算法总结如下,
算法6.1
- 初始化
μ0=(A0,B0,π0);
-
Fork=1,…,maxiteratedo
Forwardalgorithm→α(k);Backwardalgorithm→β(k);ξt(i,j)←∑i=1N∑j=1Nαt−1(i)aijbt(ot)βt(j)αt−1(i)aijbt(ot)βt(j);γt(i)←j=1∑Nξt(i,j);aij(k+1)=∑t=1T−1γt(j)∑t=1T−1ξt(i,j);bj(k+1)(n)=∑t=1Tγt(j)∑t=1Tγt(j)1(ot=vn);πi(k+1)=γ1(i);
- 算法终止,输出
μ=(A,B,π);
七、一些可用的工具
hmmlearn和seqlearn是比较好用的用于序列学习的工具包,提供了python接口,链接如下
https://pypi.org/project/hmmlearn/
https://pypi.org/project/seqlearn/
[参考资料]
[1]: 宗成庆.《统计自然语言处理》