隐马尔可夫模型的原理与实现.docx
- 文档编号:15180299
- 上传时间:2023-07-02
- 格式:DOCX
- 页数:25
- 大小:273.65KB
隐马尔可夫模型的原理与实现.docx
《隐马尔可夫模型的原理与实现.docx》由会员分享,可在线阅读,更多相关《隐马尔可夫模型的原理与实现.docx(25页珍藏版)》请在冰点文库上搜索。
隐马尔可夫模型的原理与实现
・253・
隐马尔可夫模型的原理与实现
刘河生,高小榕,杨福生
(清华大学电机工程与应用电子技术系,北京 100084)
摘要:
隐马尔可夫模型正在被愈来愈多地引入到生物医学信号的处理中。
本文旨在简述它的基本原理和实现中的问题,并且用简洁的列表形式总结它的算法步骤。
关键词:
隐马尔可夫模型;信号处理;实现算法
中图分类号:
R311;R318 文献标识码:
A 文章编号:
100121110(2002)0620253207
TheoryofhiddenMarkovmodelinganditsimplementation
LIUHe2sheng,GAOXiao2rong,YANGFu2sheng
(DepartmentofElectricalEngineeringandAppliedElectronics,TsinghuaUniversity,Beijing100084,China)Abstract:
HiddenMarkovModelisnowbeingappliedincreasinglyinbiedicalp.Thepapermakesashortreviewonitstheoryandtheproblemsencounteredinitsimon.clearlysummarizeits.algorithmicprocedures
Keywords:
hiddenMarkovmodel;signalp;ipt1 概况[1,2]
应用,它正被愈来愈多地引入到生物医学信号的处理中。
本文旨在简述它的基本原理和实现中的问题,并且用简洁的列表形式总结它的算法步骤。
隐马尔可夫模型是马尔可夫模型的进一步发展。
马尔可夫模型是马尔可夫过程的模型化,可以用图1(a)的框图形象表示。
它把一个总随机过程看成一系列状态的不断转移。
时刻t的状态用qt表示,它可以是N种状态集合S=[s1,s2,…,sN]中的任意一个。
马尔可夫模型的特性主要用“转移概率”来表示。
后一状态出现的概率决定于其前出现过的状态次序。
即:
状态qt出现的概率为Pr[qtqt-1,qt-2,…,
q1]。
如果此概率只决定于前一个状态,即Pr[qt
[3~7]
觉”、“快速眼动”、“睡眠一期”……,它们便是“状
态”;而可以观测到的则是在这些状态下的各种生理参数表现,例如在脑电图上的表现。
这些表现便构成观察。
t时刻的观察记作OtO当观察是离散型时,OtO是总的观察集合V=[u1,u2,…,uM]中的一种,如图1(b)所示。
注意M未必等于N
。
它是研究中引用得qt-1],则称为一阶马尔可夫过程。
最多的形式,即:
Pr[qtqt-1,qt-2,…,q1]=Pr[qt
qt-1]。
隐马尔可夫模型(HiddenMarkovModel,HMM)则认为模型的状态是不可观测的(这便是“隐”得名的由来)。
能观测到的只是它表现出的一些
图1 (a)马尔可夫过程 (b)隐马尔可夫过程
隐马尔可夫过程的特性可用下述参数集合来表征:
(1)转移概率aij=Pr[sjsi]即:
由状态i转移到
观测量(observations)。
例如:
睡眠的状态可分为“醒
收稿日期:
2002202221
状态j的概率(对一阶马尔可夫过程)。
由于共有N种可能的状态,因此aij共有N×N个可能的取值。
把它们用矩阵表示成
・254・
N
A=[aij] 且
∑a
j=1
ij
=1
(2)观察概率bj(k)=Pr[uk(sj],即:
在状态sj下
产生观察uk的概率。
如果共有M种可能的观察,则bj(k)组成M×N矩阵B。
N
的状态序列Q=[q1,q2,…,qT]却并不唯一,因为qt的每一个取值都能以一定概率产生给定的ot[例如,设ot=u5,则产生此u5的概率分别是b1(u5),b2(u5),…,bN(u5)]。
也就是说任何一组状态序列q1q2…qT都能依下述联合概率产生观察序列O=[o1,o2,…,
oT]:
q1bq1(o1)aq1q2b(o2)aq2q3bq3(o3)…aqt-1qTbqT(oT)Π
(1)
式中:
Πq1是初始时处在q1状态的概率;bq1(o1)是在
B=[bj(k)] 且
∑b
k=1
(k)
=1
(3)初始状态概率:
指第一个状态q1究竟取S=[s1,s2,…,sN]中哪一个的概率。
它组成1×N矢量
Π:
i=Pr[q1=si]Π
q1状态下产生观察o1的概率;aq1、aq2是由状态q1转
而Π=[Π1,Π2,…,ΠN]
以下讨论中把上述参数合起来用Κ表示:
Κ=[A,B,Π]。
它便是表征HMM的参数集合。
采用HMM进行研究工作时常遇到三类问题:
(1)评价问题:
给定模型参数Κ=[A,B,Π]及观察序列O=[o1,o2,…,oT]。
求此模型产生此观察序列的概率Pr[OΚ]。
即:
设有L,12、ΚL,且皆已知。
=o1,o2,…,oT]给予这组模型,看哪一个Pr[OΚi]最大就认为该观察属于Κ。
这也就是选择与观察最匹配的模型。
i一类
(2)解码问题:
给定模型Κ及观察序列O;问此观察序列是模型Κ中取怎样的状态次序[q1→q2→…→qT]得到的。
解决此问题的关键是采用什么作为取得结论的判据。
通常是取产生此观察序列概率最大的一组状态序列Q=[q1,q2,…,qT]作为判决。
(3)辩识(或称训练)问题:
给定HMM的结构(指状态数N,观察类数M),由给定的一组供训练
移到状态q2的概率。
由于每一组状态序列都能产生一个这样的概率,而每一个状态qi都有N种不同的可能取值。
因此类似的路径共有NT条[见图2。
图中由q1到qT的]。
把这些路径
Pr[OΚ]。
即:
图2 状态转换展开图
sN
sN
sN
q1
q1
Pr[OΚ]=
q1=s1
∑ ∑…∑Πb
q2
qT=s1
(o1)aq1q2b(o2)
…aqT-1qTbqT(oT)
简记作 =
ΞΞΞ
)Κ∑Pr(Q,O
Q
(2)
上述概念并不复杂,只是计算量过大:
(1)式包含(2T-1)次乘法,因此
(2)式包含(2T-1)NT≈2TNT
用的观察组O1,O2,O3…[3],估计该模型的最优参数。
δ=[Aδ,Bδ,Πδ]Κ
第三类问题是更基本的,因为前两类问题都建
次乘法。
为了实现Pr[OΚ]的计算,需要寻找更高效的算法。
考虑到上述计算过程中实际作了许多重复计算,因此改进算法的基本思路是:
步步为营的递推算法。
其具体策略又可分成三类:
(1)前向递推法:
其基本思路可借助于图3来说明。
令at(i)代表给定模型参数Κ下过程到t时刻为止产生部分给定观察序列o1~ot,且在t时刻状态为
si的部分概率,记作:
立在Κ已知的基础上。
但由于第三类问题算法比较复杂,因此以下讨论时把它放在前两类之后。
2 三类问题的基本解法
2.1 评价问题:
给定HMM参数Κ=[A,B,Π],同时还给定一组观察O=[o1,o2,…,oT][33],求此模型产生此观察序列的概率Pr[OΚ]。
需要着重指出的是:
尽管O已给定,但产生此O
此处O是大写字母,指状态序列
[o1,o2,…,oT];ΞΞ所谓给定“意即各ot的取值已给定;例如:
O=[o1,o2,…,oT]”o1=u2,o2=u8,…。
Ξ
・255・
at(i)=Pr[o1,o2,…,ot;qt=siΚ]
求at+1(j)=Pr[o1,o2,…,ot,ot+1;qt+1=sjΚ]与
at(j)的关系。
由图3可见:
N
at+1(j)=[
∑Α(i)Α]b(o
i=1
t
ij
j
t+1
)(3
)
图5 t时刻和t+1时刻状态si与sj的关系
N
Βt(i)=
图3 由Αt(i)推导Αt+1(j)
∑a
j=1
ij
bj(ot+1)Βt+1(i)
(4)
~N,t=T-1,T-2,…,1i=1
2所示。
初始化令ΒT(i)
=~NΒT(i)=1,i=12.递推:
N
由此可得前向递推的步骤如表1所列。
此法所
需乘法次数约为N2T,使
(2)式的计算量大为下降。
表1 前向递推法
1.初始化:
Αt(i)=Πibi(o1)2.递推:
N
t=1t=T时的赋值
Αt+1(j)=[∑Αt(i)j
(1)
i=1
~T-1,j=Nt=1
3.终止:
止于t=T
N
(3):
Βt(i)=∑Αijbj(ot+1)Βt+1(j)
j=1
~N;t=T-1,T-2,…,2,1i=1
3.终止:
止于t=1
N
即(4)式:
由右
向左逐级递推
)=Pr(OΚ
∑Α(j)
j=1
T
)=Pr(OΚ
∑Β(i)
i=1
1
(2)后向递推法:
此时改为由过程的后部分(例
如t+1,t+2,…,T
)向前递推。
其基本关系如图4所
(3)前后向递推法:
也就是把前向递推与后向递
示。
推结合起来,o~t用前向递推,T~t用后向递推(t
值可在1~T间任意选定)。
此时有:
N
)=Pr(OΚ2.2 解码问题
∑Α(i)Β(i)
i=1
t
t
(5)
此时的任务是对给定的模型参数Κ=[A,B,Π]和给定的观察序列O=[o1,o2,…,oT],求此观察序列最可能由怎样的状态序列Q=[q1,q2,…,qT]产生。
通常认为最可能产生它的状态序列就是按
(1)式
图4 由Βt+1(j)推导Βt(i)
计算得概率值最大的一个序列。
原理上说,解决此任务的概念也很简单:
按图2和
(1)式把每条可能路径的概率都求出来,然后取所得概率最大的一条路径便是所求。
这种“穷举法”的问题仍然是计算量过大。
解决之道仍然是步步为营的递推算法。
显然,定义:
令Βt(i)为在给定模型参数Κ,且qt=si条件下
产生观察o~toT的部分概率。
即:
Βt(i)=Pr[ot,ot+1,…,oTqt=si,Κ]求Βt+1(j)=
Pr[ot+1,ot+2,…,oTqt+1=sj,Κ]与Βt(i)间的关系。
由图五可见:
・256・
=Χi(t)=)Pr(OΚ
t
t
N
(6)
Q=[q1,q2,…,qT]。
3333
∑Α(i)Β(i)
i=1
~Ni=1
然后取Χi(t)最大的一个便是所求答案。
不过实际更
常采用的是逐步搜索前进的Viterbi算法。
Viterbi算法的核心是下述递推关系:
如果已经
为了方便地获得此最优状态序列,需要在递推的每一步把对应于每个∆t(j)(j=1~N)的“获胜状态”记录下来,以便到最后再通过回推得出各q3t。
因
3此,按表3,最后除得到qT外还会得到一个二维矩阵Wt(j)(t=1~T,j=1~N)。
根据W就可以逐步3反推,由q3t+1找到qt,从而把最优状态序列确定下来。
2.3 参数辩识问题
此时是用一组观察序列为训练集来优化HMM的参数Κ[A,B,Π],使Pr[OΚ]达到最大。
因此要采
确定了对给定部分观察序列o1,o2,…,ot时的最优状
态序列q1,q2,…,qt,如何找到当部分观察序列增1时[也就是o1,o2,…,ot,ot+1]的最优状态序列q1,q2,…,qt,qt+1。
定义辅助变量如下:
ax
∆t(i)=mq1,…,qt-1 Pr[q1,q2,…,qt=sio1,o2,…,ot,Κ]它的含义是:
对给定的模型参数Κ和部分观察序列~ot,满足qt处在si状态下使概率最大的状态序列o1
称此概率为∆t(i),现在以∆t(i)为出发q1,q2,…,qt-1。
点来求
∆t+i(j)=
max
q1,…,qt-1,qt
Pr[q1,q2,…,qt-1,qt-1,qt,
用一种优化调节步骤。
优化算法不止一种,下面介绍的是Baum2后向算法),它是一Welch算法(也称前、
种基于期望调节[expectationmodification]概念的算法。
Εt(i,j)Κ和观察序列O下,tsi,而t+1时刻处于sj状态的概率:
O]t(i,j)=Pr[qt=si,qt+1=sjΕΚ
qt+1=sjo1,…,ot,ot+1,Κ]
不难理解
max
∆t+i(j)={i=1~N[t(iijjt+1(7)
qts1,…,sN中取∆t(i)aij
最大的一个节点;sj产生ot+1的概率,它与qt=si的选择无关。
但值得强调的是:
随着∆t(i)中qt=si的不同选择,相应的q1,…,qt-1最优路径也将随之而变。
据此便得出∆t(i)中递推过程如表3所示。
表3 Viterbi解码过程
1.初始化:
∆1(i)=Πibi(o1)
~Ni=1
W1(i)=0 此时尚未开始记录2.递推:
ax
∆t(j)=[mi=1~N∆t-1(i)Αij]bj(ot)
~T,j=1~Nt=2
arg
t(j)=i=1~N[∆t-1(i)aij] t=2~T,j=1~N
[此步是对每个t的每个j,记录该∆t(j)的得胜节点,也就是前一步的qt-1来源,以便最后回溯]。
3.终止:
ax
最优总概率Pr3=mi=1~N[∆T(i)]
max
最终选定qT3=argi=1~N[∆r(i)]4.最优状态回溯:
3
由q3i+1反查W记录,得知qt其表示式可以借助图5来考虑。
图上突出地画出t时刻的状态si和t+1时刻的状态sj两个节点。
对于t~t+1时段以外的部分,左边(即:
时刻
不难看出:
t(i,j)=Pr[qt=s1,qt+1=sjΕO,Κ]
=
Pr[OΚ]
由图6可得:
Εt(i,j)=
N
N
t
ij
j
t+1
(8)
∑∑Α(i)Αb(o
i=1j=1N
)Βt+1(j)
上节提出的Χt(i)与此处的Εt(i,j)存在下述简单关
系:
Χt(i)=
∑Ε(i,j)
j=1
t
(9)
更进一步,如果把Χt(i)对所有t(t=1~T-1)求和,它代表着各种可能取的状态序列中si被访问的次数,它也就是由si转移出去的次数[时刻T不被考虑,因为它不再转移出去],类似地,把Εt(i,j)对
~T-1求和代表si→sj这种转移出现的次数,t=1即:
T-1
这时的相应路径才是最后判断得的最优路径
∑Χ(i)=
t=1
t
从si转移出去的期望次数
T-1
・257・
∑Ε(i,j)=
t=1
t
T-1
从si转移到sj的期望次数
因此可得aij和bj(k)的新估计值:
δij=Α
∑Ε(i,j)
t
T-1
(10)
∑Χ(i)
t
t=1
∑
T
bj(k)=
δ
s,t.∶O=uT
t=1
Χt(j)
t
Ξ
(11)(12)
图6 HMM的类型
∑Χ(j)
δ=这样我们就把原来的Κ=[A,B,Π]更新成Κδ]。
Baum证明了:
一般情况下Pr[Oδ]>Pr[Aδ,Bδ,ΠΚ
δ比Κ更接近于优化目标,因此经反复迭[OΚ],即Κ
δi=在t=1时状态si的出现概率=ΧΠ1(i)
t=1
3.2 关于HMM的状态数N和观察数M
从物理概念上看,状态数N应和待研究的任务
有关。
例如进行睡眠分析时可将状态取为目前较公认的醒觉、快速眼动及反映睡眠深度的各分期。
进行,并且将状态取为
δ代就可以得到Κ的极大似然解,而且满足,∑Πi=1,P2Q、QS、P。
不过从数学上看i
NM
它
δδ~N),∑bj(k)=1 (=1~,容许一定任意性。
Αij=1 (i=1∑我们在脑电模j=1k=1
但是,,N=3、4、5进行分析,取得大。
的最终分析效果却相差不多。
。
至于观察符号数M则决定于研究时对指定状
态提取哪些特征。
例如,在语音、脑电分析中常通过
3 实现HMM时的一些问题
建立AR模型来提取特征,此时观察便是模型参数
3.1 关于模型的类型所构成的矢量。
图6所示是几种常见类型,图中每一节点代表以上讨论中均假设观察为离散型,但实际工作
中还常遇到观察是连续信号(矢量)的情况(例如:
直一个状态。
图6(a)是转移矩阵A中各元素aij都不等
接用电生理数据作观察值)。
虽然可以通过矢量量化于0的所谓“各态遍历模型”,每一个状态都有转移
(c)则以编制码书的方式将连续变量转换成离散的码字,到任意其他状态(包括本身)的可能。
图6(b)、
但这样作却可能劣化分析结果。
因此有些情况下有是限定状态只能由低序号向高序号转移的“左—右
必要采用连续的概率密度函数(pdf)来描述连续观模型”。
其中图6(b)的约束条件是[aij≠0,当且仅当
察矢量。
最常用的办法是把观察矢量的pdf拟合成j=i或i+1];而图6(c)则容许状态有间隔一个序号
若干高斯函数的线性组合。
设O为待拟合的观察矢的跳转,其约束条件是[aij≠0,当且仅当j=i或i+
量,令1或i+2]。
对于左—右模型必有Π1为1,其余各Πi=0。
而且由于最终状态不再转移出去,因此aNN=
1。
图6(d)则是由若干并行支路合成的模型。
有时各
M
bj(O)=
并行支路之间还容许有交叉转移,如图中所示。
值得指出的是,当任何一个aij被初始约束为0后,在其后辩识过程中作参数更新估计时该aij仍必为0,因此模型结构不会被破坏。
还应该指出,模型类型的可能变种还有很多,研究人员完全可以根据自己任务的特点构作适当的模型结构。
其难点主要是参数更新公式的推导
。
(8)∑cN(O,u,∑) j=1~N
式中Λ、∑、c分别是状态j的对应观察中第m
m=1
jm
jm
jm
jm
jm
jm
个高斯分量的均值矢量、协方差阵和混合权重。
cjm满足约束条件,
M
[cjmΕ0,
∑c
m=1
jm
=1],因此bj(O)dO=1。
∫
采用连续pdf后引起的新问题是:
进行HMM
Ξ
(11)式分子的求和中有约束条件s.t.∶uk
。
・258・
参数递推更新时的算法。
由(8)式可见bj(o)决定于
cjm、Λjm、
∑
jm
三项因素,因此更新bj(o)也就是更新
这三项因素。
文[8、9]证明了它们的更新公式如表4所列。
这些公式从概念上也不难理解。
Χt(j,k)表示在
t时刻用第k个高斯分量表征ot时系统处在状态j
个概率的乘积。
由于概率值必定小于1,所以当状态
数目较多时此乘积将愈变愈小,甚至可能低于计算机数据容许动态范围的下限,使计算无从进行。
解决此困难的途径有二:
一是在计算过程的适当时刻t引入一些与状态i无关的归一化运算。
待到计算终了再行复原。
但更方便的方法是改用对数将乘法运算改成加法运算。
由于对数函数是单调函数,因此这样作不会影响最后结果。
例如,可以将表3的Viterbi算法改成对数形式如表5。
表5 对数形式的Viterbi算法
1.初始化:
下的概率,因此表中(ii)式反映此概率在全时间过程中所占百分数。
(iii)式是对(ii)式分子乘以观察值ot作为权重,因此它反映系统处在j状态下用第k个高斯分量表征O时O的均值。
表4 连续变量时的参数更新
第一步:
计算辅助变量
Χt(j,k)=[
N
][
cjkN(Ot,Λk),
M
∑Αt(i)Βt(i)
j=1
∑
m=1
cjmN(Ot,Λjm
∑
Υ1(i)=logΠi+log[bi(o1)]
])jm
(i)
2.递推:
Ωt(i)=
3.:
max
i=1~N[Ωt-1(i)+logΑij]+log[bj(ot)]
j=1~N,t=1T
3
ax~T t=1
第二步:
参数更新
δCjk=
T
TM
∑Χ(j,k)
t
ii∑∑
t=1k=1T
Χt(j,k)
t
以上对HMM作了初步但基础性的介绍。
目前,HMM不但在许多科技领域得到了广泛应用,而且
δjk=Λ
∑Χ(j,k)t
TMt=1k=1T
(iii)
∑∑Χt(j,k)∑Χ(j,k)・(o-t
t
T
t
t=1
δjk=2
Λjk)(ot-Λjk)’
(iv)
在其基本框架上也有许多发展。
其中与生物医学信号处理关系较密切的有以下方面。
(1)灵活的模型结构变种:
根据具体任务,从图6基础上发展出多种新的模型例如:
是所谓“嵌入式(EmbeddedHMM)的结构图[10,11]。
其特点是HMM”
∑Χ(j,k)
一般言之,连续观察HMM比离散观察所需训练样本数也更多。
3.3 递推辨识模型参数Κ=[A,B,Π]时初始值的
由两级状态组成,即:
超级状态及内置状态。
超级状态共有三级,而每一超级状态中又包含一组由内置状态组成的HMM。
系统的状态可在各超级状态间转移,且在每一超级状态下实际上是在该状态包含的内含状态间转移。
各超级状态的内含状态却不能直接转移。
这种结构特别适用于具有二维结构的信息分析(如图像分析等)。
又如是文[12]提出的计算生物学中用于蛋白质建模和符号序列分析的
每一位置包括三种状态,矩形块为匹配状态,HMM。
菱形块为插入状态,园形块为删除状态,各联线之间
各有一个转移概率互相联系。
其具体解释可参看文[13]。
(2)与其他方法的结合:
生物医学信号处理的发展趋势是,不能企图只靠单一方法完全解决问题,需要把多种方法分层次地结合起来,让不同方法在不同层次上发挥其优势。
文[13]介绍一种将人工神经网络与HMM相结合的处理方法,文[14]则介绍一
设定
δ只是局前已指出,前、后向算法收敛后得到的Χ
部的极大似然解。
因此如何选择参数的初始值以便使此局部极大就是全局极大是一个值得研究的问题。
应该说,目前还没有解决此问题的有理论根据的
方法。
经验证明,在给定的已知约束条件下(例如某些aij=0)随机或均匀地选择初始Π和A值是个可行办法。
但B的初始值选择则影响较大,特别是在观察是连续变量的情况下。
解决此问题只能根据具体情况来决定。
本文从略。
3.4 尺度问题(scaling)
由
(1)~(4)式都可看出,Α、Β的计算都涉及多
・259・
[7]CoastDA,StemRM,CanoGG,etal.Anapproachto
cardiacarrhythmiaanalysisusinghiddenMarkovmod2els[J].IEEETransBME,1990,BME217(9):
8262835.[8]liporaceLA.Maximumlikelihoodestimationformulti2
rateobservationofMarkovsources[J].1986,IT232
(2):
3072309.
[9]JuangBH,LevinsonSV,SondhiMM.Maximumlikeli2
hoodestimationformultivartiatemixtureobservationsofMarkovchains[J].IEEETransonInformationThe2ory,1986,IT232
(2):
3072309.
[10]KuoS,AgazziE.Keywordspottinginpoorlyprinted
documentsusingpseudo22DMarkovmodel[J].IEEETransonPatternAnalysisandMachineIntelligence,1994,PAMI216(8):
8422848.
[11]NefianAV,HayesMH.AnembeddedHMM2basedap2
proachforfac
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 隐马尔可夫 模型 原理 实现