2. 4 蒙特卡洛计算中减少方差的技巧 蒙特卡洛求积分的方差为 {}nfV= 2 σ . 其中 V 为被积函数的方差。 {}f f 公式反映出增加随机点数时，蒙特卡洛计算的精度可以得 到改善，但是精度提高非常缓慢。因此用增加蒙特卡洛计算的随 机投点数来提高精度总是耗费大量的机时。 n 另一个减少计算结果误差的途径是减少的方差V。 f {}f 重要的减少方差V的技巧。 {}f 一、 分层抽样(stratified sampling) 数学上，分层抽样是基于黎曼积分的特性： ， 0 ∫∫∫ +== 1 00 1 )()()( a a dxxfdxxfdxxfI 1<< a 将积分区域划分成小区域是在数值积分中常用的技巧。 蒙特卡洛的分层抽样技巧的抽样步骤： （1） 将积分区间（或空间）划分为不相交的子区间（或子空间）； 然后在第个子区间（或子空间）内抽取个随机点。 i i n （2） 如果将子区间长度（或子空间体积）记为{}i，我们将子区 间（或子空间）内所有点上的函数值乘上权重因子{} i ni之后迭加 起来，就得到该积分在这个子区间的积分估计值； （3） 将所有子区间的积分值迭加起来，就得到在整个区间的积 分估计值。这样得到积分的估计值的方差为： 1 {} {} (){} {} 2 2 1 2 j j j n i ij j j n j xfV n j I j σ ∑∑∑ = ? ? ? ? ? ? ? ? = = V 如果适当选择子区间{}i的大小和随机点数，就可以使计算 结果的方差得以减小。这里选择 i n {}i和的关键是要了解被积函数 在子区间内的特性。如果 i n f {}i的划分和的选择都不适当，也可 能造成更大的误差。 i n 如果我们不管被积函数的特性，而简单地将积分区域划分成相 等的子区间{}，并在各子区间内抽取相同数量的随机点数。这 种处理方法称为均匀分层抽样法。 i i n 求一维定积分的问题，比较一下用分层抽样法和用原始蒙特卡 洛方法计算得到的方差。 设所求积分为： . ∫ = 1 0 )( dxxfI 数学上可以写成 I . ∫∫ ≡= 1 0 1 1 0 )()()( dxxfxgdxxf 在[0,1]区间插入J 个点，其中10 10 =<???<<= J xxx。令 1 1 1 11 () , ()/ , , () 0, () () ,( 1,2, , ) j j j j x j x jj j j x jj x pfxdx f xpx xx fx Igxfxdxj ? ? ? ? = ? ? ≤< ?? = ?? ? ? ? = = ??? ? ? ∫ ∫ 其它 J 在上面的公式中，显然有关系式 . ∑ = = J j jj IpI 1 如果用分层抽样蒙特卡洛方法计算的积分值, 在第j个子区间上 2 以 )(xf j 分布密度函数抽取个简单子样 j n ),,2,1( Jjx ij ???= ，则积分的无偏 估计值为 j p 11 2 j σ ∫ ? j j x x ? j I 2 2 2 j j j σ 1 n i j ∑ = ? ) {}x i )( =dx dxxf 1 2 )( I J + []fIII jj )() 2 ?+? j ) 2 j I( ∑∑∑ === ? ? ? ? ? ? ? ? == J j J j n i ij j jjJ j xg n IpI 1 )( 1 . 令第j区间积分的方差为，根据方差的定义我们有关系式 {}== jijj dxxfxgxgV 1 22 )()()(σ . 则得到分层抽样计算结果的方差 { } J IV 为： {} {} 1 2 1 2 )( 1 J j ij j J j jJ n p xgV n pI ∑∑ == ==V . 如果用通常的原始蒙特卡洛方法计算，以分布密度函数抽取 N 个简单子样，则积分的无偏估计值为： )( 1 xf ∑ = = N i i xg N I 1 ( 1 它的方差为： {} N gV N I g N i 2 1 2 1 σ ≡= ∑ = V 其中又可以表示为 2 g σ [] [] ∑ ∫∫ = ? ??≡ J j x x g j j IxgxfIxg 1 1 2 1 0 2 1 )()()(σ []dxxfIIxgp j J j x x Jj j j )()( 1 2 1 ∑ ∫ = ? ??= dxxxgIIIxgp j J j x x jj j j )()()((2)(( 1 1 ∑ ∫ = ? ?+?= . ∑∑ == ?+= J j j J j jj Ipp 1 2 1 2 )σ 设分层抽样法的总抽样数为N ，我们有 J nnnN +???++= 21 . 3 比较这两种方法计算出的结果的方差，我们有 {} {} 2 1 2 1 22 )( 1 j J j j j J j jjjjJ n p IIpp N IVI σσ ∑∑ == ? ? ? ? ? ? ? ?+=?V 2 1 2 1 )( 11 IIp Nn p N p j J j jj J j j j j ?+ ? ? ? ? ? ? ? ? ? ∑∑ == σ= . 公式的右边第二项显然是大于零的量。第一项的正负则是取决于 分层抽样时子区间的划分和子区间内的抽样点数。 j n 如果上式的值大于零，则分层抽样计算积分的方差小于采用 原始蒙特卡洛方法的方差。若取 Nn p j j 1 = ，即n jj Np=，此时公式 (2.4.14)中第一项为零，公式（2.4.14）总是大于零。这就意味 着按比例的分层抽样的方差比原始蒙特卡洛方法小。这样的分层 抽样方法具有实用意义。 如果采用均匀分层抽样方法，将[0，1]区间分成J 个相等的 子区间，每个子区间内抽取的点数 J N j =n，并且这些点是均匀分 布的，即 J pxf j 1 ,1)( 1 == ，这时公式(2.4.14)中的第一项也为零，因 而(2.4.14)式的值总是正的。 由此我们也可以看出：均匀分层抽样法是一个减小方差的保 险方法。 二、 重要抽样法(importance sampling) 重要抽样法的原理起源于数学上的变量代换方法的思想，即 ∫∫ . ∫ == 1 0 1 0 1 0 )( )( )( )( )( )( )( xdG xg xf dxxg xg xf dxxf 4 此时随机点的选择不再是均匀的，而是以分布函数G分布的。 新的被积函数为乘以权重 1 。公式（2.4.15）中 )(x )(xf )(/ xg dx xdG )( ) =xg( 。 这里称为偏倚分布密度函数。该方法使原本对的抽样，变 成由另一个分布密度函数 )(xg )( xf ) ) ( * xf ( ( ) xg xf ≡ 中产生简单子样，并附带一个 权重。这种方法也称为偏倚抽样法。 )(xg 公式右边积分中被积函数的方差为{ }gf / f V。如果选择恰当， 并使它在积分域内的函数曲线形状与接近，则该方差可以变得 很小。 )(xg 函数应当满足如下条件： )(xg （1）应当是个分布密度函数。 )(xg （2）不应在积分域内起伏太大，使之尽量等于常数， 以保证方差V比V小。 )(/)( xgxf {}gf / {}f （3）分布密度函数所对应的分布函数)(xg ( )Gx能够比较方便 地解析求出。 (4) 能方便地产生在积分域内满足分布函数( )Gx分布的随机 点。 如能按上述条件找到函数,我们就可以依下列步骤求积 分： )(xg (1) 根据分布密度函数产生随机点x . 例如采用反函数法。 )(xg (2) 求出各抽样点x 的函数值，并将所有点上的该函数 值迭加起来，再除以抽样点数就得到积分结果。 )(/)( xgxf n 也可以采用作为分布密度函数，利用舍选法来舍去)(/)( xgxfw ≡ 5 或接受个随机点的x 的值。用此方法时，应至少可以事先判断 出的最大值。当然最好能从的函数中，推导出。 w )(/)( xgxf max w x N , 2 ?? {f )(xg )(xg {gf / } }g/ 上述讨论可以很容易地推广到更高维的积分中。但是要注 意如下两个方面的问题： null 在产生随机向量的所有分量后，再用舍选法往往更快，效率 更高。 x G null 在计算)(/)( xgxf GG 值之前，做随机变量 N xx ,,, 21 ???到的变 换有时是很有用的。这时需要将雅可比行列式 yyy ,,, 21 ??? ),,,( 1 NN yxxx ??? ,,(/) 21 yy ???包括在权重因子内。 重要抽样法无疑是蒙特卡洛计算中最基本和常用的技巧之 一。它无论在提高计算速度和增加数值结果的稳定性方面都有很 大的潜力。 局限性： （1） 能寻找出某分布密度函数，并能解析求出其对应的分 布函数G的情况并不多。当然我们也可以用数值计算方法求出 ，但通常这样处理不灵活，运算速度也慢，而且结果也不准 确。 )(x )(xG （2） 当所选择的在某点函数值为零或很快趋于零时（如高 斯分布），这时在该点的数值计算是十分危险的。其方差V可 能趋于无穷大。即使是在某点上函数不为零，但其值很小时， 方差 V 也可能很大。这一问题采用通常的从样本点估计方差的 方法却不一定能检查出来。这种情况会使计算结果不稳定。 )(xg 6 三、 控制变量法（相关抽样法）(control variates) 控制变量法利用数学上积分运算的线性特性： [] ∫ ∫∫ +?= dxxgdxxgxfdxxf )()()()( 选择函数时要考虑到：g在整个积分区间都是容易精确算出, 并且在上式右边第一项的运算中对( )(xg )(x )gf ?积分的方差应当要比第 二项对积分的方差小。 f 在应用这种方法时，在重要抽样法中所遇到的，当趋于 零时，被积函数趋于无穷大的困难就不再存在，因而计算出 的结果稳定性比较好。 该方法也不需要从分布密度函数， 解 析求出分布函数G。 由此我们可以看出选择所受到的限制 比重要抽样法要小些。 )(xg )(xg )( gf ? )(x )(xg 四、对偶变量法(antithetic variates) 通常在蒙特卡洛计算中采用互相独立的随机点来进行计算。 对偶变量法中却使用相关联的点来进行计算。它利用相关点 间的关系可以是正关联的，也可以是负关联的这个特点。 两个函数值和之和的方差为 1 f 2 f {}{}{} { }( ) { }( ){ }2 fEffEfEfVfVffV 22112121 ??++=+ . 如果我们选择一些点，它们使和是负关联的。这样就可以使上 式所示的方差减小。当然这需要对具体的函数和有充分的了 解。 1 f 2 f 1 f 2 f 7 例 已知是一个单调递增的函数，现求积分 )(xf . ∫ = 1 0 )( dxxfI 解 首先，按通常的方法在积分域[0，1]区间上产生均匀分布的 随机点集{。计算对应每个点的函数} i x i x [ ] 2/)1()( ii xfxf ?+的值，再将 所有点上的函数值迭加起来，除以总的随机点数，则得到(2.4.18) 式的积分值。即 [] ∑ = ?+≈ N i ii xfxf N I 1 2/)1()( 1 . 这种做法与通常的蒙特卡洛计算中将的值迭加起来不相同。 由于的单调递增性， [ )( i xf ])(xf 2/)1()( ii xfxf ?+ 的值应当比单个点的函数值 更接近于常数。因而方差也小些。 )(xf i 这实际上是采用了和)(xf )1( xf ?的积分期望值的平均值作为 结果。由于采用相同的随机数列{ } i x 1( xf ，使得和两个函数 高度负关联，因而方差比和 )(xf )1( xf ? )(xf )?两者各自积分的方差之和要 小。 8