《应用随机过程教程》教学资源（参考资料）与在算法和智能计算中的应用——第8章 Monte Carlo与Markov Chain Monte Carlo（MCMC）方法

龚光鲁,钱敏平著应用随机过程教程及其在算法与智能计算中的应用 清华大学出版社,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)<6 1用 Chebyshev不等式估计样本量 最粗略的确定样本量N的方法,是通过 Chebyshev不等式来估计.即由 P(/- f(r)dxp>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 ), 使 ò - > < b a P(| I f (x)dx | e ) d ^ . 1 用 Chebyshev 不等式估计样本量 最粗略的确定样本量 N 的方法, 是通过 Chebyshev 不等式来估计. 即由 2 ^ ^ ^ ^ (| ( ) | ) (| | ) e e e Var I P I f x dx P I E I b a - > = - > £ ò , 要求 d e e < - - = 2 2 2 2 ^ (1 ) ( ) N p p b a M Var I . 考虑到 4 1 p(1- p) £ , 我们要确定 N , 只需满 足 d e < - 2 2 2 4 1 ( ) N b a M 就够了, 也就是说只要 2 2 2 4 ( ) de b a M N - > 即可. 但是这个采样量 的估计过于粗糙. 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不等式得到的样量的估计与用中心极限定理得到的采样量的估 计作比较,可以看到只要φ。<<1,用中心极限定理得到的采样量的估计,要比用 Chebyshev不等式得到的样量的估计小得多 [注]还可以应用概率理论中的大偏差速率函数来求得采样量的估计可以证明由大偏差速率函数 能得到采样量的估计 N>-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 << , 用中心极限定理得到的采样量的估计, 要比用 Chebyshev 不等式得到的采样量的估计小得多. ［注］ 还可以应用概率理论中的大偏差速率函数来求得采样量的估计. 可以证明由大偏差速率函数 能得到采样量的估计 2 2 ln e - d N > . 当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) ，满足