龚光鲁,钱敏平著应用随机过程教程及其在算法与智能计算中的应用 清华大学出版社,2003 第8章 Monte Carlo与 Markov chain monte carlo(MCMC)方法 在许多很复杂的统计问题中,有时很难直接对各种统计方法进行理论分析.为了评估它 们的优劣,常见的实用办法是做随机模拟:即设法按问题的要求与条件去构造出一系列的模 拟样本,用它们的样本频率代替对应的概率作统计分析与推断,观察由这些摸拟样品所作出 的推断的正确率.因为在概率论初期发展时,随机模拟的原型常常来自博采,于是人们就以 博采之都 Monte carlo作为随机模拟方法的别称.久而久之, Monte carlo方法作为名称 倒比随机模拟方法更为广泛地被常用了.相仿地,人们还把组合计算中的某些随机模拟方法, 称为 Las vegas方法,这是以美国的博采城 Las vegas命名的 般随机模拟方法的长处是,计算的复杂度不依赖于计算空间的维数因此,在计算非 常高维数的积分(或多指标的求和)时, Monte Carlo方法比通常的计算方法有明显的优越 性.而很多的 Monte carlo计算问题,可归结为计算积分或庞大和数的问题 1.计算积分的 Monte carlo方法与采样量估计 通过构造独立同分布随机数,计算积分的 Monte carlo方法,也称为静态 Monte Carlo 方法,其思想可以在本节中,通过估计最简单的积分f(x)dx得到阐明.对于高维积分, 其思路与一维积分是一样的.另一方面,在甚高维情形,因为计算量太大,用静态 Monte Carlo方法处理速度太慢.我们在第2节中通过构造 Markov链的极限不变分布,来模拟 计算积分的方法(称为动态 Monte Carlo方法, Markov链 Monte carlo方法,MCMC),将 更为适用 1.1用频率估计概率来计算积分的 Monte Carlo方法 假定0≤f(x)≤M.那么由积分的面积含义有 f(x)ax=S|,(其中S|为S={(x,y):a≤x≤b,f(x)≥y}的面积) 考虑平面区域Ω=如,b×[0,M上的均匀随机变量5,则 p=P(∈S) 「/(x)k 对于N个独立的-均匀随机数;,(≤N),记 Ns={1,…5x}中落在S中的频数 于是,利用大数定律便知 (b-a)
202 龚光鲁，钱敏平著 应用随机过程教程及其在算法与智能计算中的应用 清华大学出版社,2003 第 8 章 Monte Carlo 与 Markov Chain Monte Carlo (MCMC) 方法 在许多很复杂的统计问题中,有时很难直接对各种统计方法进行理论分析．为了评估它 们的优劣,常见的实用办法是做随机模拟：即设法按问题的要求与条件去构造出一系列的模 拟样本,用它们的样本频率代替对应的概率作统计分析与推断, 观察由这些摸拟样品所作出 的推断的正确率．因为在概率论初期发展时, 随机模拟的原型常常来自博采, 于是人们就以 博采之都 Monte Carlo 作为随机模拟方法的别称. 久而久之, Monte Carlo 方法作为名称 倒比随机模拟方法更为广泛地被常用了. 相仿地, 人们还把组合计算中的某些随机模拟方法, 称为 Las Vegas 方法, 这是以美国的博采城 Las Vegas 命名的. 一般随机模拟方法的长处是，计算的复杂度不依赖于计算空间的维数. 因此, 在计算非 常高维数的积分（或多指标的求和）时, Monte Carlo 方法比通常的计算方法有明显的优越 性. 而很多的 Monte Carlo 计算问题, 可归结为计算积分或庞大和数的问题. 1. 计算积分的 Monte Carlo 方法与采样量估计 通过构造独立同分布随机数, 计算积分的 Monte Carlo 方法, 也称为静态 Monte Carlo 方法, 其思想可以在本节中,通过估计最简单的积分 ò b a f (x)dx 得到阐明. 对于高维积分, 其思路与一维积分是一样的. 另一方面，在甚高维情形, 因为计算量太大, 用静态 Monte Carlo 方法处理速度太慢. 我们在第 2 节中通过构造 Markov 链的极限不变分布，来模拟 计算积分的方法（称为动态 Monte Carlo 方法， Markov 链 Monte Carlo 方法, MCMC), 将 更为适用. 1. 1 用频率估计概率来计算积分的 Monte Carlo 方法 假定 0 £ f ( x) £ M . 那么由积分的面积含义有 f (x)dx | S | b a = ò , (其中| S | 为 S ={(x, y) : a £ x £ b, f (x) ³ y} D 的面积 ). 考虑平面区域W =[a,b]´[0,M ] D 上的均匀随机变量x , 则 b a M p P S ( ) 1 ( ) - = Î = D x ò b a f (x)dx . 对于 N 个独立的W -均匀随机数 ,(i N) xi £ , 记 { , , } NS 1 N = x L x 中落在 S 中的频数， 于是，利用大数定律便知 D = ^ I N N b a M S ( - ) （８．１）
是积分I=f(x)dx的相合估计,即对于任意的E>0,当N→>∞时,有 P(b-a)M-f(x)dxE)→0 又由于N。服从参数为(P.N)的二项分布,所以我们有E=∫/x)k,即是积分 ∫(x)dx的无偏估计.此估计的方差为 va=(b-a)2M2(- A=(b-m)2(b-m(1 p) (b-a)M (b-a)M-l]=O(-) 又因为方差代表平均平方误差,所以我们就说积分的估计的误差为O(-=) 频率法所需采样量的估计 我们知道统计估计是以区间估计来给出概率误差的,即对于允许误差E,及允许以δ的 失败概率,来确定样本量N(依赖于ε,δ),使 P(-]f(x)xE)8)=P(I-EIpe)sar/ 要求m/(b-a)2M2P(-p) <δ.考虑到p(1-p)≤,我们要确定N,只需满 交(b-0M4N6就够了,也就是说只要小、(b-2)M即可但是这个采样量 的估计过于粗糙 2用中心极限定理的近似估计样本量 在N大时(这总假定满足),用中心极限定理也可以得到采样量N的估计.中心极限定 理断言,=(-0M的近似分布为正态分布N()xmri)令满足
203 是积分 ò D = b a I f (x)dx 的相合估计, 即对于任意的 e > 0, 当 N ® ¥ 时, 有 - - ò > ® b a S f x dx N N P(| (b a)M ( ) | e ) 0． 又由于 NS 服从参数为 ( p, N) 的二项分布, 所以我们有 = ^ E I ò b a f (x)dx , 即 ^ I 是积分 ò b a f (x)dx 的无偏估计. 此估计的方差为 = ^ Var I N p p b a M (1 ) ( ) 2 2 - - N b a M I b a M I b a M ) ( ) (1 ( ) ( ) 2 2 - - - = - b a M I I N [( ) ] 1 = - - ) 1 ( N = O ． 又因为方差代表平均平方误差, 所以我们就说, 积分的估计 ^ I 的误差为 ) 1 ( N O . 频率法所需采样量的估计 我们知道统计估计是以区间估计来给出概率误差的, 即对于允许误差e , 及允许以d 的 失败概率, 来确定样本量N (依赖于e,d ), 使 ò - > = - > £ ò , 要求 d e e 即可. 但是这个采样量 的估计过于粗糙. 2 用中心极限定理的近似估计样本量 在 N 大时(这总假定满足), 用中心极限定理也可以得到采样量 N 的估计. 中心极限定 理断言， I (b a)M ^ = - N NS 的近似分布为正态分布 N( ò b a f (x)dx , ) ^ Var I . 令fd 满足
dx=1-6) 于是有 P|-]f(x)x中vamr1)≤6 选取采样量N使、Vm156,即(b-03M2-Ps52,则就有 N PQ/-∫/x)dpE)≤6 所以只需取 (b-a)M- 3采样量的估计的比较 我们把用 Chebyshev不等式得到的样量的估计与用中心极限定理得到的采样量的估 计作比较,可以看到只要φ。-hn8 当0小时,一lnδ趋于∞的速度要比φ。趋于∞的速度快.所以,以用中心极限定理得到的采样量的 1.2用样本函数的平均值估计积分的 Monte carlo方法一期望法 期望法的核心思想是把积分看成某个随机变量的期望.最自然的是看成[a,b]上均匀 随机变量的期望.设η~U[a,b]([a,b上的均匀分布)则 =「f(x)dr=(b-a)E/(m) 于是对于N个独立的[a,b]上的均匀随机数,可以用矩估计 =(b-a)f()+…+f(7) N
204 ò ¥ - = fd d p e dx x 2 2 2 1 , ( 即 ò - - = - 2 2 2 1 2 1 2 d d f f d p e dx x ). 于是有 ò - ³ £ b a P I f x dx f VarI d d (| ( ) | ) ^ 2 ^ . 选取采样量 N 使 f e d £ ^ 2 VarI , 即 2 2 2 2 2 (1 ) ( ) fd e £ - - N p p b a M , 则就有 ò - ³ £ b a P(| I f (x)dx | e ) d ^ . 所以只需取 2 2 2 2 2 4 ( ) fd e b a M N - > . 3 采样量的估计的比较 我们把用 Chebyshev 不等式得到的采样量的估计与用中心极限定理得到的采样量的估 计作比较, 可以看到只要 1 2 2 dfd . 当d 小时, - ln d 趋于 ¥ 的速度要比 2 2 fd 趋于 ¥ 的速度快. 所以, 以用中心极限定理得到的采样量的 估计为最好. １．2 用样本函数的平均值估计积分的 Monte Carlo 方法 ―期望法 期望法的核心思想是把积分看成某个随机变量的期望. 最自然的是看成 [a, b] 上均匀 随机变量的期望. 设h ~ U[a,b] ( [a, b] 上的均匀分布), 则 I f (x)dx (b a)Ef (h) b a = = - ò D . (8. 2) 于是对于 N 个独立的[a, b] 上的均匀随机数, 可以用矩估计 N f f I b a N ( ) ( ) ( ) 1 ^ ^ h + + h = - L (8. 3)
作为Ⅰ=」f(x)x的估计.显然它也是无偏估计而且 1a(1)=(6-a{a(m)=(6-a31[(m)2-E(m)] dx 1 (b-a)2f( (b-aM-l=var(n) 可见期望法比频率法更有效.但是期望法多算了N次函数值,这是需要花费计算时间的 所以在全面考虑有效性与计算时间的得失时,有时人们也愿意采用频率法 1.3减少方差的技术 用 Monte carlo方法计算积分∫f(x)时,未必一定要使用均匀随机数事实上,从 ab]上取值的任意一种随机数出发,都可以得到∫f(x)的估计量而且在f(x)20时 显见,值∫(x)大的x对于积分」f(x)dx有更大的贡献由此得到直观启发,所用的随机 数的分布密度的形状越象∫(x),则就越合理.这个思想就是重要度采样法 g-采样法 假定分布密度8(x)在f(x)非零处恒正,则积分=」f(x)h=/f(x) 8(x)8(x)dx.对于 密度为g(x)的“g-随机数”s,我们有 1= EL /(c b 于是对于N个独立的g-随机数512…5N,关于积分/=f(x)tx可取估计 A(g)A f(s1) g(51) g(s 显见它也是无偏的相合估计.利用 Schwartz不等式得到 (r)x=g(x)dx
205 作为 ò = b a I f (x)dx 的估计. 显然它也是无偏估计. 而且 N Var f Var I b a [ ( )] ( ) ( ) 2 ^ ^ h = - [ ] 2 2 2 ( ) [ ( )] 1 ( ) Ef h Ef h N = b - a - 2 2 1 2 ( ) 1 ( ) I b a N dx f x N b a b a - - = - ò 2 [( ) 1 b a MI I N £ - - ( ) ^ = Var I . 可见期望法比频率法更有效. 但是期望法多算了 N 次函数值, 这是需要花费计算时间的. 所以在全面考虑有效性与计算时间的得失时, 有时人们也愿意采用频率法. 1. 3 减少方差的技术 用 Monte Carlo 方法计算积分 ò b a f (x)dx 时, 未必一定要使用均匀随机数. 事实上, 从 [a, b]上取值的任意一种随机数出发, 都可以得到 ò b a f (x)dx 的估计量. 而且在 f (x) ³ 0 时 显见，值 f (x) 大的 x 对于积分 ò b a f (x)dx 有更大的贡献. 由此得到直观启发，所用的随机 数的分布密度的形状越象 f (x) , 则就越合理. 这个思想就是重要度采样法. 1. g - 采样法 假定分布密度g (x) 在 f (x) 非零处恒正, 则积分 I = ò b a f (x)dx ò = b a g x dx g x f x ( ) ( ) ( ) . 对于 密度为 g (x) 的 “g - 随机数” V , 我们有 ] ( ) ( ) [ V V g f I = E . (8. 4) 于是对于 N 个独立的 g - 随机数 N V , ,V 1 L , 关于积分I = ò b a f (x)dx 可取估计 ] ( ) ( ) [ 1 1 1 ( ) ^ V V g f N I g D = ] ( ) ( ) N N g f V V +L+ . (8. 5) 显见它也是无偏的相合估计. 利用 Schwartz 不等式得到 ò ò = b a b a g x dx g x dx g x f x ] ( ) ( ) ( ) ( ) [ 2 ò b a g x dx g x f x ] ( ) ( ) ( ) [ 2
> f(x) v8(r)v8()dx]=[/(x)dx] 而且上式当且仅当在√图(x)=c-(0时(即g()=(x)时)达到极小值 g(x) (x).这说明a()=可∫ ,八(了8(x)d-门]的最小值在g(x)=(x) N。g(x) 时取到.又因为g(x)为密度,故c=x1 对应的Vawr()=0.此时 f(xdx (go)1 =f(x)x是无估计误差的精确值 综上讨论可知,要使方差达到最小,就应该用g0(x)=Cf(x)作为参考密度,这时用Von Neumann取舍方法必须知道参考密度g0(x),而计算g0(x)又必须知道常数c,后者却是于 积分的倒数,这样的耦合性的存在使这种方法不能真正用于计算积分.但是由此我们至少 a(g) 得到了以下的认识即只要密度g(x)的形状与被积分的函数近似用作为∫/(x 的估计,就会降低方差.这个想法还可以推广到∫(x)未必为非负函数的情形.这就是下 面的概念 定义日,1分布密度为g(x)=(3)的g一采样,称为(关于/(x)的重要 lf(x)idx 度采样( Importance Sampling) 重要度采样不能直接通过取舍原则实现.近似地实现重要度采样可以采用下一节中的 Markov链 Monte Carlo方法,这个方法同时给出连积分f(x)|ax的近似算法 在实践中人们也往往按照重要度采样的思路,灵活地寻找常用的已知类型的密度g,使 它在峰值附近与|∫(x)较接近以便达到降低估计的方差的目的 修正的重要度采样( Modified Importance Sampling) 对于g-样.假定有非负函数h(x),满足
206 ³ 2 ( ) ] ( ) ( ) [ ò b a g x dx g x f x = 2 [ ( ) ] ò b a f x dx ． 而且上式当且仅当在 ( ) ( ) ( ) g x f x g x = c 时 ( 即 g( x) = cf (x) 时 ) 达到极小值 2 [ ( ) ] ò b a f x dx . 这说明 N Var I g 1 ( ) ( ) ^ = ] ( ) ] ( ) ( ) [ [ 2 2 g x dx I g x f x b a - ò 的最小值在 ( ) ( ) 0 g x cf x D = 时取到. 又因为 ( ) 0 g x 为密度, 故 ò = b a f x dx c ( ) 1 . 对应的 ( ) 0 ( ) ^ 0 = g Var I . 此时 ò = = b a g f x dx c I ( ) 1 ( ) ^ 0 是无估计误差的精确值． 综上讨论可知，要使方差达到最小，就应该用 ( ) ( ) 0 g x cf x D = 作为参考密度, 这时用 Von Neumann 取舍方法必须知道参考密度 ( ) 0 g x , 而计算 ( ) 0 g x 又必须知道常数c ，后者却是于 积分的倒数，这样的耦合性的存在使这种方法不能真正用于计算积分. 但是由此我们至少 得到了以下的认识，即只要密度 g (x) 的形状与被积分的函数近似, 用 ( ) ^ g I 作为 ò b a f (x)dx 的估计, 就会降低方差. 这个想法还可以推广到 f (x) 未必为非负函数的情形. 这就是下 面的概念． 定义８.１ 分布密度为 ò = b a f x dx f x g x | ( ) | | ( ) | ( ) 的 g - 采样, 称为 (关于 f (x) 的) 重要 度采样（Importance Sampling)． 重要度采样不能直接通过取舍原则实现．近似地实现重要度采样可以采用下一节中的 Markov 链 Monte Carlo 方法，这个方法同时给出连积分 ò b a ｜f (x)| dx 的近似算法． 在实践中人们也往往按照重要度采样的思路, 灵活地寻找常用的已知类型的密度g , 使 它在峰值附近与 | f (x) | 较接近, 以便达到降低估计的方差的目的. 修正的重要度采样 (Modified Importance Sampling) 对于g - 采样. 假定有非负函数h(x) ，满足
a=h(x)g(x)dx>0 而且a已知.那么我们可以采用h(x)作为修正乘积因子,显见对于g-随机数g有 f(x)dh=-8(5) Eh(s) 如果我们先放弃对于估计的无偏性要求,而只要求估计的相合性,则对于N个独立的g 随机数51…5N,我们可以通过比值构造I=f(x)dx的如下的估计量 f(s1) f(n) 1=a.8(s1)8(5N) h(s1)+…+h(sN 显见它是/=∫f(x)的相合估计在某些假定下,它还是渐进无偏的,即 E(1)= 而且了保留了重要度采样的特性,即当M=c时,是1-/,于是 只要当加(y1人(小队。会降低方差注意对于给定的分布密度 g 尽量与少近似这个修正乘积因子h(x)是用来再一次降低由于密度g(x)与被积函数 g(x) ∫(x)的倍数不够像所带来的失误而设置的.而当h(x)≡1时,就退化为是g-采样,这相 当于对g-采样不再作修正 2. Markov fE Monte Carlo(MCMC 用 Markov链的样本,来对不变分布, Gibbs分布, Gibbs场,高维分布,或样本空间非 常大的离散分布等作采样,并用以作随机模拟的方法,统称为 Markov Chain Monte CarIo McM)方法.这是动态的 Monte Carlo方法.由于这种方法的问世,使随机模拟在很多领 域的计算中,相对于决定性算法,显示出它的巨大的优越性.而有时随机模拟与决定性算法 的结合使用,会显出更多的长处 至少可以用在以下几个层面
207 ò = > b a a h(x)g(x)dx 0 ， 而且a 已知．那么我们可以采用h(x) 作为修正乘积因子．显见对于 g - 随机数V 有 ( ) ] ( ) ( ) [ ( ) V V V a Eh g f E I f x dx b a = = ò . (8. 6) 如果我们先放弃对于估计的无偏性要求, 而只要求估计的相合性, 则对于 N 个独立的 g - 随机数 N V , ,V 1 L , 我们可以通过比值, 构造 ò = b a I f (x)dx 的如下的估计量 ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 N N N h h g f g f I V V V V V V a + + + + = × Ú D L L . (8. 7) 显见它是 ò = b a I f (x)dx 的相合估计. 在某些假定下, 它还是渐进无偏的, 即 E I I N = Ú ®¥ lim ( ) . 而且 Ú I 保留了重要度采样的特性，即当 ( ) ( ) ( ) g x f x h x = c 时， Ú I 就是 ò = b a I f (x)dx ．于是， 只要当h(x) 与 ( ) ( ) g x f x 近似, 就会降低方差. 注意对于给定的分布密度 g (x) ，h(x) 应选取得 尽量与 ( ) ( ) g x f x 近似. 这个修正乘积因子h(x) 是用来再一次降低由于密度 g (x) 与被积函数 f (x) 的倍数不够像所带来的失误而设置的. 而当h(x) º 1时, 就退化为是 g - 采样, 这相 当于对 g - 采样不再作修正. 2. Markov 链 Monte Carlo (MCMC) 用 Markov 链的样本, 来对不变分布, Gibbs 分布, Gibbs 场, 高维分布, 或样本空间非 常大的离散分布等作采样, 并用以作随机模拟的方法, 统称为 Markov Chain Monte Carlo (MCMC)方法. 这是动态的 Monte Carlo 方法. 由于这种方法的问世, 使随机模拟在很多领 域的计算中, 相对于决定性算法，显示出它的巨大的优越性. 而有时随机模拟与决定性算法 的结合使用, 会显出更多的长处. MCMC 至少可以用在以下几个层面:
1.用于生成较复杂的随机数: *实现对高维分布(或高维格点分布)兀取样,得到丌随机数 是实现重要度采样的一种方法.对∫(x)的重要性采样,就是取得 兀(=丌(x)= a f(x) -)随机数 lf() dy 2.实现高维积分(或项数极多的求和)的数值计算(典例是 Gibbs分布的各种泛函的 平均值的计算).对于∫(x)≥0,作以丌=丌(x)= A f(x) 为极限分布的 Markov链X f()dy 利用遍历定理可以由这个 Markov链的一条轨道,得到分布密度x(x)的估计,记为r(x).再 用作为积分/=∫/(x女的估计 丌(x) 3.用模拟方法估计最可几轨道.例如,如果模拟了100条轨道,那么就能以大概率推 断,最可几轨道就在这些轨道的邻近.当统计量的分布未知时,可以用模拟方法从频率估 计置信限 4.用被估参数的 Bayes分布(参见第9章)的取样,来估计参数 5.求复杂样本空间上函数的极值(模拟退火) 2.1 Gibbs采样法( Gibbs sampler) 1.用MCMC方法得到 Gibbs分布的样本与估计 Gibbs分布 在第6章中,我们考虑了在d维的N一格点集上的 Ising模型的 Gibbs分布 em(m,由于所涉及的状态空间(全体组态的集合)S=+1非常大(例 如,把一幅256×256个采样点的黑白图像看成一个组态,则d=2,N=256,S中有 236256(=2>260010809个元素),这就使得分母中的求和无法实际完成而MOMC 方法就是以通过构造一个以这个Gibs分布r;= 为不变分布的离散时间的 Markov链Xn(它就是 Glauber动力学中的连续时间的 Markov链的离散时间采样),作为模 拟计算的基点的.构造的 Markov链必须易于计算,所以我们要求它的概率转移速率只容许 在组态的一个格点上变动.这样的变动方式,称为Gibs方式,这种抽样方法称为Gibs采 样法,或者 Gibbs样本生成法这个 Markov链的不变分布正是此 Gibbs分布丌,我们还要 求此 Markov链的转移矩阵满足Pn—→1π.这就是说,要求Gibs分布是 Glauber 动力学的极限分布. 于是,当n大时,X的一个样本可以近似地认为取自Gibs分布的一个样本,即按此
208 1. 用于生成较复杂的随机数： * 实现对高维分布（或高维格点分布）p 取样, 得到p 随机数. * 是实现重要度采样的一种方法. 对| f (x) |的重要性采样, 就是取得 p ( ò D = = f y dy f x x | ( )| | ( ) | p( ) )随机数. 2. 实现高维积分（或项数极多的求和）的数值计算（典例是 Gibbs 分布的各种泛函的 平均值的计算）. 对于 f (x) ³ 0 , 作以p ò D = = f y dy f x x ( ) ( ) p( ) 为极限分布的 Markov 链 Xn ， 利用遍历定理可以由这个 Markov 链的一条轨道，得到分布密度p (x) 的估计，记为 ( ) ^ p x ．再 用 ( ) ( ) ^ x f x p 作为积分 ò I = f (x)dx 的估计. 3. 用模拟方法估计最可几轨道. 例如, 如果模拟了 100 条轨道, 那么就能以大概率推 断, 最可几轨道就在这些轨道的邻近. 当统计量的分布未知时, 可以用模拟方法从频率估 计置信限. ４. 用被估参数的 Bayes 分布（参见第 9 章）的取样，来估计参数． ５. 求复杂样本空间上函数的极值(模拟退火). 2. 1 Gibbs 采样法 (Gibbs sampler) 1. 用 MCMC 方法得到 Gibbs 分布的样本与估计 Gibbs 分布 在第６章中，我们考虑了在 d 维的 N - 格点集上的 Ising 模型的 Gibbs 分布 åÎ - - = S H H e e h b h b x p x ( ) ( ) , 由于所涉及的状态空间 (全体组态的集合) N d S {1, , ) { 1,1} L = - 非常大 (例 如, 把一幅 256´ 256 个采样点的黑白图像看成一个组态, 则d = 2, N = 256 , S 中有 2 ( 2 2 10 ) 256 256 2 60000 18000 16 = > > ´ 个元素), 这就使得分母中的求和无法实际完成. 而 MCMC 方法就是以通过构造一个以这个 Gibbs 分布 åÎ - - = S H H e e h b h b x p x ( ) ( ) 为不变分布的离散时间的 Markov 链 Xn (它就是 Glauber 动力学中的连续时间的 Markov 链的离散时间采样), 作为模 拟计算的基点的． 构造的 Markov 链必须易于计算, 所以我们要求它的概率转移速率只容许 在组态的一个格点上变动．这样的变动方式, 称为 Gibbs 方式, 这种抽样方法称为 Gibbs 采 样法，或者 Gibbs 样本生成法. 这个 Markov 链的不变分布正是此 Gibbs 分布p , 我们还要 求此 Markov 链的转移矩阵满足 P n ¾n®¾¥® p T 1 . 这就是说, 要求 Gibbs 分布是 Glauber 动力学的极限分布. 于是, 当n 大时, Xn 的一个样本可以近似地认为取自 Gibbs 分布的一个样本, 即按此
Markov链沿任意一条轨道充分发展,就得到 Gibbs分布的近似取样. 再则,Gtbs分布的归一化常数(称为配分常数∑c1(),是一个巨大的求和,即 个”离散的”积分.用随机模拟法计算这个"离散的”积分的最佳随机数正服从 Gibbs分布 (即重要度采样).对于Gibs分布的取样,用通常的取舍原则常常并不可行.例如,分别取 C=1,h()=em5),而参考密度p0(5)为组态空间上的均匀分布,这时e)的值常常 小得超出计算精度,而求和变量ξ的范围是庞大的组态空间,这就导致求和无法实际计 算.所以需要用 Markov链 Monte carlo方法.用MCMC方法生成了以 Gibbs分布为极限分 布的 Markov链Xn以后,由遍历定理用 Markov链的一条轨道,可给出极限分布(Gibs分 布)的估计:对于充分大的N,可令 丌=(1(xx+1)+…l(X2) (8.8) -H(5) 再用e-(5)除以Gbs分布在ξ处的估计值 作为配分函数的估计.在理论上这个 估计应该与5的取法无关,但是,在实际计算中对多个不同的组态1分别估计此和数后, 再作平均常常能降低方差. 在第6章中,我们曾给出了用 Glauber动力学构造的两个不同的连续时间的 Markov链 (对应于两个不同的转移概率速率矩阵Q,它们都以Gibs分布π为极限分布,而且都是可 逆的 较为深入的理论研究表明,使用不可逆的,且以丌为不变分布的 Markov链作 Markov Monte carlo,会加快这个极限的收敛速度然而,在另一方面这种做法又会增加计算的复杂 程度.再则,为减少估计的方差而作的努力也常会增加计算时间.这就是说,在计算中,我 们会面临难以两全的抉择.在实际中如何采取折衷,既要看问题的性质,又要参考实践的经 验,没有统一的原则 用以完成MCMC采样操作的 Markov链,可构造如下: 在第6章中,对于d维有限格点上,由具有两个自旋的组态空间上的能量函数 H()=-22 n(x)n(y)-hEn(x) 可构造如下的转移概率速率 C(x,5),(n=5) (≠ξ的其它情形) (C(x,)的两种取法各为: (x,)=eP(H(H() (8.11)
209 Markov 链沿任意一条轨道充分发展, 就得到 Gibbs 分布的近似取样. 再则，Gibbs 分布的归一化常数(称为配分常数) b (x ) h H S e - Î å , 是一个巨大的求和，即一 个”离散的”积分． 用随机模拟法计算这个＂离散的＂积分的最佳随机数正服从 Gibbs 分布 (即重要度采样)． 对于 Gibbs 分布的取样, 用通常的取舍原则常常并不可行．例如，分别取 C =1, ( ) ( ) b x x H h e - = , 而参考密度 ( ) 0 p x 为组态空间上的均匀分布，这时 bH (x ) e - 的值常常 小得超出计算精度, 而求和变量x 的范围是庞大的组态空间, 这就导致求和无法实际计 算．所以需要用 Markov 链 Monte Carlo 方法. 用 MCMC 方法生成了以 Gibbs 分布为极限分 布的 Markov 链 Xn 以后, 由遍历定理用 Markov 链的一条轨道, 可给出极限分布（Gibbs 分 布）的估计: 对于充分大的 N ，可令 ( ( ) ( )) 1 { } 1 { } 2 ^ N X N I X I N x x p x = + +L ， （８．８） 再用 bH (x ) e - 除以 Gibbs 分布在x 处的估计值 ^ ( ) x b x p H e - , 作为配分函数的估计．在理论上这个 估计应该与x 的取法无关. 但是, 在实际计算中对多个不同的组态 i x 分别估计此和数后， 再作平均常常能降低方差． 在第 6 章中, 我们曾给出了用 Glauber 动力学构造的两个不同的连续时间的 Markov 链 (对应于两个不同的转移概率速率矩阵 Q), 它们都以 Gibbs 分布p 为极限分布, 而且都是可 逆的. 较为深入的理论研究表明, 使用不可逆的, 且以p 为不变分布的 Markov 链作 Markov Monte Carlo, 会加快这个极限的收敛速度. 然而, 在另一方面这种做法又会增加计算的复杂 程度. 再则, 为减少估计的方差而作的努力也常会增加计算时间. 这就是说，在计算中，我 们会面临难以两全的抉择. 在实际中如何采取折衷，既要看问题的性质, 又要参考实践的经 验，没有统一的原则． 用以完成 MCMC 采样操作的 Markov 链，可构造如下： 在第 6 章中, 对于d 维有限格点上，由具有两个自旋的组态空间上的能量函数 = - å - å x y相邻 x H x y h x , ( ) ( ) ( ) 2 1 (h) h h h , (8. 9) 可构造如下的转移概率速率 î í ì ¹ = = h x的其它情形） x h x xh 0 ( ( , ), ( ) x C x q (8. 10) (C(x,x ) 的两种取法各为： ( ( ) ( )) ( , ) b x x x H H C x e - - = (8. 11) 或
C(x,5)= 1+e 它们决定的 Glauber动力学,分别对应的两个不同的连续时间的 Markov链,它们都以H(n) 为能量函数的Gibs分布丌为可逆不变分布,而且丌还是转移矩阵的极限分 布:P()-y=I兀.如果考虑间隔为△t的采样时间,其中Mt充分小,我们还可以进 行如下的近似计算.假定初始组态为g.在时刻M它以概率C(x,s)△t转移到组态sx,而 以概率1-△∑C(x5)停留在原来的组态,这就近似地得到了 Markov链在△M时刻所 …N}4 处的组态51·再以类似的方式得到 Markov链在2A时刻所处的组态s2,依次下去,当样 数n充分大时,组态n在这段{51,…sn}中出现的频率,就用作zn的估计丌n Gibbs分布的一些统计平均量的模拟近似 对于上面定义的,以Gibs分布为极限的 Markov链,我们用它的一段样本s1…,sN 可以给出如下的 Gibbs统计平均量 ∑f(5)-( 的模拟近似 ∑f(s;) (8.13) 它是F的无偏估计.不难证明:存在数V(f,P,丌)<∞,使F的方差满足 Ha()=(,Px)+0)20-C(P)m 13‖ 其中|∑f(5),而C(P)为此 Markov链转移矩阵P的 Dobrushin收缩数 Gibs采样法还可以有效地用于 Bayes方法中的后验密度的采样的模拟计算(见第9章) 2.2 Metropol is采样法 (Metropol is sampler) 1. Gibbs分布的采样的 Markov链 Metropolis方法 为了突出主要思想,下面我们把组态空间(状态空间)S={-1,]}…)“简单地记为 {1,…,K}.对于组态空间上给定的分布丌, Metropolis构造了以 210
210 ( ( ) ( )) 1 1 ( , ) b x x x H H x e C x - - + = ) ． (8. 11)’ 它们决定的Glauber 动力学, 分别对应的两个不同的连续时间的Markov链, 它们都以H (h) 为能量函数 的 Gibbs 分布 p 为可逆不变分布 , 而且 p 还 是转移矩阵的 极限分 布: P (t) ¾t®¾¥® p T = 1 . 如果考虑间隔为Dt 的采样时间, 其中Dt 充分小, 我们还可以进 行如下的近似计算．假定初始组态为V . 在时刻Dt 它以概率C(x,V )Dt 转移到组态 x V , 而 以概率 å Î - D d x N t C x {1, , } 1 ( , ) L V 停留在原来的组态, 这就近似地得到了 Markov 链在Dt 时刻所 处的组态 1 V . 再以类似的方式得到 Markov 链在2Dt 时刻所处的组态 2 V ，依次下去. 当采样 数n 充分大时, 组态h 在这段{ , , } 1 n V L V 中出现的频率, 就用作ph 的估计 ^ p h . Gibbs 分布的一些统计平均量的模拟近似 对于上面定义的, 以 Gibbs 分布为极限的 Markov 链, 我们用它的一段样本 N V , ,V 1 L , 可以给出如下的 Gibbs 统计平均量 å å Î - - Î = S H H S e f e F x b x b x x x ( ) ( ) ( ) (8. 12) 的模拟近似 n f F n i å i = = 1 ^ (V ) , (8. 13) 它是 F 的无偏估计. 不难证明: 存在数V( f , P,p ) < ¥ , 使 ^ F 的方差满足 ) 1 ( ( , , ) ( ) ^ n o n V f P Var F = + p C P n f (1 ( )) 13 || ||2 - £ , (8. 14) 其中 = å x || f || f (x ), 而C(P) 为此 Markov 链的转移矩阵P 的 Dobrushin 收缩数. Gibbs 采样法还可以有效地用于 Bayes 方法中的后验密度的采样的模拟计算（见第９章）． 2. 2 Metropolis 采样法 (Metropolis sampler) １．Gibbs 分布的采样的 Markov 链 Metropolis 方法 为了突出主要思想, 下面我们把组态空间(状态空间) N d S {1, , ) { 1,1} L = - 简单地记为 {1,L,K}. 对于组态空间上给定的分布p , Metropolis 构造了以
P=Pmm(1,),(≠) (8.15) P 为转移的时齐 Markov链,其中P=(P)是一个对称的互通转移矩阵, 称为预选矩阵,或访问方案,使用它是为了减少或简化状态间的连接,以加快 Markov链的 分布向不变分布收敛的速度.显见兀是它的可逆分布(注意在 Gibbs分布情形,状态i,j体 现为组态ξ,η,于是在计算(8.15)式的转移概率时,就只需算比值 个=e(H(H(n)(2≠m),而并不需要计算配分函数, Glauber动力学的构架,也 正是用了这一点).由这个有限状态 Markov链的互通性,我们有 Pn)=1丌 因此,在时间发展充分长以后,我们可以用 Metropol is的 Markov链所处的状态,作为按分 布π取的样本.也就是说,与 Gibbs采样法一样, Metropolis方法也给出了在计算机上模 拟π一随机数的一个算法. Metropolis提出的这种采样法,称为 Metropol is采样法.它与 Gibs采样法的不同处在于,对于 Metropolis样法而言,任意两个组态ξ,n,只要预选 概率P>0就可以转移 Metropolis采样在时刻n的更新x→xm可具体采取如下的操作 (1)设当前为时刻n,取的状态为x(m=i.对它作随机扰动,即取一个分布为 (P1…,Pκ)的随机数,设为j (2)若—≥1,则将状态更新为x=j:否则进行(3) (3)独立地取一个U[0]随机数U,如果U≤—,则将状态更新为x)=j:否则状态 不更新,即令x(m+=i (请读者证明,如此由i到j的转移的可能性恰是(8.15)式规定的转移概率) P的对称性并非必要。理论分析指出,经过适当的选取(研究矩阵P的第二个特征值) 使用非对称的P可能加快收敛速度.对于非对称的预选矩阵P, Metropolis.样法所构造 的 Markov链的转移应取下式: 211
211 min( 1, ), ( ) ~ p p j i i j ij = ij ¹ p p (8. 15) ( å¹ = - j i ii ij p 1 p ) 为转移的时齐 Markov 链, 其中 ~ P ( ) ~ = pij 是一个对称的互通转移矩阵, 称为预选矩阵，或访问方案，使用它是为了减少或简化状态间的连接, 以加快 Markov 链的 分布向不变分布收敛的速度. 显见p 是它的可逆分布 （注意在 Gibbs 分布情形，状态i, j 体 现为组态 x ,h ， 于是在计算 （ ８ ． １ ５ ） 式的转移概率时 , 就只需算比值 ( ) ( ( ) ( ) x h p p b x h h x = ¹ - H -H e , 而并不需要计算配分函数． Glauber 动力学的构架，也 正是用了这一点）. 由这个有限状态 Markov 链的互通性，我们有 P n ¾n®¾¥® p T = 1 . 因此, 在时间发展充分长以后, 我们可以用 Metropolis 的 Markov 链所处的状态, 作为按分 布p 取的样本. 也就是说， 与 Gibbs 采样法一样，Metropolis 方法也给出了在计算机上模 拟p - 随机数的一个算法. Metropolis 提出的这种采样法, 称为 Metropolis 采样法. 它与 Gibbs 采样法的不同处在于, 对于 Metropolis 采样法而言，任意两个组态x ,h ， 只要预选 概率 0 ~ pxh > 就可以转移. Metropolis 采样在时刻n 的更新 (n) ® (n+1) x x 可具体采取如下的操作: (1) 设当前为时刻 n , 取的状态为 x i n = ( ) . 对它作随机扰动, 即取一个分布为 ( , , ) ~ 1 ~ i iK p L p 的随机数，设为 j ; (2) 若 ³ 1 i j p p ，则将状态更新为 x j n = ( +1) ；否则进行（３）； (3) 独立地取一个U[0,1] 随机数U , 如果 i j U p p £ , 则将状态更新为x j n = ( +1) ；否则状态 不更新，即令 x i n = ( +1) . （请读者证明，如此由i 到 j 的转移的可能性恰是（８．１５）式规定的转移概率）． ~ P的对称性并非必要. 理论分析指出, 经过适当的选取 (研究矩阵 ~ P的第二个特征值), 使用非对称的 ~ P 可能加快收敛速度. 对于非对称的预选矩阵 ~ P, Metropolis 采样法所构造 的 Markov 链的转移应取下式：
