202
龚光鲁,钱敏平著 应用随机过程教程及其在 算法与智能计算中的应用
清华大学出版社,2003
第 8 章 Monte Carlo 与 Markov Chain Monte Carlo (MCMC) 方法
在许多很复杂的统计问题中,有时很难直接对各种统计方法进行理论分析,为了评估它们的优劣,常见的实用办法是做随机模拟,即设法按问题的要求与条件去构造出一系列的模拟样本,用它们的样本频率代替对应的概率作统计分析与推断,观察由这些摸拟样品所作出的推断的正确率,因为在概率论初期发展时,随机模拟的原型常常来自博采,于是 人们就以博采之都 M onte Carlo 作为随机模拟方法的别称,久而久之,Monte Carlo 方法作为名称倒比随机模拟方法更为广泛地被常用了,相仿地,人们还把组合计算中的某些随机模拟方法,
称为 Las Vegas 方法,这是以美国的博采城 Las Vegas 命名的,
一般 随机 模拟 方法 的长处是,计算的复杂度不依赖于计算空间的维数,因此,在计算非常 高 维数的积分 （ 或多指标的求和 ） 时,Monte Carlo 方法比通常的计算方法有明显的优越性,而很多的 Monte Carlo 计算问 题,可归结为计算积分或庞大和数的问题,
1,计算积分的 Monte Carlo 方法与采样量估计
通过构造独立同分布随机数,计算积分的 Monte Carlo 方法,也称为 静态 Monte Carlo
方法,其思想可以在本节中,通过估计最简单的积分 ∫
b
a
dxxf )( 得到阐明,对于高维积分,
其思路与一维积分是一样的,另一方面,在甚高维情形,因为计算量太大,用静态 Monte
Carlo 方法处理速度太慢,我们在第 2 节中通过构造 Markov 链的极限不变分布,来 模拟计算积分的方法 （ 称为 动态 Monte Carlo 方法,Markov 链 Monte Carlo 方法,M C MC),将更为适用,
1,1 用频率估计概率来计算积分的 Monte Carlo 方法
假定 Mxf ≤≤ )(0,那么由积分的面积含义有
||)( Sdxxf
b
a
=∫,(其中 || S 为 })(,:),{( yxfbxayxS ≥≤≤=? 的面积 ),
考虑平面区域 ],0[],[ Mba ×= 上的均匀随机变量 x,则
MabSPp )(
1)(
=∈=
x
∫
b
a
dxxf )(,
对于 N 个独立的? -均匀随机数 )(,Nii ≤x,记
},,{ 1 NSN xx L= 中落在 S 中的频数,
于是,利用大数定律便知
=^I
N
NMab S)(? （ ８,１ ）
203
是积分 ∫
=
b
a
dxxfI )( 的相合估计,即对于任意的 0>e,当 ∞→N 时,有
∫ →>
b
a
S dxxf
N
NMabP 0)|)()((| e,
又由于 SN 服从参数为 ),( Np 的二项分布,所以我们有 =^IE ∫
b
a
dxxf )(,即
^I
是积分
∫
b
a
dxxf )( 的无偏估计,此估计的方差为
=^IVar N ppMab )1()( 22 N Mab
I
Mab
I
Mab
))(1()(
)( 22?
=
IIMabN ])[(1= )1( NO=,
又因为方差代表平均平方误差,所以我们就说,积分的估计
^I
的误差为 )1( NO,
频率法所需采样量的估计
我们知道统计估计是以区间估计来给出概率误差的,即对于允许误差 e,及允许以 d 的失败概率,来确定样本量 N (依赖于 de,),使
∫ <>?
b
a
dxxfIP de)|)((|
^
,
1 用 Chebyshev 不等式估计样本量
最粗略的确定样 本量 N 的方法,是通过 Chebyshev 不等式来估计,即由
2
^
^^^
)|(|)|)((| eee IVarIEIPdxxfIP
b
a
≤>?=>? ∫,
要求 dee <
= 2
22
2
^ )1()(
N
ppMab
IVar,考虑到
4
1)1( ≤? pp,我们要确定 N,只需满足 de <
2
22
4
1)(
NMab 就够了,也就是说只要
2
22
4
)(
de
MabN?> 即可,但是这个采样量的估计过于粗糙,
2 用中心极限 定理的近似估计样本量
在 N 大时 (这总假定满足 ),用中心极限定理也可以得到采样量 N 的估计,中心极限定理断言,MabI )(^?= NNS 的近似分布为正态分布 (N ∫
b
a
dxxf )( ),^IVar,令 df 满足
204
∫
∞?
=
df
dp dxe
x
2
2
2
1,( 即 ∫
=2
2
2
121 2
d
d
f
f
dp dxe
x
),
于是有
∫ ≤≥?
b
a
IVardxxfIP df d )|)((|
^
2
^
,
选取采样量 N 使 ef d ≤
^
2
IVar,即 2
2
2
22 )1()(
df
e≤
N
ppMab,则就有
∫ ≤≥?
b
a
dxxfIP de)|)((|
^
,
所以只需取
2
2
2
22
4
)(
dfe
MabN?>,
3 采样量的估计的比较
我们把用 Chebyshev 不 等式得到的采样量的估计与用中心极限定理得到的采样量的估计作比较,可以看到只要 12
2
<<ddf,用中心极限定理得到的采样量的估计,要比用
Chebyshev 不等式得到的采样量的估计小得多,
［ 注 ］ 还可以应用概率理论中的大偏差速率函数来求得采样量的估计,可以 证明由大偏差速率函数能得到采样量的估计
22
ln
e
d?>N,
当 d 小时,dln? 趋于 ∞ 的速度要比 2
2
df 趋于 ∞ 的速度快,所以,以用中心极限定理得到的采样量的估计为最好,
１,2 用样本函数的平均值估计积分的 Monte Carlo 方法 ― 期望法
期望法的核心思想是把积分看成某个随机变量的期望,最自然的是看成 ],[ ba 上均匀
随机变量的期望,设 ],[~ baUh ( ],[ ba 上的 均匀分布 ),则
)()()( hEfabdxxfI
b
a
==∫
,(8,2)
于是对于 N 个独立的 ],[ ba 上的均匀随机数,可以用矩估计
N
ffabI N )()()( 1^^ hh ++?= L (8,3)
205
作为 ∫=
b
a
dxxfI )( 的估计,显然它也是无偏估计,而且
N
fVarabIVar )]([)()( 2^^ h?= [ ]222 )]([)(1)( hh EfEf
Nab=
222 1)(1)( I
Nab
dxxf
Nab
b
a
= ∫ 2)[(1 IMIabN≤ )( ^IVar=,
可见期望法比频率法更有效,但是期望法多算了 N 次函数值,这是需要花费计算时间的,
所以在全面考虑有效性与计算时间的得失时,有时人们也愿意采用频率法,
1,3 减少方差的技术
用 Monte Carlo 方法计算积分 ∫
b
a
dxxf )( 时,未必一定要使用均匀随机数,事实上,从
],[ ba 上取值的任意一种随机数出发,都可以得到 ∫
b
a
dxxf )( 的估计量,而且在 0)( ≥xf 时
显见,值 )(xf 大的 x 对于积分 ∫
b
a
dxxf )( 有更大的贡献,由此得到直观启发,所用的随机数的分布密度的形状越象 )(xf,则 就 越合理,这个思想就是 重要度采样法,
1,?g 采样法
假定分布密度 )(xg 在 )(xf 非零处恒正,则积分 =I ∫
b
a
dxxf )( ∫=
b
a
dxxgxg xf )()( )(,对于密度为 )(xg 的,?g 随机数,V,我们有
])( )([ VVgfEI =,(8,4)
于是对于 N 个独立的?g 随机数 NVV,,1 L,关于积分 =I ∫
b
a
dxxf )( 可取估计
])( )([1
1
1
)(^
V
V
g
f
NI
g?
= ])( )(
N
N
g
f
V
V++L,(8,5)
显见它也是无偏的相合估计,利用 Schwartz 不等式得到
∫ ∫=
b
a
b
a
dxxgdxxgxg xf )()(])( )([ 2 ∫
b
a
dxxgxg xf )(])( )([ 2
206
≥
2
])()( )([∫
b
a
dxxgxg xf =
2
])([∫
b
a
dxxf,
而且上式当且仅当在 )( )()( xg xfcxg = 时 (即 )()( xcfxg = 时 ) 达到极小值
2
])([∫
b
a
dxxf,这说明 NIVar
g 1
)(
)(^
= ])(])( )([[ 22 Idxxgxg xf
b
a
∫ 的最小值在 )()(0 xcfxg?=
时取到,又因为 )(0 xg 为密度,故
∫
= b
a
dxxf
c
)(
1,对应的 0)( )(^ 0 =gIVar,此时
∫==
b
a
g
dxxfcI )(1
)(^ 0
是无估计误差的精确值,
综上讨论可知,要使方差达到最小,就应该用 )()(0 xcfxg?= 作为参考密度,这时用 Von
Neumann 取舍方法必须知道参考密度 )(0 xg,而计算 )(0 xg 又必须知道常数 c,后者却是于积分的倒数,这样的耦合性的存在使这种方法不能真正用于计算积分,但是由此我们至少得到了以下的认识,即 只要密度 )(xg 的形状与被积分的函数近似,用
)(^ g
I 作为 ∫
b
a
dxxf )(
的估计,就会降低方差,这个想法还可以推广到 )(xf 未必为非负函数的情形,这就是下面的概念,
定义 ８,１ 分布密度为
∫
= b
a
dxxf
xfxg
|)(|
|)(|)( 的?g 采样,称为 (关于 )(xf 的 ) 重要度采样 （ Importance Sampling),
重要度采样不能直接通过取舍原则实现,近似地实现重要度采样可以采用下一节中的
Markov 链 Monte Carlo 方法,这个方法同时给出连积分 ∫
b
a
dxxf |)(｜ 的近似算法,
在实践中人们也往往按照重要度采样的思路,灵活地寻找常用的已知类型的密度 g,使它在峰值附近与 | )(xf | 较接近,以便达到降低估计的方差的目的,
修正的重要度采样 (Modified Importance Sampling)
对于?g 采样,假定有非负函数 )(xh,满足
207
∫ >=
b
a
dxxgxh 0)()(a,
而且 a 已知,那么我们可以采用 )(xh 作为修正乘积因子,显见对于?g 随机数 V 有
)(
])( )([
)( VV
Va
Eh
g
fE
dxxfI
b
a
== ∫,(8,6)
如果我们先放弃对于估计的无偏性要求,而只要求估计的相合性,则对于 N 个独立的?g
随机数 NVV,,1 L,我们可以通过比值,构造 ∫=
b
a
dxxfI )( 的如下的估计量
)()(
)(
)(
)(
)(
1
1
1
N
N
N
hh
g
f
g
f
I VV V
V
V
V
a ++
++
=?∨ L
L
,(8,7)
显见它是 ∫=
b
a
dxxfI )( 的相合估计,在某些假定下,它还是渐进无偏的,即
IIEN =∨∞→ )(lim,
而且
∨I
保留了重要度采样的特性,即当 )( )()( xg xfcxh = 时,
∨I
就是 ∫=
b
a
dxxfI )(,于是,
只要当 )(xh 与 )( )(xg xf 近似,就会降低方差,注意对于给定的分布密度 )(xg,)(xh 应选取得尽量与 )( )(xg xf 近似,这个修正乘积因子 )(xh 是用来再一次降低由于密度 )(xg 与被积函数
)(xf 的倍数不够像所带来的失误而设置的,而当 1)( ≡xh 时,就退化为是?g 采样,这相当于对?g 采样不再作修正,
2,Markov 链 Monte Carlo (MCMC)
用 Markov 链的样本,来对不变分布,Gibbs 分布,Gibbs 场,高维分布,或样本空间非常大的离散分布等作采样,并用以作随机模拟的方法,统称为 Markov Chain Monte Carlo
(M C MC) 方法,这是 动态的 Monte Carlo 方法,由于这种方法的问世,使随机模拟在很多领域的计算中,相对于决定性算法,显示出它的巨大的优越性,而有时随机模拟与决定性算法的结合使用,会显出更多的长处,
MCMC 至少可以用在以下几个层面,
208
1,用于生成较复杂的随机数,
* 实现对高维分布 （ 或高维格点分布 ） p 取样,得到 p 随机数,
* 是实现重要度采样的一种方法,对 |)(| xf 的重要性采样,就是取得
p ( ∫
== dyyf xfx |)(| |)(|)(p ) 随机数,
2,实现高维积分 （ 或项数极多的求 和 ） 的数值计算 （ 典例是 Gibbs 分布的各种泛函的平均值的计算 ）,对于 0)( ≥xf,作以 p ∫
== dyyf xfx )( )()(p 为极限分布的 Markov 链 nX,
利用遍历定理可以由这个 Markov 链的一条轨道,得到分布密度 )(xp 的估计,记为 )(^ xp,再用
)(
)(
^
x
xf
p
作为 积分 ∫= dxxfI )( 的估计,
3,用模拟方法估计最可几轨道,例如,如果模拟了 100 条轨道,那么就能以大概率推断,最可几轨道就在这些轨道的邻近,当统计量的分布未知时,可以用模拟方法从频率估计置信限,
４,用被估参数的 Bayes 分布 （ 参见第 9 章 ） 的取样,来估计参数,
５,求复杂样本空间上函数的极值 ( 模拟退火 ),
2,1 Gibbs 采样法 ( Gibbs sampler )
1,用 MCMC 方法得到 Gibbs 分布的样本与估计 Gibbs 分布
在第 ６ 章中,我们考虑了在 d 维的?N 格点集上的 Ising 模型的 Gibbs 分布
∑
∈
=
S
H
H
e
e
h
hb
xb
xp )(
)(
,由于所涉及的状态空间 (全体组态的集合 ) dNS ),,1{}1,1{ L?= 非常大 (例如,把一幅 256× 256 个采样点的黑白图像看成一个组态,则 256,2 == Nd,S 中有
)1022(2 18000600002256256 16 >>=× 个元素 ),这就使得分母中的求和无法实际完成,而 MCMC
方法就是以通过构造一个以这个 Gibbs 分布 ∑
∈
=
S
H
H
e
e
h
hb
xb
xp )(
)(
为不变分布的离散时间的
Markov 链 nX (它就是 Glauber 动力学中的连续时间的 Markov 链的离散时间采样 ),作为模拟计算的基点的,构造的 Markov 链必须易于计算,所以我们要求它的概率转移速率只容许在组态的一个格点上变动,这样的变动方式,称为 Gibbs 方式,这种抽样方法称为 Gibbs 采样法,或者 Gibbs 样本生成法,这个 Markov 链的不变分布正是此 Gibbs 分布 p,我们还要求此 Markov 链的转移矩阵满足 P →? ∞→nn pT1,这就是说,要求 Gibbs 分布是 Glauber
动力学的极限分布,
于是,当 n 大时,nX 的一个样本可以近似地认为取自 Gibbs 分布的一个样本,即按此
209
Markov 链沿任意一条轨道充分发展,就得到 Gibbs 分布的近似取样,
再则,Gibbs 分布的归一化常数 (称为 配分常数 ) )(xb
h
H
S
e?
∈
∑,是一个巨大的求和,即一个,离散的,积分,用随机模拟法计算这个 ＂ 离散的 ＂ 积分的最佳随机数正服从 Gibbs 分布
(即重要度采样 ),对于 Gibbs 分布的取样,用通常的取舍原则常常并不可行,例如,分别取
1=C,)()( xbx Heh?=,而参考密度 )(0 xp 为组态空间上的均匀分布,这时 )(xbHe? 的值常常小得超出计算精度,而求和变量 x 的范围是庞大的组态空间,这就导致求和无法实际计算,所以需要用 Markov 链 Monte Carlo 方法,用 MCMC 方法生成了以 Gibbs 分布为极限分布的 Markov 链 nX 以后,由遍历定理用 Markov 链的一条轨道,可给出极限分布 （ Gibbs 分布 ） 的估计,对于充分大的 N,可令
))()((1 2}{1}{^ NN XIXIN xxxp L+= +,（ ８,８ ）
再用 )(xbHe? 除以 Gibbs 分布在 x 处的估计值 ^
)(
x
xb
p
He?
,作为配分函数的估计,在理论上这个估计应该与 x 的取法无关,但是,在实际计算中对多个不同的组态 ix 分别估计此和数后,
再作平均常常能降低方差,
在第 6 章中,我们曾给出了用 Glauber 动力学构造的两个不同的连续时间的 Markov 链
(对应于两个不同的转移概率速率矩阵 Q),它们都以 Gibbs 分布 p 为极限分布,而且都是可逆的,
较为深入的理论研究表明,使用不可逆的,且以 p 为不变分布的 Markov 链作 Markov
Monte Carlo,会加快这个极限的收敛速度,然而,在另一方面这种做法又会增加计算的复杂程度,再则,为减少估计的方差而作的努力也常会增加计算时间,这就是说,在计算中,我们会面临难以两全的抉择,在实际中如何采取折衷,既要看问题的性质,又要参考实践的经验,没有统一的原则,
用以完成 MCMC 采样操作的 Markov 链,可构造如下,
在第 6 章中,对于 d 维有限格点上,由具有两个自旋的组态空间上的能量函数
∑ ∑=
相邻yx x
xhyxH
,
)()()(21)( hhhh,(8,9)
可构造 如下的转移概率速率
≠ == 的其它情形 ）xh xhxxh (0 )(),,(
xxC
q (8,10)
( ),( xxC 的两种取法各为,
))()((),( xxbx HHexC= (8,11)
或
210
))()((1
1),(
xxbx HH xexC+= ),(8,11) ’
它们决定的 Glauber 动力学,分别对应的两个不同的连续时间的 Markov链,它们都 以 )(hH
为能量函数 的 Gibbs 分布 p 为可逆不变分布,而且 p 还 是转移矩阵的 极限分布,P?→? ∞→tt)( pT1=,如果考虑间隔为 t? 的采样时间,其中 t? 充分小,我们还可以进行如下的近似 计算,假定初始组态为 V,在时刻 t? 它以概率 txC?),( V 转移到组态 xV,而以概率 ∑
∈
dNx
xCt
},,1{
),(1
L
V 停留在原来的组态,这就近似地得到了 Markov 链在 t? 时刻所处的组态 1V,再以类似的方式得到 Markov 链在 t?2 时刻所处的组态 2V,依次下去,当采样数 n 充分大时,组态 h 在这段 },,{ 1 nVV L 中出现的频率,就用作 hp 的估计
^
hp,
Gibbs 分布的一些统计平均量的模拟近似
对于上面定义的,以 Gibbs 分布为极限的 Markov 链,我们用它的 一段样本 NVV,,1 L,
可以给出如下的 Gibbs 统计平均量
∑
∑
∈
∈=
S
H
H
S
e
ef
F
x
xb
xb
x
x
)(
)()(
(8,12)
的模拟近似
n
f
F
n
i
i∑
== 1^
)(V
,(8,13)
它是 F 的无偏估计,不难证明,存在数 ∞<),,( pPfV,使
^F
的方差满足
)1(),,()( ^ nonPfVFVar += p nPC f ))(1( ||||13
2
≤,(8,14)
其中 ∑=
x
x )(|||| ff,而 )(PC 为此 Markov 链的转移矩阵 P的 Dobrushin 收缩 数,
Gibbs 采样法还可以有效地用于 Bayes 方法中的后验密度的采样的模拟计算 （ 见第 ９ 章 ）,
2,2 M etr opolis 采样法 ( Metropolis sampler )
１,Gibbs 分布的采样的 Markov 链 Metropolis 方法
为了突出主要思想,下面我们把组态空间 ( 状态空间 ) dNS ),,1{}1,1{ L?= 简单地记为
},,1{ KL,对于组态空间上给定的分布 p,Metropolis 构造了以
211
)(),,1min(
~
ijpp
i
j
ijij ≠= p
p
(8,15)
( ∑
≠
=
ij
ijii pp 1 ) 为转移的时齐 Markov 链,其中
~P )( ~
ijp= 是一个对称的互通 转移矩阵,
称为 预选矩阵,或访问方案,使用它是为了减少 或简化 状态间的连接,以加快 Markov 链的分布向不变分布收敛的速度,显见 p 是它的可逆分布 （ 注意在 Gibbs 分布情形,状态 ji,体现为组态 hx,,于是在计算 （ ８,１５ ） 式的转移概率时,就只需算比值
)()()(( hxpp hxb
h
x ≠= HHe,而并不需要计算配分函数,Glauber 动力学的构架,也正是用了这一点 ）,由这个有限状态 Markov 链的互通性,我们有
P →? ∞→nn pT1=,
因此,在时间发展充分长以后,我们可以 用 Metropolis 的 Markov 链所处的状态,作为按分布 p 取的样本,也就是说,与 Gibbs 采样法一样,Metropolis 方法也给出了在计算机上模拟 p - 随机数的一个算法,Metropolis 提出的这种采样法,称为 Metropolis 采样法,它与
Gibbs 采样法的不同处在于,对于 Metropolis 采样法而言,任意两个组态 hx,,只要预选概率 0
~
>xhp 就可以转移,
Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作,
(1) 设当前为时刻 n,取的状态为 ix n =)(,对它作随机扰动,即取一个分布为
),,( ~1~ iKi pp L 的随机数,设为 j ;
(2) 若 1≥
i
j
p
p
,则将状态更新为 jx n =+ )1( ； 否则 进行 （ ３ ）；
(3) 独立地取一个 ]1,0[U 随机数 U,如果
i
jU
p
p≤
,则将状态更新为 jx n =+ )1( ； 否则状态不更新,即令 ix n =+ )1(,
（ 请读者证明,如此由 i 到 j 的转移的可能性恰是 （ ８,１５ ） 式规定的转移概率 ）,
~P
的对称性并非必要,理论分析指出,经过适当的选取 (研究矩阵
~P
的第二个特征值 ),
使用非对称的
~P
可能加快收敛速度,对于 非对称的预选矩阵
~P
,Metropolis 采样法所 构造的 Markov 链的转移 应取 下 式,
212
)(),,1min( ~
~
~
ij
p
ppp
iji
jij
ijij ≠=
p
p
,(8,16)
此时所构造的 Markov 链仍然以 p 为可 逆分布 （ 请读者验证,并给出这种 Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作步骤 ）．,
具有非对称的预选矩阵的 Metropolis 采样法 的一个特殊情形就是 独立的 Metropolis 采样,即 预选转移矩阵与出发点无关 的 Metropolis 采样,此时预选矩阵 各行是一样的,因此 可以写为
1P =~ rT,
记
j
j
jw r
p=
,结合拒收原则,对独立的 Metropolis 采样实际 操作 为
（ 1） 按分布 r 采样,得到 j ；
（ 2） 若 1
)(
≥
nx
j
w
w,则取 jx n =+ )1( ; 否则进行 (3);
（ 3） 取 ]1,0[U 随机数 U,修改 )(nx 为 )1( +nx,
≤=+
其他情形 ）(
),1min((
)(
)1(
)(
n
x
j
n
x
w
wUj
x n 。
注 １ 一般地 Markov 链的两个状态之间都可能转移,因此 Markov 链相当于连续取
值情形中的 ＂ 高维情形 ＂,预选矩阵的设计是尽可能地减少,或者简化 状态间的连接,以减少计算量,有时 它相当于局部化,
注 ２ 在独立的 Metropolis 采样法,且 p 是一维密度的情形,我们还可以把 Metropolis
算法视为 Von Neuman 拒收原则在某种意义下的推广,此时 拒收原则中的 )(0 xp 就相当于预选矩阵,为了方便 做 比较,我们把 p 改记为 )( xp,此时 独立 的 Metropolis 算法可以设计如下,取一个与 分布密度 )( xp 取值范围差不多的随机数的独立样本 L,,10 hh,记其分布密度为 )(0 xp,再取与它独立的 ]1,0[U 独立随机数列 L,,21 UU,令 00 hx =,然后归纳地定义
≤=
+
+
+
+
)
p
p
p
pU
n
n
n
n
n
nn
n
其它情形(
)1,)( )()( )(min((
10
01
1
1
x
h
x
x
hh
x,
于是当 n 大时,nx 就近似地是 p 随机数,(这里 )(0 xp 起了 ijp 的作用 ),把 它与拒收原则相
213
比较可以看出,Metropolis 算法不需要假定 )( )(
0 xp
xp 的有界性,为此付出的代价是,不能从是否有 )( )(
0 h
h
Cp
pU ≤ 来决定是否 将 x 取 为 h,代 替一个 更新的 随机数 h 以一个 随机数的 近似序列 nh,并且从是否有
≤nU 序列 )( )(
10
1
+
+
n
n
p
p
h
h 与其前面的值
)(
)(
0 n
n
p
p
x
x 的比
来决定是否把 nx,更新 为,1+nh,所以 可以说,独立的 Metropolis 算法主要用在找不到
)(0 xp 使 )( )(
0 xp
xp 有界的情形,特别是在 )( xp 只 能 得到测量的数据值,而没有显式表示的情形,然而,在实际计算中 一般地 只在 ＂ 高维情形 ＂ 才有必要使用 Metropolis 算法,
注 ３ 比较这 Gibbs 采样法 与 Metropolis 采样法 的收敛速度,是一个重要的理论问题,
注 ４ 使用 非对称的
~P
的 另一种 Markov 转移的 取法,是 如下的 Hastings 推广了的 Metropolis 采样
)(,
1 ~
~
~
ij
p
p
spp
iji
jij
ij
ijij ≠
+
=
p
p
,(8,17)
其中
~P
为非对称的转移矩阵,满足,ijp
~
与 jip
~
要么同时为零,要么同时为正,而 )( ijs 是一个非负的对称矩阵,满足
)(,1
1 ~
~ ij
p
p
s
iji
jij
ij ≠≤
+
p
p
,(8,18)
以上几种采样法都需要知道 p 的解析表达式,在并不知道它的表达式时,下面的方法对于采集以它为分布的样本,能起很好的作用,
2,3 通过条件分布对分布 p 作 随机采样的 Gibbs 方 法
Gibbs 采样法 中 Markov 链的状态的更新 按 Gibbs 方式 进行,即更新时 只变动组态的一个分量,它 实行起来较为简单,所以这种思想可以用来得到高维总体的样本,记 高维总体的分布为 p （ 未必是 Gibbs 分布 ）,Gibbs 方式 更新 的思想,在用于甚高维总体或复杂总体的取样时,主要是 通过 p 的 条件分布族,来构造一个不可约正常返的马氏链 nX,使它以 p 为不变分布,由于 n 充分大时,nX 的分布渐近于 p,nX 就可以近似地作为 分布 p 的 总体
214
的样本,下面我们以可数状态空间的 分布 p 为例来进行 讨论,其处理问题的思想,也适用于生成有密度的总体的样本,
1 由条件分布族得到以 p 为 不变分布的 Markov 转移的思想
Gibbs 样本生成法的要义,在于把高维总体的取样,化为一系列一维分布的 取 样,而后者是容易做到的,甚至可以利用 SAS 等软件得到,但是,在具体运作时,对于较大的 m 及一个给定 的 m 维分布 p,事情就远非如此简单,设 p ),,( 1 mxx Lp=,记 ),,( 1 mxxx L=,
),,( 1 myyy L=,再记在第 1 至 k - 1 个分量固定为 11,,?kyy L,同时将第 1+k 至第 m
个分量固定为 mk x,,x L1+ 的条件下,第 k 个分量在 ky 的条件概率为
|( kyp 11,,?kyy L,)x,,x mk L1+,
我们定义如下的转移概率
∏
=
==
m
k
xy yxpp
1
),( |( kyp 11,,?kyy L,)x,,x mk L1+,(8,19)
容易验证 )( xyp 是一个转移矩阵,我们 证 明 p 是它的不变分布,事实上
∑
x
yxpx ),()(p ∑=
x
x)(p ∏
=
m
k 1
|( kyp 11,,?kyy L,),,1 mk xx L+
∑ ∑∑=
mxx
x
m
m
m
x xx
xxyxx
,,1
21
1
2
1
1 ),,(
),,,()],,([
L L
LL
p
pp ∏
=
m
k 2
|( kyp 11,,?kyy L,),,1 mk xx L+
∑=
mxx
mxxy
,,
21
2
),,,(
L
Lp ∏
=
m
k 2
|( kyp 11,,?kyy L,),,1 mk xx L+
∑ ∑∑=
mxx
x
m
m
m
x xxy
xxyyxxy
,,21
321
21
3
2
2 ),,,(
),,,,()],,,([
L L
LL
p
pp ∏
=
m
k 3
|( kyp 11,,?kyy L,),,1 mk xx L+
∑=
mxx
mxxyy
,,
321
3
),,,,(
L
Lp ∏
=
m
k 3
|( kyp 11,,?kyy L,),,1 mk xx L+
)(),,( 1 yyy m pp === LL,
[ 注 ] 事实上,我们并不需要知道 )(xp 的表达式,而只需知道在固定其它一切分量的条件下,余下的某一个分量的条件分布,就可以得到以 ),( yxp 为转移的时齐马氏链的样本,
2,构造 Markov 的方 法
从已得到的 Markov 链在时刻 n 的样本 nX,去求转移到时刻 1+n 的样本 1+nX 时,转移 ),( yxp 是一维 的 条件 分布 的乘积,而 这些一维条件分布的样本是较为容易生成的,即可
215
按以下程序逐个地得到 1+nX 的样本的各个分量 ),,( 1 myy L,
（ １ ） 先 得到服从分布 }},,1{),,,|({ 121 NSyxxy m LL?=∈p 的随机变量 （ 记为
1,1+nX ） 的一个样本 1y ( 在用 Von Neumann 取舍原则取样时,对于
),,,(
),,,(),,|(
2
21
21
m
m
m xx
xxyxxy
L
LL
+∞= p
pp,只需假定对固定的 ),,(
2 mxx L,存在 )( 10 yp 及常数 C,满足 )(),,,( 1021 yCpxxy m ≤Lp,并对 ]1,0[U 均匀随机数 U,看 )( ),,,(
0
2
UCp
xxU mLp
是否不小于 U,来决定取舍样本 U ),
（ ２ ） 用同样的方法,再得到服从分布 }},,1{),,,,|({ 2312 NSyxxyy m LL?=∈p 的随
机变量 （ 记为 2,1+nX ） 的一个样本 2y,
（ ３ ） 依此下 去,得到服从分布 }),,,,,,|({ 111 Syxxyyy kmkkk ∈+? LLp
的随机变量 （ 记为 knX,1+ ） 的一个样本 ky,)1(?≤ mk,
（ ４ ） 最后得到服从分布 }}),,,|({ 1 Syyyy mmm ∈Lp 的随机变量 （ 记为
mnX,1+ ） 的一个样本 my,于是 ),,( 1 myyy L= 就是 的一个样 本,
现在我们任取 )0(0 yX =,再按上面方法得到 1X 的一个样本 )1(y,对 n 归纳地我们可用上面方法 分别 得到 nXX,,2 L 的样本 )(,),2( nyy L,当 n 充分大时,由于 Markov
链 nX 的分布近似于 )(?p,我们就可以认为 )(ny 近似地是 分布 )(?p 的一个样本,
［ 注 ］ 上面的 Gibbs 采样也可以看成为如下的算法
（ 1） 按循环次序更新状态的各个分量,在更新了第 d 个分量后,就接着更新第 － 个分量 ；
（ 2） 按条件分布 ),,,,,|( )()( 1)( 1)(1 ndnjnjn xxxx LL +p 选取 )1( +njx,再把状态 ),,( )()(1 ndn xx L 更新为
+++
=),,,( )1()1()1(1 ndnjn xxx LL ),,,,,,( )()( 1)1()( 1)(1 ndnjnjnjn xxxxx LL ++? 。
[注 ] 回顾 Gibbs 采样,还可以从理论角度得到以下认识,如果状态空间是实数,我们把
Gibbs 采样中第 k 个分量的条件更新看成 kx 到 kkk xy x+= 的随机 (加法 )变换,kk xTy
kx
=,
那么 Gibbs 采样的更新就是 说,,只要 kx 服从条件分布 |( kyp 11,,?kyy L,),,1 nk xx L+,则目标函数 (即不变密度 )(xp )在这类随机变换下是不变的,这种看法有助于进一步设计有效
216
的 MCMC 采样方法,
３ Gibbs 采样法与 Metropolis 采样法的变种应用
第 １ 种 变种 （ 随机变种 ）
这时,状态 ),,( )()(1 ndn xx L 的更新,不按指标次序轮番更替,代之以
（ 1） 每次按某个 ＂ 整值 ＂ 随机变量 h 的采样值 j,来随机选取指标 ；
（ 2） 再按条件分布 ),,,,,|( )()( 1)( 1)(1 ndnjnjn xxxx LL +p 选取 )1( +njx ；
（ 3） 然后把 ),,( )()(1 ndn xx L 更新为,
+++
=),,,( )1()1()1(1 ndnjn xxx LL ),,,,,,( )()( 1)1()( 1)(1 ndnjnjnjn xxxxx LL ++? 。
第 ２ 种 变种 （ 分裂变种 ）
设状态变量 x 可以分为两部分 ),( 21 xxx =,其中 2x 是丢失 数据 的部分,分别记这两部分随机变量相互间的条件分布为 )|( 121|2 xxp,)|( 212|1 xxp,如果我们可能在 1x 已知时对 )|( 121|2 xxp 进行采样,也可能在 2x 已知时对 )|( 212|1 xxp 进行采样,那么,我们可以 按下面步骤进行,
(1) 任意设置一个初值,例如
)0(
1x ;
(2) 由条件分布 )|( )(22|1 nx?p 中抽取
)1(
1
+nx;
(3) 再由条件分布 )|( )1(11|2 +? nxp 中抽取
)1(
2
+nx
,这样依次交替进行,
只要 )|( 121|2 xxp,)|( 212|1 xxp 恒正,那么由 上 面 的 Gibbs 采样的一般理论可知,),( )(2)(1 nn xx 就是一个 Markov 链的采样值,而且 ),( 21 xx 的联合分布是这个链的不变分布,于是当 n 充分大的时候,可以认为 ),( )(2)(1 nn xx 是 近似地 采自 ),( 21 xx 的联合分布的一个样值 。
第 ３ 种 变种 （ Metropolis 采样 型 变种 ）
（ 1） 每次按某个随机变量 h 的采样值 j,来随机选取指标 ；
（ 2） 再按条件分布 ),,,,,,|( )()( 1)()( 1)(1 ndnjnjjnjn xxxxxx LL +? ≠?p 选取一个样本 z ；
（ 3） 取一个 与之独立的 ]1,0[U 样本 U,令
≤
= +?
+?
+
其他情形 ）(
)),,,,,|( ),,,,,|(,1min((
)(
)()(
1
)(
1
)(
1
)()(
1
)(
1
)(
1
)(
)1(
n
j
n
d
n
j
n
j
n
n
d
n
j
n
j
nn
j
n
j
x
xxxxz
xxxxxUz
x LL
LL
p
p；
（ 4） 然后把 ),,( )()(1 ndn xx L 更新为
217
+++
=),,,( )1()1()1(1 ndnjn xxx LL ),,,,,,( )()( 1)1()( 1)(1 ndnjnjnjn xxxxx LL ++? 。
４,利用可逆链防止多峰情形的偏取
在 Gibbs 样本生成过程中,只有当马氏链 { nX } 发展的,时间,n 充分大时,它的分布才近似于 p,在实践中怎样才能判断 n 是否已经足够大了呢？ 通常的判断标准一般是看在计算中的 nX 的分布是否比较稳定了,但是,在事实上这种判断方法 不尽正确,下面的示意图给出了 p )( xp= 具有两峰的情形,从不同初值出发的轨道,可能分别稳定在 1E 与 2E
两个区域,而事实上 nX 还远未达到,近似于,不变分布的要求,
E 1 E 2
也就是说,按上一段的方法得到的 Markov 链,从不同的初值出发,在计算达到看起来稳定时,它们还远不能近似 π (x) 的样本,原因在于这时 的 平稳分布 π (x) 有两个峰,且在两峰中间的转移概率非常小,以使从 E 1 中的点出发能到达 E 2 和从 E 2 中的点出发能到达 E 1 的概率都极为微小,从而在一条轨道上要从 E 1 到达 E 2 就要化费十分长的时间,有时甚至在现实生活中可能不会在有意义的时间内达到,发现与使用 Gibbs 样本生成法的先驱们往往希望通过选取多个初始值,得到多条轨道以解决这一问题,但是这并不能完全解决问题,正如在当前的例子中的两条轨道各自只局限在不同的区域 E 1 和 E 2 中,它们的分布,根本反映不出分别限制在 E 1 和 E 2 上的 p )( xp= 如何拼接的,即以多大的比例 1C,2C,使
)( 21 CC + p 1C≈ p
1
|E + 2C p
2
| E,
只有找到了这个比例 1C,2C,才能从分别限制在 E 1 和 E 2 上的 p
1
|E 和 p
2
| E,拼接出整体的 π (x),
为了找 出 比例 1C,2C,我们建议利用可逆 Markov 链来生成 Gibbs 样本,并 从可逆
Markov 链的特性得到不变分布在它的各个峰区上的比例,下面我们以双峰时两条轨道的情
218
形,来阐述确定与估计这个比例的方法,
设分别稳定在 E 1 和 E 2 上两条轨道的样本分布对应于分布 )(1 xp 与 )(2 xp,我们近似地把 )(xp 在 E 1 和 E 2 外为 0,于是由可逆 Markov 链的性质可知,它们应该分别与 { π (x):x
∈ E 1} 及 { π (x):x ∈ E 2 } 成比例,即
),(,)( )()( )( i
i
i Eyx
y
x
y
x ∈?=
p
p
p
p,
于是
)(
)(
^
x
xi
p
p ≈
i
i C
x
x?=
)(
)(
p
p,
因此
)(
1
i
i EC p≈,
具体确定 iC 可 如下 进行,对 ii Ex ∈,利用可逆性我们有
)()(
)()(
122
211
2
1
xx
xx
C
C
pp
pp=
),()(
),()(
1222
2111
xxpx
xxpx
p
p=,
于是对 iEyx ∈,,有
∑∑
∈∈
==
ii Ex
i
Ex
i yxp
xypyx
),(
),()()(1 pp ]
),(
1[)(
1 yXp
Ey yip=,
对于 1X 的 N 次取样 )(1)1(1,,NXX L,我们可以得到
1
1
]),( 1[)(
=
yXpEy yip 的估计
1
1
)(
1
^
),(
11)(?
=
= ∑N
k
ki yXpNyp,
这里的近似只需对从 y 出发的 Markov 链,作独立的 N 次一步转移的模拟,它比 Markov 链的分布近似于不变分布要求的步数要少得多,
余下的问题是从不变分布 π (x) 具体构造可逆 Markov 链,这只需把前面构造的马氏转移
),( yxp 修正为
+=
)(
),()(),(
2
1),(^
x
xypyyxpyxp
p
p,(8,20)
容易检查 它 是以 )(xp 为可逆分布的转移阵,这恰是我们所希望要的,
219
为了便于具体操作,我们还希望 ),(^ yxp 也按 Gibbs 方式转移,即也能归结于计算每次只变动一个分量的 m 个,条件 转移概率,的乘积,为此我们注意,
)(
),()(
x
xypy
p
p
),,(
),,|(),,,|(),,|(),,(
1
11312211
m
mmmmm
xx
xxxyyxxyyxyy
L
LLLLL
p
pppp?=
),,(
),,(
1
1
m
m
xx
yy
L
L
p
p=
dyyyy
yyx
m
m
),,,(
),,,(
2
21
L
L
p
p
∫ dyyyxx
yxx
mm
mm
),,,,(
),,,(
21
11
∫ L
LL
p
p
dyyxx
xx
m
m
),,,(
),,(
11
1
∫
× LLp p
).,|( 21 myyy Lp= ),.,,|( 111 mmm yxxy LLp ).,|( 11?mm xxy Lp,
即 )( ),()( x xypypp 具有类似于 ),( yxp 的形式,不同处只在于所有条件分布中的变量次序倒了一下,而 ),(^ yxp 的含意为,分别以 21 的概率取 ),( yxp 与取 )( ),()( x xypypp,所以,从模拟程序的实现来看,这仅仅增加了一倍计算量,
2,4 MCMC 应用于 Bayes 参数估计
１ 关于 Bayes 统计
我们将在第 ９ 章中较多地介绍与评述 Bayes 统计,Bayes 统计方法的基本想法是,把概率函数中的 未知参数 q 当作随机变量,
由于在 Bayes 方法中不再区分参数与随机变量,所以未知参数的分布的确认是至关重要的,在抽取样本以前,就只能根据先验知识,设置未知参数的分布,称为 先验分布 。 在取样以后,根据对于抽 取 到的样本的概率规律的了解,使用 Beyes 公式,把参数的先验分布改进为 后验分布,后验分布 也称 Bayes 分布 。
设随机变量 x 的分布为 ),( qxp （ 未知参数 q 可以是向量 ）,把 q 看成随机变量 （ 或随机向量 ） 后,),( qxp 就是,在 q 已知的条件下,x 的条件分布,即 ),()|(| qqqx xpxp =,假定 q 的先验分布的分布密度为 )(Jf,于是 ),( qx 的联合分布为 ),()( qqf xp,
在 J 固定条件下,设 nxx,,1 L 为采自 ),( qxp 的独立随机样本,把在 nxx,,1 L 已知的条件下 q 的条件分布,记为 ),,|( 1 nxxJj L,那么,它等于 q 与 ( nxx,,1 L )的联合密度
220
),()(
1
qxqf i
n
i
p∏
=
,除以 ( nxx,,1 L )的边缘密度,即
),,|( 1 nxxJj L
∫ ∏
∏
=
=
n
i
i
n
i
i
dp
p
1
1
),()(
),()(
JJxJf
JxJf
＝,
这就是 q 的 Bayes 分布,所以至关重要的是如何选取参数的先验分布,这个问题将在第 ９ 章中作一些较为浅显的介绍,
2 用 MCMC 方法对 Bayes 分布参数作估计
对于参数 q 的 Bayes 分布,用 MCMC 方法 做 以它为极限分布 的 Markov 链 nJ,于是,
由于 nJ 的分布近似于此 Bayes 分布,我们可以取
^J
N
Nmm ++ ++= JJ L1,(
0nm ≥,0n 充分大 ) 作为参数 J 估计,
３,用 MCMC 方法得到较为复杂的参数估计的大意
较复杂的参数估计的典型情形是 混合总体的参数 Bayes 估计,设样本 nXX,,1 L 来自权重参数未知的混合总体,这些未知权重也就作为另加的待估计参数,在对混合分布作参数估计时,可以引进一些新的随机变量 (在 Bayes 方法中,原则上参数与随机变量不再严格区分,可以把它们看成,补偿参数,) Z,然后用 Gibbs 采样法,交替使用,在? 已知的条件下,Z 的条件分布 )|(1|2?zp,与在 Z 已知的条件下,
的条件分布 )|(2|1 zp? 。 由此可得到 ),( Z? 的 Bayes 分布的采样,其步骤如下,
（ 1） 任给初值 0? ；
（ 2） 按条件分布 )|(1|2 nzp? 采样,得样本 nz ；
（ 3） 按条件分布 )|(2|1 nzp? 采样,得更新样本 1+?n 。
于是只要 )|(1|2?zp 与 )|(2|1 zp? 恒正,那么就有,),( nn z? 是 Markov 链,而且以 ),( Zθ 的联合分布为极限分布 。 从而当 m 大时,可以认为 ),( mm z? 近似地是出自 ),( Zθ 的联合分布的采样,同时可以用取
^J
N
Nmm ++ ++= JJ L1 作为? 的估计量 。
以上的关键在于知道 )|(1|2?zp 与 )|(2|1 zp?,而且易于抽取样本,它们还不能以正概率取 0,以保证所作的 Markov 链的转移函数的收敛性,关于 )|(1|2?zp 与 )|(2|1 zp? 的取法与获得,在现代 一些 统计论文中,有一些典型的做法,本书不再叙述,
221
４ MCMC 方法用于缺失数据情形下的参数估计的大意
如果 ),( YZX = 的分布密度有未知参数?,),,(?yzp,其中 Y 是可以通过采样观测到的,但是,Z 却观测不到 （ 或者是由数据的缺失而引起的,它称为缺失数据 ）,用 Bayes 统计的思想,在观测样本 Y 条件分布给定的条件下,记条件分布为
)|(1|2?zp ),,(?= Yzp,)|(2|1 zp?,
只要它们恒正,而且有一种方法可以近似地取得服从这两个分布的样本,那么就可以用第 ３ 段中的方法得到 缺失数据 Z 的估计,
3,模拟退火 ( simulated annealing )
３,１ 模拟退火方法的基本想法
一般的优化问题,都可以归结为求某个目标函数 f 在一定范围内的最小值,但是许多常用的方法,往往都是运用某种局部的比较去决定优化过程的发展,以逐步达到最小值,例如梯度法,就是每次将目标函数在当时所考察的点与附近的点的取值加以比较,根据梯度方向就是目标函数值下降最快的方向,导出了梯度下 降法,即在被考察的点沿函数下降最快的方向移动,希望经过多次这样的移动,最终达到目标函数 f 最小的点,但是这个算法的实际结果是停止在梯度为零的某个点上,这就是说,在计算进行到到达目标函数的某个局部极小值点时,在考察点的目标函数值就小于它附近所有点的目标函数值,于是搜索过程就停止了,我们知道目标函数的任意一个局部极小值点都满足上述条件,然而局部极小值点处的函数值,有可能比函数的实际最小值大得多,于是求得最小值的努力就完全寄托于初值的选取,
为了克服上述算法的这个缺 点,Metropolis 等人在 1953 年提出的一种有效地模拟固体达到热平衡发展的算法,而 Kh asi ’minski 在 1965 年对决定性的梯度算法,引入小的随机噪声干扰,使之在达到平衡时,以高概率进入目标函数的最小值附近,但是当时没有受到重视,Kirkpatrick 等人于 1983 年对组合优化问题,重新提出模拟退火算法,其基本想法正是 Kh asi ’minski 的构想,即引入人工噪声使得当算法达到局部极小值时,能有一个小概率
,逸出,局部极小值的陷井,这个想法源自冶金退火过程的物理思想,即将金属加温至充分高,再使其慢慢地退火冷却,在加温时,金属内部粒子随温度增加变为无序,能量增大 ; 而在退火时,在慢慢冷却的过程中粒子渐趋于有序,在每个温度时都达到平衡态,最后在常温时达到能量最大的基态,而用退火模拟优化问题时,则将能量模拟为目标函数,温度演化成代表噪声水平的控制参数 b,这样就得到优化问题的模拟退火算法,
简单而实用的模拟退火模型是有限状态模型,典型情形是组合优化问题,这时假定自变量只能取自有限个值,而加入随机噪声的算法,就是把原来的计算过程变成随机过程 — 也就是一个 Markov 链,并使得此 Markov 链是互通的,非周期的,正常返的,且具有如下形式的不变分布,
222
))((
))((
)( *
*
1 fif
j
fjfi ee
∑=
b
b
bp )(
)(
1 if
j
jf ee
∑=
b
b,（ ８,20 ）
其中 f 是我们要优化的目标函数,)(min* iff i=,而 b 是一个控制噪声水平 ( 随机性水平 )
的参数,这个分布有一个重要的性质,当 ∞→b 时,其极限分布集中在目标函数 f 取最小的集合上,即在 })(:{ *fifi = 上 ( f 取最小值的点可能不唯一 ),于是我们可以利用这个
Markov 链的发展,使其到达不变分布,然后令 β →＋ ∞,以达到 f 取最小值的点集上,而固定 b 时的 Markov 链的发展的取样 过程,则可以通过上一节所叙述的各种 M arkov M onte
C arlo 方法来实现,例如可用 Gibbs 采样法,Metropolis 采样法,或条件概率的采样法,
(8,20) 中的不变分布的形式,很像在物理上的 Boltzmann 分布,而 b 则是倒温度,即 b ＝
T
1 （Ｔ 是绝对温度 ）,若把无噪声 ( 即 b ＝＋ ∞ 情形 ) 的决定性算法,看成是绝对 零度 ( Ｔ
＝ 01 =b ) 的情形,那么,β ＜ ∞ 的情形就是将温度升为正温度,于是上面所述的达到目标函数取最小值的集合的过程,即对固定的 b,构造 Markov 链的样本,再让 ∞→b 的过程,
恰似先将温度由绝对零度升高 ( 即 b ＜ ∞ ),然后再让温度降回到绝对零度 ( 即令 b = ＋ ∞ )
的过程,这个升温再降温的过程,就是,退火,过程的模拟,因而这个算法称为 模拟退火算法,或称 随机松驰法,简言之,即 升温以得 到 不变分布,再降温 ( b →＋ ∞ ) 使不变分布集中于目标函数 f 的取最小值的集合上,
为了加快计算速度,人们在理论上发现,可以采取受人工噪声干扰的非时齐的 Markov
链,把 Markov 链的时间发展与降温这两步,归并成一步作模拟退火,即 使所选的 Markov
链的第 n 步转移矩阵随 nb 而改变,当时间 n 趋于无穷时,只要 nb (称为 冷却进度 （ Cooling
Schedule） ) 以适当缓慢的速度趋于 ∞,则此非时齐的 Markov 链,就依概率收敛到目标函数达到最小值的点集上,这个非时齐的 退火的过程完全由冷却进度所控制,如果 nb 减小过快,则在直观上就相当于在每个温度段还未达到平衡,退火过程就还会陷入能量函数的局部极小 ; 反之,如果 nb 减小过慢,相当于几乎是常值,退火过程就会按平衡态 (即不变分布 )发展,这两种情形都不能达到整体极小,
更确切地说,存在一个由目标函数 f 确定的正常数 a,只要取
nn log?= gb ( )g < a,(8,21)
223
( 即取温度
n
T b1= ),则当 n →∞ 的进程就达到了,适当缓慢地降低温度,这个要求,当然 g 取得越接近 a 效果越好,计算接近目标函数取最小值的点集的进程越快,
在实际计算时,当 n 大时,按 (8,21) 式取的 nb 的变化非常小,以致在实际上几乎看不出变化,所以有时人们就简便地使用一个固定的很小的 b,构造具有 （ ８,２０ ） 式的不变分布的 Markov 链,用其发展以代替模拟退火,显然,这种 Markov 发展的终极实际上是一个平衡分布,但是只要这个平衡分布高度集中于目标函数全局极小的集合附近,那么这样的计算误差也就可以容忍了,在具体操作时,到底如何掌握分寸,是一个实际经验问题,至今还缺少中肯的理论分析,
在许多情况下,由于所处理的系统非常复杂,或者实际维数太高,还有的时候目标函数没有明确的分析表达式,这时人们就引入局部改变的邻域,记状态 i 的邻域为 )(iN （ 这种邻域是根据问题的要求自然地确定的 ）,作矩阵
∈==?
))((0
))(()( ),( 1
,
~
,
~~
iNj
iNjiNppP
jiji,（ ８,２２ ）
它规定了只 许在近邻之间转移,即当系统处于状态 i 时,它下一时刻只能转移为 )(iN 中的状态,再规定所要的 Markov 链第 n 步转移概率矩阵由如下的 Metropolis 采样,
P )(n ))(( nijp b
D
=,（ ８,２３ ）
=
≠
=
∑ ≠ ijp
ijp
p
i
k
ikik
i
j
ij
ji
},1min{1
},1min{
)(
)(
)(~
)(
)(~
,
b
b
b
b
p
p
p
p
b （ ８,２４ ）
所给定,那么,运 转这一族转移矩阵 { P )(n }所决定的非时齐 Markov 链,只要温度
n
nT b
1=
以合适的速度下降,就收敛到使 f 取最小值的点集,
上面指出的算法的局部化思想,是模拟退火算法中为了克服因为状态太多而引起的计算爆炸的有力手段,这里的矩阵
~P
可以一般化地取一个互通的概率转移矩阵,显见它越稀疏 ( 即其中 0 越多 ) 则进程 可能 越快,它称为 局部化预选矩阵,简称 预选矩阵,它相当于规定了可以转移的,邻域,,把整体搜索变为局部搜索,以减少计算量,在研究第 9 章中的
Gibbs 场的 MCMC 及模拟退火时,就以邻域转移代替这个预选矩阵,
3.2 有关模拟退火算法的非时齐马氏链的理论背景
224
我们的仍旧假定状态集 ( 或称为组态空间 ) 是有限集,即 },,2,1{ MS L=,设 f 为 Ｓ 上的一个实值函数,令 L(S) 为 S 上 全体实值函数的集合,对于任意 ))(,),1(( Mfff L=,
))(,),1(( Mggg L= ∈L(S),定义
范数 ∑
=
=
M
i
iff
1
|)(|||||
D
,距离 ||||),( gfgfd?=,
那么 L(S) 就成为完备的距离空间,显见 S 上的概率分布都是 L(S) 中的元素,而且 S 上的任意概率分布的序列在此距离空间中的极限仍 是一个概率分布,( 但是,如果 S 不是有限集而是具有可数个元素的集合时,
此结论不真 ),
由定义我们立即推得,对于任意一个转移矩阵 P,恒有 f|| P ||||| f≤,其中 f P指 f 与 P作为矩阵相乘,
(1) 时齐的 Markov 链优化的 模拟退火 方法
定理 ８,２ 设以 P ))(()( bb ijp= 为转移矩阵的 Markov 链记为 }{ bx n,则
1)(limlim * ==∞→∞→ ffP
nn bxb
,
证明
==∞→∞→ )(limlim *ffP
nn bxb
}))(:{(limlim *fifiP nn =∈∞→∞→ bb x
1lim )(
})(:{ *
== ∑
=
∞→
b
b p i
fifi
,
( 2 ) 非时齐 Markov 链的模 拟退火
非时齐 Markov 链的分布的收敛性有
定理 ８,３ （ Dobrushin - Isaacson - Madsen 定理 ） 设 nx 是 状态空间 S 上的非时齐 Markov 链,其第 n 步转移矩阵为 P )(n,如果它们满足,
(A.1) P )(n 作为时齐的转移矩阵时有不变分布 )(np ;
(A.2) ∑ ∞<? +
n
nn |||| )1()( pp ;
(A.3) 或者满足 Isaacson - Madsen 条件,对于任意概率分布向量 nm,,及正整数 j,恒有
P)(|| nm? L)j( P 0|| n)n( →? ∞→ ;
或者满足 Dobrushin 条件,对于任意 j,收缩系数满足
(C )()( nj PP L ) 0 →? ∞→n,
那么,
(1) )()( lim nn pp ∞→?∞ = 存在 ;
225
(2) 无论什么初分布 0m 下,此非时齐的 Mar kov 链 }{ nx 在时刻 n 的绝对分布 nm 恒有极限
)(∞→ pm
n,
证明 首先,由 (A.2) 利用三角形不等式立得 }{ )(np 是 L(S) 中的 Cauchy 列,故而 (1) 正确,
往证 (2),用三角形不等式与 Dobrushin 不等式,我们得到
|||| )()()1(0 ∞?pm nPP L
||)(|| )()()()()1(0 njn PPPP LL ∞?≤ pm |||| )()()()( ∞∞?+ pp nj PP L
≤2 (C )()( nj PP L ) |||| )()()()( ∞∞?+ pp nj PP L,
在满足 Isaacson - Madsen 条件时,上式的第一项趋于 0,( 而由 Do brushin 不等式推出
||)(|| )()()()()1(0 njn PPPP LL ∞? pm ≤2 (C )()( nj PP L ),
即在 Dobrushin 条件成立时,Isaacson - Madsen 条件必然成立 ),对第二项利用不变分布条件
)()()( mmm Ppp = 推得对于 Mj ≥ 有
|||| )()()()( ∞∞?pp nj PP L
)()()()( )(|| njj PP Lpp?= ∞ ||)( )()()()( ∞?+ pp njj PP L
)()()()( )(|| njj PP Lpp?= ∞ ||)( )()()1()( ∞+?+ pp njj PP L
)()()()( )(|| njj PP Lpp?= ∞ ||)()( )()()()()()1(
1
∞++?+
=
+?+ ∑ pppp nnkjkjkj
jn
k
PP L
||)(||sup )()( nMn pp?≤ ∞≥ ||||sup|||| )()()1()( ∞≥+
≥
+?+ ∑ pppp nMnnn
Mn
0→,
[ 注 ] 在乘积 )()( nj PP L Mkinjikp ≤
=,),( )( 中,)|( 1)1,( ikPp jnnjik === ++ xx 为第 n 次转移矩阵为 P )(n 的非时齐 Markov 链 }{ nx 的多步转移概率,由于我们研究的是有限状态的情形,于是
Isaacson - Madsen 条件可以改写为如下的等价条件 ( 称为 弱遍历条件 ),
(A,2) ’ 对于任意状态 Skli ∈,,,任意正整数 j,恒有
0)(lim ),(),( =?∞→ njlknjikn pp,
对于模拟退火的收敛性,利用较为深入的,转移概率为指数型衰减的时齐的 Markov 链序列的大偏差理论 ( 粗略地说,它研究此链越出给定小概率状态集合的概率的指数型衰减 ),可以证明 ( 这需要相当深入的 叙述与 证 明,本书无法列入 ) 下面的定理 8.4,
设 )(if 为定义在 },,1{ MS L= 上的函数,令
226
p ),,( )()(1)( nMnn pp L=,)(nip )()(1 if
j
jf
n
n
ee∑= bb,（ ８,２５ ）
P )(n )( )(nijp
D
=,
=
≠
=
∑ ≠ ijp
ijp
p
n
i
n
k
ikik
n
i
n
j
ij
n
ij
},1min{1
},1min{
)(
)(~
)(
)(~
)(
p
p
p
p
,（ ８,２６ ）
定理 ８,４ （ G emam - Gemam ）
假定
(1) ∞↑< nb0,满足,存在 0n,使 0nn ≥ 后,有
nMn ln1?≤b,（ ８,２７ ）
其中 D 是目标函数 f 的振幅,
)(min)(max ifif ii?=DD ; （ ８,２８ ）
(2) 预选矩阵
~P )( ~
ijp= 是互通的概率转移矩阵,
那么,以 P )(n 为第 n 步转移矩阵列的非时齐的 Markov 链 }{ nx 满足弱遍历条件,因而 (A,2) 满足,
定理 ８,５ （ 模拟退火的收敛性定理 ）
在定理 8.4 的条件满足时,无论什么初 始分布 0m 下,以 P )(n 为第 n 步转移矩阵列的非时齐的 Markov
链 }{ nx,在时刻 n 的绝对分布 nm 有极限 )(∞→ pmn,其中 )(∞p 为集合 min})(:{ =ifi 上的离散均匀分布,从而有
1)( * = →? ∞→ ffP nnx,
证明 我们来验证 Dobrushin - Isaa cson - Madsen 定理的条件成立,由 G emam - Gemam 定理,我们只需验证条件 (A,1),为此,首先注意当 n 充分大时 ( 例如 0nn > 时 ),利用 nb 的单调不减性,我们有,对于固定的 i,)(nip )()(1 if
j
jf
n
n
ee∑= bb 在 0)( >if 时,关于 n 是单调不减的 ( 对 nb 的微商非负 );
而在时 0)( ≤if 时,关于 n 是单调不增的 ( 对 nb 的微商非正 ),再注意 ( 这是在第 5 章 Dobrushin 不等式中应用过的事实 )
|||| )()1( nn pp?+ ∑?= +
i
n
i
n
i ||||
)()1( pp ∑ ++?=
i
n
i
n
i )(
)()1( pp ∑?+?+
i
n
i
n
i )(
)()1( pp
∑ ++?=
i
n
i
n
i )(2
)()1( pp ∑?+?=
i
n
i
n
i )(2
)()1( pp,
227
于是我们得到
∑
=
N
nn 0
|||| )()1( nn pp?+ =2 ∑
=
N
nn 0
∑ ++?
i
n
i
n
i )(
)()1( pp ∑∑
=
++
>
=
N
nn
n
i
n
i
if 0
)(2 )()1(
0)(
pp
∑∑
=
+
>
=
N
nn
n
i
n
i
if 0
)(2 )()1(
0)(
pp ][2 )()1(
0)(
0n
i
N
i
if
pp?= +
>
∑ 2≤,
令 ∞→N,便得 ∑
∞
= 0nn
|||| )()1( nn pp?+ 2≤,所以 (A,1) 成立,定理证毕,
习题 ８
1,证明对于 Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 采取的操作得到的由 i 到 j 的转移,就 是
（ ８,１５ ） 式 定义 的转移概率,
2,证明由 （ ８,１６ ） 给出的转移矩阵 所对应的 Markov 链以 p 为可逆分布,再给出这种 Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作步骤,
3,证明由 Hastings 推广了的 Metropolis 采样所对应的 Markov 链以 p 为可逆分布,再给出这种 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作步骤,
龚光鲁,钱敏平著 应用随机过程教程及其在 算法与智能计算中的应用
清华大学出版社,2003
第 8 章 Monte Carlo 与 Markov Chain Monte Carlo (MCMC) 方法
在许多很复杂的统计问题中,有时很难直接对各种统计方法进行理论分析,为了评估它们的优劣,常见的实用办法是做随机模拟,即设法按问题的要求与条件去构造出一系列的模拟样本,用它们的样本频率代替对应的概率作统计分析与推断,观察由这些摸拟样品所作出的推断的正确率,因为在概率论初期发展时,随机模拟的原型常常来自博采,于是 人们就以博采之都 M onte Carlo 作为随机模拟方法的别称,久而久之,Monte Carlo 方法作为名称倒比随机模拟方法更为广泛地被常用了,相仿地,人们还把组合计算中的某些随机模拟方法,
称为 Las Vegas 方法,这是以美国的博采城 Las Vegas 命名的,
一般 随机 模拟 方法 的长处是,计算的复杂度不依赖于计算空间的维数,因此,在计算非常 高 维数的积分 （ 或多指标的求和 ） 时,Monte Carlo 方法比通常的计算方法有明显的优越性,而很多的 Monte Carlo 计算问 题,可归结为计算积分或庞大和数的问题,
1,计算积分的 Monte Carlo 方法与采样量估计
通过构造独立同分布随机数,计算积分的 Monte Carlo 方法,也称为 静态 Monte Carlo
方法,其思想可以在本节中,通过估计最简单的积分 ∫
b
a
dxxf )( 得到阐明,对于高维积分,
其思路与一维积分是一样的,另一方面,在甚高维情形,因为计算量太大,用静态 Monte
Carlo 方法处理速度太慢,我们在第 2 节中通过构造 Markov 链的极限不变分布,来 模拟计算积分的方法 （ 称为 动态 Monte Carlo 方法,Markov 链 Monte Carlo 方法,M C MC),将更为适用,
1,1 用频率估计概率来计算积分的 Monte Carlo 方法
假定 Mxf ≤≤ )(0,那么由积分的面积含义有
||)( Sdxxf
b
a
=∫,(其中 || S 为 })(,:),{( yxfbxayxS ≥≤≤=? 的面积 ),
考虑平面区域 ],0[],[ Mba ×= 上的均匀随机变量 x,则
MabSPp )(
1)(
=∈=
x
∫
b
a
dxxf )(,
对于 N 个独立的? -均匀随机数 )(,Nii ≤x,记
},,{ 1 NSN xx L= 中落在 S 中的频数,
于是,利用大数定律便知
=^I
N
NMab S)(? （ ８,１ ）
203
是积分 ∫
=
b
a
dxxfI )( 的相合估计,即对于任意的 0>e,当 ∞→N 时,有
∫ →>
b
a
S dxxf
N
NMabP 0)|)()((| e,
又由于 SN 服从参数为 ),( Np 的二项分布,所以我们有 =^IE ∫
b
a
dxxf )(,即
^I
是积分
∫
b
a
dxxf )( 的无偏估计,此估计的方差为
=^IVar N ppMab )1()( 22 N Mab
I
Mab
I
Mab
))(1()(
)( 22?
=
IIMabN ])[(1= )1( NO=,
又因为方差代表平均平方误差,所以我们就说,积分的估计
^I
的误差为 )1( NO,
频率法所需采样量的估计
我们知道统计估计是以区间估计来给出概率误差的,即对于允许误差 e,及允许以 d 的失败概率,来确定样本量 N (依赖于 de,),使
∫ <>?
b
a
dxxfIP de)|)((|
^
,
1 用 Chebyshev 不等式估计样本量
最粗略的确定样 本量 N 的方法,是通过 Chebyshev 不等式来估计,即由
2
^
^^^
)|(|)|)((| eee IVarIEIPdxxfIP
b
a
≤>?=>? ∫,
要求 dee <
= 2
22
2
^ )1()(
N
ppMab
IVar,考虑到
4
1)1( ≤? pp,我们要确定 N,只需满足 de <
2
22
4
1)(
NMab 就够了,也就是说只要
2
22
4
)(
de
MabN?> 即可,但是这个采样量的估计过于粗糙,
2 用中心极限 定理的近似估计样本量
在 N 大时 (这总假定满足 ),用中心极限定理也可以得到采样量 N 的估计,中心极限定理断言,MabI )(^?= NNS 的近似分布为正态分布 (N ∫
b
a
dxxf )( ),^IVar,令 df 满足
204
∫
∞?
=
df
dp dxe
x
2
2
2
1,( 即 ∫
=2
2
2
121 2
d
d
f
f
dp dxe
x
),
于是有
∫ ≤≥?
b
a
IVardxxfIP df d )|)((|
^
2
^
,
选取采样量 N 使 ef d ≤
^
2
IVar,即 2
2
2
22 )1()(
df
e≤
N
ppMab,则就有
∫ ≤≥?
b
a
dxxfIP de)|)((|
^
,
所以只需取
2
2
2
22
4
)(
dfe
MabN?>,
3 采样量的估计的比较
我们把用 Chebyshev 不 等式得到的采样量的估计与用中心极限定理得到的采样量的估计作比较,可以看到只要 12
2
<<ddf,用中心极限定理得到的采样量的估计,要比用
Chebyshev 不等式得到的采样量的估计小得多,
［ 注 ］ 还可以应用概率理论中的大偏差速率函数来求得采样量的估计,可以 证明由大偏差速率函数能得到采样量的估计
22
ln
e
d?>N,
当 d 小时,dln? 趋于 ∞ 的速度要比 2
2
df 趋于 ∞ 的速度快,所以,以用中心极限定理得到的采样量的估计为最好,
１,2 用样本函数的平均值估计积分的 Monte Carlo 方法 ― 期望法
期望法的核心思想是把积分看成某个随机变量的期望,最自然的是看成 ],[ ba 上均匀
随机变量的期望,设 ],[~ baUh ( ],[ ba 上的 均匀分布 ),则
)()()( hEfabdxxfI
b
a
==∫
,(8,2)
于是对于 N 个独立的 ],[ ba 上的均匀随机数,可以用矩估计
N
ffabI N )()()( 1^^ hh ++?= L (8,3)
205
作为 ∫=
b
a
dxxfI )( 的估计,显然它也是无偏估计,而且
N
fVarabIVar )]([)()( 2^^ h?= [ ]222 )]([)(1)( hh EfEf
Nab=
222 1)(1)( I
Nab
dxxf
Nab
b
a
= ∫ 2)[(1 IMIabN≤ )( ^IVar=,
可见期望法比频率法更有效,但是期望法多算了 N 次函数值,这是需要花费计算时间的,
所以在全面考虑有效性与计算时间的得失时,有时人们也愿意采用频率法,
1,3 减少方差的技术
用 Monte Carlo 方法计算积分 ∫
b
a
dxxf )( 时,未必一定要使用均匀随机数,事实上,从
],[ ba 上取值的任意一种随机数出发,都可以得到 ∫
b
a
dxxf )( 的估计量,而且在 0)( ≥xf 时
显见,值 )(xf 大的 x 对于积分 ∫
b
a
dxxf )( 有更大的贡献,由此得到直观启发,所用的随机数的分布密度的形状越象 )(xf,则 就 越合理,这个思想就是 重要度采样法,
1,?g 采样法
假定分布密度 )(xg 在 )(xf 非零处恒正,则积分 =I ∫
b
a
dxxf )( ∫=
b
a
dxxgxg xf )()( )(,对于密度为 )(xg 的,?g 随机数,V,我们有
])( )([ VVgfEI =,(8,4)
于是对于 N 个独立的?g 随机数 NVV,,1 L,关于积分 =I ∫
b
a
dxxf )( 可取估计
])( )([1
1
1
)(^
V
V
g
f
NI
g?
= ])( )(
N
N
g
f
V
V++L,(8,5)
显见它也是无偏的相合估计,利用 Schwartz 不等式得到
∫ ∫=
b
a
b
a
dxxgdxxgxg xf )()(])( )([ 2 ∫
b
a
dxxgxg xf )(])( )([ 2
206
≥
2
])()( )([∫
b
a
dxxgxg xf =
2
])([∫
b
a
dxxf,
而且上式当且仅当在 )( )()( xg xfcxg = 时 (即 )()( xcfxg = 时 ) 达到极小值
2
])([∫
b
a
dxxf,这说明 NIVar
g 1
)(
)(^
= ])(])( )([[ 22 Idxxgxg xf
b
a
∫ 的最小值在 )()(0 xcfxg?=
时取到,又因为 )(0 xg 为密度,故
∫
= b
a
dxxf
c
)(
1,对应的 0)( )(^ 0 =gIVar,此时
∫==
b
a
g
dxxfcI )(1
)(^ 0
是无估计误差的精确值,
综上讨论可知,要使方差达到最小,就应该用 )()(0 xcfxg?= 作为参考密度,这时用 Von
Neumann 取舍方法必须知道参考密度 )(0 xg,而计算 )(0 xg 又必须知道常数 c,后者却是于积分的倒数,这样的耦合性的存在使这种方法不能真正用于计算积分,但是由此我们至少得到了以下的认识,即 只要密度 )(xg 的形状与被积分的函数近似,用
)(^ g
I 作为 ∫
b
a
dxxf )(
的估计,就会降低方差,这个想法还可以推广到 )(xf 未必为非负函数的情形,这就是下面的概念,
定义 ８,１ 分布密度为
∫
= b
a
dxxf
xfxg
|)(|
|)(|)( 的?g 采样,称为 (关于 )(xf 的 ) 重要度采样 （ Importance Sampling),
重要度采样不能直接通过取舍原则实现,近似地实现重要度采样可以采用下一节中的
Markov 链 Monte Carlo 方法,这个方法同时给出连积分 ∫
b
a
dxxf |)(｜ 的近似算法,
在实践中人们也往往按照重要度采样的思路,灵活地寻找常用的已知类型的密度 g,使它在峰值附近与 | )(xf | 较接近,以便达到降低估计的方差的目的,
修正的重要度采样 (Modified Importance Sampling)
对于?g 采样,假定有非负函数 )(xh,满足
207
∫ >=
b
a
dxxgxh 0)()(a,
而且 a 已知,那么我们可以采用 )(xh 作为修正乘积因子,显见对于?g 随机数 V 有
)(
])( )([
)( VV
Va
Eh
g
fE
dxxfI
b
a
== ∫,(8,6)
如果我们先放弃对于估计的无偏性要求,而只要求估计的相合性,则对于 N 个独立的?g
随机数 NVV,,1 L,我们可以通过比值,构造 ∫=
b
a
dxxfI )( 的如下的估计量
)()(
)(
)(
)(
)(
1
1
1
N
N
N
hh
g
f
g
f
I VV V
V
V
V
a ++
++
=?∨ L
L
,(8,7)
显见它是 ∫=
b
a
dxxfI )( 的相合估计,在某些假定下,它还是渐进无偏的,即
IIEN =∨∞→ )(lim,
而且
∨I
保留了重要度采样的特性,即当 )( )()( xg xfcxh = 时,
∨I
就是 ∫=
b
a
dxxfI )(,于是,
只要当 )(xh 与 )( )(xg xf 近似,就会降低方差,注意对于给定的分布密度 )(xg,)(xh 应选取得尽量与 )( )(xg xf 近似,这个修正乘积因子 )(xh 是用来再一次降低由于密度 )(xg 与被积函数
)(xf 的倍数不够像所带来的失误而设置的,而当 1)( ≡xh 时,就退化为是?g 采样,这相当于对?g 采样不再作修正,
2,Markov 链 Monte Carlo (MCMC)
用 Markov 链的样本,来对不变分布,Gibbs 分布,Gibbs 场,高维分布,或样本空间非常大的离散分布等作采样,并用以作随机模拟的方法,统称为 Markov Chain Monte Carlo
(M C MC) 方法,这是 动态的 Monte Carlo 方法,由于这种方法的问世,使随机模拟在很多领域的计算中,相对于决定性算法,显示出它的巨大的优越性,而有时随机模拟与决定性算法的结合使用,会显出更多的长处,
MCMC 至少可以用在以下几个层面,
208
1,用于生成较复杂的随机数,
* 实现对高维分布 （ 或高维格点分布 ） p 取样,得到 p 随机数,
* 是实现重要度采样的一种方法,对 |)(| xf 的重要性采样,就是取得
p ( ∫
== dyyf xfx |)(| |)(|)(p ) 随机数,
2,实现高维积分 （ 或项数极多的求 和 ） 的数值计算 （ 典例是 Gibbs 分布的各种泛函的平均值的计算 ）,对于 0)( ≥xf,作以 p ∫
== dyyf xfx )( )()(p 为极限分布的 Markov 链 nX,
利用遍历定理可以由这个 Markov 链的一条轨道,得到分布密度 )(xp 的估计,记为 )(^ xp,再用
)(
)(
^
x
xf
p
作为 积分 ∫= dxxfI )( 的估计,
3,用模拟方法估计最可几轨道,例如,如果模拟了 100 条轨道,那么就能以大概率推断,最可几轨道就在这些轨道的邻近,当统计量的分布未知时,可以用模拟方法从频率估计置信限,
４,用被估参数的 Bayes 分布 （ 参见第 9 章 ） 的取样,来估计参数,
５,求复杂样本空间上函数的极值 ( 模拟退火 ),
2,1 Gibbs 采样法 ( Gibbs sampler )
1,用 MCMC 方法得到 Gibbs 分布的样本与估计 Gibbs 分布
在第 ６ 章中,我们考虑了在 d 维的?N 格点集上的 Ising 模型的 Gibbs 分布
∑
∈
=
S
H
H
e
e
h
hb
xb
xp )(
)(
,由于所涉及的状态空间 (全体组态的集合 ) dNS ),,1{}1,1{ L?= 非常大 (例如,把一幅 256× 256 个采样点的黑白图像看成一个组态,则 256,2 == Nd,S 中有
)1022(2 18000600002256256 16 >>=× 个元素 ),这就使得分母中的求和无法实际完成,而 MCMC
方法就是以通过构造一个以这个 Gibbs 分布 ∑
∈
=
S
H
H
e
e
h
hb
xb
xp )(
)(
为不变分布的离散时间的
Markov 链 nX (它就是 Glauber 动力学中的连续时间的 Markov 链的离散时间采样 ),作为模拟计算的基点的,构造的 Markov 链必须易于计算,所以我们要求它的概率转移速率只容许在组态的一个格点上变动,这样的变动方式,称为 Gibbs 方式,这种抽样方法称为 Gibbs 采样法,或者 Gibbs 样本生成法,这个 Markov 链的不变分布正是此 Gibbs 分布 p,我们还要求此 Markov 链的转移矩阵满足 P →? ∞→nn pT1,这就是说,要求 Gibbs 分布是 Glauber
动力学的极限分布,
于是,当 n 大时,nX 的一个样本可以近似地认为取自 Gibbs 分布的一个样本,即按此
209
Markov 链沿任意一条轨道充分发展,就得到 Gibbs 分布的近似取样,
再则,Gibbs 分布的归一化常数 (称为 配分常数 ) )(xb
h
H
S
e?
∈
∑,是一个巨大的求和,即一个,离散的,积分,用随机模拟法计算这个 ＂ 离散的 ＂ 积分的最佳随机数正服从 Gibbs 分布
(即重要度采样 ),对于 Gibbs 分布的取样,用通常的取舍原则常常并不可行,例如,分别取
1=C,)()( xbx Heh?=,而参考密度 )(0 xp 为组态空间上的均匀分布,这时 )(xbHe? 的值常常小得超出计算精度,而求和变量 x 的范围是庞大的组态空间,这就导致求和无法实际计算,所以需要用 Markov 链 Monte Carlo 方法,用 MCMC 方法生成了以 Gibbs 分布为极限分布的 Markov 链 nX 以后,由遍历定理用 Markov 链的一条轨道,可给出极限分布 （ Gibbs 分布 ） 的估计,对于充分大的 N,可令
))()((1 2}{1}{^ NN XIXIN xxxp L+= +,（ ８,８ ）
再用 )(xbHe? 除以 Gibbs 分布在 x 处的估计值 ^
)(
x
xb
p
He?
,作为配分函数的估计,在理论上这个估计应该与 x 的取法无关,但是,在实际计算中对多个不同的组态 ix 分别估计此和数后,
再作平均常常能降低方差,
在第 6 章中,我们曾给出了用 Glauber 动力学构造的两个不同的连续时间的 Markov 链
(对应于两个不同的转移概率速率矩阵 Q),它们都以 Gibbs 分布 p 为极限分布,而且都是可逆的,
较为深入的理论研究表明,使用不可逆的,且以 p 为不变分布的 Markov 链作 Markov
Monte Carlo,会加快这个极限的收敛速度,然而,在另一方面这种做法又会增加计算的复杂程度,再则,为减少估计的方差而作的努力也常会增加计算时间,这就是说,在计算中,我们会面临难以两全的抉择,在实际中如何采取折衷,既要看问题的性质,又要参考实践的经验,没有统一的原则,
用以完成 MCMC 采样操作的 Markov 链,可构造如下,
在第 6 章中,对于 d 维有限格点上,由具有两个自旋的组态空间上的能量函数
∑ ∑=
相邻yx x
xhyxH
,
)()()(21)( hhhh,(8,9)
可构造 如下的转移概率速率
≠ == 的其它情形 ）xh xhxxh (0 )(),,(
xxC
q (8,10)
( ),( xxC 的两种取法各为,
))()((),( xxbx HHexC= (8,11)
或
210
))()((1
1),(
xxbx HH xexC+= ),(8,11) ’
它们决定的 Glauber 动力学,分别对应的两个不同的连续时间的 Markov链,它们都 以 )(hH
为能量函数 的 Gibbs 分布 p 为可逆不变分布,而且 p 还 是转移矩阵的 极限分布,P?→? ∞→tt)( pT1=,如果考虑间隔为 t? 的采样时间,其中 t? 充分小,我们还可以进行如下的近似 计算,假定初始组态为 V,在时刻 t? 它以概率 txC?),( V 转移到组态 xV,而以概率 ∑
∈
dNx
xCt
},,1{
),(1
L
V 停留在原来的组态,这就近似地得到了 Markov 链在 t? 时刻所处的组态 1V,再以类似的方式得到 Markov 链在 t?2 时刻所处的组态 2V,依次下去,当采样数 n 充分大时,组态 h 在这段 },,{ 1 nVV L 中出现的频率,就用作 hp 的估计
^
hp,
Gibbs 分布的一些统计平均量的模拟近似
对于上面定义的,以 Gibbs 分布为极限的 Markov 链,我们用它的 一段样本 NVV,,1 L,
可以给出如下的 Gibbs 统计平均量
∑
∑
∈
∈=
S
H
H
S
e
ef
F
x
xb
xb
x
x
)(
)()(
(8,12)
的模拟近似
n
f
F
n
i
i∑
== 1^
)(V
,(8,13)
它是 F 的无偏估计,不难证明,存在数 ∞<),,( pPfV,使
^F
的方差满足
)1(),,()( ^ nonPfVFVar += p nPC f ))(1( ||||13
2
≤,(8,14)
其中 ∑=
x
x )(|||| ff,而 )(PC 为此 Markov 链的转移矩阵 P的 Dobrushin 收缩 数,
Gibbs 采样法还可以有效地用于 Bayes 方法中的后验密度的采样的模拟计算 （ 见第 ９ 章 ）,
2,2 M etr opolis 采样法 ( Metropolis sampler )
１,Gibbs 分布的采样的 Markov 链 Metropolis 方法
为了突出主要思想,下面我们把组态空间 ( 状态空间 ) dNS ),,1{}1,1{ L?= 简单地记为
},,1{ KL,对于组态空间上给定的分布 p,Metropolis 构造了以
211
)(),,1min(
~
ijpp
i
j
ijij ≠= p
p
(8,15)
( ∑
≠
=
ij
ijii pp 1 ) 为转移的时齐 Markov 链,其中
~P )( ~
ijp= 是一个对称的互通 转移矩阵,
称为 预选矩阵,或访问方案,使用它是为了减少 或简化 状态间的连接,以加快 Markov 链的分布向不变分布收敛的速度,显见 p 是它的可逆分布 （ 注意在 Gibbs 分布情形,状态 ji,体现为组态 hx,,于是在计算 （ ８,１５ ） 式的转移概率时,就只需算比值
)()()(( hxpp hxb
h
x ≠= HHe,而并不需要计算配分函数,Glauber 动力学的构架,也正是用了这一点 ）,由这个有限状态 Markov 链的互通性,我们有
P →? ∞→nn pT1=,
因此,在时间发展充分长以后,我们可以 用 Metropolis 的 Markov 链所处的状态,作为按分布 p 取的样本,也就是说,与 Gibbs 采样法一样,Metropolis 方法也给出了在计算机上模拟 p - 随机数的一个算法,Metropolis 提出的这种采样法,称为 Metropolis 采样法,它与
Gibbs 采样法的不同处在于,对于 Metropolis 采样法而言,任意两个组态 hx,,只要预选概率 0
~
>xhp 就可以转移,
Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作,
(1) 设当前为时刻 n,取的状态为 ix n =)(,对它作随机扰动,即取一个分布为
),,( ~1~ iKi pp L 的随机数,设为 j ;
(2) 若 1≥
i
j
p
p
,则将状态更新为 jx n =+ )1( ； 否则 进行 （ ３ ）；
(3) 独立地取一个 ]1,0[U 随机数 U,如果
i
jU
p
p≤
,则将状态更新为 jx n =+ )1( ； 否则状态不更新,即令 ix n =+ )1(,
（ 请读者证明,如此由 i 到 j 的转移的可能性恰是 （ ８,１５ ） 式规定的转移概率 ）,
~P
的对称性并非必要,理论分析指出,经过适当的选取 (研究矩阵
~P
的第二个特征值 ),
使用非对称的
~P
可能加快收敛速度,对于 非对称的预选矩阵
~P
,Metropolis 采样法所 构造的 Markov 链的转移 应取 下 式,
212
)(),,1min( ~
~
~
ij
p
ppp
iji
jij
ijij ≠=
p
p
,(8,16)
此时所构造的 Markov 链仍然以 p 为可 逆分布 （ 请读者验证,并给出这种 Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作步骤 ）．,
具有非对称的预选矩阵的 Metropolis 采样法 的一个特殊情形就是 独立的 Metropolis 采样,即 预选转移矩阵与出发点无关 的 Metropolis 采样,此时预选矩阵 各行是一样的,因此 可以写为
1P =~ rT,
记
j
j
jw r
p=
,结合拒收原则,对独立的 Metropolis 采样实际 操作 为
（ 1） 按分布 r 采样,得到 j ；
（ 2） 若 1
)(
≥
nx
j
w
w,则取 jx n =+ )1( ; 否则进行 (3);
（ 3） 取 ]1,0[U 随机数 U,修改 )(nx 为 )1( +nx,
≤=+
其他情形 ）(
),1min((
)(
)1(
)(
n
x
j
n
x
w
wUj
x n 。
注 １ 一般地 Markov 链的两个状态之间都可能转移,因此 Markov 链相当于连续取
值情形中的 ＂ 高维情形 ＂,预选矩阵的设计是尽可能地减少,或者简化 状态间的连接,以减少计算量,有时 它相当于局部化,
注 ２ 在独立的 Metropolis 采样法,且 p 是一维密度的情形,我们还可以把 Metropolis
算法视为 Von Neuman 拒收原则在某种意义下的推广,此时 拒收原则中的 )(0 xp 就相当于预选矩阵,为了方便 做 比较,我们把 p 改记为 )( xp,此时 独立 的 Metropolis 算法可以设计如下,取一个与 分布密度 )( xp 取值范围差不多的随机数的独立样本 L,,10 hh,记其分布密度为 )(0 xp,再取与它独立的 ]1,0[U 独立随机数列 L,,21 UU,令 00 hx =,然后归纳地定义
≤=
+
+
+
+
)
p
p
p
pU
n
n
n
n
n
nn
n
其它情形(
)1,)( )()( )(min((
10
01
1
1
x
h
x
x
hh
x,
于是当 n 大时,nx 就近似地是 p 随机数,(这里 )(0 xp 起了 ijp 的作用 ),把 它与拒收原则相
213
比较可以看出,Metropolis 算法不需要假定 )( )(
0 xp
xp 的有界性,为此付出的代价是,不能从是否有 )( )(
0 h
h
Cp
pU ≤ 来决定是否 将 x 取 为 h,代 替一个 更新的 随机数 h 以一个 随机数的 近似序列 nh,并且从是否有
≤nU 序列 )( )(
10
1
+
+
n
n
p
p
h
h 与其前面的值
)(
)(
0 n
n
p
p
x
x 的比
来决定是否把 nx,更新 为,1+nh,所以 可以说,独立的 Metropolis 算法主要用在找不到
)(0 xp 使 )( )(
0 xp
xp 有界的情形,特别是在 )( xp 只 能 得到测量的数据值,而没有显式表示的情形,然而,在实际计算中 一般地 只在 ＂ 高维情形 ＂ 才有必要使用 Metropolis 算法,
注 ３ 比较这 Gibbs 采样法 与 Metropolis 采样法 的收敛速度,是一个重要的理论问题,
注 ４ 使用 非对称的
~P
的 另一种 Markov 转移的 取法,是 如下的 Hastings 推广了的 Metropolis 采样
)(,
1 ~
~
~
ij
p
p
spp
iji
jij
ij
ijij ≠
+
=
p
p
,(8,17)
其中
~P
为非对称的转移矩阵,满足,ijp
~
与 jip
~
要么同时为零,要么同时为正,而 )( ijs 是一个非负的对称矩阵,满足
)(,1
1 ~
~ ij
p
p
s
iji
jij
ij ≠≤
+
p
p
,(8,18)
以上几种采样法都需要知道 p 的解析表达式,在并不知道它的表达式时,下面的方法对于采集以它为分布的样本,能起很好的作用,
2,3 通过条件分布对分布 p 作 随机采样的 Gibbs 方 法
Gibbs 采样法 中 Markov 链的状态的更新 按 Gibbs 方式 进行,即更新时 只变动组态的一个分量,它 实行起来较为简单,所以这种思想可以用来得到高维总体的样本,记 高维总体的分布为 p （ 未必是 Gibbs 分布 ）,Gibbs 方式 更新 的思想,在用于甚高维总体或复杂总体的取样时,主要是 通过 p 的 条件分布族,来构造一个不可约正常返的马氏链 nX,使它以 p 为不变分布,由于 n 充分大时,nX 的分布渐近于 p,nX 就可以近似地作为 分布 p 的 总体
214
的样本,下面我们以可数状态空间的 分布 p 为例来进行 讨论,其处理问题的思想,也适用于生成有密度的总体的样本,
1 由条件分布族得到以 p 为 不变分布的 Markov 转移的思想
Gibbs 样本生成法的要义,在于把高维总体的取样,化为一系列一维分布的 取 样,而后者是容易做到的,甚至可以利用 SAS 等软件得到,但是,在具体运作时,对于较大的 m 及一个给定 的 m 维分布 p,事情就远非如此简单,设 p ),,( 1 mxx Lp=,记 ),,( 1 mxxx L=,
),,( 1 myyy L=,再记在第 1 至 k - 1 个分量固定为 11,,?kyy L,同时将第 1+k 至第 m
个分量固定为 mk x,,x L1+ 的条件下,第 k 个分量在 ky 的条件概率为
|( kyp 11,,?kyy L,)x,,x mk L1+,
我们定义如下的转移概率
∏
=
==
m
k
xy yxpp
1
),( |( kyp 11,,?kyy L,)x,,x mk L1+,(8,19)
容易验证 )( xyp 是一个转移矩阵,我们 证 明 p 是它的不变分布,事实上
∑
x
yxpx ),()(p ∑=
x
x)(p ∏
=
m
k 1
|( kyp 11,,?kyy L,),,1 mk xx L+
∑ ∑∑=
mxx
x
m
m
m
x xx
xxyxx
,,1
21
1
2
1
1 ),,(
),,,()],,([
L L
LL
p
pp ∏
=
m
k 2
|( kyp 11,,?kyy L,),,1 mk xx L+
∑=
mxx
mxxy
,,
21
2
),,,(
L
Lp ∏
=
m
k 2
|( kyp 11,,?kyy L,),,1 mk xx L+
∑ ∑∑=
mxx
x
m
m
m
x xxy
xxyyxxy
,,21
321
21
3
2
2 ),,,(
),,,,()],,,([
L L
LL
p
pp ∏
=
m
k 3
|( kyp 11,,?kyy L,),,1 mk xx L+
∑=
mxx
mxxyy
,,
321
3
),,,,(
L
Lp ∏
=
m
k 3
|( kyp 11,,?kyy L,),,1 mk xx L+
)(),,( 1 yyy m pp === LL,
[ 注 ] 事实上,我们并不需要知道 )(xp 的表达式,而只需知道在固定其它一切分量的条件下,余下的某一个分量的条件分布,就可以得到以 ),( yxp 为转移的时齐马氏链的样本,
2,构造 Markov 的方 法
从已得到的 Markov 链在时刻 n 的样本 nX,去求转移到时刻 1+n 的样本 1+nX 时,转移 ),( yxp 是一维 的 条件 分布 的乘积,而 这些一维条件分布的样本是较为容易生成的,即可
215
按以下程序逐个地得到 1+nX 的样本的各个分量 ),,( 1 myy L,
（ １ ） 先 得到服从分布 }},,1{),,,|({ 121 NSyxxy m LL?=∈p 的随机变量 （ 记为
1,1+nX ） 的一个样本 1y ( 在用 Von Neumann 取舍原则取样时,对于
),,,(
),,,(),,|(
2
21
21
m
m
m xx
xxyxxy
L
LL
+∞= p
pp,只需假定对固定的 ),,(
2 mxx L,存在 )( 10 yp 及常数 C,满足 )(),,,( 1021 yCpxxy m ≤Lp,并对 ]1,0[U 均匀随机数 U,看 )( ),,,(
0
2
UCp
xxU mLp
是否不小于 U,来决定取舍样本 U ),
（ ２ ） 用同样的方法,再得到服从分布 }},,1{),,,,|({ 2312 NSyxxyy m LL?=∈p 的随
机变量 （ 记为 2,1+nX ） 的一个样本 2y,
（ ３ ） 依此下 去,得到服从分布 }),,,,,,|({ 111 Syxxyyy kmkkk ∈+? LLp
的随机变量 （ 记为 knX,1+ ） 的一个样本 ky,)1(?≤ mk,
（ ４ ） 最后得到服从分布 }}),,,|({ 1 Syyyy mmm ∈Lp 的随机变量 （ 记为
mnX,1+ ） 的一个样本 my,于是 ),,( 1 myyy L= 就是 的一个样 本,
现在我们任取 )0(0 yX =,再按上面方法得到 1X 的一个样本 )1(y,对 n 归纳地我们可用上面方法 分别 得到 nXX,,2 L 的样本 )(,),2( nyy L,当 n 充分大时,由于 Markov
链 nX 的分布近似于 )(?p,我们就可以认为 )(ny 近似地是 分布 )(?p 的一个样本,
［ 注 ］ 上面的 Gibbs 采样也可以看成为如下的算法
（ 1） 按循环次序更新状态的各个分量,在更新了第 d 个分量后,就接着更新第 － 个分量 ；
（ 2） 按条件分布 ),,,,,|( )()( 1)( 1)(1 ndnjnjn xxxx LL +p 选取 )1( +njx,再把状态 ),,( )()(1 ndn xx L 更新为
+++
=),,,( )1()1()1(1 ndnjn xxx LL ),,,,,,( )()( 1)1()( 1)(1 ndnjnjnjn xxxxx LL ++? 。
[注 ] 回顾 Gibbs 采样,还可以从理论角度得到以下认识,如果状态空间是实数,我们把
Gibbs 采样中第 k 个分量的条件更新看成 kx 到 kkk xy x+= 的随机 (加法 )变换,kk xTy
kx
=,
那么 Gibbs 采样的更新就是 说,,只要 kx 服从条件分布 |( kyp 11,,?kyy L,),,1 nk xx L+,则目标函数 (即不变密度 )(xp )在这类随机变换下是不变的,这种看法有助于进一步设计有效
216
的 MCMC 采样方法,
３ Gibbs 采样法与 Metropolis 采样法的变种应用
第 １ 种 变种 （ 随机变种 ）
这时,状态 ),,( )()(1 ndn xx L 的更新,不按指标次序轮番更替,代之以
（ 1） 每次按某个 ＂ 整值 ＂ 随机变量 h 的采样值 j,来随机选取指标 ；
（ 2） 再按条件分布 ),,,,,|( )()( 1)( 1)(1 ndnjnjn xxxx LL +p 选取 )1( +njx ；
（ 3） 然后把 ),,( )()(1 ndn xx L 更新为,
+++
=),,,( )1()1()1(1 ndnjn xxx LL ),,,,,,( )()( 1)1()( 1)(1 ndnjnjnjn xxxxx LL ++? 。
第 ２ 种 变种 （ 分裂变种 ）
设状态变量 x 可以分为两部分 ),( 21 xxx =,其中 2x 是丢失 数据 的部分,分别记这两部分随机变量相互间的条件分布为 )|( 121|2 xxp,)|( 212|1 xxp,如果我们可能在 1x 已知时对 )|( 121|2 xxp 进行采样,也可能在 2x 已知时对 )|( 212|1 xxp 进行采样,那么,我们可以 按下面步骤进行,
(1) 任意设置一个初值,例如
)0(
1x ;
(2) 由条件分布 )|( )(22|1 nx?p 中抽取
)1(
1
+nx;
(3) 再由条件分布 )|( )1(11|2 +? nxp 中抽取
)1(
2
+nx
,这样依次交替进行,
只要 )|( 121|2 xxp,)|( 212|1 xxp 恒正,那么由 上 面 的 Gibbs 采样的一般理论可知,),( )(2)(1 nn xx 就是一个 Markov 链的采样值,而且 ),( 21 xx 的联合分布是这个链的不变分布,于是当 n 充分大的时候,可以认为 ),( )(2)(1 nn xx 是 近似地 采自 ),( 21 xx 的联合分布的一个样值 。
第 ３ 种 变种 （ Metropolis 采样 型 变种 ）
（ 1） 每次按某个随机变量 h 的采样值 j,来随机选取指标 ；
（ 2） 再按条件分布 ),,,,,,|( )()( 1)()( 1)(1 ndnjnjjnjn xxxxxx LL +? ≠?p 选取一个样本 z ；
（ 3） 取一个 与之独立的 ]1,0[U 样本 U,令
≤
= +?
+?
+
其他情形 ）(
)),,,,,|( ),,,,,|(,1min((
)(
)()(
1
)(
1
)(
1
)()(
1
)(
1
)(
1
)(
)1(
n
j
n
d
n
j
n
j
n
n
d
n
j
n
j
nn
j
n
j
x
xxxxz
xxxxxUz
x LL
LL
p
p；
（ 4） 然后把 ),,( )()(1 ndn xx L 更新为
217
+++
=),,,( )1()1()1(1 ndnjn xxx LL ),,,,,,( )()( 1)1()( 1)(1 ndnjnjnjn xxxxx LL ++? 。
４,利用可逆链防止多峰情形的偏取
在 Gibbs 样本生成过程中,只有当马氏链 { nX } 发展的,时间,n 充分大时,它的分布才近似于 p,在实践中怎样才能判断 n 是否已经足够大了呢？ 通常的判断标准一般是看在计算中的 nX 的分布是否比较稳定了,但是,在事实上这种判断方法 不尽正确,下面的示意图给出了 p )( xp= 具有两峰的情形,从不同初值出发的轨道,可能分别稳定在 1E 与 2E
两个区域,而事实上 nX 还远未达到,近似于,不变分布的要求,
E 1 E 2
也就是说,按上一段的方法得到的 Markov 链,从不同的初值出发,在计算达到看起来稳定时,它们还远不能近似 π (x) 的样本,原因在于这时 的 平稳分布 π (x) 有两个峰,且在两峰中间的转移概率非常小,以使从 E 1 中的点出发能到达 E 2 和从 E 2 中的点出发能到达 E 1 的概率都极为微小,从而在一条轨道上要从 E 1 到达 E 2 就要化费十分长的时间,有时甚至在现实生活中可能不会在有意义的时间内达到,发现与使用 Gibbs 样本生成法的先驱们往往希望通过选取多个初始值,得到多条轨道以解决这一问题,但是这并不能完全解决问题,正如在当前的例子中的两条轨道各自只局限在不同的区域 E 1 和 E 2 中,它们的分布,根本反映不出分别限制在 E 1 和 E 2 上的 p )( xp= 如何拼接的,即以多大的比例 1C,2C,使
)( 21 CC + p 1C≈ p
1
|E + 2C p
2
| E,
只有找到了这个比例 1C,2C,才能从分别限制在 E 1 和 E 2 上的 p
1
|E 和 p
2
| E,拼接出整体的 π (x),
为了找 出 比例 1C,2C,我们建议利用可逆 Markov 链来生成 Gibbs 样本,并 从可逆
Markov 链的特性得到不变分布在它的各个峰区上的比例,下面我们以双峰时两条轨道的情
218
形,来阐述确定与估计这个比例的方法,
设分别稳定在 E 1 和 E 2 上两条轨道的样本分布对应于分布 )(1 xp 与 )(2 xp,我们近似地把 )(xp 在 E 1 和 E 2 外为 0,于是由可逆 Markov 链的性质可知,它们应该分别与 { π (x):x
∈ E 1} 及 { π (x):x ∈ E 2 } 成比例,即
),(,)( )()( )( i
i
i Eyx
y
x
y
x ∈?=
p
p
p
p,
于是
)(
)(
^
x
xi
p
p ≈
i
i C
x
x?=
)(
)(
p
p,
因此
)(
1
i
i EC p≈,
具体确定 iC 可 如下 进行,对 ii Ex ∈,利用可逆性我们有
)()(
)()(
122
211
2
1
xx
xx
C
C
pp
pp=
),()(
),()(
1222
2111
xxpx
xxpx
p
p=,
于是对 iEyx ∈,,有
∑∑
∈∈
==
ii Ex
i
Ex
i yxp
xypyx
),(
),()()(1 pp ]
),(
1[)(
1 yXp
Ey yip=,
对于 1X 的 N 次取样 )(1)1(1,,NXX L,我们可以得到
1
1
]),( 1[)(
=
yXpEy yip 的估计
1
1
)(
1
^
),(
11)(?
=
= ∑N
k
ki yXpNyp,
这里的近似只需对从 y 出发的 Markov 链,作独立的 N 次一步转移的模拟,它比 Markov 链的分布近似于不变分布要求的步数要少得多,
余下的问题是从不变分布 π (x) 具体构造可逆 Markov 链,这只需把前面构造的马氏转移
),( yxp 修正为
+=
)(
),()(),(
2
1),(^
x
xypyyxpyxp
p
p,(8,20)
容易检查 它 是以 )(xp 为可逆分布的转移阵,这恰是我们所希望要的,
219
为了便于具体操作,我们还希望 ),(^ yxp 也按 Gibbs 方式转移,即也能归结于计算每次只变动一个分量的 m 个,条件 转移概率,的乘积,为此我们注意,
)(
),()(
x
xypy
p
p
),,(
),,|(),,,|(),,|(),,(
1
11312211
m
mmmmm
xx
xxxyyxxyyxyy
L
LLLLL
p
pppp?=
),,(
),,(
1
1
m
m
xx
yy
L
L
p
p=
dyyyy
yyx
m
m
),,,(
),,,(
2
21
L
L
p
p
∫ dyyyxx
yxx
mm
mm
),,,,(
),,,(
21
11
∫ L
LL
p
p
dyyxx
xx
m
m
),,,(
),,(
11
1
∫
× LLp p
).,|( 21 myyy Lp= ),.,,|( 111 mmm yxxy LLp ).,|( 11?mm xxy Lp,
即 )( ),()( x xypypp 具有类似于 ),( yxp 的形式,不同处只在于所有条件分布中的变量次序倒了一下,而 ),(^ yxp 的含意为,分别以 21 的概率取 ),( yxp 与取 )( ),()( x xypypp,所以,从模拟程序的实现来看,这仅仅增加了一倍计算量,
2,4 MCMC 应用于 Bayes 参数估计
１ 关于 Bayes 统计
我们将在第 ９ 章中较多地介绍与评述 Bayes 统计,Bayes 统计方法的基本想法是,把概率函数中的 未知参数 q 当作随机变量,
由于在 Bayes 方法中不再区分参数与随机变量,所以未知参数的分布的确认是至关重要的,在抽取样本以前,就只能根据先验知识,设置未知参数的分布,称为 先验分布 。 在取样以后,根据对于抽 取 到的样本的概率规律的了解,使用 Beyes 公式,把参数的先验分布改进为 后验分布,后验分布 也称 Bayes 分布 。
设随机变量 x 的分布为 ),( qxp （ 未知参数 q 可以是向量 ）,把 q 看成随机变量 （ 或随机向量 ） 后,),( qxp 就是,在 q 已知的条件下,x 的条件分布,即 ),()|(| qqqx xpxp =,假定 q 的先验分布的分布密度为 )(Jf,于是 ),( qx 的联合分布为 ),()( qqf xp,
在 J 固定条件下,设 nxx,,1 L 为采自 ),( qxp 的独立随机样本,把在 nxx,,1 L 已知的条件下 q 的条件分布,记为 ),,|( 1 nxxJj L,那么,它等于 q 与 ( nxx,,1 L )的联合密度
220
),()(
1
qxqf i
n
i
p∏
=
,除以 ( nxx,,1 L )的边缘密度,即
),,|( 1 nxxJj L
∫ ∏
∏
=
=
n
i
i
n
i
i
dp
p
1
1
),()(
),()(
JJxJf
JxJf
＝,
这就是 q 的 Bayes 分布,所以至关重要的是如何选取参数的先验分布,这个问题将在第 ９ 章中作一些较为浅显的介绍,
2 用 MCMC 方法对 Bayes 分布参数作估计
对于参数 q 的 Bayes 分布,用 MCMC 方法 做 以它为极限分布 的 Markov 链 nJ,于是,
由于 nJ 的分布近似于此 Bayes 分布,我们可以取
^J
N
Nmm ++ ++= JJ L1,(
0nm ≥,0n 充分大 ) 作为参数 J 估计,
３,用 MCMC 方法得到较为复杂的参数估计的大意
较复杂的参数估计的典型情形是 混合总体的参数 Bayes 估计,设样本 nXX,,1 L 来自权重参数未知的混合总体,这些未知权重也就作为另加的待估计参数,在对混合分布作参数估计时,可以引进一些新的随机变量 (在 Bayes 方法中,原则上参数与随机变量不再严格区分,可以把它们看成,补偿参数,) Z,然后用 Gibbs 采样法,交替使用,在? 已知的条件下,Z 的条件分布 )|(1|2?zp,与在 Z 已知的条件下,
的条件分布 )|(2|1 zp? 。 由此可得到 ),( Z? 的 Bayes 分布的采样,其步骤如下,
（ 1） 任给初值 0? ；
（ 2） 按条件分布 )|(1|2 nzp? 采样,得样本 nz ；
（ 3） 按条件分布 )|(2|1 nzp? 采样,得更新样本 1+?n 。
于是只要 )|(1|2?zp 与 )|(2|1 zp? 恒正,那么就有,),( nn z? 是 Markov 链,而且以 ),( Zθ 的联合分布为极限分布 。 从而当 m 大时,可以认为 ),( mm z? 近似地是出自 ),( Zθ 的联合分布的采样,同时可以用取
^J
N
Nmm ++ ++= JJ L1 作为? 的估计量 。
以上的关键在于知道 )|(1|2?zp 与 )|(2|1 zp?,而且易于抽取样本,它们还不能以正概率取 0,以保证所作的 Markov 链的转移函数的收敛性,关于 )|(1|2?zp 与 )|(2|1 zp? 的取法与获得,在现代 一些 统计论文中,有一些典型的做法,本书不再叙述,
221
４ MCMC 方法用于缺失数据情形下的参数估计的大意
如果 ),( YZX = 的分布密度有未知参数?,),,(?yzp,其中 Y 是可以通过采样观测到的,但是,Z 却观测不到 （ 或者是由数据的缺失而引起的,它称为缺失数据 ）,用 Bayes 统计的思想,在观测样本 Y 条件分布给定的条件下,记条件分布为
)|(1|2?zp ),,(?= Yzp,)|(2|1 zp?,
只要它们恒正,而且有一种方法可以近似地取得服从这两个分布的样本,那么就可以用第 ３ 段中的方法得到 缺失数据 Z 的估计,
3,模拟退火 ( simulated annealing )
３,１ 模拟退火方法的基本想法
一般的优化问题,都可以归结为求某个目标函数 f 在一定范围内的最小值,但是许多常用的方法,往往都是运用某种局部的比较去决定优化过程的发展,以逐步达到最小值,例如梯度法,就是每次将目标函数在当时所考察的点与附近的点的取值加以比较,根据梯度方向就是目标函数值下降最快的方向,导出了梯度下 降法,即在被考察的点沿函数下降最快的方向移动,希望经过多次这样的移动,最终达到目标函数 f 最小的点,但是这个算法的实际结果是停止在梯度为零的某个点上,这就是说,在计算进行到到达目标函数的某个局部极小值点时,在考察点的目标函数值就小于它附近所有点的目标函数值,于是搜索过程就停止了,我们知道目标函数的任意一个局部极小值点都满足上述条件,然而局部极小值点处的函数值,有可能比函数的实际最小值大得多,于是求得最小值的努力就完全寄托于初值的选取,
为了克服上述算法的这个缺 点,Metropolis 等人在 1953 年提出的一种有效地模拟固体达到热平衡发展的算法,而 Kh asi ’minski 在 1965 年对决定性的梯度算法,引入小的随机噪声干扰,使之在达到平衡时,以高概率进入目标函数的最小值附近,但是当时没有受到重视,Kirkpatrick 等人于 1983 年对组合优化问题,重新提出模拟退火算法,其基本想法正是 Kh asi ’minski 的构想,即引入人工噪声使得当算法达到局部极小值时,能有一个小概率
,逸出,局部极小值的陷井,这个想法源自冶金退火过程的物理思想,即将金属加温至充分高,再使其慢慢地退火冷却,在加温时,金属内部粒子随温度增加变为无序,能量增大 ; 而在退火时,在慢慢冷却的过程中粒子渐趋于有序,在每个温度时都达到平衡态,最后在常温时达到能量最大的基态,而用退火模拟优化问题时,则将能量模拟为目标函数,温度演化成代表噪声水平的控制参数 b,这样就得到优化问题的模拟退火算法,
简单而实用的模拟退火模型是有限状态模型,典型情形是组合优化问题,这时假定自变量只能取自有限个值,而加入随机噪声的算法,就是把原来的计算过程变成随机过程 — 也就是一个 Markov 链,并使得此 Markov 链是互通的,非周期的,正常返的,且具有如下形式的不变分布,
222
))((
))((
)( *
*
1 fif
j
fjfi ee
∑=
b
b
bp )(
)(
1 if
j
jf ee
∑=
b
b,（ ８,20 ）
其中 f 是我们要优化的目标函数,)(min* iff i=,而 b 是一个控制噪声水平 ( 随机性水平 )
的参数,这个分布有一个重要的性质,当 ∞→b 时,其极限分布集中在目标函数 f 取最小的集合上,即在 })(:{ *fifi = 上 ( f 取最小值的点可能不唯一 ),于是我们可以利用这个
Markov 链的发展,使其到达不变分布,然后令 β →＋ ∞,以达到 f 取最小值的点集上,而固定 b 时的 Markov 链的发展的取样 过程,则可以通过上一节所叙述的各种 M arkov M onte
C arlo 方法来实现,例如可用 Gibbs 采样法,Metropolis 采样法,或条件概率的采样法,
(8,20) 中的不变分布的形式,很像在物理上的 Boltzmann 分布,而 b 则是倒温度,即 b ＝
T
1 （Ｔ 是绝对温度 ）,若把无噪声 ( 即 b ＝＋ ∞ 情形 ) 的决定性算法,看成是绝对 零度 ( Ｔ
＝ 01 =b ) 的情形,那么,β ＜ ∞ 的情形就是将温度升为正温度,于是上面所述的达到目标函数取最小值的集合的过程,即对固定的 b,构造 Markov 链的样本,再让 ∞→b 的过程,
恰似先将温度由绝对零度升高 ( 即 b ＜ ∞ ),然后再让温度降回到绝对零度 ( 即令 b = ＋ ∞ )
的过程,这个升温再降温的过程,就是,退火,过程的模拟,因而这个算法称为 模拟退火算法,或称 随机松驰法,简言之,即 升温以得 到 不变分布,再降温 ( b →＋ ∞ ) 使不变分布集中于目标函数 f 的取最小值的集合上,
为了加快计算速度,人们在理论上发现,可以采取受人工噪声干扰的非时齐的 Markov
链,把 Markov 链的时间发展与降温这两步,归并成一步作模拟退火,即 使所选的 Markov
链的第 n 步转移矩阵随 nb 而改变,当时间 n 趋于无穷时,只要 nb (称为 冷却进度 （ Cooling
Schedule） ) 以适当缓慢的速度趋于 ∞,则此非时齐的 Markov 链,就依概率收敛到目标函数达到最小值的点集上,这个非时齐的 退火的过程完全由冷却进度所控制,如果 nb 减小过快,则在直观上就相当于在每个温度段还未达到平衡,退火过程就还会陷入能量函数的局部极小 ; 反之,如果 nb 减小过慢,相当于几乎是常值,退火过程就会按平衡态 (即不变分布 )发展,这两种情形都不能达到整体极小,
更确切地说,存在一个由目标函数 f 确定的正常数 a,只要取
nn log?= gb ( )g < a,(8,21)
223
( 即取温度
n
T b1= ),则当 n →∞ 的进程就达到了,适当缓慢地降低温度,这个要求,当然 g 取得越接近 a 效果越好,计算接近目标函数取最小值的点集的进程越快,
在实际计算时,当 n 大时,按 (8,21) 式取的 nb 的变化非常小,以致在实际上几乎看不出变化,所以有时人们就简便地使用一个固定的很小的 b,构造具有 （ ８,２０ ） 式的不变分布的 Markov 链,用其发展以代替模拟退火,显然,这种 Markov 发展的终极实际上是一个平衡分布,但是只要这个平衡分布高度集中于目标函数全局极小的集合附近,那么这样的计算误差也就可以容忍了,在具体操作时,到底如何掌握分寸,是一个实际经验问题,至今还缺少中肯的理论分析,
在许多情况下,由于所处理的系统非常复杂,或者实际维数太高,还有的时候目标函数没有明确的分析表达式,这时人们就引入局部改变的邻域,记状态 i 的邻域为 )(iN （ 这种邻域是根据问题的要求自然地确定的 ）,作矩阵
∈==?
))((0
))(()( ),( 1
,
~
,
~~
iNj
iNjiNppP
jiji,（ ８,２２ ）
它规定了只 许在近邻之间转移,即当系统处于状态 i 时,它下一时刻只能转移为 )(iN 中的状态,再规定所要的 Markov 链第 n 步转移概率矩阵由如下的 Metropolis 采样,
P )(n ))(( nijp b
D
=,（ ８,２３ ）
=
≠
=
∑ ≠ ijp
ijp
p
i
k
ikik
i
j
ij
ji
},1min{1
},1min{
)(
)(
)(~
)(
)(~
,
b
b
b
b
p
p
p
p
b （ ８,２４ ）
所给定,那么,运 转这一族转移矩阵 { P )(n }所决定的非时齐 Markov 链,只要温度
n
nT b
1=
以合适的速度下降,就收敛到使 f 取最小值的点集,
上面指出的算法的局部化思想,是模拟退火算法中为了克服因为状态太多而引起的计算爆炸的有力手段,这里的矩阵
~P
可以一般化地取一个互通的概率转移矩阵,显见它越稀疏 ( 即其中 0 越多 ) 则进程 可能 越快,它称为 局部化预选矩阵,简称 预选矩阵,它相当于规定了可以转移的,邻域,,把整体搜索变为局部搜索,以减少计算量,在研究第 9 章中的
Gibbs 场的 MCMC 及模拟退火时,就以邻域转移代替这个预选矩阵,
3.2 有关模拟退火算法的非时齐马氏链的理论背景
224
我们的仍旧假定状态集 ( 或称为组态空间 ) 是有限集,即 },,2,1{ MS L=,设 f 为 Ｓ 上的一个实值函数,令 L(S) 为 S 上 全体实值函数的集合,对于任意 ))(,),1(( Mfff L=,
))(,),1(( Mggg L= ∈L(S),定义
范数 ∑
=
=
M
i
iff
1
|)(|||||
D
,距离 ||||),( gfgfd?=,
那么 L(S) 就成为完备的距离空间,显见 S 上的概率分布都是 L(S) 中的元素,而且 S 上的任意概率分布的序列在此距离空间中的极限仍 是一个概率分布,( 但是,如果 S 不是有限集而是具有可数个元素的集合时,
此结论不真 ),
由定义我们立即推得,对于任意一个转移矩阵 P,恒有 f|| P ||||| f≤,其中 f P指 f 与 P作为矩阵相乘,
(1) 时齐的 Markov 链优化的 模拟退火 方法
定理 ８,２ 设以 P ))(()( bb ijp= 为转移矩阵的 Markov 链记为 }{ bx n,则
1)(limlim * ==∞→∞→ ffP
nn bxb
,
证明
==∞→∞→ )(limlim *ffP
nn bxb
}))(:{(limlim *fifiP nn =∈∞→∞→ bb x
1lim )(
})(:{ *
== ∑
=
∞→
b
b p i
fifi
,
( 2 ) 非时齐 Markov 链的模 拟退火
非时齐 Markov 链的分布的收敛性有
定理 ８,３ （ Dobrushin - Isaacson - Madsen 定理 ） 设 nx 是 状态空间 S 上的非时齐 Markov 链,其第 n 步转移矩阵为 P )(n,如果它们满足,
(A.1) P )(n 作为时齐的转移矩阵时有不变分布 )(np ;
(A.2) ∑ ∞<? +
n
nn |||| )1()( pp ;
(A.3) 或者满足 Isaacson - Madsen 条件,对于任意概率分布向量 nm,,及正整数 j,恒有
P)(|| nm? L)j( P 0|| n)n( →? ∞→ ;
或者满足 Dobrushin 条件,对于任意 j,收缩系数满足
(C )()( nj PP L ) 0 →? ∞→n,
那么,
(1) )()( lim nn pp ∞→?∞ = 存在 ;
225
(2) 无论什么初分布 0m 下,此非时齐的 Mar kov 链 }{ nx 在时刻 n 的绝对分布 nm 恒有极限
)(∞→ pm
n,
证明 首先,由 (A.2) 利用三角形不等式立得 }{ )(np 是 L(S) 中的 Cauchy 列,故而 (1) 正确,
往证 (2),用三角形不等式与 Dobrushin 不等式,我们得到
|||| )()()1(0 ∞?pm nPP L
||)(|| )()()()()1(0 njn PPPP LL ∞?≤ pm |||| )()()()( ∞∞?+ pp nj PP L
≤2 (C )()( nj PP L ) |||| )()()()( ∞∞?+ pp nj PP L,
在满足 Isaacson - Madsen 条件时,上式的第一项趋于 0,( 而由 Do brushin 不等式推出
||)(|| )()()()()1(0 njn PPPP LL ∞? pm ≤2 (C )()( nj PP L ),
即在 Dobrushin 条件成立时,Isaacson - Madsen 条件必然成立 ),对第二项利用不变分布条件
)()()( mmm Ppp = 推得对于 Mj ≥ 有
|||| )()()()( ∞∞?pp nj PP L
)()()()( )(|| njj PP Lpp?= ∞ ||)( )()()()( ∞?+ pp njj PP L
)()()()( )(|| njj PP Lpp?= ∞ ||)( )()()1()( ∞+?+ pp njj PP L
)()()()( )(|| njj PP Lpp?= ∞ ||)()( )()()()()()1(
1
∞++?+
=
+?+ ∑ pppp nnkjkjkj
jn
k
PP L
||)(||sup )()( nMn pp?≤ ∞≥ ||||sup|||| )()()1()( ∞≥+
≥
+?+ ∑ pppp nMnnn
Mn
0→,
[ 注 ] 在乘积 )()( nj PP L Mkinjikp ≤
=,),( )( 中,)|( 1)1,( ikPp jnnjik === ++ xx 为第 n 次转移矩阵为 P )(n 的非时齐 Markov 链 }{ nx 的多步转移概率,由于我们研究的是有限状态的情形,于是
Isaacson - Madsen 条件可以改写为如下的等价条件 ( 称为 弱遍历条件 ),
(A,2) ’ 对于任意状态 Skli ∈,,,任意正整数 j,恒有
0)(lim ),(),( =?∞→ njlknjikn pp,
对于模拟退火的收敛性,利用较为深入的,转移概率为指数型衰减的时齐的 Markov 链序列的大偏差理论 ( 粗略地说,它研究此链越出给定小概率状态集合的概率的指数型衰减 ),可以证明 ( 这需要相当深入的 叙述与 证 明,本书无法列入 ) 下面的定理 8.4,
设 )(if 为定义在 },,1{ MS L= 上的函数,令
226
p ),,( )()(1)( nMnn pp L=,)(nip )()(1 if
j
jf
n
n
ee∑= bb,（ ８,２５ ）
P )(n )( )(nijp
D
=,
=
≠
=
∑ ≠ ijp
ijp
p
n
i
n
k
ikik
n
i
n
j
ij
n
ij
},1min{1
},1min{
)(
)(~
)(
)(~
)(
p
p
p
p
,（ ８,２６ ）
定理 ８,４ （ G emam - Gemam ）
假定
(1) ∞↑< nb0,满足,存在 0n,使 0nn ≥ 后,有
nMn ln1?≤b,（ ８,２７ ）
其中 D 是目标函数 f 的振幅,
)(min)(max ifif ii?=DD ; （ ８,２８ ）
(2) 预选矩阵
~P )( ~
ijp= 是互通的概率转移矩阵,
那么,以 P )(n 为第 n 步转移矩阵列的非时齐的 Markov 链 }{ nx 满足弱遍历条件,因而 (A,2) 满足,
定理 ８,５ （ 模拟退火的收敛性定理 ）
在定理 8.4 的条件满足时,无论什么初 始分布 0m 下,以 P )(n 为第 n 步转移矩阵列的非时齐的 Markov
链 }{ nx,在时刻 n 的绝对分布 nm 有极限 )(∞→ pmn,其中 )(∞p 为集合 min})(:{ =ifi 上的离散均匀分布,从而有
1)( * = →? ∞→ ffP nnx,
证明 我们来验证 Dobrushin - Isaa cson - Madsen 定理的条件成立,由 G emam - Gemam 定理,我们只需验证条件 (A,1),为此,首先注意当 n 充分大时 ( 例如 0nn > 时 ),利用 nb 的单调不减性,我们有,对于固定的 i,)(nip )()(1 if
j
jf
n
n
ee∑= bb 在 0)( >if 时,关于 n 是单调不减的 ( 对 nb 的微商非负 );
而在时 0)( ≤if 时,关于 n 是单调不增的 ( 对 nb 的微商非正 ),再注意 ( 这是在第 5 章 Dobrushin 不等式中应用过的事实 )
|||| )()1( nn pp?+ ∑?= +
i
n
i
n
i ||||
)()1( pp ∑ ++?=
i
n
i
n
i )(
)()1( pp ∑?+?+
i
n
i
n
i )(
)()1( pp
∑ ++?=
i
n
i
n
i )(2
)()1( pp ∑?+?=
i
n
i
n
i )(2
)()1( pp,
227
于是我们得到
∑
=
N
nn 0
|||| )()1( nn pp?+ =2 ∑
=
N
nn 0
∑ ++?
i
n
i
n
i )(
)()1( pp ∑∑
=
++
>
=
N
nn
n
i
n
i
if 0
)(2 )()1(
0)(
pp
∑∑
=
+
>
=
N
nn
n
i
n
i
if 0
)(2 )()1(
0)(
pp ][2 )()1(
0)(
0n
i
N
i
if
pp?= +
>
∑ 2≤,
令 ∞→N,便得 ∑
∞
= 0nn
|||| )()1( nn pp?+ 2≤,所以 (A,1) 成立,定理证毕,
习题 ８
1,证明对于 Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 采取的操作得到的由 i 到 j 的转移,就 是
（ ８,１５ ） 式 定义 的转移概率,
2,证明由 （ ８,１６ ） 给出的转移矩阵 所对应的 Markov 链以 p 为可逆分布,再给出这种 Metropolis 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作步骤,
3,证明由 Hastings 推广了的 Metropolis 采样所对应的 Markov 链以 p 为可逆分布,再给出这种 采样在时刻 n 的更新 )1()( +→ nn xx 可具体采取如下的操作步骤,
