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 构造了以