265
龚光鲁,钱敏平著 应用随机过程教程 – 与在算法和智能计算中的应用
清华大学出版社,2003
第 10 章 隐马氏模型 （ Hidden Markov Model,HMM ） 及其应用
１,熵与相对熵
１,1 离散分布的熵与相对熵
熵的概念出自 C.Shannon,引进这个指标的 目的在于刻 画 一个离散分布 ( 一个离散随机变量 ) 或一个分布密度 ( 一个连续型随机变量 ) 的不确定性的大小,也就是说知道了此随机变量的取值所获得的信息的大小,
定义 １０,１ 对于离散分布 p ),,,( 1 LL npp=,我们定义它的 熵 为
H( p )= ∑?
i
ii pp ln,
又定义分布 p 关于分布 q ),,,( 1 LL nqq= 的 Kullback - Leibler 相对熵 为
h （ p,q ） = ∑
i i
i
i q
pp ln,
命题 １０,２ 相对熵 h （ p,q ） ≥ 0,而且 h （ p,q ） =0,当且仅当 p = q 时成立等号,
证明 ),0[ ∞ 上 函数 tttg ln1)(=
D
在 1≠t 时恒正 ( 这 一结论可 由 g 的导数直接可以看出 ),且 0)1( =g,取
i
i
p
qt =,于是
i
i
p
qln ≤
i
i
p
q - 1,即
i
i
q
pln
i
i
p
q?≥ 1,而且等号仅当
1=
i
i
q
p 时成立,从而有
∑
i i
i
i q
pp ln ∑ =?≥
i i
i
i p
qp 0)1(,
所以结论成立,
这个命题表明,相对熵在相当程度上表达了 p 与 q 的差别,当 p = q 时,h （ p,q ） =0,而当所有的 ip 都与 iq 接近时,h （ p,q ） 就很小,从而 h （ p,q ） 可以看成 p 与 q 之间的一种
,准距离,,这里我们之所以称它为准距离,是因为它既不对称 ( 即 h （ p,q ） ≠ h （ q,p ） ),
也不满足三角形不等式 。 所以不满足 第 9 章 中的距离公理,
例 １０,３ ( 有限个值的分布 的熵 )
分布 p ),,( 1 Npp L= 的 熵 满足
H( p )= Npp
i
ii lnln ≤? ∑,
且等号当且仅当 分布 在 N 个 值均匀时成立,即以相同概率取 N 个值的分布 ( 称为离散均匀分布 ) 的熵最大,
证明 记分布 n = )1,,1( NN L,于是 本结论是 相对熵 不等式 h （ p,n ） ≥ 0 的变形,
例 １０,４ ( 数学期望固定条件下的离散 的 最大熵分布 )
假定 存在实数 a,使 )1(,≥≥ ixi a,对于固定的 ),,,( 1 LL nxx 与 m,在满足
mxx === EpxP ii,)(
266
的离散分布 p ),,,( 1 LL npp= 中,具有最大熵的分布为 p *,且
∑ ==
i
xx
i
ii eCCep ))((,1* 00 ll,
其中 0l 满足 条件
∑ ∞<?
i
xie 0l,∑ ∞<?
i
x
i
iex 0l,ii
x
i i
i
x exe 00 llm∑ ∑=,
这个最大熵 为 Cln?０ml,
证明 对于任意 0≠l,只要 ∑ ∞<?
i
xie l,相对熵不等式
≤0 h （ p,p * ） ∑?=
i
x
i
i iCe
pp
lln
∑∑?+=
i
ii
i
ii Cxppp lnln l
就是
H( p )= Cpp
i
ii lnln?≤? ∑ lm,
所以
H( p * ) ≥?=?= ∑ CCeCe
i
xi ln]ln[
0
00 mlll H( p ),
例 １０,５ （ 数学期望和方差都固定 时 的最大熵 ）
在 均值为 μ,方差为 σ 2,且 取值 为 kx ),2,1( L=k 的 离散分布中,分布 p *,
1
)()(
* )(,2
2
2
2
∑==
k
xx
k
kk
eCCep b
a
b
a
的熵最大,此 最大熵 为 Cln)( 2
22
+ b ams,其中 ba,满足
)2,1,0(,)( 2
2)(
=∞<∑
lex
k
x
l
k
k
b
a
,∑
k
xk
e 2
2)(
b
a
m,2
2)(
∑
=
k
x
k
k
ex b
a
=∑
2
)(
2 )( 2
2
k
xk
e b
a
s ][ 2
2)(
∑
k
xk
e b
a
][ 2
2)(
2∑
k
x
k
k
ex b
a
∑
k
x
k
k
ex 2
)(
)( 2
2
b
a
,
( 证明类似 ),
１,２ 分布密度 的熵与相对熵
定义 １０,６ 对于概率分布密度,我们可以仿效前面的思想,把 分布密度 )(xp 的熵 ( 其实是微分熵 ) 定义为
H ∫?= dxxpxpp )(ln)()(,
而把 分布密度 p(x) 对 分布密度 q(x) 的相对熵 定义为
∫= dxxq xpxpqph )( )(ln)(),(,
这个定义可以推广到多维密度 ),,( 1 dxxp L 与 ),,( 1 dxxq L,
267
H ∫= d
d
d
d dxdxxxq
xxpxxpp L
L
LL
1
1
1
1 ),,(
),,(ln),,()( ),
∫= d
d
d
d dxdxxxq
xxpxxpqph L
L
LL
1
1
1
1 ),,(
),,(ln),,(),(,
同样 也有 相对熵是非负的,从而由熵的定义,再利用相对熵非负性质 可 推得下述结论,
例 １０,７ ( 半直线上均值已知时的最大熵分布密度 )
在 ），［ ∞0 上具有同均值 m 的密度中,指数分布
m
1Exp 的熵最大,这个 最大 熵为 mlog1+,
例 １０,８ ( 均值与方差已知的 分布密度的 最大熵 )
在均值 为 m,方差为 2s 的 分布 密度中,正态分布的熵最大,
例 １０,９ ( 均值向量与 ( 协 ) 方差矩阵已知的有 联合密度的 最大熵 )
在均值向量为 m,方差矩阵为 Σ 的多维密度中,Gauss 分布 ),( ΣmN 的熵最大,请读者自己写出这个最大熵的表达式,
例 １０,１０ ( 有限区间上的最大上密度 )
在有限区间 ],[ ba 上 取值 的 分布 密度中,均匀分布的熵最大,此最大 熵为 )ln( ab?,
例 １０,１０ 也可以推广到多维情形,
例 １０,１１ ( 在 取值于 混合概率向量 的分布密度 中,关于某个 概率向量 的平均相对熵给定条件下,求 具最大相对 熵 的密度 )
令样本空间为
{=? x = }1),(,10:),,( 11 =++≤<< did xxdixxx LL,
中的 点可以解释为概率向量,而取值于? 的连续型随机变量的分布密度可以解释为,广义的混合分布,,给定 W∈= ),,( 1 dzzz L,把? 中的点看成概率向量,就有相对熵 ),( xzh,
它 是 取值于? 的 x函数,对于一个取值? 的分布密度 xp( ),,() dx xxp L
= 定义对于这个分布密度的平均相对熵,
(m z xdxpxzh )(),() ∫=
D
,
我们要在平均相对熵 )( zm 相等的 分布密度 )(xp 中,选取 一个 )(* xp,使 其 熵
∫? xdxpxp )(ln)( ** 达到 最大,
为此我们 断言,如下的 Gibbs 型分布密度
)()(* zCxp = ),()( xzhze l? )(( zC= ))()(1
ln)(
1 di
ii zz
d
zz
zzz
xxe ll
l
L∑
,
其中
1),()( )()(∫= xdezC xzhzl,
)(zl 满足
xdexzhzCz xzhz ),()(),()()( lm?∫=,
就是我们要的具有最大熵的密度,证明如下,对于任意一个满足 (m z xdxpxzh )(),() ∫=D
的分布密度 )(xp,由 相对熵不等式
268
),(0 *pph≤ xd
ezC
xpxp
xzhz∫?
=
),()()(
)(ln)(
l
= H( )(ln)()() zCzzp?+ ml
得到
H( )p )(ln)()( zCzz?≤ ml,
于是 我们有
H( *p ) = xdezCezC
xzhzxzhz ])(ln[) ),()(),()( ll∫? （
)(ln)()( zCzz?ml＝ ≥ H( p ),
[ 注 1] 上面我们用相对熵 ),( xzh 作为准距离,而不是用更为合理的 [21 ),( xzh )],( zxh+,
也可以用 与它相应的 ),( zxh,只 要 易于计算,不管用 那一个 都行,
[ 注 ２ ] Diriclet 分布的密度函数为
)1{
11
11 1
1),,(
=++
=
d
d
xxdd IxCxxxf LLL
aa,
其中
∫
∞
=Γ
++Γ
ΓΓ=
0
1
1
1 )(,
)(
)()( dxexC x
d
d aa
aa
aa
L
L,
例 １０,１１ 中的最大熵分布恰是 Dirichlet 分布,在生物信息论中,常常需要估计概率向量组成的空间上的分布密度,例 １０,１１ 说明,在平均相对熵给定的条件下,Dirichlet 分布是 ＂ 最吃亏的 ＂ 分布,用它作为先验分布 ＂ 看起来更为保险 ＂,这就解释了在生物信息论中,人们常常喜欢 用 Dirichlet 分布作为先验分布的原因,
相对熵应用的一个实例 － 在各个测试特征 （ 统计量 ） 中选择几个最为有效的相对熵方法
（ 特征量选取 的相对熵方法 ）
假定我们有两组性质完全不同的群体,例如一组是健康人,另一组是 SARS（ 非典型性肺炎 ） 病毒持带人,又假定人们已经提出了 N 种区别健康人与 SARS 病毒持带人的不同特征,我们要在其中选取区别效果最好的 M 个特征,相对熵方法是较为有效的一种方法,其实际操作为,
用一个给定的区 分特征,对上面的一组健康人测定了一组数据,简称为甲数据组 ； 又对上面的
SARS (Severe Acute Respiratory Syndrome) 病毒持带人测定了一组数据,简称为乙数据组,再进行如下的步骤,
（ １ ） 分别找出此两个数据组的近似分布,从数据组出发,应用统计中的核估计的思想,分别得到其近似分布密度曲线 ；
（ ２ ） 计算此两个分布密度的相对熵 0h ；
（ 实践证明,用数值计算求此两个分布密度的相对熵,对于计算格点大小的划分并不太敏感,相反地,如果直接将两 个直方图作为离散分布求相对熵,则对于直方图的计算格点大小的划分十分敏感,
以致得到的计算结果很不稳定 ）,
（ ３ ） 将甲乙两组数据合并,再随机地重组为和以前个数相同的两组 (随机地重排 （ Random
Sorting） 后,再按原来的各组的个数顺序分成两组 ),用 （ １ ），（ ２ ） 步骤计算其相对熵 ；
（ ４ ） 重复地作 （ ３ ） 多次 （ 例如 １ 万次 ）,计算其中相对熵大于 0h 的次数所占的频率,记为 p,
269
于是对应于一个区分特征就对应地有一个 p 值,它是一个客观的标准,和各种测试方法的结果
（ 统计量 ） 的绝对大小,量纲都无关,如果要在 N 个不同的用以区分是否感染了 SARS 的特征中选取 M 个相对地最为有效的特征,我们只需选取对应的 p 值最小的前 M 个特征即可,
２,隐 Markov 模型 (HMM)
２,１ 一个实例
先举一个例子以 直观地 理解这个 模型 的实质,
例 １０,１２ 设 某人 在 三个 装有红白两种颜色的球的 盒子 中,任取一个盒子,然后在此盒子中每次抽取一个球,连续地抽取 m 次,假定各个盒子的内容与抽取 方式 分别为,
红球数 白球数 每次抽取方式
盒 1 90 10 随机取一球,记下颜色后不放回,而放进一个与它颜色不同的球,
盒 2 50 50 随机取一球,记下颜色后放回
盒 3 40 60 随机取一球,记下颜色后不放回,并放进一个红球,
现在如果 某人 用上述方法得到了一个记录 (红,红,红,红,白 )（ 即 5=m ）,但是不告诉我们 球出自 哪个盒 子,我们应如何推测他是从哪 个盒子中 抽取的观测样本呢?
令
=)(knS 在第 k 个盒子 )3,2,1( =k 中第 n 次抽取完成后在各盒中的红球数,
那么,在 k 分别 固定 为 １,２,３ 时,}0:{ )( ≥nS kn 分别为 Markov 链,且 其 概率 转移 矩 阵分别为,
+=?
=
=
)(0
)1(1001
)1(100
)1(
其它情形
iji
iji
pij,（ 逢红,红减 １ ； 逢白,红加 １ ）,
≠
==
)(0
)(1)2(
ij
ijp
ij,（ 内容总是不变 ）
+=?
=
=
)(0
)1(1001
)(100
)3(
其它情形
iji
iji
pij,（ 逢红不变 ； 逢白,红加 １ ）
而且初值 分别为,40,50,90 )3(0)2(0)1(0 === SSS,于是这 ３ 个盒子就分别对应于 ３ 个不同的
Markov链 模型,把这 ３ 个模型分别简记为 321,,lll,并 把 某人 观测到的样本序列中的第 n个记为 nO,即令
=nO 抽到的记录列中第 n个记录中的白球数 （ 只能 为 1 或 0）,
此例是一个玩具模型,从这个例子可以看出,在观测出自哪个盒子已知时,状态随机变
270
量序列 }{ nS 与某人提供的观测随机变量序列 }{ nO 之间的条件概率计算的关系可以直观地 写为,
mmmim SOSOSOS,,,,,,1110L,
其中在前面的一段随机变 量序列取定值的条件下,继后的那个随机变量取值的条件概率就完全确定了,而且 在这 ３ 个模型下分别都有
)|()|(),,,|( 121111 SSPSSPSOSOSP nnnnn == ++ L,
)|()|(),,,,,|( 1211111 SOPSOPSOSOSOOP nnmmmimn == ++ L,
于是
),,,,,,( 1110 mmmim SOSOSOSPL
)|()|()|()|()|()( 111201010= mmmm SSPSOPSOPSSPSOPSP L,
具体地,我们有
在 模型 1l 下 （ 把取到的球换色 ）
)|)1,0,0,0,0(),,,,(( 154321 l=OOOOOP =0.9× 0.89× 0.88× 0.87× 0.13≈ 0.08,
在 模型 2l 下 （ 球的内容不变 ）,nO 是独立同分布的随机变量序列,且
2
1
2
1
10
~nO,
所以
)|)1,0,0,0,0(),,,,(( 254321 l=OOOOOP 5)21(= ≈ 0.03,
在 模型 3l 下 （ 取到红则不变,取到 白则换为红球 ）
)|)1,0,0,0,0(),,,,(( 354321 l=OOOOOP ≈×= 6.0)4.0( 4 0.015,
再 用 Bayes 公 式 分别 得到 （ 即取上面 3 个概率的归一化值 ）
== ))1,0,0,0,0(),,,,(( 54321 OOOOOP ｜１l 0.64，
== ))1,0,0,0,0(),,,,(( 54321 OOOOOP ｜２l 0.24
== ))1,0,0,0,0(),,,,(|( 543213 OOOOOP l 0.12,
可见从第一盒抽出样本 (红,红,红,红,白 )的概率要比从其他两盒中抽出该样本的概率要大得多,
这个模型的意义绝对不在于游戏,从中可以 抽象出一种能广泛应用的数学模型 － 隐 Markov
模型 （ Hidden Markov Model,简记为 HMM ）,在这个例子中的不同盒的抽取方式可以抽象为不同的编码 方式,而这个例子 正 启示了用隐 Markov模型作模式识别的一种粗略梗概,
２,２ 隐 Markov 模型的描述
定义 １０,１３ 设 }1:{ ≥nX n 是取值于有限状态 },,1{ LL 的随机变量列,称为 状态链,1X 的分布称为其初始分布,记为 m,假定 nX 是时齐的 Markov 链,记
)|( 1 iXjXPa nnij === +,LjiijaA ≤=,)(,(10,1)
为它的转移矩阵,假定 nX 的取值,初始分布 m 与转移矩阵 A,都不能测量得到,而能测量到的是另一个与它有联系的,且可以观测到的一个取值于与有限集 },,{ 1 Mnn L 的随机变量序列 nY,称为 观测链,它们合起来还要求 满 足 如下的 隐 Markov 条件 （ HMM 条件 ）,
271
NNNNNN yiiiyiiiyii
bababxXyYP
11121111
),(
=== Lm,(10,2)
其中 N 为样本观测的时间长度,而
),,( 1 NXXX L=,),,( 1 NYYY L=,),,( 1 Nyyy L=,),,( 1 Niix L=,
NnLiy nMn ≤≤≤≤∈ 1,1},,{ 1 nn L,初始概率向量 （ 初始分布 ） 为 ),,( 1 Lmmm L=,由未知状态链与测量到的观测链一起 ),( nn YX,就构成了 隐 Markov 模型,这里,隐,的含义是说状态链是隐藏起来的,
隐 Markov 模型的基本 假定 是,
参数向量 m,参数矩阵 A 与 MLilbB ×
= )( （ },,{},,,1{ 1 MlLi nn LL ∈∈ ） 都是未知的,我们 将 它们 合 记为 参数组 ),,( BAml
=,参数组 ),,( BAml
= 完全确定了状态链与观测连的联合统计规律,所以,我们通常用一个 l 表示一个隐 Markov 模型,并 称 之为 隐 Markov 模型
l ( 更确切地为 隐 Markov 链 ),在例
10.12 中,３ 个不同的模型就对应了 ３ 个不同的参数组,只要令 1,+== nnnn OYSX,它们
满足 HMM 条件,因而纳入了隐 Markov 模型的框架,
(10.2) 是 ),( YX 的 联合分布 通过 参数表达 的形式,它是 计算各种边缘概率与条件概率 的出发点,而 HMM 条件的含义是,状态链与观测链的联合分布是由一系列简单转移与条件概率的乘积表达的,
２,３ 隐 Markov 模型的等价表述
由例 １０,１２ 的启示,可以想到下面的等价条件
命题 １０,１４
HMM 条件等价于,对于任意的 },,1( Li L∈ 以及 },,{ 1 Ml nn L∈,有
,)|(),,,,,|(
11 111111 ilnlnllnnnnln
biXvYPXiYXiYiXvYP
n
Dnn ==========
L
(10,3)
)|(),,,,,|( 1111111
1
iXjXPYiXiXYiXjXP nnlnnlnnn
n
========= ++ nn L ija?=,
(10,4)
这两个等式的证明只需利用条件概率的定义,所以我们把它留给读者,它们的直观含义就是,
nY 与 1+nX 相对于历史条件的 统计规律只与时间上最接近的 nX 有关,而与其它更早的历史无关,在实际问题中,这是比较容易判断的,
２,４ 非线性滤波作为 隐 Markov 模型 的特例
设 nX 是取值于状态空间 },,2,1{ LS L= 的 Markov 链,其 样本不能被实际测量 得到,而能测量到的是 如下的 nY 的样本
nnn wXgY += )(,(10,5)
其中 g 是一个未知函数,}{ nw 是独立同分布的 随机 干扰,只取有限个值,且 与随机过程
}{ nX 独立,那么 },{ nn YX 就是一个 （ 二维的 ） 隐 Markov 模型,这就把 （ １０,５ ） 的非线性滤波放进了 HMM 的框架,
２,５ 在应用中研究 隐 Markov 模型 的 主要 方面
(1) 从一段观测序列 },{ mkYk ≤ 及 已知 模型 ),,( BAml = 出发,估计 nX 的最佳值,
称为 解码 问题,这是状态估计问题,
272
(2) 从一段观测序列 },{ mkYk ≤ 出发,估计模型 参数组 ),,( BAml =,称为 学习问题,就 是参数估计问题,
(3) 对于一个特定的观测链 },{ mkYk ≤,已知它 可能 是由已经学习好的若干模型之一所得的观测,要决定 此观测 究竟 是得自其中哪一个模型,这 称为 识别问题,就 是分类问题,
在实际问题中,这 3 个问题 有时并不完全 能分开,有时也并不 需要 解全 3 个问题,例如,
在语音识别或手写体汉字 或数字的 脱机识别中,我们只需要作 (2) 与 (3),这相当于将一个,标准的,语音 音素 或一个,标准的,手 写 体 汉字,学习成,一个隐 Markov 模型,即把它与 一个或几个 特定的隐 Markov 模型 ( 更确切地说,是特定的一组参数 ) 相 对应起来,以便 把该模型作为这个音 素 或 手写 字的代表 模板,这是 学习相位 ; 而在进一步用这些模板中的最合适者,作为对于一 个需要识别的 音 素 或 手写 字 的分类归属,这是 运转相位,至于 归入那个模板最合适,就要用合理的距离,或准距离 （ 常用的是相对熵 ）,在此意义下优化,
[ 注 1] nX 也可以推广为,取值于平面有限格点的 Markov 场,而 nY 与它的关系仍如上述,
这就定义了一个 隐 Markov 场,
[ 注 2] 一般地,nY 还可以是连续型随机变量,如果 记 在 xX n = 的条件下,nY 的条件分布密度为 f,则 ),,( xfAml = 就也称为一个 （ 连续的 ） 隐 Markov 模型,这时 f 常是分布类型已知而带有未知参数的密度,状态是连续的 隐 Markov 模型至今还 未见 有人使用,
３ 解码问题 -- 已知模型 l 与观测 yY = 时状态 X 的估计
３,１ 出现当前的观测的概率 )|( lyYP = 的计算
我们仍旧 沿用 记号
),,( 1 NXXX L=,),,( 1 NYYY L=,),,( 1 Nyyy L=,),,( 1 Niix L=,
NnLiy nMn ≤≤≤≤∈ 1,1},,{ 1 nn L,初始概率向量 ),,( 1 Nmmm L=,
由 (10,2),利用条件概率的性质 容易算出
NN yiyiN
bbiiXyYP LL
11
)),,,(|( 1 === l,
NNNNNN yiiiyyiiyii
bababxXyYP
1121111
)|,(
=== Lml,
NNNN
N
yiiiiiyii
ii
baabyYP
121111
1,,
)|(
∑
== L
L
ml,(10,6)
对于 Nn ≤≤1 及观测样值 yY =,记 （ 因为观测样值 y 是固定的,所以下面我们将在足标中把它略去 ）
)|,,,()( 11 la iXyYyYPi nnnn ==== L ( 依赖 y ),(10,7)
则 在模型 l 给定下,关于 观测资料 ),,( 1 nyy L 的长度 n )( Nn <,我们有递推公式 ( 称为 向前递推公式 或 向前算法 )
1
)()(1
+∑
=+
niyj jinn
bji aaa,(10,8)
由此得到
(1) 基于 )(ina 的向前递推公式 计算 )|( lyYP = 的 步骤,
计算初值
1
)(1 iyibi ma =,
273
用递推公式
1
)()(1
+∑
=+
niyj jinn
bji aaa,(10,9 )
最后 得到结论 )()|( iyYP N
i
al ∑==,(10,10 )
计算 )|( lyYP = 也可以通过另一种途径,为此 记 （ 观测样值 y 也是是固定的,我们也把它在足标中略去 ）
),|,,(()( 11 lb iXyYyYPi nNNnnn ==== ++ L,(10,11)
则在模型 l 给定下,关于观测资料 ),,( 1 nyy L 的长度 n )( Nn <,我们 还有递推公式 ( 称为 向后递推公式 或 向后 算法 )
1
)()( 1
+∑ +
=
njyj ijnn
bji abb
这就得到计算 )|( lyYP = 的另一个算法,即
(2) 基于 )(inb 的向后递推公式 计算 )|( lyYP = 的步骤
定义 1)( =iNb,
用递推公式
1
)()( 1
+∑ +
=
njyj ijnn
bji abb ( 10,12)
最后 得到结论 ∑==
i
iyibiyYP 1)()|( 0 mbl (10,13)
３,２ 解码问题 -- 已知模型 l 与观测 yY = 时状态 X 的估计
令
),|()( lg yYiXPi nn ===,(10,14)
那么
∑∑
=
==
===
i
nn
nn
n
i
n
n
ii
ii
iXyYP
iXyYPi
)()(
)()(
)|,(
)|,()(
ba
ba
l
lg,(10,15)
( 推导的思路 如下,令 F n ),,( nmYX mm ≤= s ( 右 方 表示括号内随机变量 所 生成的 s 代数 ),
},,{ 11 nnn yYyYY ===? L,},,{ 11 NNnnn yYyYY === +++ L,因为 l 是固定的,我们在下面记号中 也把它 （ 为了记号的方便,我们在下面 有时并 不区分取概率的符号与取期望的符号,取概率 也 理解为 取期望 ）,于是 我们有
)|,( liXyYP n == |,([( iXyYPE n === F n )],|(,,[ iXYPiXYE nnnn === +? F n )]
)]|(,,[ iXYPiXYE nnnn === +? )|(),( iXYPiXYP nnnn === +? )()( ii nn ba= ),
对于解码问题的解决,我们用如下两种考虑 方法,
1,对 单个 时刻的 状态的最大概率估计
如果 *i 满足 )(max)( * ii nNin gg ≤=,则取 *
^
iX n =,
单个 时刻的 状态的最大概率估计的长处在于简单易算,但是它忽视了状态链在不同 时刻之 间的联系,也 就 不能充分地利用已知的信息,下面的 路径估计法,是在考虑不同时间的联系后的改进 方法,
274
2 路径的 ( 轨道 ) 最大概率估计 -- Viterbi 算法
这是基于动态规划的思想的一 种整体性考虑的方法 （ 所谓动态规划实质上 就 是一种从后向前 逐步取优的优化方法 ）,令在时刻 n 取到状态 i 的 各个 路径 的最大 概率为
)|,,,,,,(max)( 111111,,
11
ld yYyYiXiXiXPi nnnnniin
n
======
LLL,(10,16)
那么我们有
))((max)(
11 jinjiyn
ajbi
n
dd
+
=+,(10,17)
用来 计算 它的实际 程式 是 如下的 Viterbi 算法,假定最佳路径为 ),,(
^
1
^
NXX L,今设置一
个从 1
^
+nX 得到 nX
^
的对应 )(,1
^
1
^
1 +++ = nnnn XX yy,计算 )(ind 与 ny 的递推算法 的 具体 设计 如下,
设置初值 1)(1 iyi bi md = ;
用递推公式 ))((max)(
11 jinjiyn
ajbi
n
dd
+
=+ ;
再 定义 *1 )( jin =+y,其中 *j 只要 满足
ijniyn ajbi n *1 )()(
*
1 dd +=+,(10,18)
最终采取如下的向后递推算法 ( 相当于动态规划,动态规划的核心思想是,在一切最佳
路径上的任意节点出发,其后也应该是最佳路径,因此就能 从后 向 前 地逐步求得最佳路径 )
得到 如下的,最佳,轨道估计 ),,(
^
1
^
NXX L,
计算终值 )}(max)(:{*^ jiiiX NjNN dd =∈=,(10,19)
反向归纳地取 )(
^
11
^
++= nnn XX y,(10.20)
４ 学习 问题 – 由观测 yY = 估计模型参数 l
４,1 状态链样本已知时的参数的频率估计
设 隐 Markov 模型的 状态链 有 充分长 的样本,把 其中从状态 i 到 下一个时刻为 状态 j 转移的频数 记 为 ijA,那么,粗略地 可取
∑
=
= N
j
ij
ij
ij
A
Aa
1
^
为 ija 的估计,同样 设此观测链样本充分长,记其中 从状态 i 到同一时刻的观测 l 转移的频数为 ijB,则可取
∑
=
= M
l
il
ilil
B
Bb
1
^
为 ilb 的估计,不幸的是状态链是经过估计得到的,所以不修正地使用频 率估计会增加误差,再则,这种估计 也 过于敏感,
４,2 模型参数估计的 EM 算法的思想
我们取 l 的最大似然估计
^
l,即它 满足
)|(max)|(
^
ll l yYPyYP === (10.21)
275
但是,对 于似然函数
NNNN
N
yiiiiiyii
ii
baabyYP
121111
1,,
)|(
∑
== L
L
ml 而言,这个最大值问题 的 计算量过大,在实际中是不可能被采用的,为此,人们灵活地 利用 使似然函数尽量地大的思想,
采取折衷的方案,构造一个递推算法,使之 能相当合理地给出模型参数 l 的粗略估计,其核心思想是,并不要求备选估计 l 使 )|( lyYP = （ 关于 l ） 达 到最大或局部极大,而只要求使
)|( lyYP = 相当大,从而 计算变为实际可能,
为此,我们 定义一个 描述模型,趋势,的量 )|( * llQ 代替似然函数 )|( lyYP =,其定义 为
)|,(ln)|,()|( ** llll yYxXPyYxXPQ
x
=====∑?
)|,(ln)|,( *∑=
x
yxPyxP ll
简记为
(10,22)
( 这个量并不是相对熵,但是在形式上很类似 ),利用在 10 << x 时的不等式 1ln?≤ xx
( 证明,令 xxxf ln)1()(=,则 0)1( =f 且 01)(" 2 >= xxf,即 )( xf 为凸函数,从而 0)( ≥xf ) 得到
∑=?
x yxP
yxPyxPQQ
)|,(
)|,(ln)|,()|()|( **
l
llllll
)|()|()1
)|,(
)|,()(|,( ** ll
l
ll yYPyYP
yxP
yxPyxP
x
=?==?≤ ∑,
由此 可见对于固定的 l,只要 )|()|( * llll QQ >,就有 )|()|( * ll yYPyYP =>=,
于是 想 把模型 nl 修改为更好的模型 1+nl （ 即找 1+nl 使 )|()|( 1 nn yYPyYP ll =>= + ）,
只需找 1+nl 使
lll sup)|( 1 =+ nnQ )|( nQ ll,(10,23)
即只要把 )|( nQ ll 关于 l 的最大 值 处取成 1+nl,就 有 )|()|( 1 nn yYPyYP ll =>= +,
这样得到的模型序列 }{ nl 能 保了 )|( nyYP l= 关于 n 是 严格递增 的,虽然 在这里 还不能在理论上证明 )|( nyYP l= 收敛到 )|(max l
l
yYP =,但是当 n充分大时,nl 也还能提供在实际中较为满意的粗略近似,
综上论述,我们 把如 上得到近似模型列 nl 的方法 归结为两个步骤,
1 E 步骤 （ 求期望 ）,计算
=)|( * llQ )|,(ln)|,( *∑
x
yxPyxP ll ；
2 M 步骤 （ 求最大 ）,求 1+nl 使
lll sup)|( 1 =+ nnQ )|( nQ ll 。
276
这两个步骤合起来构成的算法,称为 EM 算法,EM 算法是针对 在 测量数据不完全时,求 参数的一种近似于最大似然估计的统计方法,隐 Markov 模型 的 建 模 （ 模型 参数 的 估计 ）,是 EM
算法的一个最常见,且 极有用的 一种 典型例子,下面将看到,隐 Markov 模型 中的 M – 步骤的解可以有显式表示,这就 是 一组把模型参数修改为新的模型参数的递推公式,这组公式 正好是在 隐 Markov 模型 中普遍应用的著名的 Welch - Baum 公式,下面我们给出 得到这组递推公式的 梗概,首先我们 有下面的引理,
引理 １０,１５ 设 )(0 Nizi ≤≥,则在约束条件 ∑
=
=
N
i
ix
1
1,)(0 Nixi ≤≥
下,函数 ∑
=
N
i
ii xz
1
ln 在
∑
=
= N
i
i
i
i
z
zx
1
)( Ni ≤ 处取得最大,
证明 注意 ),,(1 1
1
NN
j
j
zz
z
L
∑
=
对 ),,( 1 Nxx L 的相对熵 0
)(
ln)( 1
1
1
≥
∑
∑ ∑ =
=
=
i
K
j
j
i
K
i
K
j
j
i
x
z
z
z
z,
４,3 隐 Markov 模型 中 M – 步骤的 求 解
现在我们用引理 10,1,来求 )|( nQ ll 关于 l 的最大值的位置 1+nl,注意
NNNN yiiiiiyii
baabxXyYPxyP
121111
)|,()|,(
====? Lmll,
∑ ∑
=
=
+
++=
N
n
N
n
iiyii nnnn abxyP
1
1
1
11
lnlnln)|,(ln ml,
由此得
=)|( * llQ )|,( lxXyYP
x
==∑ )lnln(ln
1
1
1
***
11 ∑ ∑
=
=
+
++
N
n
N
n
iiyii nnnn abm
)(),(ln)|,(ln }{
1
*
1
1 nln
Nn
il
M
lii
i yIiXyYPbiXyYP ==+=== ∑∑∑∑
≤≤=
lm
∑∑?
=
+ ===+
ji
N
n
nnij jXiXyYPa
,
1
1
1
* )|,,(ln l ∑ ∑++=?
i i
iaib aQbQQ
ii
)|()|()|( *** lllmm,
此处 带 * 的参数都表示修改 后 的参数,而
**,
ii ba 分别表示待修改的参数 BA,的第 i列列向量,显见 ∑ =
i
i,1
*m ∑ =
j
ija,1
* ∑ =
l
ilb 1
*,这样可以通过分别最大化上面等式右边各项得到 将模型 ll =n 更新 为
*
1 ll =+n 的修正 模型 1+nl,而其中每一项都可以用引理 １０ 。 １５
最大化,这样从最大化 就得到 如下的结论,
Baum － Welch 算法 （ 递推关系 ）,
277
)(
)|(
)|,(
1
1
^
)1( i
yYP
iXyYP
n
nn
i gl
lm =
=
===+,(10,24)
∑
∑
∑
∑
=
=
=
=
+
+ =
==
===
= 1
1
1
1
1
1
1
1
1^
)1(
)(
),(
),|(
),|,(
N
n
n
N
n
n
N
n
nn
N
n
nnn
n
ij
i
ji
yYiXP
yYjXiXP
a
g
V
l
l 记为
,(10,25)
∑
∑
∑
∑
=
==
=
=+ =
==
==
= N
n
n
N
lyn
n
N
n
nn
N
n
nlnn
n
il
i
i
iXyYP
yIiXyYP
b n
1
,1
1
1
}{^
)1(
)(
)(
)|,(
)()|,(
g
g
l
l 记为
,(10,26)
其中
∑
=
i
nn
nn
n
ii
iii
)()(
)()()(
ba
bag 如 (10,14) 与 (10,15),而
∑
+
+
+=====
i
nn
njyijn
nnn
ii
jbaiyYjXiXPji
n
)()(
)()(),|,(),( 1
1
1
ba
balV
,(10,27)
（ (10.27) 的推导需要先算 )|,,( 1 ljXiXyYP nn === +,它与计算 )|,( liXyYP n == 的方法几乎完全一样 ),此时还有
∑=
j
nn jii ),()( Vg,(10,28)
在上面所述的算法中,初始 值 0l 的设置会直接影响到估计的好坏,为此常用的 一种方法是,根据先验知识设置一条较长的,标准,的虚拟状态链 },,{ )0()0(1 NXX L,用 ４,０ 段 中的频率估计粗估一个 l,并 把它取为 0l,
对于隐 Markov 模型,甚至对于第 15 章中的一般 EM 模型,观测链提供了所有的已知信息,
但是,与一般的 EM 方法相比,HMM 远为简单,因为 如果我们 仔细考察 B aum - Welch 算法,就可以发现它并不依赖于状态 链 的估计,这就是说,隐马氏模型中的参数估计与状态链的估计 之 间是 解 耦的,而一般 EM 算法 常常 都是 有 耦合的,即在估计参数时需要知道状态的估计,而在估计状态是又需要知道参数的估计,
在向前与向后两个递推公式基础上的 B aum - Welch 算法,给出了估计模型参数的合理方法,
虽然 在理论上 还没有 证明 其极限 就是最大似然估计,
[ 注 1] 可以对 HMM 条件作一点小的推广,假定 对 },,{},,,1{ 1 MlNi nn LL ∈∈ 有
),(),|(
),,,,,|(
1
111111 11
lbiXjXvYP
YiXYiXiXvYP
ijnnln
llnnnnln n
=====
======
nn L
(10,29)
)|(),,,,,|( 1111111
1
iXjXPYiXiXYiXjXP nnlnnlnnn
n
========= ++ nn L ija记为 =
(10,30)
与 n无关,在实际问题中常常会遇到这一种情形,
[ 注 2] 在 HMM 中,nn YX,都还可以是多维的 向量,这时的矩阵 BA,就 相当 复杂,
278
[ 注 3] 如果模型中的时间参数 n 相应地变为地点参数,那么就导致 二维的隐 随机 场 模型,特别对于给定的 －? 领域系统,就有隐 －? Markov 场模型,此时 B aum - Welch 算法就无法进行了,这时的 学习问题可以,例如用如下的 交错估计 方法 ( 这个方法是 耦合的,在这点上它与 B aum - Welch 算法不同,我们仍称这个算法为 推广的 V i t er b i 算法 ),
(1) 先在某个
)(nll =
固定下,用 推广的 V it er bi 算法得到优化的状态轨道 ;
(2) 再 求
^
l,使 ll sup),|(
^*
=== xXyYP ),|( lxXyYP ==,其中
),,( **1* nxxx L= 是随 n增加时轨道的最佳估计,也可以使用轨道频率,
(3) 令
^)1(
ll =+n,把
)(nl
更新为
)1( +nl;,
(4) 如此反复 迭代,
与之联系的有以下一些问题,
问题 1,用推广的 Vertabi 算法频率算法代替 Baum - Welch 算法后,解决了隐 Markov 场时,
Baum - Welsh 算法,在二维难于递推的困难,但是如何给出递推公式? 能否在某些条件下,得到递推更新中每次修正都往好的方向发展的结论,例如,使相对熵增加?
问 题 2,如果问题 1 中的估计结果较好,那么用一个随机取的 x 代替 ),,( **1* nxxx L=,
是否还能接受?
问题 3 隐 Markov模型中,如何引进 Bootstrap 式思维,
[ 注 4] 在生物信息论 中,还常用 隐半 Markov 模型,这时的做法可以是,
(1) 假定一个初始模型 (0)l ( 包括初分布,停留时间,转移参数 ),由观测资料通过
Viterbi 算法估计轨道 ;
(2) 由此轨道通过计算频率来更新初分布,停留时间,转移参数,得到更新模型
(1)
l ;
(3) 重复步骤 (1) 与 (2),直至较为稳定,
在语音识别中,也有人使用 隐半 Markov模型,但是与用隐 Markov模型相比较,在评估是否得到显著的改进这一点上,还颇多争议 。
５ 关于 隐 Markov模型 的评注
５,１ 隐 Markov模型 包容度大有非常宽的应用面
用隐 Markov模型 建模的优点
（ 1 ） 隐 Markov 模型 是一种简单的数学模型,它的包容性很大,内涵很广,在应用中的弹性相当大,可以很好地利用我们对于被建模的对象所了解的先验知识,因而有很广的适应 性,
（ 2 ） 隐 Markov 模型 是一种不完全数据的统计模型 （ 包含了非线性滤波模型 ）,它 之所以被广泛采用,在于这种模型既能反映对象的随机性,又能反映对象的潜在结构,便于利用 对象的结构与局部联系性质等方面的知识,以及 对研究对象的直观的与先验的了解,
（ 3 ） 隐 Markov 模型是 很少的几个既从物理模型出发,又与数据拟合相直接联系的算法之一,
（ 4 ） 隐 Markov 模型 虽然也具有某些黑箱的特点,但是与纯黑箱操作的人工神经网络算法相比,有明显的优点,它的参数常常有较为实质的含义,而线性模型,时间序列模型 中的参数一般只是作为被拟合了的参数而出现的,缺少 较为 真实含义,
（ 5 ） 隐 Markov 模型 有快速而有效的学习算法,
以上诸点使 隐 Markov 模型称为 随机建模与拟合的重要基石,在应用中有其普适性,
279
在具体应用 隐 Markov模型 建模时,首先要设定 Markov 链的状态集及其规模,即总状态数
L,这也有相当大的弹性,然后,确定相应的观测过程,为了使得观测链 ( 或场 ) 与马氏链 ( 或场 )
之间有隐 Markov 关系,我们常常需要对实际提供的观测链 ( 场 ) 加以改造,例如,一个脱机手写体汉字是一张黑白象素图 ( 一个 取值 0 或 1 矩阵 ),它和刻划笔道的 －? Markov 场之间的依赖关系并不能满足隐 Markov 关系的要求,需要对它进行 预 加工,也就要 先 进行,特征提取,,
以尽量简化得使隐 Markov 关系能够成立,
在确定了模型规模以后,对于隐 Markov模型的处理,又分为两个不同的相位,学习相位与运转相位,学习相位是由给定的观测过程的一组样本出发去估计隐 Markov 模型 的参数
),,( BAml =，,从而完全确定模型 ; 而运转相 位则是用学习相位所 确定的模型参数,对给定的某个观测过程的样本,来估计相应隐 Markov模型的状态,或计算出它在各个隐 Markov 模型下出现的概率,
语音识别和脱机手写体汉字识别最终需要知道的是给定的语音或脱机手写体汉字样本应出于哪个模型 (哪个 参数组 ),这是一个识别过程,通常我们可以按照 Bayes 方法 （ 后验概率 方法 ）
去决定模型,即选取使得给定样本出现的概率最大的模型作为它的分类归属,这就是识别的结果,而在 DNA碱基序列的排序的问题中,考虑的方式则与之不同,在 DNA碱基序列中,需要知道的是如何排序,也就是需要知道给定的某个观测过程的样本所对应的隐 Markov 链的状态,
在学习相位中,我们面临的是一个不完全数据的参数估计问题,除了 Baum –Welsh算法以外,还 可以进行交错估计,即可先设置一组初值参数,由此出发配合观测列按 Bayes 方法 去估计马氏链的状态列,再用此状态列配合观测列用最大似然估计去重新估计模型参数,交替地不断重复此二步骤,…,以达到比较稳定的最终结果,也可以从设置马氏链的一个状态列作为开始,去估计参数组,再由此参数组对应的模型去估计状态,如此交替进行,
在应用 隐 Markov 模型 于语音识别或手写体 汉字识别时,往往把一个音,一个词,一个字对应于一个隐 Markov模型 l,这时当然也可以把一个标准音或标准字的状态设置为马氏链的初值,再开始上述交替估计的过程,
值得指出的是,利用已知的先验知识,来得到简单的递推计算公式,是使得 隐 Markov 模型的计算能够 化简 的关键,在实际应用中,参数组的维数常常非常高,这就使得计算无法实际进行,
所以在交替估计的过程中,我们还可 以借助 Gibbs 采样的思想,一次只变化一个状态分量或一个参数,以简化计算,也常常使用 Markov 链 Monte Carlo 方法,
对于隐 Markov模型的状态是否取得合理,它们是否包含了所研究的问题的必要的信息,观测过程是否包含了认识模型的足够信息等问题,也是使用中必需了解的,怎样从我们能够得到的信息 (通常主要是样本数据 )来判断上述问题,在统计上也是 至关 重要的,
５,3 隐 Markov 模型 的更为一般的形式
隐 Markov 模型 的观测链 }{ nY 还可以是连续型随机变量,其典型情形有正态分布,Gamma
分布 ),( lkΓ,椭圆对称分布 ( 密度具有 ))()((
||
1 1
2
1 mm?Σ?
Σ
yyf T 形式 ),或上述各种分布的混合分布,例如 ),()1(),( 22221 smlsml NN?+,21 exp)1(exp mm ll?+ 等等,
更一般地,隐 Markov 模型 中的一个状态甚至还可以是一个整个的随机过程,
［ 注 ］ 隐 Markov 模型在应用方面的参考文献量十分大,而它在理论方面的参考文献还不多,例 如关于大数 定 律 的有
Brain G Leroox,Maximum - liklihood estimation for HMM,Stoch astic processes
and their Appl,40(1992),127 - 143,
关于中心极限定理 的有
(1) P,J,Bickel and Y,Ratof,Inference in hidden Markov models I,(local
280
asymptotic normality in the stationary cas e ) Bernoulli,Vol.2,No.3,199 - 228,
1996,
(2) P,J,Bickel,Y,Ratof and T,Ryden,A symptotic normality for the maximum
likelihood estimates for general hidden Markov models,Annals of Statistics,Vol,
26,No.4,1614 - 1635,1998,
６,隐 Markov 模型的应用例 子梗概
６,１ 语音的机器识别
语音识别与下面 ６,２ 段的 脱机手写体汉字识别,这两个问题都有一个共同的特点,它们都有随机性,但 是 也都有一个潜在的基本结构,例如,一个语音也有其基本结构,而 发音时口型的大小,长短,强弱,在口腔中的位置等,却因人而会有随机的变化,即使是同一个人,在不同时间发同一个音也是有随机差异的,而一个汉字的基本笔画及相对关系是确定的,但是笔画的长短,间隔与倾斜度却因人而异,用隐 Markov 模型来描述这种具有结构的随机性相当有效,
这里 以简单化了 的 实际问题为例,这样 可以 使我们摆脱许多细节而 突出 思路,一个词的语音可由许多音素组成,如前所述,隐 Markov 模型在语音识别中,首先是建立一个对应关系,
使一个音对应于一个隐 Markov 模型 ( 或几个 隐 Markov 模型 ),这里的隐 Markov 模型的状态,
一般取为这个音所含的全部可能的音素 ( 或其细分,或其组合 ),例如一个汉语,高,的发音可能由
‘ggg(ge)(ga)(ga)(ga) aaa(ao)(ao)ooo’ (10.31)
( 后面可加上足够的 □ ( 空 )) 等音素组成,而且每个音素究竟维持多少时间单位,对不同的人,
对同一人在不同的时刻 发 这个音,都可能有很大的不同,于是对应于,高,这个音,是 一列相继的状态组成的一个现实 ( 在上面序列中一个字母代表相应汉语拼音,而一个括号中有两个字母则表示过渡 ),由于声转移有内在统计规律,我们 就可以自然地把它考虑为一个 Markov 链
)( nX 的状态序列,但是,)( nX 实际上并不可观测,即我们能测量到的并非这个音对应的
Markov 链本身,而是 每个字母对应的声波 信号的振幅 的采样 ( 或经过某种 预 处理的一些特征量 ),这个可以 观测 到 的现实就是 对于一个音的观测链样本列,状态链与观测链在一起就组成了一个隐 Marlov 模型,其中的状态序列,就是该隐马氏模型的马氏链的轨道 ( 包括,□
（ 空 ),),所以说一个音对应于一个隐 Markov 模型 ( 由于发音人音质的差异,一个音也可以对应与多个隐 Markov 模型 ),一般在语音识别中,把音频信号直接作为观测值常会使计算时的数据量过大,更为实际的是把音频信号的一些特征量,作为压缩了的观测数据,这称为 预 处理,
这些特征可以是音频信号 的过零率,线性预测系数 (LPC,即 Linear Predict Coefficient,
它是把采样值作为 AR 模型 （ 参见第 １１ 章 ） 拟合 后,得到的模型的系数 ),倒频系数 ( CEP,记信号 )(ns 的离散 Fourier 系数 为 ))(( nsF,则 |)))((|(ln1 nsFFCEP?
= ) 等等,一般地常取 12
至 16 个 LPC 作为压缩了的观测数据,也就是把原始观测数据 用 AR( p ) 模型 ( 其中 p = 12
至 16 ) 建模,再 把这样拟合的参数作为压缩了的观测数据,( 为了保证状态到观测的转移只依赖于前一个状态,有时还需要把邻近的状态,堆,起来,例如说,把相继出现的 k 个状态看成为一个状态向量,并把与此 k 个状态向后错一位的相继出现的 k 个状态,看成为其后时刻的一个状态向量,…,使一个时刻的,状态,或,观测,成为一个向量 ),再则,有时也常把在状态已知的条件下,观测向量的条 件分布假定为 ( 一维的或 多维 的 ) 正态分布 ( 一般假定方差矩阵为对角矩阵 ),总之,我们可以把一些特征量作为观测链 { nY,n ≥ 1},其中 nY 可以是多维的 ( 向量 ),序列 ( 10.31 ) 并不能观测到,它正是我们想要知道的状态链,为了解它,就要
281
由对该音的一组观测样本 ( 若干个发这 一 音的声信号或提取的一些特征量 ),去作相应的 隐
Markov 模型的 参数 ),,( 0 BAml = 估计,这就是学习相位,于是语 音的识别 对应于 设置一些模板,一个音对应于一个或数个模板 l ( 有时也可以用最佳状态列的估计作为音或词的代表 ),
这时 模板的代表性,选多少个状态,选多少个观测,对于计算量的控制与,都至关重要,为了能 有效地学习,就需要有足够多的样本,也就是要有这个音或词的重复发音,例如把 m 个不同的人的对同一个音或词的发音的 预 处理加工,作为 m条观测链样本,学习的结果是建立了标准人的 发音与隐 Markov 模型间的对应 ( 也可以让一个音或一个词对应 于 几个隐马氏模型 ),
以 得到一个模板库,
在学习得到模板后,就可以作语音的机器识别问题,即对于一条新给的观测链,在一族 (可能是上千个 )已建立好的隐 Markov 模型的模板中,选取一个与之,最近的,(相对熵最小的 ) 隐马氏链模型 l,它所代表的音或词,就可以判定为这个观测所对应的音或词,这就是 识别 问题,
隐 Markov 模型的基本假定是,状态链是时齐的 Markov 链,而且状态到观 测的转移也是时齐的,但是,在实际语音识别中,常会遇到非时齐的情形,这时如果 照 搬时齐情形的做法,就会引起振荡,然而,用非时齐的隐马氏模型来建模,无论在理论上或 对 实际计算,都远非可取,所以运用隐 Markov 模型 时 极需要灵活处理,实际上,在实践中遇到的问题,绝非只是理论的直接应用,
其复杂的程度 常常 不 比理论问题 差,
如果状态链在一个状态上的停留时间分布,显著 地 不是指数分布 ( 离散采样时它对应于几何分布 ),那么在理论上可以想象,应该改用隐半马氏模型代替隐马氏模型,可惜的是,这种修改后 的模型 的 计算 远为复杂,而在实际使用中,这两个模型的差别也并不如我们想象的那样大,因之这种修改也并没有被广泛采用,
６,２ 脱机手写体汉字识别
汉字的模板称为码本,在线手写体汉字或数字可以看成一维的,用 隐 Markov 模型 的框架建模 就较为简单,而脱机手写体汉字是二维的,即在隐 Markov 模型中的时间参数 n 变成为二维的空间参数 ),( mn,而二维的 Markov 性就远较一维复杂,因为它没有很自然的序,在 确 定相应的观测过程时,为了使得观测链 ( 场 ) 与状态链 (Markov 场 ) 之间有象 ( 10,3) 及 ( 10,4 ) 那样的简单关系,我们常常需要对实际提供的观测链 ( 场 ) 加以改造,例如,在实际上我们得到的手写体汉字是一张黑白象素图 ( 一个取值 0 或 1 矩阵 ),它和笔画的隐 Markov 场之间的依赖关系不能满足 ( 10,3 ) 或 ( 10,4 ) 的要求,我们要对它进行加工,也就是,特征提取,,尽量简化 得 使
( 10,3 ) 或 ( 10,4 ) 能够成立,为了使脱机手写体汉字纳入隐 Markov 场的框架,选取什么作为状态与观测才能较好地体现 Mar kov 性,这是具体操作中的 最 重要的,粗略地说,一般可以选取一个格点上八个方向的笔道作为观测的特征,
此外,数字识别,文字识别,印刷中的 错别字自动检测,航海中的 动态手势识别等等,
也可以用隐 Markov 模型处理,只是有的 情形的 "时间参数 "n 应该 相应地改为 空间参数,
６,３ DNA 序列片断装配及启动子识别
在遗传学中,生物体属性的遗传是由细胞核中的染色体中所含的基因所指导,基因是由
DNA (脱氧核糖核酸 ) 碱基的序列组成,如果把基因看成语句,DNA 的 A,T,G,C 4 种碱基,就象语句中的字母,一条 DNA 链就是一个 A,T,G,C 这 ４ 个 字母组成的序列,读出 DNA 碱基序列对于基因工程,遗传病及许多癌症的发现与治疗都具有重要意义,然而,人类基因组计划目的并不是解决,人类由甚么组成的,,而是要回答,是甚么使得一种有机体区别于另一种有机体,,一个 人 的 基因组的全序列,是一个 包含约 30 亿个由 A,T,C,G组成的序列,它们分为 23
对染色体,与一条性染色体,生物体中最小的 DNA链,约含一百万碱基对 (bp,base pair),即由
A,T,C,G 4 个字母组成的长度为一百万的一个序列,基因是染色体 DNA序列的一些子序列,它们给出了 有关 蛋白质的信息,至今,人们知道在基因之间还存在 DNA 序列,可是还不清楚它们的作用,一个人的基因组的全序列可以说是他的生命的设计书,它能为了解遗传与发展等生命
282
奥秘奠定基础,而 得到基因组的全序列,只是我们读懂这本生命之书的起始的一步,生物信息论的第一步是 DNA 序列的分析,包括 DNA 序列的比对 (alig n ment),拼接 (assembly); 寻找基因 ( gene findi ng ),调控研究 (regulation),基因功能 (function)研究等,
(1999 年 12 月 1 日,人类基因组科学家完成了人类 22 号染色体的测序,这是在 23对染色体中第二小的染色体,有 3350 万个碱基对,由长短不同的 12个 DNA 片断组成,最长的片断有 2300 万碱基对,其它的片断,长度在
1000 碱基对到 58.3 碱基对之间,在 22 染色体上发现了 679个基因,其中 545 个是功能基因,与先天心脏病,低免疫性,精神分裂,智力低,许多恶性肿瘤如白血病等有关,另外的 134 个不是 功能基因,2000 年 5月 8 日,人类基因组又基本上给出了第 21号染色体的基因图谱,这是最小的染色体,有 225 个基因,与白血病,痴呆,肌肉萎缩,
燥狂性抑郁,部分癌症等有关,人类基因的总数比水稻要少,说明人类的基因可能更为有效 ),
然而,DNA 序列 是 非常长的链,用 目前的设备是无法读出的,为了读取一条链,就 需要 线把它 打断成一段一段,使每一段都可以用设备读出,而在解读 DNA序列时,遇到的困难是,
(1)目前快速读出 DNA序列的设备只能读出几百到一千个碱基对的链,
(2)快速打断 DNA序列的方法是随机的,不能控制断点的位置及段数,
(3)打断后的片断混为一体无法区分次序先后,
解决困难的一个办法是,同时对同一个 DNA序列的多个拷贝 (copy) 作随机打断及作随机片断读出,利用断点几乎概率为 1 地不会在同样位置重复,并利用不同的拷贝片断的相互核对,将片断连接起来 (Assembly)得到整体的排序,但是,在 读出时可能有错读,漏读及误增字符等错误,
且在原序列中还有重复的片断,这就使拼接问题十分复杂,上世纪末的人类基因测序中,人类基因组的科学家对染色体 DNA 序列,利用一些只出现一次的序列 ( 称为 marker),将整条序列分段 ( 段长为 2M - 3M) 读出并拼接起来,以降低拼接的复杂性,而 Cerela 公司的科学家则更着重于技术与算法,他们利用最新的毛细管测序机,可以从一段序列的两头同时测得两个
500 - 700bp 小段 DNA 序列的这个特点,将人类基因组科学家拼接思路中忽视的匹配 ( mate ) 信息利用起来,为不用 marker 直接测序开辟了新路,
读出了全部基因组的 DNA 序列后,基因检测是后基因处理中的首要问题,即要知道序列中各部分的含义 ( 如,哪里是基因,哪里是基因的蛋白编码部分,哪里是其调控部分,…),不知道这些,DNA 测序的结果不会有用,每个基因有几部分构成,启动子 ( Promoter ) – 像是基因的引言 ; 外显子 ( Exons ) – 编码部分,它在基因 DNA 序列中也可能是被内含子 ( introns ) 分隔的子序列 ; 内含子 ( Introns ) – 外显子间的序列,它的功能至今还不清楚 ; 终止子
( Terminator ),它们可示意为,
启动子 外显子 1 内含子 1 外显子 2 内含子 2 外显子 3 终止子
对于 已 读出的 DNA 序列,Karlin 和 Berger 有一个基因检测的软件,称为 Genescan,利用基因编码区的统计特征,即在没有漏读 ( 也称 删除,记为 D)与没有 增 读 (也称插入,记为 I)时,
利用 DNA 编码序列具有很强的 3-连子特性 (A,T,C,G 4 个核酸组成 64 个 3 – 连子,即 3 个 相连的 字母 组成 的 一个 编码,在此 ６４ 个编码 中除了 3个编码代表终止符号以外,其它的 ６１ 个编码 分别代表 20个氨基酸 (amino) 的编码 ( 一个氨基酸可以有几个 不同的 编码 ),在编码区无漏读或无 增 读 的条件下,Karlin 和 Berger 采用了动态规划方法,编制了基因检测软件 Genescan,
达到了很高的识别率,但是在编码区有漏读或有 增 读时,就会发生 3 - 连子的错误,且 随后的所有 3 - 连子就全 被 搞乱,接着 就 会 出现一系列的错误,如果还用 Genescan 就会发生严重的错误,所以,在用这个软件前,必需先在观测序列中,标识错读,漏读 ( 删除 ) 或 增 读 ( 插入 ) 的碱基对,以便避免出现大量的编码错误,并尽可能的校正错读的碱基对,
对于标识错读,漏读或增读碱基对,G.Churchill 首次提出了 用隐 Markov 模型解决 DNA
片断装配问题,他的基本思想是将由随机枪打断长链中的每个碱基对,附以读出的属性,即读对 R,错读 M,漏读 D 或增读 I ( 其中最误事的是漏读 D 或增读 I,因为它们引起一系列后继错误 ),把这些属性看成 一个 Markov 链的 状态的 样本,而把读出的碱基序列片段看成为相应于这个 状态 链的观测链,即 设观测过程为可能带有小概率的 漏读 D与 增读 I （ 在实际解读中,漏读与增读度不可能很多 ） 的 DNA码区序列,这时,隐 Markov链的状态 空间就是 },,,{ IDMR,
其中 R 表示读对的碱基对,M 表示读错的碱基对,而隐马氏链的样本列,就是对观测序列上的
283
碱基对标出的属性 MDR,,与 I 的序列,这样的考虑 是,基于在编码区中的序列上,氨基酸编码的转移则相当集中,因而,隐 Markov模型能很好地刻画与鉴别有 无删除,错读与插入 （ 增读 ）,这样就 把问题纳入了隐 Markov 链的框架,再利用统计性质,就可以 纠正在使用 寻找基因的软件 Genescan 时可能导致的错误,在学习相位 ( 即学习问题 ),由给定的观测过程的一组样本找出隐 Markov 模型的参数族,可以用 Baum - Welch 算法,而在运转相位 ( 即解码问题 ),对给定的观测样本找出相应的 Markov 模型的状态,可以用 Viterbi 算法,( 相仿的思路也出自借鉴 Gibbs 场的思想,即在学习时先设定初值,再轮番修改参数和状态 ( 见第 15 章中的 EM 算法 ),而在运转时用 Bayes 决策找出 最佳 状态 ),
在 DNA 序列中,启动子是一个特殊位点,它 反映了在 DNA 转录过程中发挥定位功能的
DNA 片断表现在序列上的规律性,所以,启动子识别是生物信息论中首先要解决的问题,隐
Markov 模型也可以用来识别启动子 (Promotor ) 等特殊位点及编码区,( 这些特殊位点还包括增强子
( enhancer ),抑制子 ( silencer )),转录起始位置 ( TSS,transcription start ing site),外显子 (exon),内含子 (i ntron ),受子 ( accepter ),给子 ( donor ) 等等 ),
利用已有的生物学知识及数据库,构造合适的数学模型,例如,合适的隐 Markov 模型,再与数据进行反复的交叉验证与修改,是得到成功的识别方法的关键,
习题 10
１,证明 在均值为 m,方差为 2s 的分布密度中,正态分布的熵最大 ； 而 在均值向量为 m,方差矩阵为
Σ 的多维密度中,Gauss 分布 ),( ΣmN 的熵最大,请读者自己写出最大熵的表达式,
２,证明在有限区间 ],[ ba 上取值的分布密度中,均匀分布的熵最大,这个熵为 )log( ab?,再把这个
结论推广到多维情形,
３,设 )(,),(1 xgxg mL 都是非负有界函数,那么,在满足条件
)(,)()( midxxpxg ii ≤=∫ a
的分布密度 )(xp 中,当且只当
))()((
0
11)( xgxg mmCexp ll ++?= L
时熵为最大,其中
∫ ++?= dxeC xgxg mm )()(( 11
1
ll L,
而 mll,,1 L 是由
dxexg xgxgi mm ))()(( 11)( ll ++?∫ L )(,))()(( 11 midxe xgxgi mm ≤= ++?∫ lla L
确定的正数,
４,能不能把 3 推广到 )(,),(1 xgxg mL 不全为非负,也不全为有界的情形?
龚光鲁,钱敏平著 应用随机过程教程 – 与在算法和智能计算中的应用
清华大学出版社,2003
第 10 章 隐马氏模型 （ Hidden Markov Model,HMM ） 及其应用
１,熵与相对熵
１,1 离散分布的熵与相对熵
熵的概念出自 C.Shannon,引进这个指标的 目的在于刻 画 一个离散分布 ( 一个离散随机变量 ) 或一个分布密度 ( 一个连续型随机变量 ) 的不确定性的大小,也就是说知道了此随机变量的取值所获得的信息的大小,
定义 １０,１ 对于离散分布 p ),,,( 1 LL npp=,我们定义它的 熵 为
H( p )= ∑?
i
ii pp ln,
又定义分布 p 关于分布 q ),,,( 1 LL nqq= 的 Kullback - Leibler 相对熵 为
h （ p,q ） = ∑
i i
i
i q
pp ln,
命题 １０,２ 相对熵 h （ p,q ） ≥ 0,而且 h （ p,q ） =0,当且仅当 p = q 时成立等号,
证明 ),0[ ∞ 上 函数 tttg ln1)(=
D
在 1≠t 时恒正 ( 这 一结论可 由 g 的导数直接可以看出 ),且 0)1( =g,取
i
i
p
qt =,于是
i
i
p
qln ≤
i
i
p
q - 1,即
i
i
q
pln
i
i
p
q?≥ 1,而且等号仅当
1=
i
i
q
p 时成立,从而有
∑
i i
i
i q
pp ln ∑ =?≥
i i
i
i p
qp 0)1(,
所以结论成立,
这个命题表明,相对熵在相当程度上表达了 p 与 q 的差别,当 p = q 时,h （ p,q ） =0,而当所有的 ip 都与 iq 接近时,h （ p,q ） 就很小,从而 h （ p,q ） 可以看成 p 与 q 之间的一种
,准距离,,这里我们之所以称它为准距离,是因为它既不对称 ( 即 h （ p,q ） ≠ h （ q,p ） ),
也不满足三角形不等式 。 所以不满足 第 9 章 中的距离公理,
例 １０,３ ( 有限个值的分布 的熵 )
分布 p ),,( 1 Npp L= 的 熵 满足
H( p )= Npp
i
ii lnln ≤? ∑,
且等号当且仅当 分布 在 N 个 值均匀时成立,即以相同概率取 N 个值的分布 ( 称为离散均匀分布 ) 的熵最大,
证明 记分布 n = )1,,1( NN L,于是 本结论是 相对熵 不等式 h （ p,n ） ≥ 0 的变形,
例 １０,４ ( 数学期望固定条件下的离散 的 最大熵分布 )
假定 存在实数 a,使 )1(,≥≥ ixi a,对于固定的 ),,,( 1 LL nxx 与 m,在满足
mxx === EpxP ii,)(
266
的离散分布 p ),,,( 1 LL npp= 中,具有最大熵的分布为 p *,且
∑ ==
i
xx
i
ii eCCep ))((,1* 00 ll,
其中 0l 满足 条件
∑ ∞<?
i
xie 0l,∑ ∞<?
i
x
i
iex 0l,ii
x
i i
i
x exe 00 llm∑ ∑=,
这个最大熵 为 Cln?０ml,
证明 对于任意 0≠l,只要 ∑ ∞<?
i
xie l,相对熵不等式
≤0 h （ p,p * ） ∑?=
i
x
i
i iCe
pp
lln
∑∑?+=
i
ii
i
ii Cxppp lnln l
就是
H( p )= Cpp
i
ii lnln?≤? ∑ lm,
所以
H( p * ) ≥?=?= ∑ CCeCe
i
xi ln]ln[
0
00 mlll H( p ),
例 １０,５ （ 数学期望和方差都固定 时 的最大熵 ）
在 均值为 μ,方差为 σ 2,且 取值 为 kx ),2,1( L=k 的 离散分布中,分布 p *,
1
)()(
* )(,2
2
2
2
∑==
k
xx
k
kk
eCCep b
a
b
a
的熵最大,此 最大熵 为 Cln)( 2
22
+ b ams,其中 ba,满足
)2,1,0(,)( 2
2)(
=∞<∑
lex
k
x
l
k
k
b
a
,∑
k
xk
e 2
2)(
b
a
m,2
2)(
∑
=
k
x
k
k
ex b
a
=∑
2
)(
2 )( 2
2
k
xk
e b
a
s ][ 2
2)(
∑
k
xk
e b
a
][ 2
2)(
2∑
k
x
k
k
ex b
a
∑
k
x
k
k
ex 2
)(
)( 2
2
b
a
,
( 证明类似 ),
１,２ 分布密度 的熵与相对熵
定义 １０,６ 对于概率分布密度,我们可以仿效前面的思想,把 分布密度 )(xp 的熵 ( 其实是微分熵 ) 定义为
H ∫?= dxxpxpp )(ln)()(,
而把 分布密度 p(x) 对 分布密度 q(x) 的相对熵 定义为
∫= dxxq xpxpqph )( )(ln)(),(,
这个定义可以推广到多维密度 ),,( 1 dxxp L 与 ),,( 1 dxxq L,
267
H ∫= d
d
d
d dxdxxxq
xxpxxpp L
L
LL
1
1
1
1 ),,(
),,(ln),,()( ),
∫= d
d
d
d dxdxxxq
xxpxxpqph L
L
LL
1
1
1
1 ),,(
),,(ln),,(),(,
同样 也有 相对熵是非负的,从而由熵的定义,再利用相对熵非负性质 可 推得下述结论,
例 １０,７ ( 半直线上均值已知时的最大熵分布密度 )
在 ），［ ∞0 上具有同均值 m 的密度中,指数分布
m
1Exp 的熵最大,这个 最大 熵为 mlog1+,
例 １０,８ ( 均值与方差已知的 分布密度的 最大熵 )
在均值 为 m,方差为 2s 的 分布 密度中,正态分布的熵最大,
例 １０,９ ( 均值向量与 ( 协 ) 方差矩阵已知的有 联合密度的 最大熵 )
在均值向量为 m,方差矩阵为 Σ 的多维密度中,Gauss 分布 ),( ΣmN 的熵最大,请读者自己写出这个最大熵的表达式,
例 １０,１０ ( 有限区间上的最大上密度 )
在有限区间 ],[ ba 上 取值 的 分布 密度中,均匀分布的熵最大,此最大 熵为 )ln( ab?,
例 １０,１０ 也可以推广到多维情形,
例 １０,１１ ( 在 取值于 混合概率向量 的分布密度 中,关于某个 概率向量 的平均相对熵给定条件下,求 具最大相对 熵 的密度 )
令样本空间为
{=? x = }1),(,10:),,( 11 =++≤<< did xxdixxx LL,
中的 点可以解释为概率向量,而取值于? 的连续型随机变量的分布密度可以解释为,广义的混合分布,,给定 W∈= ),,( 1 dzzz L,把? 中的点看成概率向量,就有相对熵 ),( xzh,
它 是 取值于? 的 x函数,对于一个取值? 的分布密度 xp( ),,() dx xxp L
= 定义对于这个分布密度的平均相对熵,
(m z xdxpxzh )(),() ∫=
D
,
我们要在平均相对熵 )( zm 相等的 分布密度 )(xp 中,选取 一个 )(* xp,使 其 熵
∫? xdxpxp )(ln)( ** 达到 最大,
为此我们 断言,如下的 Gibbs 型分布密度
)()(* zCxp = ),()( xzhze l? )(( zC= ))()(1
ln)(
1 di
ii zz
d
zz
zzz
xxe ll
l
L∑
,
其中
1),()( )()(∫= xdezC xzhzl,
)(zl 满足
xdexzhzCz xzhz ),()(),()()( lm?∫=,
就是我们要的具有最大熵的密度,证明如下,对于任意一个满足 (m z xdxpxzh )(),() ∫=D
的分布密度 )(xp,由 相对熵不等式
268
),(0 *pph≤ xd
ezC
xpxp
xzhz∫?
=
),()()(
)(ln)(
l
= H( )(ln)()() zCzzp?+ ml
得到
H( )p )(ln)()( zCzz?≤ ml,
于是 我们有
H( *p ) = xdezCezC
xzhzxzhz ])(ln[) ),()(),()( ll∫? （
)(ln)()( zCzz?ml＝ ≥ H( p ),
[ 注 1] 上面我们用相对熵 ),( xzh 作为准距离,而不是用更为合理的 [21 ),( xzh )],( zxh+,
也可以用 与它相应的 ),( zxh,只 要 易于计算,不管用 那一个 都行,
[ 注 ２ ] Diriclet 分布的密度函数为
)1{
11
11 1
1),,(
=++
=
d
d
xxdd IxCxxxf LLL
aa,
其中
∫
∞
=Γ
++Γ
ΓΓ=
0
1
1
1 )(,
)(
)()( dxexC x
d
d aa
aa
aa
L
L,
例 １０,１１ 中的最大熵分布恰是 Dirichlet 分布,在生物信息论中,常常需要估计概率向量组成的空间上的分布密度,例 １０,１１ 说明,在平均相对熵给定的条件下,Dirichlet 分布是 ＂ 最吃亏的 ＂ 分布,用它作为先验分布 ＂ 看起来更为保险 ＂,这就解释了在生物信息论中,人们常常喜欢 用 Dirichlet 分布作为先验分布的原因,
相对熵应用的一个实例 － 在各个测试特征 （ 统计量 ） 中选择几个最为有效的相对熵方法
（ 特征量选取 的相对熵方法 ）
假定我们有两组性质完全不同的群体,例如一组是健康人,另一组是 SARS（ 非典型性肺炎 ） 病毒持带人,又假定人们已经提出了 N 种区别健康人与 SARS 病毒持带人的不同特征,我们要在其中选取区别效果最好的 M 个特征,相对熵方法是较为有效的一种方法,其实际操作为,
用一个给定的区 分特征,对上面的一组健康人测定了一组数据,简称为甲数据组 ； 又对上面的
SARS (Severe Acute Respiratory Syndrome) 病毒持带人测定了一组数据,简称为乙数据组,再进行如下的步骤,
（ １ ） 分别找出此两个数据组的近似分布,从数据组出发,应用统计中的核估计的思想,分别得到其近似分布密度曲线 ；
（ ２ ） 计算此两个分布密度的相对熵 0h ；
（ 实践证明,用数值计算求此两个分布密度的相对熵,对于计算格点大小的划分并不太敏感,相反地,如果直接将两 个直方图作为离散分布求相对熵,则对于直方图的计算格点大小的划分十分敏感,
以致得到的计算结果很不稳定 ）,
（ ３ ） 将甲乙两组数据合并,再随机地重组为和以前个数相同的两组 (随机地重排 （ Random
Sorting） 后,再按原来的各组的个数顺序分成两组 ),用 （ １ ），（ ２ ） 步骤计算其相对熵 ；
（ ４ ） 重复地作 （ ３ ） 多次 （ 例如 １ 万次 ）,计算其中相对熵大于 0h 的次数所占的频率,记为 p,
269
于是对应于一个区分特征就对应地有一个 p 值,它是一个客观的标准,和各种测试方法的结果
（ 统计量 ） 的绝对大小,量纲都无关,如果要在 N 个不同的用以区分是否感染了 SARS 的特征中选取 M 个相对地最为有效的特征,我们只需选取对应的 p 值最小的前 M 个特征即可,
２,隐 Markov 模型 (HMM)
２,１ 一个实例
先举一个例子以 直观地 理解这个 模型 的实质,
例 １０,１２ 设 某人 在 三个 装有红白两种颜色的球的 盒子 中,任取一个盒子,然后在此盒子中每次抽取一个球,连续地抽取 m 次,假定各个盒子的内容与抽取 方式 分别为,
红球数 白球数 每次抽取方式
盒 1 90 10 随机取一球,记下颜色后不放回,而放进一个与它颜色不同的球,
盒 2 50 50 随机取一球,记下颜色后放回
盒 3 40 60 随机取一球,记下颜色后不放回,并放进一个红球,
现在如果 某人 用上述方法得到了一个记录 (红,红,红,红,白 )（ 即 5=m ）,但是不告诉我们 球出自 哪个盒 子,我们应如何推测他是从哪 个盒子中 抽取的观测样本呢?
令
=)(knS 在第 k 个盒子 )3,2,1( =k 中第 n 次抽取完成后在各盒中的红球数,
那么,在 k 分别 固定 为 １,２,３ 时,}0:{ )( ≥nS kn 分别为 Markov 链,且 其 概率 转移 矩 阵分别为,
+=?
=
=
)(0
)1(1001
)1(100
)1(
其它情形
iji
iji
pij,（ 逢红,红减 １ ； 逢白,红加 １ ）,
≠
==
)(0
)(1)2(
ij
ijp
ij,（ 内容总是不变 ）
+=?
=
=
)(0
)1(1001
)(100
)3(
其它情形
iji
iji
pij,（ 逢红不变 ； 逢白,红加 １ ）
而且初值 分别为,40,50,90 )3(0)2(0)1(0 === SSS,于是这 ３ 个盒子就分别对应于 ３ 个不同的
Markov链 模型,把这 ３ 个模型分别简记为 321,,lll,并 把 某人 观测到的样本序列中的第 n个记为 nO,即令
=nO 抽到的记录列中第 n个记录中的白球数 （ 只能 为 1 或 0）,
此例是一个玩具模型,从这个例子可以看出,在观测出自哪个盒子已知时,状态随机变
270
量序列 }{ nS 与某人提供的观测随机变量序列 }{ nO 之间的条件概率计算的关系可以直观地 写为,
mmmim SOSOSOS,,,,,,1110L,
其中在前面的一段随机变 量序列取定值的条件下,继后的那个随机变量取值的条件概率就完全确定了,而且 在这 ３ 个模型下分别都有
)|()|(),,,|( 121111 SSPSSPSOSOSP nnnnn == ++ L,
)|()|(),,,,,|( 1211111 SOPSOPSOSOSOOP nnmmmimn == ++ L,
于是
),,,,,,( 1110 mmmim SOSOSOSPL
)|()|()|()|()|()( 111201010= mmmm SSPSOPSOPSSPSOPSP L,
具体地,我们有
在 模型 1l 下 （ 把取到的球换色 ）
)|)1,0,0,0,0(),,,,(( 154321 l=OOOOOP =0.9× 0.89× 0.88× 0.87× 0.13≈ 0.08,
在 模型 2l 下 （ 球的内容不变 ）,nO 是独立同分布的随机变量序列,且
2
1
2
1
10
~nO,
所以
)|)1,0,0,0,0(),,,,(( 254321 l=OOOOOP 5)21(= ≈ 0.03,
在 模型 3l 下 （ 取到红则不变,取到 白则换为红球 ）
)|)1,0,0,0,0(),,,,(( 354321 l=OOOOOP ≈×= 6.0)4.0( 4 0.015,
再 用 Bayes 公 式 分别 得到 （ 即取上面 3 个概率的归一化值 ）
== ))1,0,0,0,0(),,,,(( 54321 OOOOOP ｜１l 0.64，
== ))1,0,0,0,0(),,,,(( 54321 OOOOOP ｜２l 0.24
== ))1,0,0,0,0(),,,,(|( 543213 OOOOOP l 0.12,
可见从第一盒抽出样本 (红,红,红,红,白 )的概率要比从其他两盒中抽出该样本的概率要大得多,
这个模型的意义绝对不在于游戏,从中可以 抽象出一种能广泛应用的数学模型 － 隐 Markov
模型 （ Hidden Markov Model,简记为 HMM ）,在这个例子中的不同盒的抽取方式可以抽象为不同的编码 方式,而这个例子 正 启示了用隐 Markov模型作模式识别的一种粗略梗概,
２,２ 隐 Markov 模型的描述
定义 １０,１３ 设 }1:{ ≥nX n 是取值于有限状态 },,1{ LL 的随机变量列,称为 状态链,1X 的分布称为其初始分布,记为 m,假定 nX 是时齐的 Markov 链,记
)|( 1 iXjXPa nnij === +,LjiijaA ≤=,)(,(10,1)
为它的转移矩阵,假定 nX 的取值,初始分布 m 与转移矩阵 A,都不能测量得到,而能测量到的是另一个与它有联系的,且可以观测到的一个取值于与有限集 },,{ 1 Mnn L 的随机变量序列 nY,称为 观测链,它们合起来还要求 满 足 如下的 隐 Markov 条件 （ HMM 条件 ）,
271
NNNNNN yiiiyiiiyii
bababxXyYP
11121111
),(
=== Lm,(10,2)
其中 N 为样本观测的时间长度,而
),,( 1 NXXX L=,),,( 1 NYYY L=,),,( 1 Nyyy L=,),,( 1 Niix L=,
NnLiy nMn ≤≤≤≤∈ 1,1},,{ 1 nn L,初始概率向量 （ 初始分布 ） 为 ),,( 1 Lmmm L=,由未知状态链与测量到的观测链一起 ),( nn YX,就构成了 隐 Markov 模型,这里,隐,的含义是说状态链是隐藏起来的,
隐 Markov 模型的基本 假定 是,
参数向量 m,参数矩阵 A 与 MLilbB ×
= )( （ },,{},,,1{ 1 MlLi nn LL ∈∈ ） 都是未知的,我们 将 它们 合 记为 参数组 ),,( BAml
=,参数组 ),,( BAml
= 完全确定了状态链与观测连的联合统计规律,所以,我们通常用一个 l 表示一个隐 Markov 模型,并 称 之为 隐 Markov 模型
l ( 更确切地为 隐 Markov 链 ),在例
10.12 中,３ 个不同的模型就对应了 ３ 个不同的参数组,只要令 1,+== nnnn OYSX,它们
满足 HMM 条件,因而纳入了隐 Markov 模型的框架,
(10.2) 是 ),( YX 的 联合分布 通过 参数表达 的形式,它是 计算各种边缘概率与条件概率 的出发点,而 HMM 条件的含义是,状态链与观测链的联合分布是由一系列简单转移与条件概率的乘积表达的,
２,３ 隐 Markov 模型的等价表述
由例 １０,１２ 的启示,可以想到下面的等价条件
命题 １０,１４
HMM 条件等价于,对于任意的 },,1( Li L∈ 以及 },,{ 1 Ml nn L∈,有
,)|(),,,,,|(
11 111111 ilnlnllnnnnln
biXvYPXiYXiYiXvYP
n
Dnn ==========
L
(10,3)
)|(),,,,,|( 1111111
1
iXjXPYiXiXYiXjXP nnlnnlnnn
n
========= ++ nn L ija?=,
(10,4)
这两个等式的证明只需利用条件概率的定义,所以我们把它留给读者,它们的直观含义就是,
nY 与 1+nX 相对于历史条件的 统计规律只与时间上最接近的 nX 有关,而与其它更早的历史无关,在实际问题中,这是比较容易判断的,
２,４ 非线性滤波作为 隐 Markov 模型 的特例
设 nX 是取值于状态空间 },,2,1{ LS L= 的 Markov 链,其 样本不能被实际测量 得到,而能测量到的是 如下的 nY 的样本
nnn wXgY += )(,(10,5)
其中 g 是一个未知函数,}{ nw 是独立同分布的 随机 干扰,只取有限个值,且 与随机过程
}{ nX 独立,那么 },{ nn YX 就是一个 （ 二维的 ） 隐 Markov 模型,这就把 （ １０,５ ） 的非线性滤波放进了 HMM 的框架,
２,５ 在应用中研究 隐 Markov 模型 的 主要 方面
(1) 从一段观测序列 },{ mkYk ≤ 及 已知 模型 ),,( BAml = 出发,估计 nX 的最佳值,
称为 解码 问题,这是状态估计问题,
272
(2) 从一段观测序列 },{ mkYk ≤ 出发,估计模型 参数组 ),,( BAml =,称为 学习问题,就 是参数估计问题,
(3) 对于一个特定的观测链 },{ mkYk ≤,已知它 可能 是由已经学习好的若干模型之一所得的观测,要决定 此观测 究竟 是得自其中哪一个模型,这 称为 识别问题,就 是分类问题,
在实际问题中,这 3 个问题 有时并不完全 能分开,有时也并不 需要 解全 3 个问题,例如,
在语音识别或手写体汉字 或数字的 脱机识别中,我们只需要作 (2) 与 (3),这相当于将一个,标准的,语音 音素 或一个,标准的,手 写 体 汉字,学习成,一个隐 Markov 模型,即把它与 一个或几个 特定的隐 Markov 模型 ( 更确切地说,是特定的一组参数 ) 相 对应起来,以便 把该模型作为这个音 素 或 手写 字的代表 模板,这是 学习相位 ; 而在进一步用这些模板中的最合适者,作为对于一 个需要识别的 音 素 或 手写 字 的分类归属,这是 运转相位,至于 归入那个模板最合适,就要用合理的距离,或准距离 （ 常用的是相对熵 ）,在此意义下优化,
[ 注 1] nX 也可以推广为,取值于平面有限格点的 Markov 场,而 nY 与它的关系仍如上述,
这就定义了一个 隐 Markov 场,
[ 注 2] 一般地,nY 还可以是连续型随机变量,如果 记 在 xX n = 的条件下,nY 的条件分布密度为 f,则 ),,( xfAml = 就也称为一个 （ 连续的 ） 隐 Markov 模型,这时 f 常是分布类型已知而带有未知参数的密度,状态是连续的 隐 Markov 模型至今还 未见 有人使用,
３ 解码问题 -- 已知模型 l 与观测 yY = 时状态 X 的估计
３,１ 出现当前的观测的概率 )|( lyYP = 的计算
我们仍旧 沿用 记号
),,( 1 NXXX L=,),,( 1 NYYY L=,),,( 1 Nyyy L=,),,( 1 Niix L=,
NnLiy nMn ≤≤≤≤∈ 1,1},,{ 1 nn L,初始概率向量 ),,( 1 Nmmm L=,
由 (10,2),利用条件概率的性质 容易算出
NN yiyiN
bbiiXyYP LL
11
)),,,(|( 1 === l,
NNNNNN yiiiyyiiyii
bababxXyYP
1121111
)|,(
=== Lml,
NNNN
N
yiiiiiyii
ii
baabyYP
121111
1,,
)|(
∑
== L
L
ml,(10,6)
对于 Nn ≤≤1 及观测样值 yY =,记 （ 因为观测样值 y 是固定的,所以下面我们将在足标中把它略去 ）
)|,,,()( 11 la iXyYyYPi nnnn ==== L ( 依赖 y ),(10,7)
则 在模型 l 给定下,关于 观测资料 ),,( 1 nyy L 的长度 n )( Nn <,我们有递推公式 ( 称为 向前递推公式 或 向前算法 )
1
)()(1
+∑
=+
niyj jinn
bji aaa,(10,8)
由此得到
(1) 基于 )(ina 的向前递推公式 计算 )|( lyYP = 的 步骤,
计算初值
1
)(1 iyibi ma =,
273
用递推公式
1
)()(1
+∑
=+
niyj jinn
bji aaa,(10,9 )
最后 得到结论 )()|( iyYP N
i
al ∑==,(10,10 )
计算 )|( lyYP = 也可以通过另一种途径,为此 记 （ 观测样值 y 也是是固定的,我们也把它在足标中略去 ）
),|,,(()( 11 lb iXyYyYPi nNNnnn ==== ++ L,(10,11)
则在模型 l 给定下,关于观测资料 ),,( 1 nyy L 的长度 n )( Nn <,我们 还有递推公式 ( 称为 向后递推公式 或 向后 算法 )
1
)()( 1
+∑ +
=
njyj ijnn
bji abb
这就得到计算 )|( lyYP = 的另一个算法,即
(2) 基于 )(inb 的向后递推公式 计算 )|( lyYP = 的步骤
定义 1)( =iNb,
用递推公式
1
)()( 1
+∑ +
=
njyj ijnn
bji abb ( 10,12)
最后 得到结论 ∑==
i
iyibiyYP 1)()|( 0 mbl (10,13)
３,２ 解码问题 -- 已知模型 l 与观测 yY = 时状态 X 的估计
令
),|()( lg yYiXPi nn ===,(10,14)
那么
∑∑
=
==
===
i
nn
nn
n
i
n
n
ii
ii
iXyYP
iXyYPi
)()(
)()(
)|,(
)|,()(
ba
ba
l
lg,(10,15)
( 推导的思路 如下,令 F n ),,( nmYX mm ≤= s ( 右 方 表示括号内随机变量 所 生成的 s 代数 ),
},,{ 11 nnn yYyYY ===? L,},,{ 11 NNnnn yYyYY === +++ L,因为 l 是固定的,我们在下面记号中 也把它 （ 为了记号的方便,我们在下面 有时并 不区分取概率的符号与取期望的符号,取概率 也 理解为 取期望 ）,于是 我们有
)|,( liXyYP n == |,([( iXyYPE n === F n )],|(,,[ iXYPiXYE nnnn === +? F n )]
)]|(,,[ iXYPiXYE nnnn === +? )|(),( iXYPiXYP nnnn === +? )()( ii nn ba= ),
对于解码问题的解决,我们用如下两种考虑 方法,
1,对 单个 时刻的 状态的最大概率估计
如果 *i 满足 )(max)( * ii nNin gg ≤=,则取 *
^
iX n =,
单个 时刻的 状态的最大概率估计的长处在于简单易算,但是它忽视了状态链在不同 时刻之 间的联系,也 就 不能充分地利用已知的信息,下面的 路径估计法,是在考虑不同时间的联系后的改进 方法,
274
2 路径的 ( 轨道 ) 最大概率估计 -- Viterbi 算法
这是基于动态规划的思想的一 种整体性考虑的方法 （ 所谓动态规划实质上 就 是一种从后向前 逐步取优的优化方法 ）,令在时刻 n 取到状态 i 的 各个 路径 的最大 概率为
)|,,,,,,(max)( 111111,,
11
ld yYyYiXiXiXPi nnnnniin
n
======
LLL,(10,16)
那么我们有
))((max)(
11 jinjiyn
ajbi
n
dd
+
=+,(10,17)
用来 计算 它的实际 程式 是 如下的 Viterbi 算法,假定最佳路径为 ),,(
^
1
^
NXX L,今设置一
个从 1
^
+nX 得到 nX
^
的对应 )(,1
^
1
^
1 +++ = nnnn XX yy,计算 )(ind 与 ny 的递推算法 的 具体 设计 如下,
设置初值 1)(1 iyi bi md = ;
用递推公式 ))((max)(
11 jinjiyn
ajbi
n
dd
+
=+ ;
再 定义 *1 )( jin =+y,其中 *j 只要 满足
ijniyn ajbi n *1 )()(
*
1 dd +=+,(10,18)
最终采取如下的向后递推算法 ( 相当于动态规划,动态规划的核心思想是,在一切最佳
路径上的任意节点出发,其后也应该是最佳路径,因此就能 从后 向 前 地逐步求得最佳路径 )
得到 如下的,最佳,轨道估计 ),,(
^
1
^
NXX L,
计算终值 )}(max)(:{*^ jiiiX NjNN dd =∈=,(10,19)
反向归纳地取 )(
^
11
^
++= nnn XX y,(10.20)
４ 学习 问题 – 由观测 yY = 估计模型参数 l
４,1 状态链样本已知时的参数的频率估计
设 隐 Markov 模型的 状态链 有 充分长 的样本,把 其中从状态 i 到 下一个时刻为 状态 j 转移的频数 记 为 ijA,那么,粗略地 可取
∑
=
= N
j
ij
ij
ij
A
Aa
1
^
为 ija 的估计,同样 设此观测链样本充分长,记其中 从状态 i 到同一时刻的观测 l 转移的频数为 ijB,则可取
∑
=
= M
l
il
ilil
B
Bb
1
^
为 ilb 的估计,不幸的是状态链是经过估计得到的,所以不修正地使用频 率估计会增加误差,再则,这种估计 也 过于敏感,
４,2 模型参数估计的 EM 算法的思想
我们取 l 的最大似然估计
^
l,即它 满足
)|(max)|(
^
ll l yYPyYP === (10.21)
275
但是,对 于似然函数
NNNN
N
yiiiiiyii
ii
baabyYP
121111
1,,
)|(
∑
== L
L
ml 而言,这个最大值问题 的 计算量过大,在实际中是不可能被采用的,为此,人们灵活地 利用 使似然函数尽量地大的思想,
采取折衷的方案,构造一个递推算法,使之 能相当合理地给出模型参数 l 的粗略估计,其核心思想是,并不要求备选估计 l 使 )|( lyYP = （ 关于 l ） 达 到最大或局部极大,而只要求使
)|( lyYP = 相当大,从而 计算变为实际可能,
为此,我们 定义一个 描述模型,趋势,的量 )|( * llQ 代替似然函数 )|( lyYP =,其定义 为
)|,(ln)|,()|( ** llll yYxXPyYxXPQ
x
=====∑?
)|,(ln)|,( *∑=
x
yxPyxP ll
简记为
(10,22)
( 这个量并不是相对熵,但是在形式上很类似 ),利用在 10 << x 时的不等式 1ln?≤ xx
( 证明,令 xxxf ln)1()(=,则 0)1( =f 且 01)(" 2 >= xxf,即 )( xf 为凸函数,从而 0)( ≥xf ) 得到
∑=?
x yxP
yxPyxPQQ
)|,(
)|,(ln)|,()|()|( **
l
llllll
)|()|()1
)|,(
)|,()(|,( ** ll
l
ll yYPyYP
yxP
yxPyxP
x
=?==?≤ ∑,
由此 可见对于固定的 l,只要 )|()|( * llll QQ >,就有 )|()|( * ll yYPyYP =>=,
于是 想 把模型 nl 修改为更好的模型 1+nl （ 即找 1+nl 使 )|()|( 1 nn yYPyYP ll =>= + ）,
只需找 1+nl 使
lll sup)|( 1 =+ nnQ )|( nQ ll,(10,23)
即只要把 )|( nQ ll 关于 l 的最大 值 处取成 1+nl,就 有 )|()|( 1 nn yYPyYP ll =>= +,
这样得到的模型序列 }{ nl 能 保了 )|( nyYP l= 关于 n 是 严格递增 的,虽然 在这里 还不能在理论上证明 )|( nyYP l= 收敛到 )|(max l
l
yYP =,但是当 n充分大时,nl 也还能提供在实际中较为满意的粗略近似,
综上论述,我们 把如 上得到近似模型列 nl 的方法 归结为两个步骤,
1 E 步骤 （ 求期望 ）,计算
=)|( * llQ )|,(ln)|,( *∑
x
yxPyxP ll ；
2 M 步骤 （ 求最大 ）,求 1+nl 使
lll sup)|( 1 =+ nnQ )|( nQ ll 。
276
这两个步骤合起来构成的算法,称为 EM 算法,EM 算法是针对 在 测量数据不完全时,求 参数的一种近似于最大似然估计的统计方法,隐 Markov 模型 的 建 模 （ 模型 参数 的 估计 ）,是 EM
算法的一个最常见,且 极有用的 一种 典型例子,下面将看到,隐 Markov 模型 中的 M – 步骤的解可以有显式表示,这就 是 一组把模型参数修改为新的模型参数的递推公式,这组公式 正好是在 隐 Markov 模型 中普遍应用的著名的 Welch - Baum 公式,下面我们给出 得到这组递推公式的 梗概,首先我们 有下面的引理,
引理 １０,１５ 设 )(0 Nizi ≤≥,则在约束条件 ∑
=
=
N
i
ix
1
1,)(0 Nixi ≤≥
下,函数 ∑
=
N
i
ii xz
1
ln 在
∑
=
= N
i
i
i
i
z
zx
1
)( Ni ≤ 处取得最大,
证明 注意 ),,(1 1
1
NN
j
j
zz
z
L
∑
=
对 ),,( 1 Nxx L 的相对熵 0
)(
ln)( 1
1
1
≥
∑
∑ ∑ =
=
=
i
K
j
j
i
K
i
K
j
j
i
x
z
z
z
z,
４,3 隐 Markov 模型 中 M – 步骤的 求 解
现在我们用引理 10,1,来求 )|( nQ ll 关于 l 的最大值的位置 1+nl,注意
NNNN yiiiiiyii
baabxXyYPxyP
121111
)|,()|,(
====? Lmll,
∑ ∑
=
=
+
++=
N
n
N
n
iiyii nnnn abxyP
1
1
1
11
lnlnln)|,(ln ml,
由此得
=)|( * llQ )|,( lxXyYP
x
==∑ )lnln(ln
1
1
1
***
11 ∑ ∑
=
=
+
++
N
n
N
n
iiyii nnnn abm
)(),(ln)|,(ln }{
1
*
1
1 nln
Nn
il
M
lii
i yIiXyYPbiXyYP ==+=== ∑∑∑∑
≤≤=
lm
∑∑?
=
+ ===+
ji
N
n
nnij jXiXyYPa
,
1
1
1
* )|,,(ln l ∑ ∑++=?
i i
iaib aQbQQ
ii
)|()|()|( *** lllmm,
此处 带 * 的参数都表示修改 后 的参数,而
**,
ii ba 分别表示待修改的参数 BA,的第 i列列向量,显见 ∑ =
i
i,1
*m ∑ =
j
ija,1
* ∑ =
l
ilb 1
*,这样可以通过分别最大化上面等式右边各项得到 将模型 ll =n 更新 为
*
1 ll =+n 的修正 模型 1+nl,而其中每一项都可以用引理 １０ 。 １５
最大化,这样从最大化 就得到 如下的结论,
Baum － Welch 算法 （ 递推关系 ）,
277
)(
)|(
)|,(
1
1
^
)1( i
yYP
iXyYP
n
nn
i gl
lm =
=
===+,(10,24)
∑
∑
∑
∑
=
=
=
=
+
+ =
==
===
= 1
1
1
1
1
1
1
1
1^
)1(
)(
),(
),|(
),|,(
N
n
n
N
n
n
N
n
nn
N
n
nnn
n
ij
i
ji
yYiXP
yYjXiXP
a
g
V
l
l 记为
,(10,25)
∑
∑
∑
∑
=
==
=
=+ =
==
==
= N
n
n
N
lyn
n
N
n
nn
N
n
nlnn
n
il
i
i
iXyYP
yIiXyYP
b n
1
,1
1
1
}{^
)1(
)(
)(
)|,(
)()|,(
g
g
l
l 记为
,(10,26)
其中
∑
=
i
nn
nn
n
ii
iii
)()(
)()()(
ba
bag 如 (10,14) 与 (10,15),而
∑
+
+
+=====
i
nn
njyijn
nnn
ii
jbaiyYjXiXPji
n
)()(
)()(),|,(),( 1
1
1
ba
balV
,(10,27)
（ (10.27) 的推导需要先算 )|,,( 1 ljXiXyYP nn === +,它与计算 )|,( liXyYP n == 的方法几乎完全一样 ),此时还有
∑=
j
nn jii ),()( Vg,(10,28)
在上面所述的算法中,初始 值 0l 的设置会直接影响到估计的好坏,为此常用的 一种方法是,根据先验知识设置一条较长的,标准,的虚拟状态链 },,{ )0()0(1 NXX L,用 ４,０ 段 中的频率估计粗估一个 l,并 把它取为 0l,
对于隐 Markov 模型,甚至对于第 15 章中的一般 EM 模型,观测链提供了所有的已知信息,
但是,与一般的 EM 方法相比,HMM 远为简单,因为 如果我们 仔细考察 B aum - Welch 算法,就可以发现它并不依赖于状态 链 的估计,这就是说,隐马氏模型中的参数估计与状态链的估计 之 间是 解 耦的,而一般 EM 算法 常常 都是 有 耦合的,即在估计参数时需要知道状态的估计,而在估计状态是又需要知道参数的估计,
在向前与向后两个递推公式基础上的 B aum - Welch 算法,给出了估计模型参数的合理方法,
虽然 在理论上 还没有 证明 其极限 就是最大似然估计,
[ 注 1] 可以对 HMM 条件作一点小的推广,假定 对 },,{},,,1{ 1 MlNi nn LL ∈∈ 有
),(),|(
),,,,,|(
1
111111 11
lbiXjXvYP
YiXYiXiXvYP
ijnnln
llnnnnln n
=====
======
nn L
(10,29)
)|(),,,,,|( 1111111
1
iXjXPYiXiXYiXjXP nnlnnlnnn
n
========= ++ nn L ija记为 =
(10,30)
与 n无关,在实际问题中常常会遇到这一种情形,
[ 注 2] 在 HMM 中,nn YX,都还可以是多维的 向量,这时的矩阵 BA,就 相当 复杂,
278
[ 注 3] 如果模型中的时间参数 n 相应地变为地点参数,那么就导致 二维的隐 随机 场 模型,特别对于给定的 －? 领域系统,就有隐 －? Markov 场模型,此时 B aum - Welch 算法就无法进行了,这时的 学习问题可以,例如用如下的 交错估计 方法 ( 这个方法是 耦合的,在这点上它与 B aum - Welch 算法不同,我们仍称这个算法为 推广的 V i t er b i 算法 ),
(1) 先在某个
)(nll =
固定下,用 推广的 V it er bi 算法得到优化的状态轨道 ;
(2) 再 求
^
l,使 ll sup),|(
^*
=== xXyYP ),|( lxXyYP ==,其中
),,( **1* nxxx L= 是随 n增加时轨道的最佳估计,也可以使用轨道频率,
(3) 令
^)1(
ll =+n,把
)(nl
更新为
)1( +nl;,
(4) 如此反复 迭代,
与之联系的有以下一些问题,
问题 1,用推广的 Vertabi 算法频率算法代替 Baum - Welch 算法后,解决了隐 Markov 场时,
Baum - Welsh 算法,在二维难于递推的困难,但是如何给出递推公式? 能否在某些条件下,得到递推更新中每次修正都往好的方向发展的结论,例如,使相对熵增加?
问 题 2,如果问题 1 中的估计结果较好,那么用一个随机取的 x 代替 ),,( **1* nxxx L=,
是否还能接受?
问题 3 隐 Markov模型中,如何引进 Bootstrap 式思维,
[ 注 4] 在生物信息论 中,还常用 隐半 Markov 模型,这时的做法可以是,
(1) 假定一个初始模型 (0)l ( 包括初分布,停留时间,转移参数 ),由观测资料通过
Viterbi 算法估计轨道 ;
(2) 由此轨道通过计算频率来更新初分布,停留时间,转移参数,得到更新模型
(1)
l ;
(3) 重复步骤 (1) 与 (2),直至较为稳定,
在语音识别中,也有人使用 隐半 Markov模型,但是与用隐 Markov模型相比较,在评估是否得到显著的改进这一点上,还颇多争议 。
５ 关于 隐 Markov模型 的评注
５,１ 隐 Markov模型 包容度大有非常宽的应用面
用隐 Markov模型 建模的优点
（ 1 ） 隐 Markov 模型 是一种简单的数学模型,它的包容性很大,内涵很广,在应用中的弹性相当大,可以很好地利用我们对于被建模的对象所了解的先验知识,因而有很广的适应 性,
（ 2 ） 隐 Markov 模型 是一种不完全数据的统计模型 （ 包含了非线性滤波模型 ）,它 之所以被广泛采用,在于这种模型既能反映对象的随机性,又能反映对象的潜在结构,便于利用 对象的结构与局部联系性质等方面的知识,以及 对研究对象的直观的与先验的了解,
（ 3 ） 隐 Markov 模型是 很少的几个既从物理模型出发,又与数据拟合相直接联系的算法之一,
（ 4 ） 隐 Markov 模型 虽然也具有某些黑箱的特点,但是与纯黑箱操作的人工神经网络算法相比,有明显的优点,它的参数常常有较为实质的含义,而线性模型,时间序列模型 中的参数一般只是作为被拟合了的参数而出现的,缺少 较为 真实含义,
（ 5 ） 隐 Markov 模型 有快速而有效的学习算法,
以上诸点使 隐 Markov 模型称为 随机建模与拟合的重要基石,在应用中有其普适性,
279
在具体应用 隐 Markov模型 建模时,首先要设定 Markov 链的状态集及其规模,即总状态数
L,这也有相当大的弹性,然后,确定相应的观测过程,为了使得观测链 ( 或场 ) 与马氏链 ( 或场 )
之间有隐 Markov 关系,我们常常需要对实际提供的观测链 ( 场 ) 加以改造,例如,一个脱机手写体汉字是一张黑白象素图 ( 一个 取值 0 或 1 矩阵 ),它和刻划笔道的 －? Markov 场之间的依赖关系并不能满足隐 Markov 关系的要求,需要对它进行 预 加工,也就要 先 进行,特征提取,,
以尽量简化得使隐 Markov 关系能够成立,
在确定了模型规模以后,对于隐 Markov模型的处理,又分为两个不同的相位,学习相位与运转相位,学习相位是由给定的观测过程的一组样本出发去估计隐 Markov 模型 的参数
),,( BAml =，,从而完全确定模型 ; 而运转相 位则是用学习相位所 确定的模型参数,对给定的某个观测过程的样本,来估计相应隐 Markov模型的状态,或计算出它在各个隐 Markov 模型下出现的概率,
语音识别和脱机手写体汉字识别最终需要知道的是给定的语音或脱机手写体汉字样本应出于哪个模型 (哪个 参数组 ),这是一个识别过程,通常我们可以按照 Bayes 方法 （ 后验概率 方法 ）
去决定模型,即选取使得给定样本出现的概率最大的模型作为它的分类归属,这就是识别的结果,而在 DNA碱基序列的排序的问题中,考虑的方式则与之不同,在 DNA碱基序列中,需要知道的是如何排序,也就是需要知道给定的某个观测过程的样本所对应的隐 Markov 链的状态,
在学习相位中,我们面临的是一个不完全数据的参数估计问题,除了 Baum –Welsh算法以外,还 可以进行交错估计,即可先设置一组初值参数,由此出发配合观测列按 Bayes 方法 去估计马氏链的状态列,再用此状态列配合观测列用最大似然估计去重新估计模型参数,交替地不断重复此二步骤,…,以达到比较稳定的最终结果,也可以从设置马氏链的一个状态列作为开始,去估计参数组,再由此参数组对应的模型去估计状态,如此交替进行,
在应用 隐 Markov 模型 于语音识别或手写体 汉字识别时,往往把一个音,一个词,一个字对应于一个隐 Markov模型 l,这时当然也可以把一个标准音或标准字的状态设置为马氏链的初值,再开始上述交替估计的过程,
值得指出的是,利用已知的先验知识,来得到简单的递推计算公式,是使得 隐 Markov 模型的计算能够 化简 的关键,在实际应用中,参数组的维数常常非常高,这就使得计算无法实际进行,
所以在交替估计的过程中,我们还可 以借助 Gibbs 采样的思想,一次只变化一个状态分量或一个参数,以简化计算,也常常使用 Markov 链 Monte Carlo 方法,
对于隐 Markov模型的状态是否取得合理,它们是否包含了所研究的问题的必要的信息,观测过程是否包含了认识模型的足够信息等问题,也是使用中必需了解的,怎样从我们能够得到的信息 (通常主要是样本数据 )来判断上述问题,在统计上也是 至关 重要的,
５,3 隐 Markov 模型 的更为一般的形式
隐 Markov 模型 的观测链 }{ nY 还可以是连续型随机变量,其典型情形有正态分布,Gamma
分布 ),( lkΓ,椭圆对称分布 ( 密度具有 ))()((
||
1 1
2
1 mm?Σ?
Σ
yyf T 形式 ),或上述各种分布的混合分布,例如 ),()1(),( 22221 smlsml NN?+,21 exp)1(exp mm ll?+ 等等,
更一般地,隐 Markov 模型 中的一个状态甚至还可以是一个整个的随机过程,
［ 注 ］ 隐 Markov 模型在应用方面的参考文献量十分大,而它在理论方面的参考文献还不多,例 如关于大数 定 律 的有
Brain G Leroox,Maximum - liklihood estimation for HMM,Stoch astic processes
and their Appl,40(1992),127 - 143,
关于中心极限定理 的有
(1) P,J,Bickel and Y,Ratof,Inference in hidden Markov models I,(local
280
asymptotic normality in the stationary cas e ) Bernoulli,Vol.2,No.3,199 - 228,
1996,
(2) P,J,Bickel,Y,Ratof and T,Ryden,A symptotic normality for the maximum
likelihood estimates for general hidden Markov models,Annals of Statistics,Vol,
26,No.4,1614 - 1635,1998,
６,隐 Markov 模型的应用例 子梗概
６,１ 语音的机器识别
语音识别与下面 ６,２ 段的 脱机手写体汉字识别,这两个问题都有一个共同的特点,它们都有随机性,但 是 也都有一个潜在的基本结构,例如,一个语音也有其基本结构,而 发音时口型的大小,长短,强弱,在口腔中的位置等,却因人而会有随机的变化,即使是同一个人,在不同时间发同一个音也是有随机差异的,而一个汉字的基本笔画及相对关系是确定的,但是笔画的长短,间隔与倾斜度却因人而异,用隐 Markov 模型来描述这种具有结构的随机性相当有效,
这里 以简单化了 的 实际问题为例,这样 可以 使我们摆脱许多细节而 突出 思路,一个词的语音可由许多音素组成,如前所述,隐 Markov 模型在语音识别中,首先是建立一个对应关系,
使一个音对应于一个隐 Markov 模型 ( 或几个 隐 Markov 模型 ),这里的隐 Markov 模型的状态,
一般取为这个音所含的全部可能的音素 ( 或其细分,或其组合 ),例如一个汉语,高,的发音可能由
‘ggg(ge)(ga)(ga)(ga) aaa(ao)(ao)ooo’ (10.31)
( 后面可加上足够的 □ ( 空 )) 等音素组成,而且每个音素究竟维持多少时间单位,对不同的人,
对同一人在不同的时刻 发 这个音,都可能有很大的不同,于是对应于,高,这个音,是 一列相继的状态组成的一个现实 ( 在上面序列中一个字母代表相应汉语拼音,而一个括号中有两个字母则表示过渡 ),由于声转移有内在统计规律,我们 就可以自然地把它考虑为一个 Markov 链
)( nX 的状态序列,但是,)( nX 实际上并不可观测,即我们能测量到的并非这个音对应的
Markov 链本身,而是 每个字母对应的声波 信号的振幅 的采样 ( 或经过某种 预 处理的一些特征量 ),这个可以 观测 到 的现实就是 对于一个音的观测链样本列,状态链与观测链在一起就组成了一个隐 Marlov 模型,其中的状态序列,就是该隐马氏模型的马氏链的轨道 ( 包括,□
（ 空 ),),所以说一个音对应于一个隐 Markov 模型 ( 由于发音人音质的差异,一个音也可以对应与多个隐 Markov 模型 ),一般在语音识别中,把音频信号直接作为观测值常会使计算时的数据量过大,更为实际的是把音频信号的一些特征量,作为压缩了的观测数据,这称为 预 处理,
这些特征可以是音频信号 的过零率,线性预测系数 (LPC,即 Linear Predict Coefficient,
它是把采样值作为 AR 模型 （ 参见第 １１ 章 ） 拟合 后,得到的模型的系数 ),倒频系数 ( CEP,记信号 )(ns 的离散 Fourier 系数 为 ))(( nsF,则 |)))((|(ln1 nsFFCEP?
= ) 等等,一般地常取 12
至 16 个 LPC 作为压缩了的观测数据,也就是把原始观测数据 用 AR( p ) 模型 ( 其中 p = 12
至 16 ) 建模,再 把这样拟合的参数作为压缩了的观测数据,( 为了保证状态到观测的转移只依赖于前一个状态,有时还需要把邻近的状态,堆,起来,例如说,把相继出现的 k 个状态看成为一个状态向量,并把与此 k 个状态向后错一位的相继出现的 k 个状态,看成为其后时刻的一个状态向量,…,使一个时刻的,状态,或,观测,成为一个向量 ),再则,有时也常把在状态已知的条件下,观测向量的条 件分布假定为 ( 一维的或 多维 的 ) 正态分布 ( 一般假定方差矩阵为对角矩阵 ),总之,我们可以把一些特征量作为观测链 { nY,n ≥ 1},其中 nY 可以是多维的 ( 向量 ),序列 ( 10.31 ) 并不能观测到,它正是我们想要知道的状态链,为了解它,就要
281
由对该音的一组观测样本 ( 若干个发这 一 音的声信号或提取的一些特征量 ),去作相应的 隐
Markov 模型的 参数 ),,( 0 BAml = 估计,这就是学习相位,于是语 音的识别 对应于 设置一些模板,一个音对应于一个或数个模板 l ( 有时也可以用最佳状态列的估计作为音或词的代表 ),
这时 模板的代表性,选多少个状态,选多少个观测,对于计算量的控制与,都至关重要,为了能 有效地学习,就需要有足够多的样本,也就是要有这个音或词的重复发音,例如把 m 个不同的人的对同一个音或词的发音的 预 处理加工,作为 m条观测链样本,学习的结果是建立了标准人的 发音与隐 Markov 模型间的对应 ( 也可以让一个音或一个词对应 于 几个隐马氏模型 ),
以 得到一个模板库,
在学习得到模板后,就可以作语音的机器识别问题,即对于一条新给的观测链,在一族 (可能是上千个 )已建立好的隐 Markov 模型的模板中,选取一个与之,最近的,(相对熵最小的 ) 隐马氏链模型 l,它所代表的音或词,就可以判定为这个观测所对应的音或词,这就是 识别 问题,
隐 Markov 模型的基本假定是,状态链是时齐的 Markov 链,而且状态到观 测的转移也是时齐的,但是,在实际语音识别中,常会遇到非时齐的情形,这时如果 照 搬时齐情形的做法,就会引起振荡,然而,用非时齐的隐马氏模型来建模,无论在理论上或 对 实际计算,都远非可取,所以运用隐 Markov 模型 时 极需要灵活处理,实际上,在实践中遇到的问题,绝非只是理论的直接应用,
其复杂的程度 常常 不 比理论问题 差,
如果状态链在一个状态上的停留时间分布,显著 地 不是指数分布 ( 离散采样时它对应于几何分布 ),那么在理论上可以想象,应该改用隐半马氏模型代替隐马氏模型,可惜的是,这种修改后 的模型 的 计算 远为复杂,而在实际使用中,这两个模型的差别也并不如我们想象的那样大,因之这种修改也并没有被广泛采用,
６,２ 脱机手写体汉字识别
汉字的模板称为码本,在线手写体汉字或数字可以看成一维的,用 隐 Markov 模型 的框架建模 就较为简单,而脱机手写体汉字是二维的,即在隐 Markov 模型中的时间参数 n 变成为二维的空间参数 ),( mn,而二维的 Markov 性就远较一维复杂,因为它没有很自然的序,在 确 定相应的观测过程时,为了使得观测链 ( 场 ) 与状态链 (Markov 场 ) 之间有象 ( 10,3) 及 ( 10,4 ) 那样的简单关系,我们常常需要对实际提供的观测链 ( 场 ) 加以改造,例如,在实际上我们得到的手写体汉字是一张黑白象素图 ( 一个取值 0 或 1 矩阵 ),它和笔画的隐 Markov 场之间的依赖关系不能满足 ( 10,3 ) 或 ( 10,4 ) 的要求,我们要对它进行加工,也就是,特征提取,,尽量简化 得 使
( 10,3 ) 或 ( 10,4 ) 能够成立,为了使脱机手写体汉字纳入隐 Markov 场的框架,选取什么作为状态与观测才能较好地体现 Mar kov 性,这是具体操作中的 最 重要的,粗略地说,一般可以选取一个格点上八个方向的笔道作为观测的特征,
此外,数字识别,文字识别,印刷中的 错别字自动检测,航海中的 动态手势识别等等,
也可以用隐 Markov 模型处理,只是有的 情形的 "时间参数 "n 应该 相应地改为 空间参数,
６,３ DNA 序列片断装配及启动子识别
在遗传学中,生物体属性的遗传是由细胞核中的染色体中所含的基因所指导,基因是由
DNA (脱氧核糖核酸 ) 碱基的序列组成,如果把基因看成语句,DNA 的 A,T,G,C 4 种碱基,就象语句中的字母,一条 DNA 链就是一个 A,T,G,C 这 ４ 个 字母组成的序列,读出 DNA 碱基序列对于基因工程,遗传病及许多癌症的发现与治疗都具有重要意义,然而,人类基因组计划目的并不是解决,人类由甚么组成的,,而是要回答,是甚么使得一种有机体区别于另一种有机体,,一个 人 的 基因组的全序列,是一个 包含约 30 亿个由 A,T,C,G组成的序列,它们分为 23
对染色体,与一条性染色体,生物体中最小的 DNA链,约含一百万碱基对 (bp,base pair),即由
A,T,C,G 4 个字母组成的长度为一百万的一个序列,基因是染色体 DNA序列的一些子序列,它们给出了 有关 蛋白质的信息,至今,人们知道在基因之间还存在 DNA 序列,可是还不清楚它们的作用,一个人的基因组的全序列可以说是他的生命的设计书,它能为了解遗传与发展等生命
282
奥秘奠定基础,而 得到基因组的全序列,只是我们读懂这本生命之书的起始的一步,生物信息论的第一步是 DNA 序列的分析,包括 DNA 序列的比对 (alig n ment),拼接 (assembly); 寻找基因 ( gene findi ng ),调控研究 (regulation),基因功能 (function)研究等,
(1999 年 12 月 1 日,人类基因组科学家完成了人类 22 号染色体的测序,这是在 23对染色体中第二小的染色体,有 3350 万个碱基对,由长短不同的 12个 DNA 片断组成,最长的片断有 2300 万碱基对,其它的片断,长度在
1000 碱基对到 58.3 碱基对之间,在 22 染色体上发现了 679个基因,其中 545 个是功能基因,与先天心脏病,低免疫性,精神分裂,智力低,许多恶性肿瘤如白血病等有关,另外的 134 个不是 功能基因,2000 年 5月 8 日,人类基因组又基本上给出了第 21号染色体的基因图谱,这是最小的染色体,有 225 个基因,与白血病,痴呆,肌肉萎缩,
燥狂性抑郁,部分癌症等有关,人类基因的总数比水稻要少,说明人类的基因可能更为有效 ),
然而,DNA 序列 是 非常长的链,用 目前的设备是无法读出的,为了读取一条链,就 需要 线把它 打断成一段一段,使每一段都可以用设备读出,而在解读 DNA序列时,遇到的困难是,
(1)目前快速读出 DNA序列的设备只能读出几百到一千个碱基对的链,
(2)快速打断 DNA序列的方法是随机的,不能控制断点的位置及段数,
(3)打断后的片断混为一体无法区分次序先后,
解决困难的一个办法是,同时对同一个 DNA序列的多个拷贝 (copy) 作随机打断及作随机片断读出,利用断点几乎概率为 1 地不会在同样位置重复,并利用不同的拷贝片断的相互核对,将片断连接起来 (Assembly)得到整体的排序,但是,在 读出时可能有错读,漏读及误增字符等错误,
且在原序列中还有重复的片断,这就使拼接问题十分复杂,上世纪末的人类基因测序中,人类基因组的科学家对染色体 DNA 序列,利用一些只出现一次的序列 ( 称为 marker),将整条序列分段 ( 段长为 2M - 3M) 读出并拼接起来,以降低拼接的复杂性,而 Cerela 公司的科学家则更着重于技术与算法,他们利用最新的毛细管测序机,可以从一段序列的两头同时测得两个
500 - 700bp 小段 DNA 序列的这个特点,将人类基因组科学家拼接思路中忽视的匹配 ( mate ) 信息利用起来,为不用 marker 直接测序开辟了新路,
读出了全部基因组的 DNA 序列后,基因检测是后基因处理中的首要问题,即要知道序列中各部分的含义 ( 如,哪里是基因,哪里是基因的蛋白编码部分,哪里是其调控部分,…),不知道这些,DNA 测序的结果不会有用,每个基因有几部分构成,启动子 ( Promoter ) – 像是基因的引言 ; 外显子 ( Exons ) – 编码部分,它在基因 DNA 序列中也可能是被内含子 ( introns ) 分隔的子序列 ; 内含子 ( Introns ) – 外显子间的序列,它的功能至今还不清楚 ; 终止子
( Terminator ),它们可示意为,
启动子 外显子 1 内含子 1 外显子 2 内含子 2 外显子 3 终止子
对于 已 读出的 DNA 序列,Karlin 和 Berger 有一个基因检测的软件,称为 Genescan,利用基因编码区的统计特征,即在没有漏读 ( 也称 删除,记为 D)与没有 增 读 (也称插入,记为 I)时,
利用 DNA 编码序列具有很强的 3-连子特性 (A,T,C,G 4 个核酸组成 64 个 3 – 连子,即 3 个 相连的 字母 组成 的 一个 编码,在此 ６４ 个编码 中除了 3个编码代表终止符号以外,其它的 ６１ 个编码 分别代表 20个氨基酸 (amino) 的编码 ( 一个氨基酸可以有几个 不同的 编码 ),在编码区无漏读或无 增 读 的条件下,Karlin 和 Berger 采用了动态规划方法,编制了基因检测软件 Genescan,
达到了很高的识别率,但是在编码区有漏读或有 增 读时,就会发生 3 - 连子的错误,且 随后的所有 3 - 连子就全 被 搞乱,接着 就 会 出现一系列的错误,如果还用 Genescan 就会发生严重的错误,所以,在用这个软件前,必需先在观测序列中,标识错读,漏读 ( 删除 ) 或 增 读 ( 插入 ) 的碱基对,以便避免出现大量的编码错误,并尽可能的校正错读的碱基对,
对于标识错读,漏读或增读碱基对,G.Churchill 首次提出了 用隐 Markov 模型解决 DNA
片断装配问题,他的基本思想是将由随机枪打断长链中的每个碱基对,附以读出的属性,即读对 R,错读 M,漏读 D 或增读 I ( 其中最误事的是漏读 D 或增读 I,因为它们引起一系列后继错误 ),把这些属性看成 一个 Markov 链的 状态的 样本,而把读出的碱基序列片段看成为相应于这个 状态 链的观测链,即 设观测过程为可能带有小概率的 漏读 D与 增读 I （ 在实际解读中,漏读与增读度不可能很多 ） 的 DNA码区序列,这时,隐 Markov链的状态 空间就是 },,,{ IDMR,
其中 R 表示读对的碱基对,M 表示读错的碱基对,而隐马氏链的样本列,就是对观测序列上的
283
碱基对标出的属性 MDR,,与 I 的序列,这样的考虑 是,基于在编码区中的序列上,氨基酸编码的转移则相当集中,因而,隐 Markov模型能很好地刻画与鉴别有 无删除,错读与插入 （ 增读 ）,这样就 把问题纳入了隐 Markov 链的框架,再利用统计性质,就可以 纠正在使用 寻找基因的软件 Genescan 时可能导致的错误,在学习相位 ( 即学习问题 ),由给定的观测过程的一组样本找出隐 Markov 模型的参数族,可以用 Baum - Welch 算法,而在运转相位 ( 即解码问题 ),对给定的观测样本找出相应的 Markov 模型的状态,可以用 Viterbi 算法,( 相仿的思路也出自借鉴 Gibbs 场的思想,即在学习时先设定初值,再轮番修改参数和状态 ( 见第 15 章中的 EM 算法 ),而在运转时用 Bayes 决策找出 最佳 状态 ),
在 DNA 序列中,启动子是一个特殊位点,它 反映了在 DNA 转录过程中发挥定位功能的
DNA 片断表现在序列上的规律性,所以,启动子识别是生物信息论中首先要解决的问题,隐
Markov 模型也可以用来识别启动子 (Promotor ) 等特殊位点及编码区,( 这些特殊位点还包括增强子
( enhancer ),抑制子 ( silencer )),转录起始位置 ( TSS,transcription start ing site),外显子 (exon),内含子 (i ntron ),受子 ( accepter ),给子 ( donor ) 等等 ),
利用已有的生物学知识及数据库,构造合适的数学模型,例如,合适的隐 Markov 模型,再与数据进行反复的交叉验证与修改,是得到成功的识别方法的关键,
习题 10
１,证明 在均值为 m,方差为 2s 的分布密度中,正态分布的熵最大 ； 而 在均值向量为 m,方差矩阵为
Σ 的多维密度中,Gauss 分布 ),( ΣmN 的熵最大,请读者自己写出最大熵的表达式,
２,证明在有限区间 ],[ ba 上取值的分布密度中,均匀分布的熵最大,这个熵为 )log( ab?,再把这个
结论推广到多维情形,
３,设 )(,),(1 xgxg mL 都是非负有界函数,那么,在满足条件
)(,)()( midxxpxg ii ≤=∫ a
的分布密度 )(xp 中,当且只当
))()((
0
11)( xgxg mmCexp ll ++?= L
时熵为最大,其中
∫ ++?= dxeC xgxg mm )()(( 11
1
ll L,
而 mll,,1 L 是由
dxexg xgxgi mm ))()(( 11)( ll ++?∫ L )(,))()(( 11 midxe xgxgi mm ≤= ++?∫ lla L
确定的正数,
４,能不能把 3 推广到 )(,),(1 xgxg mL 不全为非负,也不全为有界的情形?
