蒙特卡罗方法初步—椭圆内均匀分布的抽样技巧

2023年12月20日发(作者：将喜欢的人写进数学试卷)

蒙特卡罗方法初步—椭圆内均匀分布的抽样技巧

q一l

1996年12爆轰波与冲击波第4期

蒙特卡罗方法初步——

椭圆内均匀分布的抽样技巧

胡熙静刘桂贤

-

—,———～

A耍以简明通俗的语言.避开统计学中抽象的数学公式和符号,介绍了作为计算数学一

十独立分支的蘩特卡罗方法的基本思想和概念,奠特卡罗方法在粒子{宣运问题中的应用以爰

蒙特卡罗(MonteCarlo)是意大利的一个城市名,以赌博闻名于世界,可以想象,蒙特

卡罗与赌博有某种相似之处.没有接触过蒙特卡罗方法的人,对它有一种神秘感,其实它

的基本思想并不新颖.早在十七世纪,人们就知道用频数来决定概率的方法.随着现代科

学技术的发展和电子计算机的发明,应用愈来愈广,方法本身也得到不断的发展和完善.

成为计算数学的一个独立分支.

蒙特卡罗方法又称随机抽样技巧或统计试验方法.本文的目的是一般介绍蒙特卡罗

方法的重要概念,基本思想,在粒子输运问题中的应用,常见的几种抽佯方法.最后,我们

介绍在椭圆内均匀分布的一种效率较高的抽样技巧.

2蒙特卡罗方法的基本思想

i,,

』蒙特卡罗方法的基本思想是,把所要求解的问题转换成某种事件出现的概率,然后通

过对事件的模拟”试验”得到该事件出现的频率,并以此近似代替该事件出现的概率,从而

获得问题的解.用统计学的语言来表述就是,把所求问题转换成某个随机变量的期望值,

通过模拟”试验的方法,获得这个随机变量的平均值.

有一个大家都熟悉的饲子,一批产品的质量好坏用合格率来表示.合格率P兰n/Ⅳ,

其中Ⅳ为随意确定的抽样数,n为合格数,近似号表示Ⅳ不可能等于全部产品数.这里把

求解的问题(产品质量),转换为概率(合格率),通过抽样捡验的方法,求

得-/Ⅳ作为合格

率的近似值.再举一个例子:实验确定圆周率.设想地面上有一组平行直线,相距为.

手里拿一根长为2f的针(<a),作Ⅳ次任意的投针试验,与平行线相交者记为1,否则记

为0.这样得到针与平行线的相交概率P:/Ⅳ.从统计学的观点可以证明Ⅲ此概率的正

确结果为21/~a甩P近似代替正确概率2Z/~a,就得到的值为

Ⅱ兰(2~/a)(Ⅳ)

这就是古典概率论中着名的蒲丰氏问题.

可以把蒙特卡罗方法归结为三个主要步骤:构造所求问题的概率;实现从已知分布中

抽样;建立估计量.在上述的蒲丰氏问题中就是把求的问题转换为投针相交率,进行投

爆轰波与冲击波1995年12月

针试验(抽样),求得相交概率的估计量.现在当然不会有人真正做这种费事的投针试验.

在计算机上可以用蒙特卡罗方法进行模拟,做到投针次数雎非常大,就获得非常精确的

Ⅱ值..

上述的例子属于本来不是随机性质的确定性闻蘧.这类伺题酌例子还

有计算定积分,

解线性方程组,偏微分方程的边值问题等.这类问题用蒙特卡罗方法时须要事先构造一个

人为的概率过程.即将不具有随机性质的问题转换为随机性质的问题,使它的某些参量正

好是所要求的同题的解.由于这类问题往往有其它方便的计算方法求解.蒙特卡罗方法使

用不多.

实际问题中有一类本身具有随机性质,如粒子输运伺题.蒙特卡罗方法的主要任务是

正确地模拟这个随机过程.随着计算机技术的飞速发展,在粒子输运领域内,蒙特卡罗方

法得到广泛的应用.’

3粒子输运问题

在反应堆设计,核物理实验,辐射(中子或光子)屏蔽,皂子与物质相互作用问题中,都_’●

涉及大量的中子,光子,电子的输运问题.这是蒙特卡罗方法应用的主要领域.粒子在介质

中的作用过程本身具有随机性质.一个由源发出的粒子,在它的运动方向上,在啷一点碰

撞是偶然的,但有一定的概率分布I与原子(或康子核)发生碰撞又有各种概率不同的碰撞

类型涠撞后的能量和运动方向,遵从一定的

概率分布.粒子可能被介质吸收或从系统中

逃脱,这时粒子运动过程结束,否则继续下一

次类似的运动.一个粒子在介质中的运动情

况,可由它所经历的碰撞反映出来.所以只要

粒子与原子或核的碰撞规律在物理上是清楚

的,那么这种过程完全能够用蒙特卡罗方法

正确地模拟,然后对所求的物理量进行统计

平均,就求得问题的解.例如对于屏蔽问题,

把模拟粒子数记为Ⅳ,穿透屏蔽层的粒子效

记为n,则穿透率P兰/Ⅳ.蒙特卡罗方法计

算粒子输运问题的大致步骤如图1所示.

图t中:源分布抽样包括源粒子能量,位

置,发射方向的抽样;碰撞过程需根据具体同

题的实际过程来进行;历史终止包括粒子被

吸收,从系统逃脱,能量或权重下降到某一规

定的最小值等等;统计处理表示对该计算所

要获得的结果量进行统计平均;给定历史数

表示事先确定的源粒子抽样总教.

图1蒙特卡罗方法计算程序结构框图

碰撞过程模拟中的主要任务是从具有一定概率分布的碰撞过程中随机抽样,从而确

第{期胡熙静等,蒙特卡罗方法初步——椭腰内均匀分布的抽样技巧

定该粒子确_切的碰撞结果.1逮机抽样就需要用到随机数.用随机数对已知概率分布的抽

样方法将在后面介绍,下面先介绍随机数.

4伪随机数

随机数分为物理随机数和伪随机数两种.物理随机数的来源有两种,一种是根据放射

性物质的放射性,另一种是利用电子计算机的固有噪音.这种方法的缺点是价格昂贵并且

不能重复,使程序无法进行重复运算.在计算机上最常用的,是用数学方法0在区间[O,1]

上产生的均匀分布的数列,称为伪随机数.不同的数学方法产生的0艟机数列品质是不同

的,在使用中要注意挑选品质好的随机数列.品质的好坏包括随机数的均匀性,数列周期

长短.用数学方法产生鹊伪随机数列,在经过一定的个数后井始重复前面的数列,这种性

质叫随机数的周期.我们希望周期愈长愈好.如何检验品质的好坏,我们这里不作介绍.

现在已有不少比较成熟的产生伪随机数的数学递推关系,可供挑选使甩由于经常使用的

是伪随机数,在下面的叙述中我们简称为随机数.:

5抽样方法

从已知分布概率F”)中,用已知的[O,1]区间内均匀分布的随机数,确定直,并保

证茬次运干乍后.,获得的十z值曲分布概率与,()近似相同,这个过程Ⅱ崎-做抽谇.根

据不同的分布概率F(z),有不同的抽佯方法.我们只介绍两种最简单的抽伴方法以说明

抽样的基本概念,详细的知识可参看文献[1].一一:～一lll|_

5.1离散型分布概率的抽样方法

我以光子与物质的相互作用为倒光子与原子碰撞时发生三种反应光电吸收,康

普顿散射,电子对生成,它们的截面分别为,以,,总截面为以++.酋先构造一

个归一化分布概率p

.

-

?

?

.

.

P.=(1/唧)(1,2,3)-.

并规定P.0o然后拙取[0,1]内均匀分布的随机数乱,如果满足不等式

PJ-l≤{l<P.

则J值为抽样确定的反应道.倒如;2表示抽样确定该粒子在此点发生康普顿散射.

5.2连续型分布概率的抽样方法’

在上述倒子中,确定了碰撞为康普顿散射反应,然后需要确定散射方『句.’我们知道康

普顿散射微分截面,0)与散射角口有关,这是随连续变化的一个函数.同样,首先构造

一

个归一化的分布概率P(日)

P()=I,()sinOdO/I.,(日)sjnd日.,.

抽取随机敬,如果P()一{,则即为抽样值,即确定了散射方向需要说明.这个简单

的抽样方法在实际使用中不方便,必须加以改进.以上两种都是直接抽样的基本思想.

5.$半径为R的圆内均匀分布的粒子坐标点抽样

设一束电子的束半径为,粒子密度在R内均匀分布,抽样确定某一粒子在圃内的

爆轰波与冲击波1995年l2月

坐标,值.我们知道圆方程是+,..设玉=/R.;/R,则圆方程为

磊+器:l

在单位圊外作外切正方形(见图2).抽样的基本思想是首先在单位正方形内作均匀分布

抽样,这可以从随机数,获得,然后保留满足+;≤l的情况,不满足的舍去,抽样

框图如图3所示.

(0

‘‘

匿2匮内均匀分布的抽样冒3抽样撮圈

这种根据条件决定取舍的抽样叫做挑选抽样.挑选抽样存在抽样效率,在这个例子中

抽拌效率是圆和正方形面积之比=n/4.

B椭圆内均匀分布蜉抽样方法

我们在进行lOMeVLIA.闪光机的x光转换靶上的蒙特卡罗方法计算x光源的能

谱,剂量角分布时需要对椭圆内均匀分布的粒子抽样.我们知道加速器出射的电子柬,在

四维相空间,,,,,)中是个捕球,其中,为坐标,,,,为运动方向与轴的夹角.

如果假定在该球内的分布为K-V分布.则可证明在二维相空间,,)投影的椭圆内也

是均匀分布的.这里介绍一种满足这种要求的效率较高的新的抽样方法井证明其正确性.

为了作比较,我们先筒述一种比较直观的抽样方法.设,,平面内有一椭圆,粒子在

椭圆内均匀分布.捕四方程为

+7-”=l(I)

其中a,b分别为椭圆的长,短轴.首先以长轴2a为边作一个外接正方形(见图4),把椭圆

包围在正方形内抽取随机数,如,在正方形内,抽样粒子的坐标是

‘:(2)=(2岛一1)d

-,岛是在[0,1]区间内均匀分布的随机数,2{一l是在[一l,门区间内均匀分布的随机

数.然后把,值代人下面不等式

().+(})≤l

决定取舍.显然.随机抽得的,在正方形内均匀分布,在椭圆内也一定是均匀分布的.

第4期胡熙静等:荣特卡罗方法初步——椭四内均匀分布的抽样技巧

这种抽样方法比较直观,其抽样效率为正方

形与椭圆的面积比

=丽nab=bqt,>6(3)西’.

下面讨论一种新的抽样方法.

把椭圆方程(1)作变换

(+:l(4)(三)+;一=l(4)

‘一

詈’圈4椭四内的抽样方法图形

则(4)式变为圆方程

(旦).+(旦)=1

新的抽样方法是,在,)相空间内,用第5节中介绍的圆内均匀抽样方法,抽得,的坐

标值

f”‘磕一h(6)

l(2一1)d

取满足不等式

+≤l(7)

dd

的粒子为样品粒子坐标,由变换(5)得到在椭圆内的样品粒子坐标

==

(2{l一1)d

:

bT:(2一1)6(8)

==Lz一l

(8)式的意义是在圆内均匀抽样得的样品粒子坐标,压缩到椭圆内的坐标,.并且这

种压缩保证了在椭圆内仍然是均匀分布的

(见图5).

需要指出的是,在变换式(4)中并不要求

必需>6.在上述的情况中,如果a>b,则变

换式(4)是把=.的圆包围了椭圆.反之如

果a<b,则R=的圆是在椭圆的内部,此时

变换式(8)的意义是在圆内均匀抽样得的样

品粒子坐标,扩散到椭圆内的坐标,.

我们需要证明的是这种扩散或”压缩仍保证

在椭圆内均匀分布.

命题;在半径为的圆内,均匀抽样的

,经(5)式的线性变换,在椭圆内的样品坐

标,也是均匀分布的.

图5椭圆和圆之间的关系圈

爆轰波与;中击渡1995年12月

首先,我们考查方向.由(5)式知与是线性关系即/一a.这就是说对任一

值,在方向的粒子密度都以比倒a邝均匀压缩(或扩散).再考查密度P随方阿的

变化.在轴任取一点和小间隔d,作垂直于轴的窄条面(如图5阴影部份).此窄

条面交椭圆于z.,o,交圆于.,..当d很小时,椭国内的阴影面积—d?o,在圆内

的面积sA,.设圆内密度为.(为一常数),椭圆内密度(压缩后)为p.在压缩

过程

中粒子数守恒.则有一,即

=Po

爱:㈤p

由圆方程

0一~/Ⅱ一;

由椭圆方程

.=

÷一z5

注意到o,褥

一

(Io)

Oo

从而

Ⅱ

P—Poi

可见,压缩后的密度p与无关,且为一常数.所以压缩后在方向也是均匀的.

此方法的抽样效率是圆的直径为边的正井彩面积与圆面积之比,即—/与长短

辅a’b无关,可见<2.

7结束语

蒙特卡罗方法作为一门计算数学的独立分支,内容是非常丰富的.随着现代电子计算

机的高速发展,应用也越来越广泛.各种实际问题中的模拟技巧及特殊抽样方法也在_不断

发展.文献113比较全面,系统地阐述了蒙特卡罗方法的基本思想方法,误差和改进等问

题,是一本比较好的参考书.我们在10MeVL1A闪光机的x光源转换靶的x光能谱,角

分布,剂量分布的计算中使用了蒙特卡罗方法,获得了较好的结果.

参考文献

1裴鹿成,张李泽蒙特卡罗方法及其在粒子箱运问题中的应用?北京:科学出版社,1980

2Forsythe0EC~neTationandtvstingoFTandomdigitsNBSApp]MalhSerie~,195】,12}34