2.5实用蒙特卡洛计算复合技术 在象特卡洛方法应用中澂小方些的基本技术:量要抽样法, 分层抽桦法,控侧变量法和对偶变量法。嶽而,单独伩用四种 减小方鎏的技巧仍有其局限性。 人仉发展了一些复合象特卡洛计箕技术,如应性Ω卡洛 方法和多熴毀兮卡洛抽祥方法亭。这些象卡涪技巧对于被积函 數在积分觉国內具有多个尖峰的情况,特别具有实用价值 适应性蒙特卡洛方法( adaptive Monte Carlo method) 适应性Ω守卡洛方法是一种在执行过程中邋过试探,了解被 积函数习性,然后有针对性地采用求特卡洛技万来减少方鑾的算 采用此方法的子程序有利帕格(G.P. Lepage)的 VEGAS[5]。它 是用于计算多量积分的子程序,广泛地庇用在高能物理领域。 VEGAS编程的基本思知是将量要抽样法和分层抽样法结合到進代 算法之中,该算法能够散自动调,将对被积函数的计算中到 被积函数值最大的区间。 以一錐定积分为例, VEGAS程序一开始处于试探阶段,即将积 分区间[]划分为正交子区间,并在每一个子区间中进行积分;然 后按脆音个子区间积分得到的绪果来调蕘子区间大小以备下 次迭代计调蓬子区间大小的原则是按服该子区间对总积分贡
2. 5 实用蒙特卡洛计算复合技术 在蒙特卡洛方法应用中减小方差的基本技术：重要抽样法， 分层抽样法，控制变量法和对偶变量法。然而，单独使用这四种 减小方差的技巧仍然有其局限性。 人们发展了一些复合蒙特卡洛计算技术，如适应性蒙特卡洛 方法和多道蒙特卡洛抽样方法等。这些蒙特卡洛技巧对于被积函 数在积分范围内具有多个尖峰的情况，特别具有实用价值。 一、 适应性蒙特卡洛方法(adaptive Monte Carlo method) 适应性蒙特卡洛方法是一种在执行过程中通过试探，了解被 积函数习性，然后有针对性地采用蒙特卡洛技巧来减少方差的算 法。 采用此方法的子程序有利帕格(G.P. Lepage)的 VEGAS[5]。它 是用于计算多重积分的子程序，广泛地应用在高能物理领域。 VEGAS 编程的基本思想是将重要抽样法和分层抽样法结合到迭代 算法之中，该算法能够做自动调整，将对被积函数的计算集中到 被积函数值最大的区间。 以一维定积分为例，VEGAS 程序一开始处于试探阶段，即将积 分区间 划分为正交子区间，并在每一个子区间中进行积分；然 后按照各个子区间积分得到的结果来调整子区间大小以备下一 次迭代计算，调整子区间大小的原则是按照该子区间对总积分贡 [0,1] 1
的大小来确定。贡大的子区间调恐更小一些。贡小的子 区间调蓬孙更大一些。 VEGAS程序用这种方式实际上就是采用了量要抽样法。它采用 阶梯函数來近似子区间的最佳几率寧度函数。该最能几率密度函 教定义 f( Po 「f(x)dk 对于高维(d錐)积分间题,由于存贮的要,我们必须采用分 高变量的分布度函教 p(l4,2,u)=p()p{u2)…p(u) 最后通过考干次选代,达到在要求的度下,各子区间(或子空 间)的积分估计值都相普,则我仉就找劁了优化的子区间(或子 空间)。在调蓬子区间(或子空间)过程中为了免子区间(或 子空间)剧烈的变化,子区间(子空间)大小的调通常有一个 衰减项。在程序的第二阶段,子区间(子空间)网格就固定下来 了,并过通常的蒙管卡洛方法得到在这些优化子区间(子空间) 中选代计算高精度的积分结果。每次选代积分部得到一个佔计值 E和一个方差P}: p(, ) h 其中N示在第二阶段第j次选代(=1,2,,m)积分的随机个。 在该阶段第j次积分值对总积分的觥权量应当为
献的大小来确定，贡献大的子区间调整得更小一些，贡献小的子 区间调整得更大一些。 VEGAS 程序用这种方式实际上就是采用了重要抽样法。它采用 阶梯函数来近似子区间的最佳几率密度函数。该最佳几率密度函 数定义 ( ) ( ) ( ) 0 1 0 f x p x f x dx = ∫ . 对于高维（d 维）积分问题，由于存贮的需要，我们必须采用分 离变量的分布密度函数 ( ) ( ) ( ) ( ) p u u ud p u p u p ud , ,..., = ⋅ ⋅...⋅ 1 2 1 2 . 最后通过若干次迭代，达到在要求的精度下，各子区间（或子空 间）的积分估计值都相等，则我们就找到了优化的子区间（或子 空间）。在调整子区间（或子空间）过程中为了避免子区间（或 子空间）剧烈的变化，子区间（子空间）大小的调整通常有一个 衰减项。在程序的第二阶段，子区间（子空间）网格就固定下来 了，并通过通常的蒙特卡洛方法得到在这些优化子区间（子空间） 中迭代计算高精度的积分结果。每次迭代积分都得到一个估计值 E{I j}和一个方差V{I j}： { } ( ) ( ), 1 1 ∑= = N j n n n j j p x f x N E I { } ( ) ( ) { j} N n n n j j j E I p x f x N I j 2 2 1 1 − = ∑= V . 其中N j表示在第二阶段第 j 次迭代( j = 1,2,...,m)积分的随机点个数。 在该阶段第 j 次积分值对总积分的贡献权重应当为 { j} j j V I N w = 2
将每次选代的积分值乘上与投点数N和方鑾相尝的权量,蝼后尿 加起来求平均。则总的积分佑计值为 E 此外, VEGAS程序巡回时出每个自由度( per degree of freedom) 的z2为 I PENJ-E(P , 这个结果可以作为检验各种估计值是否一政。我们期望x2/lf不 要比1大很多。 换照上面的思, VEGAS程序计算高维积分的步骤可以概括 如下 (1)将积分区域(域空间)划分为大量不相交的子区间(域子 空间)。原则上可以任意划分,但为了方便起见,住佳采用均勻 划分的办法。 (2)用原始蒙怜卡洛方法计每个子区间(域子空间)上的积 分值,开将个积分值选加起来作为个积分域上的佔计值 (3)调董子区间(或子空间)的边界,使得被积函数在个子区 间(域子空间)内的积分估计值大政相普。 (4)量复(1)-(3)的过程,利用原始蒙卡洛方法计算每次 选代的积分估计值,直到在要求达到的度下,备子区间(或子 空间)的积分估计值都相普。此时才将得到的子区间(或子空间)
将每次迭代的积分值乘上与投点数 和方差相关的权重，然后累 加起来求平均。则总的积分估计值为 N j { } { } { } { } { } = ∑ ∑ = ∑ ∑ = = = = m j j j m j j j j m j j m j j j V I N V I N E I E I w E I w 1 1 1 1 . 此外，VEGAS 程序返回时给出每个自由度(per degree of freedom) 的 为 2 χ ( { } { }) { } ∑= − − = m j j j V I E I E I m dof 1 2 2 1 1 χ / . 这个结果可以作为检验各种估计值是否一致。我们期望 不 要比 1 大很多。 / dof 2 χ 按照上面的思想，VEGAS 程序计算高维积分的步骤可以概括 如下： （1） 将积分区域（或空间）划分为大量不相交的子区间（或子 空间）。原则上可以任意划分，但为了方便起见，往往采用均匀 划分的办法。 （2） 用原始蒙特卡洛方法估计每个子区间（或子空间）上的积 分值，再将各个积分值迭加起来作为整个积分域上的估计值。 （3） 调整子区间（或子空间）的边界，使得被积函数在个子区 间（或子空间）内的积分估计值大致相等。 （4） 重复（1）-（3）的过程，利用原始蒙特卡洛方法计算每次 迭代的积分估计值，直到在要求达到的精度下，各子区间（或子 空间）的积分估计值都相等。此时才将得到的子区间（或子空间） 3
固定下来。以上为程序计算的第一阶段。在这一阶段,批京个数 可以少一些,并不记这个阶段的积分结果(因为一般方差都很 大)。 (5)最后,采用蒙特卡洛方法,换胤公式(2.5.3)计算各子区间 (域子空间)积分值和方塾,瘘后利用公式(2.5.5)将每次选代 计算的积分值加权累加平均得到该积分在总区间的积分佔计值, 用公式(2.5.6)计算每个自由度的z2。这就得到该积分的数億计 犷結鼎。这是程序计算的第二阶段。 二、多觉蒙管卡洛抽禅方法( multi- channel Monte Carlo method) 我们必熟以前面一维定积分为例,如景被积函数∫(x)在被积 区间有尖峰,则采用原始的象管卡洛积分的结果误鑾背定是很大 的。当然,我们可以采用变量变换,换照量耍抽禅方法,将被积 函數的峰篁兮性吸收到随机点的分布函数中。这种力法多少会 积分結果度得到欧菩。但是,可能会有下面的情况:被积函数 ∫(x)在被积区间的不同区域有多个不同的尖峰在这神情况下, 往往不可能找到一个变量代换,它既吸收了被积函教所有的峰值 性。又校容忌选行换物定分布的随机教抽样。多效象镎卡洛 方法就是针对被积函数∫(x)具有多个尖峰情况下的计算方法 它的基本思剋是瓴于泉兮卡洛方法的加原则加上量聂抽 样法。该方法疝用的前提是每个峰绪袍变换成的近似抽样分布鲁
固定下来。以上为程序计算的第一阶段。在这一阶段，投点个数 可以少一些，并不记这个阶段的积分结果（因为一般方差都很 大）。 （5） 最后，采用蒙特卡洛方法，按照公式(2.5.3)计算各子区间 （或子空间）积分值和方差，然后利用公式(2.5.5)将每次迭代 计算的积分值加权累加平均得到该积分在总区间的积分估计值， 用公式(2.5.6)计算每个自由度的 。这就得到该积分的数值计 算结果。这是程序计算的第二阶段。 2 χ 二、多道蒙特卡洛抽样方法(multi-channel Monte Carlo method) 我们仍然以前面一维定积分为例，如果被积函数 f ( x)在被积 区间有尖峰，则采用原始的蒙特卡洛积分的结果误差肯定是很大 的。当然，我们可以采用变量变换，按照重要抽样方法，将被积 函数的峰值特性吸收到随机点的分布函数中。这种方法多少会使 积分结果精度得到改善。但是，可能会有下面的情况：被积函数 f ( x)在被积区间的不同区域有多个不同的尖峰。在这种情况下， 往往不可能找到一个变量代换，它既吸收了被积函数所有的峰值 特性，又比较容易进行按特定分布的随机数抽样。多道蒙特卡洛 方法就是针对被积函数 f ( x)具有多个尖峰情况下的计算方法。 它的基本思想是源于蒙特卡洛方法的迭加原则加上重要抽 样法。该方法应用的前提是每个峰结构变换成的近似抽样分布密 4
度函教形式巳经知熴。每个峰的变换称为一个觉。如果对应第i个 馗抽样,它所选择的分布鲁度函数为h(x)。据分布鲁度函数的 正定和归一化,我们有:「(x)=1,(=1,m),美中为道数。 令a为非负的实教,并且满足 由于们定义 h(x)=∑ah(x) 这就衰明:对分布鲁度函数h(x)抽样时,可以分别对h(x)抽样 但是选择对第i个熴的分布密度函数抽祥的几率为a。明显地被积 函数∫(x)与分布鲁度函数h(x)具有同样的多峰值结袍。利用关系 式(2.5.8),我们将积分如做如下形式推导: 1(k=-0 M)/(x)=∑[(x) 5)()a f(x h(r) dH, (x) 此时被积函数∫(x)/h(x)中已经没有原来所有的峰值特性了,这些 峰值管性已经被分布鲁度函数h(x)抵消。这就是我们多觉蒙管卡 洛方法计箕积分的基本公式。从公式中可以看出:我们可以還过 对各个道,按分布密度函数h(x)(对分布函数为H(x)产生 随机数x。创如具体抽祥时可以用反函数法: 实践中,我们换照离散型变量抽样法以a为取分布密度函数h(x)
度函数形式已经知道。每个峰的变换称为一个道。如果对应第 个 道抽样，它所选择的分布密度函数为 i hi ( x)。根据分布密度函数的 正定和归一化，我们有： ( ) 1 0 1 i dxh x = ∫ ， (i =1,...,m)，其中i为道数。 令αi为非负的实数，并且满足 ， 1 1 m i i α = ∑ = 由于我们定义 ( ) ( ) 1 m i i i h x α h x = = ∑ 这就表明：对分布密度函数h x( ) 抽样时，可以分别对 抽样。 但是选择对第 个道的分布密度函数抽样的几率为 hi ( x) i αi。明显地被积 函数 f ( x)与分布密度函数h x( ) 具有同样的多峰值结构。利用关系 式(2.5.8)，我们将积分I 如做如下形式推导： ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 0 0 0 1 m i i i f x f x I f x dx h x dx h x dx h x h x α = = = = ∫ ∫ ∑∫ ( ) ( ) ( ) 1 0 1 m i i i f x dH x h x α = = ∑ ∫ . 此时被积函数 f ( ) x h( x)中已经没有原来所有的峰值特性了，这些 峰值特性已经被分布密度函数h x( ) 抵消。这就是我们多道蒙特卡 洛方法计算积分的基本公式。从公式中可以看出：我们可以通过 对各个道i，按分布密度函数hi ( x)（对应分布函数为 ）产生 随机数 Hi ( x) i n x 。例如具体做抽样时可以用反函数法： ( ) 1 i i n i n x H ξ − = 实践中，我们按照离散型变量抽样法，以αi为取分布密度函数h x i ( ) 5
抽样的概率。若固定总投京数为N,计算第i时所投京次数大约 是N≈a1N。采用这样的方法,可以孙到积分的象卡洛佔计为 该象卡洛积分的误塾期望值为 (a)-/2 N 其中W(a)定义为 W(a)=∑a H1(x) 该方法计算的关之一是确定参数a。实中佳佳漉过调參数 a,使W(a)达到最小值来优化参数a的选择。理论上积分计算值 应与参数a值的选择无於,因此在积分过程中敦们可以改变a的 夜,而不食影响积分寶的佔计当这是食影响计算结果的误差 议首先采用一组物始参数{},在换照上面介绍的步骤选 行几百次物卡治计算后,再计 H(a)=/(x) (1)/(x)女 最后蛰下面公式量新决定参数a,。 ∑o(r(a) 根据经验,议上面公式中的参数B值取在[/41/2]之间
抽样的概率。若固定总投点数为 ，计算第i道时所投点次数大约 是 N Ni ≈αiN 。采用这样的方法，可以得到积分的蒙特卡洛估计值为 (α ) 2 H (α ( ) i x α { } ( ) 1 1 ( ) 1 i i i i m N n i n n f x E I N = = h x = ∑∑ . 该蒙特卡洛积分的误差期望值为 ( ) 2 W I N α − 其中W 定义为 ( ) ( ) ( ) ( ) 1 0 1 m i i i f x h x α α= = W d ∑ ∫ x . 该方法计算的关键之一是确定参数αi。实践中往往通过调整参数 αi，使W )达到最小值来优化参数αi的选择。理论上积分计算值I 应与参数αi值的选择无关，因此在积分过程中我们可以改变αi的 值，而不会影响积分值的估计。当然这是会影响计算结果的误差 的。 建议首先采用一组初始参数{αi ′}，在按照上面介绍的步骤进 行几百次蒙特卡洛计算后，再计算 ( ) ( ) ( ) 2 1 0 i i f x h x α ′ = W h ∫ x d 最后按照下面公式重新决定参数 。 ( ) ( ) ( ) ( ) 1 i i i m j j j W W β β α α α α α = ′ ′ = ∑ ′ ′ 根据经验，建议上面公式中的参数β 值取在[1/ 4,1/ 2]之间。 6
