2.4蒙特卡洛计算中嫩少方差的技巧 象特卡洛求积分的方塾为 n 其中v为被积函数/的方差。 公式风映出增加随机点教n时,象譬卡涪计犷的崭度可以得 到改善,但是精度提高非常缳幔。因此用增加象卡涪讣犷的 杋投氚数来提高贛度总是耗费大量的机时。 另一个藏少计算结果误的滤径是减少∫的方v}o 要的藏少方}的抆巧。 分暴抽祥( stratified sampling) 数学上,分屡抽样是基于黎曼积分的特性: =x)减k=(x)+J(x,0a 将积分区域划分威小区域是在数值积分中常用的技巧。 象特卡洛的分层抽样技巧的抽样步骤: (1)将积分区间(或空间)划分为不相交的子区间(或子空间); 然后在第个子区间(或子空间)内抽取n个随机点。 (2)如果将子区间长度(或子空间体积)记为!},我们将子区 间(或子空间)内所有点上的函数值乘上权量因子}n之后加 起来,就得到该积分在这个子区间的积分估计值; (3)将所有子区间的积分值选加起來,就得到在个区间的积 分佔计值。这祥讣到积分的佑计值的方为:
2. 4 蒙特卡洛计算中减少方差的技巧 蒙特卡洛求积分的方差为 = V{ }f n 2 σ . 其中V{f }为被积函数 f 的方差。 公式反映出增加随机点数 时，蒙特卡洛计算的精度可以得 到改善，但是精度提高非常缓慢。因此用增加蒙特卡洛计算的随 机投点数来提高精度总是耗费大量的机时。 n 另一个减少计算结果误差的途径是减少 f 的方差V{f }。 重要的减少方差V{f }的技巧。 一、 分层抽样(stratified sampling) 数学上，分层抽样是基于黎曼积分的特性： ， 0 ∫ ∫ ∫ = = + 1 0 0 1 ( ) ( ) ( ) a a I f x dx f x dx f x dx < a <1 将积分区域划分成小区域是在数值积分中常用的技巧。 蒙特卡洛的分层抽样技巧的抽样步骤： （1） 将积分区间（或空间）划分为不相交的子区间（或子空间）； 然后在第i个子区间（或子空间）内抽取ni个随机点。 （2） 如果将子区间长度（或子空间体积）记为{i}，我们将子区 间（或子空间）内所有点上的函数值乘上权重因子{ } ni i 之后迭加 起来，就得到该积分在这个子区间的积分估计值； （3） 将所有子区间的积分值迭加起来，就得到在整个区间的积 分估计值。这样得到积分的估计值的方差为： 1
如果远当选择子区间的大小和随机点数n,就可以计算 结果的方離得以藏小。这里选择!和n的关是要了解被积函数 ∫在子区间内的特性。如果的划分和n的选择都不适当,也可 能造成大的误。 如果们不被积函教的兮性。而筒地将积分区域灲分相 普的子区间引,外在各子区间内抽取相同数量的随机点数n。这 种处理方法欷为均訇分层抽祥法。 求一定积分的问题。比校一下用分层抽祥法和用原始象兮卡 洛方法计犷到的方。 设所求积分为 I=f(x)dx 教学上可以写成 =f(x)dtx≡g(x)f(x)dr 在[0,1区间播入J个泉,其中0=x0<x1<…<x=1。令 p,='f()dr f()=()lo<x 其它 l=g(x)(x)d(j=12…,J) 在上面的公式中,显幾有关系式 =∑p 如果用分层抽样象管卡洛方法计算的积分值,在第j个子区间上
{ } { } { } ( ) { } 2 2 1 2 j j j n i ij j j n j V f x n j I j ∑ ∑ = ∑ σ = = V 如果适当选择子区间{i}的大小和随机点数 ，就可以使计算 结果的方差得以减小。这里选择 ni {i}和 的关键是要了解被积函数 在子区间内的特性。如果 ni f {i}的划分和 的选择都不适当，也可 能造成更大的误差。 ni 如果我们不管被积函数的特性，而简单地将积分区域划分成相 等的子区间{ }，并在各子区间内抽取相同数量的随机点数 。这 种处理方法称为均匀分层抽样法。 i ni 求一维定积分的问题，比较一下用分层抽样法和用原始蒙特卡 洛方法计算得到的方差。 设所求积分为： . ∫ = 1 0 I f (x)dx 数学上可以写成 I . ∫ ∫ = ≡ 1 0 1 1 0 f (x)dx g(x) f (x)dx 在[0,1]区间插入 J 个点，其中0 1 = x0 < x1 < ⋅ ⋅⋅ < xJ = 。令 1 1 1 1 1 ( ) , ( )/ , , ( ) 0, ( ) ( ) ,( 1, 2, , ) j j j j x j x j j j j x j j x p f x dx f x p x x x f x I g x f x dx j − − − = ≤ < = = = ⋅⋅⋅ ∫ ∫ 其它 J 在上面的公式中，显然有关系式 ∑ . = = J j j j I p I 1 如果用分层抽样蒙特卡洛方法计算的积分值, 在第 j 个子区间上 2
以(x)分布鲁度函数抽取n个简单子样x(=12…/),则积分的无 佑计值为 =∑Pl=∑ g(,) 令第j区间积分的方为d,根据方的定义我们有关系式 a;=vig(x,))=a(x)f, (x)cdr-1; 则得到分层抽祥计算结果的方鎏为 }=∑p1∑v(x}=∑p 如果用遢常的原給靈特卡洛方法计算,以分布密度函数f(x)抽取 N个筒单子样,则积分的无偏估计值为 g(x,) 它的方塾为 vig(x 其中2又可以表示为 :=C(x)-1(x)d=∑区(x) f(x)da ∑P∫"[区(x)-b+1-1f(x ∑P∫”kgx)-1)2+201-/Xg(x)-1)+(,-1) P P 设分层抽样法的总抽样教为N。我们有 N
以 f (x) j 分布密度函数抽取n j 个简单子样x ( j 1,2, , J ) ij = ⋅⋅⋅ ，则积分的无偏 估计值为 p j 1 1 2 σj ∫ − j j x x − j I 2 2 2 j j j σ 1 n i j ∑= ⋅ ) { } xi ( ) dx = f x dx 1 2 ( ) I J + [ ] I I I f j j ) ( ) 2 j ) − + − 2 j (I ∑ ∑ ∑ = = = = = J j J j n i ij j J j j j g x n I p I 1 ( ) 1 . 令第 j 区间积分的方差为 ，根据方差的定义我们有关系式 j = { } ij = j V g x g x f x dx 1 2 2 σ ( ) ( ) ( ) . 则得到分层抽样计算结果的方差V {I J }为： { } { } 1 2 1 2 ( ) 1 J j ij j J j J j n p V g x n I ∑ p ∑ = = V = = . 如果用通常的原始蒙特卡洛方法计算，以分布密度函数 抽取 N 个简单子样，则积分的无偏估计值为： ( ) 1f x ∑= = N i i g x N I 1 ( 1 它的方差为： { } N V g N I g N i 2 1 2 1 σ = ∑ ≡ = V 其中 又可以表示为 2 σ g [ ] ∑ [ ] ∫ ∫ = − ≡ − − J j x x g j j g x I f x g x I 1 1 2 1 0 2 1 σ ( ) ( ) ( ) p [ ] g x I I f x dx j J j x x j J j j ( ) ( ) 1 2 ∑ 1 ∫ = − = − − p g x I I I g x x dx j J j x x j j j j ( ( ) 2( )( ( ) ( ) 1 ∑ 1 ∫ = − = − + − ∑ ∑ . = = = + − J j j J j j j p p I 1 2 1 2 σ ) 设分层抽样法的总抽样数为 N ，我们有 J N = n + n + ⋅⋅⋅ + n 1 2 . 3
比較这丌神方法计算出的结果的方差,我们有 -}=∑p+p(-0)∑ N2P( 公式的右第二项显然是大于导的量。第一项的正负则是取决子 分层抽样时子区间的划分和子区间内的抽禅点数n 如果上式的值大于零,则分层抽样计算积分的方差小于采用 原始蒙特卡洛方法的方。考取P=1,即n=N,此时公式 (2.4.14)中算一项为粵,公式(2.4.14)总是大于粵。这就意味 着换比例的分层抽祥的方鎏比原娢象管卡洛方法小。逭样的分层 抽禅方法具有实用定义。 如果釆用均訇分层抽禅方法,将[0,1区间分成J个相奇的 子区间,每个子区间内抽取的点数n=N,并且这些点是均勻分 布的,即(x)=1p=1,这时公式(2.4.14)中的第一项也为零,因 而(2.4.14)式的值总是正的。 由此们也可以看出:均訇分层抽样法是一个减小方塾的保 险方法。 二、量要抽桦法( importance sampling) 量抽法的原理起淑于数学上的变量代换方法的思丸,即 f(x)=(x) dG(x)
比较这两种方法计算出的结果的方差，我们有 { } { } 2 1 2 1 2 2 ( ) 1 j J j j j J j J j j j j n p p p I I N I V I ∑ σ ∑ σ = = − V − = + − 2 1 2 1 ( ) 1 1 p I I n N p N p j J j j j J j j j j + − ∑ − ∑ = = = σ . 公式的右边第二项显然是大于零的量。第一项的正负则是取决于 分层抽样时子区间的划分和子区间内的抽样点数n j 。 如果上式的值大于零，则分层抽样计算积分的方差小于采用 原始蒙特卡洛方法的方差。若取 n N p j j 1 = ，即 n j = Np j ，此时公式 (2.4.14)中第一项为零，公式（2.4.14）总是大于零。这就意味 着按比例的分层抽样的方差比原始蒙特卡洛方法小。这样的分层 抽样方法具有实用意义。 如果采用均匀分层抽样方法，将[0，1]区间分成 J 个相等的 子区间，每个子区间内抽取的点数 J N n j = ，并且这些点是均匀分 布的，即 J f x p j 1 ( ) 1, 1 = = ，这时公式(2.4.14)中的第一项也为零，因 而(2.4.14)式的值总是正的。 由此我们也可以看出：均匀分层抽样法是一个减小方差的保 险方法。 二、 重要抽样法(importance sampling) 重要抽样法的原理起源于数学上的变量代换方法的思想，即 ∫ ∫ . ∫ = = 1 0 1 0 1 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) dG x g x f x g x dx g x f x f x dx 4
此时隨机点的选择不开是均句的,而是以分布函数G(x)分布的。 新的被积函数为f(x)乘以权量1g(。谷式(2.4.15)中g(=o 这里g(x)称为偏倚分布度函数。该方法使原本对f(x)的抽样,变 成由另一个分布密度函数广(x)=中产生简单于样,并附带一个 权量g(x)。这种方法也称为偏抽法。 公式右边积分中被积函数的方差为v/g}o如景g(x)选择恰当, 并伩也在积分城内的函教曲能形状与∫接抚,则该方差可以变得 函数g(x)砬当满足如下永件: (1)g(x)应当是个分布密度函数。 、(2)()g()不应在积分城内爬伏太大,傀之尽量等子常瓢, 保证才些/g)比小。 3)分布密度函数g(x)所对应的分布函数G(x)能够比较方良 地解析求出 (4)能方便地产生在积分域内滴足分布函数G(x)分布的随机 点 如能换上述永件找到函数g(x),我们就可以依下列步求积 (1)根据分布響度函数g(x)广生隨机痕x.例如釆用汊函数法。 (2)求出各抽桦点x的函数值()/g(x),并将所有点上的该函数 寶逸加起来,再除以抽样点数n就得到积分绪果。 也可以采用≡f(x)/g(x)作为分布密度函数,利用會选法来會去
此时随机点的选择不再是均匀的，而是以分布函数G 分布的。 新的被积函数为 乘以权重1 。公式（2.4.15）中 (x) f (x) / g(x) dx dG(x) g(x) = 。 这里 称为偏倚分布密度函数。该方法使原本对 的抽样，变 成由另一个分布密度函数 g(x) f ( x ) ) ) ( * f x ( ( ) g x f x ≡ 中产生简单子样，并附带一个 权重g(x)。这种方法也称为偏倚抽样法。 公式右边积分中被积函数的方差为 {f / g} f V 。如果 选择恰当， 并使它在积分域内的函数曲线形状与 接近，则该方差可以变得 很小。 g(x) 函数g(x)应当满足如下条件： （1）g(x)应当是个分布密度函数。 （2） 不应在积分域内起伏太大，使之尽量等于常数， 以保证方差V 比V 小。 f (x)/ g(x) { } f / g {f } （3）分布密度函数 g(x)所对应的分布函数G x( )能够比较方便 地解析求出。 (4) 能方便地产生在积分域内满足分布函数G x( )分布的随机 点。 如能按上述条件找到函数 ,我们就可以依下列步骤求积 分： g(x) (1) 根据分布密度函数g(x)产生随机点 x . 例如采用反函数法。 (2) 求出各抽样点 x 的函数值 ，并将所有点上的该函数 值迭加起来，再除以抽样点数 就得到积分结果。 f (x)/ g(x) n 也可以采用w ≡ f (x)/ g(x)作为分布密度函数，利用舍选法来舍去 5
或接受个随机点的x的值。用此方法时,应至少可以寡先判断 出w的最大寶。当辚最好能从f(x)/g(x)的函数中,推导出wno 上述讨论可以很率易地推广到夏商錐的积分中。但是滤 危如下两个方面的同题: 在产生随机向量的所有分量后,再用合选法佳往萸快,效率 頁高。 ●在计算f(x)/g()篁之前,散随机变量x1,x2…x到y1y2y的变 换有时是很有用的。逭时猾要将雅可比行列式 (x1,x2…x)/On,y2;yx包括在权置因子内。 要抽样法无疑是Ω卡洛计犷中最基本和常用的抆巧之 它无论在提高计箕遠度和增加数值结录的定性方面都有很 大的潜力。 局限性 (1)能寻找出分布鲁度函数g(x),并能解析求出其对应的分 布函数G(x)的情况外不多。当燚我们也可以用教值讣算方法求出 G(x),但通常这样处理不灵活,還遠度也慢,而且結果也不准 (2)当所选择的g(x)在点函数值为零或很快于零时(如高 斯分布),这时在该点的数值计算是十分危险的。其方/g)可 能趋于无穷大。即良是在某点上函数g(x)不为零,但其值很小时, 方差v{1g}也可能很大。这一问题釆用通常的从样本京估计方的 方法卸不一定能检查出来。这种情况会良计算绪景不獍定
或接受个随机点的 x 的值。用此方法时，应至少可以事先判断 出w的最大值。当然最好能从 f (x)/ g(x)的函数中，推导出wmax。 x N , 2 ⋅⋅ {f g(x) g(x) {f / g} / g} 上述讨论可以很容易地推广到更高维的积分中。但是要注 意如下两个方面的问题： z 在产生随机向量 的所有分量后，再用舍选法往往更快，效率 更高。 x G z 在计算 f (x)/ g(x) G G 值之前，做随机变量 N , x , , x 1 2 ⋅⋅⋅ 到 的变 换有时是很有用的。这时需要将雅可比行列式 y , y , , y 1 2 ⋅⋅⋅ ( , , , ) 1 N N ∂ x x ⋅ x )/ ( , , ⋅ y 1 2 ∂ y y ⋅⋅ 包括在权重因子内。 重要抽样法无疑是蒙特卡洛计算中最基本和常用的技巧之 一。它无论在提高计算速度和增加数值结果的稳定性方面都有很 大的潜力。 局限性： （1） 能寻找出某分布密度函数 ，并能解析求出其对应的分 布函数G 的情况并不多。当然我们也可以用数值计算方法求出 ，但通常这样处理不灵活，运算速度也慢，而且结果也不准 确。 (x) G(x) （2） 当所选择的 在某点函数值为零或很快趋于零时（如高 斯分布），这时在该点的数值计算是十分危险的。其方差V 可 能趋于无穷大。即使是在某点上函数 不为零，但其值很小时， 方差V 也可能很大。这一问题采用通常的从样本点估计方差的 方法却不一定能检查出来。这种情况会使计算结果不稳定。 g(x) 6
三、控制变量法(相关抽禅法)( control variates 控制变量法利用數学上积分算的性特性: ∫f(x)k=jU(x)-g()女+g(x 选择函数g(x)时考虎到:g(x)在蓬个积分区间部是容易确算出, 并且在上式右边第一项的還算中对(-g)积分的方差应当比第 项对j积分的方鑾小。 在应用这种方法时,在量要抽样法中所遇到的,当g(x)急于 粵时,被积函数-g)子无穷大的困难就不再亭在,因而计算出 的结果稳定性比較好。该方法也不卿耍从分布鲁度函数g(x),解 析求出分布函数G(x)。由此们可以看出选择g(x)所受到的限制 比量县抽样法要小些。 四、对偶变量法( antithetic variates 遢常在蒙特卡洛讣犷中采用互相独立的随机点来选行计算。 对偶变量法中却使用相关联的点来选行讣算。它利用相关点 间的关系可以是正头联的。也可以是负关联的这个兮点 两个函数值/和之和的方姜为 r+f2}={}+rV2}+2E-EX2-E2 如录们选择一些点,它们良/和是负关联的。这样就可以良上 式所示的方些藏小。当然这鲁要对具体的函数/和/有充分的了
三、 控制变量法（相关抽样法）(control variates) 控制变量法利用数学上积分运算的线性特性： [ ] ∫ ∫ ∫ f (x)dx = f (x) − g(x) dx + g(x)dx 选择函数 时要考虑到：g 在整个积分区间都是容易精确算出, 并且在上式右边第一项的运算中对( g(x) (x) f − g)积分的方差应当要比第 二项对 f 积分的方差小。 在应用这种方法时，在重要抽样法中所遇到的，当 趋于 零时，被积函数 趋于无穷大的困难就不再存在，因而计算出 的结果稳定性比较好。 该方法也不需要从分布密度函数 ， 解 析求出分布函数G 。 由此我们可以看出选择 所受到的限制 比重要抽样法要小些。 g(x) g(x) ( f − g) (x) g(x) 四、对偶变量法(antithetic variates) 通常在蒙特卡洛计算中采用互相独立的随机点来进行计算。 对偶变量法中却使用相关联的点来进行计算。它利用相关点 间的关系可以是正关联的，也可以是负关联的这个特点。 两个函数值 f1和 f 2之和的方差为 V{ } f1 + f 2 =V{f1}+V{f 2}+ 2E{( f1 − E{f1})( f 2 − E{f 2 })} . 如果我们选择一些点，它们使 和 是负关联的。这样就可以使上 式所示的方差减小。当然这需要对具体的函数 和 有充分的了 解。 1f 2 f 1f 2 f 7
例已知f(x)是一个单调增的函数,现求积分 I=f(x)dx 解首先,换逦常的方法在积分域[0,1区间上产生均勻分布的 随机点集{}。计算对应每个x的函数[f(x)+f(1-x)/2的值,开将 所有点上的函数迭加炙來,除以总的隨机点数,则得到(2.4.18) 式的积分值。即 ≈∑[(x)+f(-x, 这种散法与通常的象特卡洛计算中将f(x)的寶选加起來不相同。 由于f(x)的单调递增性,[(x)+f(-x)2的值应当比单个点的函数值 f(x)接近于常数。因而方也小些。 这实际上是采用了f(x)和f(-x)的积分期望值的平均值作为 结果。由于采用相同的隨机数列x},傯得f(x)和f(1-x)刘个函数 高度负关联,因而方差比f(x)和f(1-x)两者各自积分的方差之和要
例 已知 f (x)是一个单调递增的函数，现求积分 = ∫ . 1 0 I f (x)dx 解 首先，按通常的方法在积分域[0，1]区间上产生均匀分布的 随机点集{xi}。计算对应每个xi点的函数[ ( ) (1 )]/ 2 i i f x + f − x 的值，再将 所有点上的函数值迭加起来，除以总的随机点数，则得到(2.4.18) 式的积分值。即 ∑[ ] = ≈ + − N i i i f x f x N I 1 ( ) (1 ) / 2 1 . 这种做法与通常的蒙特卡洛计算中将 的值迭加起来不相同。 由于 的单调递增性，[ ( )i f x f (x) f (xi) + f (1− xi)]/ 2的值应当比单个点的函数值 f (xi)更接近于常数。因而方差也小些。 这实际上是采用了 f (x) 和 f (1− x) 的积分期望值的平均值作为 结果。由于采用相同的随机数列{xi} f (1 x ，使得 和 两个函数 高度负关联，因而方差比 和 f (x) f (1− x) f (x) − )两者各自积分的方差之和要 小。 8