2. 5 实用蒙特卡洛计算复合技术 在蒙特卡洛方法应用中减小方差的基本技术：重要抽样法， 分层抽样法，控制变量法和对偶变量法。然而，单独使用这四种 减小方差的技巧仍然有其局限性。 人们发展了一些复合蒙特卡洛计算技术，如适应性蒙特卡洛 方法和多道蒙特卡洛抽样方法等。这些蒙特卡洛技巧对于被积函 数在积分范围内具有多个尖峰的情况，特别具有实用价值。 一、 适应性蒙特卡洛方法(adaptive Monte Carlo method) 适应性蒙特卡洛方法是一种在执行过程中通过试探，了解被 积函数习性，然后有针对性地采用蒙特卡洛技巧来减少方差的算 法。 采用此方法的子程序有利帕格(G.P. Lepage)的VEGAS[5]。它 是用于计算多重积分的子程序，广泛地应用在高能物理领域。 VEGAS编程的基本思想是将重要抽样法和分层抽样法结合到迭代 算法之中，该算法能够做自动调整，将对被积函数的计算集中到 被积函数值最大的区间。 以一维定积分为例，VEGAS程序一开始处于试探阶段，即将积 分区间划分为正交子区间，并在每一个子区间中进行积分；然 后按照各个子区间积分得到的结果来调整子区间大小以备下一 次迭代计算，调整子区间大小的原则是按照该子区间对总积分贡 [1,0 ] 1 献的大小来确定，贡献大的子区间调整得更小一些，贡献小的子 区间调整得更大一些。 VEGAS程序用这种方式实际上就是采用了重要抽样法。它采用 阶梯函数来近似子区间的最佳几率密度函数。该最佳几率密度函 数定义 () ( ) () 0 1 0 fx px f xdx = ∫ . 对于高维（d维）积分问题，由于存贮的需要，我们必须采用分 离变量的分布密度函数 ()( ) ( ) ( ) dd upupupuuup ???= ...,...,, 2121 . 最后通过若干次迭代，达到在要求的精度下，各子区间（或子空 间）的积分估计值都相等，则我们就找到了优化的子区间（或子 空间）。在调整子区间（或子空间）过程中为了避免子区间（或 子空间）剧烈的变化，子区间（子空间）大小的调整通常有一个 衰减项。在程序的第二阶段，子区间（子空间）网格就固定下来 了，并通过通常的蒙特卡洛方法得到在这些优化子区间（子空间） 中迭代计算高精度的积分结果。每次迭代积分都得到一个估计值 { } j IE和一个方差{ } j IV： {} () () , 1 1 ∑ = = j N n n n j j xp xf N IE {} () () {} j N n n n j jj IE xp xf N I j 2 2 1 1 ? ? ? ? ? ? ? ? ? = ∑ = V . 其中表示在第二阶段第 j N j次迭代( )mj ,...,2,1=积分的随机点个数。 在该阶段第j次积分值对总积分的贡献权重应当为 {} j j j IV N w = 2 将每次迭代的积分值乘上与投点数和方差相关的权重，然后累 加起来求平均。则总的积分估计值为 j N {} {} { } {} {} ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? == ∑∑∑∑ ==== m j j j m j j jj m j j m j jj IV N IV IEN wIEwIE 1111 . 此外，VEGAS程序返回时给出每个自由度(per degree of freedom) 的为 2 χ { } { }( ) {} ∑ = ? ? = m j j j IV IEIE 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)计算每个自由度的。这就得到该积分的数值计 算结果。这是程序计算的第二阶段。 2 χ 二、多道蒙特卡洛抽样方法(multi-channel Monte Carlo method) 我们仍然以前面一维定积分为例，如果被积函数()f x在被积 区间有尖峰，则采用原始的蒙特卡洛积分的结果误差肯定是很大 的。当然，我们可以采用变量变换，按照重要抽样方法，将被积 函数的峰值特性吸收到随机点的分布函数中。这种方法多少会使 积分结果精度得到改善。但是，可能会有下面的情况：被积函数 ()f x在被积区间的不同区域有多个不同的尖峰。在这种情况下， 往往不可能找到一个变量代换，它既吸收了被积函数所有的峰值 特性，又比较容易进行按特定分布的随机数抽样。多道蒙特卡洛 方法就是针对被积函数( )f x具有多个尖峰情况下的计算方法。 它的基本思想是源于蒙特卡洛方法的迭加原则加上重要抽 样法。该方法应用的前提是每个峰结构变换成的近似抽样分布密 4 度函数形式已经知道。每个峰的变换称为一个道。如果对应第个 道抽样，它所选择的分布密度函数为 i ( ) i hx。根据分布密度函数的 正定和归一化，我们有：() 1 0 1 i dxh x = ∫ ， ( )1,...,im=，其中i为道数。 令 i α为非负的实数，并且满足 ， 1 1 m i i α = = ∑ 由于我们定义 () () 1 m ii i hx h xα = = ∑ 这就表明：对分布密度函数( )hx抽样时，可以分别对抽样。 但是选择对第个道的分布密度函数抽样的几率为 () i hx i i α。明显地被积 函数()f x与分布密度函数( )hx具有同样的多峰值结构。利用关系 式(2.5.8)，我们将积分I如做如下形式推导： () ( ) () () ( ) () () 11 1 00 0 1 m ii i fx fx I f x dx h x dx h x dx hx hx α = ?? ?? == ?? ?? ?? ?? ∑ ∫∫ ∫ () () () 1 0 1 m ii i fx dH x hx α = ?? = ?? ?? ∑ ∫ . 此时被积函数 () ()f xhx中已经没有原来所有的峰值特性了，这些 峰值特性已经被分布密度函数( )hx抵消。这就是我们多道蒙特卡 洛方法计算积分的基本公式。从公式中可以看出：我们可以通过 对各个道i，按分布密度函数( ) i hx（对应分布函数为）产生 随机数 () i Hx i n x。例如具体做抽样时可以用反函数法： () 1 ii nin xHξ ? = 实践中，我们按照离散型变量抽样法，以 i α为取分布密度函数 () i hx 5 抽样的概率。若固定总投点数为，计算第i道时所投点次数大约 是 N ii NNα≈。采用这样的方法，可以得到积分的蒙特卡洛估计值为 () α 2 H (α () i x α {} ( ) () 11 1 i i i i N m n in n fx EI N hx== ?? ?? = ?? ∑∑ . 该蒙特卡洛积分的误差期望值为 () 2 WI N α ? 其中W定义为 () () () () 1 0 1 m ii i fx hx αα = ?? = ?? ?? ∑ ∫ Wd. x 该方法计算的关键之一是确定参数 i α。实践中往往通过调整参数 i α，使W )达到最小值来优化参数 i α的选择。理论上积分计算值I 应与参数 i α值的选择无关，因此在积分过程中我们可以改变 i α的 值，而不会影响积分值的估计。当然这是会影响计算结果的误差 的。 建议首先采用一组初始参数{ } i α′，在按照上面介绍的步骤进 行几百次蒙特卡洛计算后，再计算 () () () 2 1 0 ii fx hx α ?? ′ = ?? ?? ∫ Wh d 最后按照下面公式重新决定参数。 ()() () () 1 i i i m j j j W W β β αα α αα = ′′ = ′′ ∑ 根据经验，建议上面公式中的参数β值取在 [ ] 1/4,1/2之间。 6