蒙特卡 法及具回用 1993~1997 裴鹿成主编 王仲奇 ogO
目录 理论研究部分 蒙特卡罗方法与随机性问题…………………………………裴鹿成(3) 彼得堡悖论与小概率大贡献问题…………………………………………裴鹿成(10) 二维随机几何模型的蒙特卡罗研究……………………… 王仲奇(18) 智能”型深穿透辐射输运蒙特卡罗模拟 杜凤英等(21) MC计算柱通量的指向概率方法… 王瑞宏等(24) 非归一分布随机抽样方法研究 程锦荣等(28) 几种新的伪随机数发生器………………………………裴鹿成(34) 离数型非归一分布的舍选抽样方法……锦荣等(38) Metropolis抽样中随机游动步长因子的确定…… ……程锦荣等(40) 应用及软件部分 非均匀系统的临界计算… ………沈雷生等(45) 光子一电子一模一核耦合输运问题的蒙特卡罗模拟………………许淑艳等(50) 外部噪声法实现超混沌同步…… 方锦清等(5 就地HP℃e谱仪探测器角响应校正因子的MC计算 肖雪夫等(60) 闪光照相1:1静态样品的数值模拟………………………………邹志高等(71) GL谱仪氚响应函数MC计算………………………………………吴建华等(78) 重水球慢化3C刻度装置的谱特征……………………………宁静等(83) 二维随机摆放过程的计算机模拟… …………王仲奇(86) 中子-Y耦合输运 Monte carlo程序MCNP在PM下的并行化… 邓力(88) 略论MCNP程序中的蒙特卡罗技巧及其不足……………………………裴鹿成(93) 核燃料后处理临界安全程序MCFR10介绍…………… ………课题组(102) KENO程序在临界分析中的应用……………… …薛小刚等(1065) OCTOPUS燃耗计算程序系统简介… 张宝成(13)
MORSE程序中的 Monte carlo技巧扩充… ………杨锦安等(117) FAMS-MC程序中测量修正中的应用…………… ·沈冠仁等(122) CCMC程序及其中反应堆控制棒均匀化参数计算中的应用 姚栋等(130) Monte carlo修正程序中使用数据的处理方法……………………毛孝勇等(13) CHMCK-Ⅲ:解任意几何临界问题的蒙特卡罗程序… 裴鹿成(137)
22222》282》2》28882别2》888》》》》》动 理论研究部分
蒙特卡罗方法与随机性问题 裴鹿成 (中国原子能科学研究院) 摘要随着科学技术的迅速发展,越来越复杂的随机性问题被提了出来。一般来讲,对 于现代科学技术中所提出来的随机性问題除极少数情况外,要想给出它的严格解是根本不 可能的,用确定性方法给出其近似解也常常是非常困难的,有时甚至也是不可能的。 蒙特卡罗方法以对随机性问题进行仿真为其基本特征1-5,这就决定了蒙特卡罗方法 对于解决随机性问题具有很强的能力。本文给出了4个用蒙特卡罗方法解决随机性问题的 实例,目的是想通过这些有趣的实例说明,对于许多用确定性方法所难以解决的随机性问 题,用蒙特卡罗方法可以比较方便地解决。 1Caus问题 us于1812年写信给 Laplace,提出了如下一个著名题:在(0,1)中任取一数,将它表 示成简单连分数,试问其第m个完全商的小数部分小于x(0<x<1)的概率为多少? 上述Gas问题可以改述成:在(0,1)中,任取一数a1,按如下办法确定其后的a2,a3 当a≠0时 1) 当a;=0时 其中{*}表示取数的小数部分,试问an小于x的概率为多少? Gas问题提出后,直到1928年,经过长达16年之久,苏联数学家KyM才给出了 一个渐近结果0。用Pm(x)表示an小于x的概率KyMm的结果是 P(x)=In(i+s (2) 公式(2)中a为一正常数。到1948年,他又特其中的渐近误差阶改进成O(am),0<a<1,并 且指出,这个渐近误差阶不可能再改进了。 为了解决Gas问题,很明显,用蒙特卡罗方法解不仅不会遇到任何困难,而且还非常 简单,其主要步骤如下: (1)为计算作准备 令n=0,P=0 (2)确定初值 令n=n+1;i=1,a1=(为在(0,1)上均匀分布的随机数) (3)确定下一个值 a;+1由公式(1)确定;令i=i+1 (4)是否巳经确定了
当i<m时,an尚未确定,转至步骤(3);否则,an已经确定。 (5)记录贡献 P+no am< r 其中n(*)表示条件函数当条件*成立时为1否则为零 (6)抽样是否结束 当n<N时(N为样本总数),抽样尚未结束转至步骤(2);否则,抽样结束,进入下一步 骤 (7)给出计算结果与误差 计算结果为 误差则为240 P(IPm. M(*)-Pa(a)1)<1. 96/ -m, M(x)(1-p(s2-095 (5) 其中P(*)表示事件*发生的概率。 进一步比较上述KysM方法和蒙特卡罗方法的优缺点。由于前者所给出的实际上 只是一个渐近公式,若根据式(2)用ln(1+x)/hn2作为所求Pn(x)的近似,当m较小时,相 差可能较大;当m较大时,相差虽然可能较小,可是,在KyMm的结果中只给出了渐近误 差的阶,无法确定误差,因此,要想用ln(1+x)/ln2确定Pn(x)对于任意m的值,除非m足 够大,是靠不住的。蒙特卡罗方法的情况则完全不是这样,它可以给Pn(x)对于任意m的 近似值及其误差估计,因此如果目的仅限于给出关于Rn(x)的数值结果,用蒙特卡罗方法 计算出来的Pnx(x),要比用渐近公式n(1+x)1n2计算出来的更具有实际意义 2随机徘徊问题 处在S维空间格点上的质点,每步向2S个相邻格点中的任一个移动的机会是均等的。 所谓随机徘徊问题是问,质点返回初给位置的概率P(S)(简称回返概率)为多少? 早在191年, Polya曾用分析的方法确定,P(1)=1,P(2)=1,P(3)≠1。可是,大概 是由于计算量太大的原因,直到1940年,事隔19年之久,才由另外两位学者给出了P(3)的 近似值:P(3)≈0.358。下面我们来考虑一般的S维空间的随机徘徊问题。 用Un表示质点随机徘徊于第m步返回初始位置的概率。根据Un的定义,不难看出 当m为奇数时,Um=0;当m为偶数时 40.(4!)2 另一方面,若简单地用P表示回返概率P(S),则根据回返概率和Um的定义,一定有如下等 式成立
P mp°(1-P) 1-P 由此方程立即得到关于回返概率P的解析表达式如下: /(1 (8) 计算回返概率的确定性方法就是,用下式 (9) 作为P的近似估计,其中M是一个足够大的数 按照随机徘徊问题的基本假设,与上述计算回返概率的确定性方法(9)相对应的,蒙特 卡罗方法计算回返概率的详细步骤如下 (1)为计算作准备 令n=1,R=0 (2)开始第n个质点的随机徘徊 令m=1,l1=0,i=1,2,…,5 (3)确定质点的新位置 令[],其中[]表示取大于数*的最小整数 l2+1当}>1/2 (4)第n个质点的随机徘徊是否结束 当n<M时,或 l;|≠0 时,第n个质点的随机徘徊尚未结束,令m=m+1,转至步骤(3);否则第n个质点的随机 徘徊结束,记录贡献 R=R+n(∑1l21=0) (12) 进入下一步骤。 (5)N个质点的随机徘徊是否已完成 当n<N时,N个质点的随机徘徊尚未完成,令n=n+1,转至步骤(2);否则,N个质 点的随机徘徊已完成,进入下一步骤。 (6)给出计算结果 Pa≈Py=R 比较解决随机徘徊问题的上述两种不同方法,以S=3和M=2000.例,在平均每秒可 完成250万次四则运算的 CYBER70/825机上计算,确定性方法的情况是,需要计算机CP 时间26081分钟(其中所有阶乘计算都采用了节省机时的办法:事先算好N!存放在数组 元素F(N)中,需要N!时,直接调用F(N),P≈0.341;蒙特卡罗方法的情况是,抽样总 数N=1000,要计算机CPU时间74.43分钟,P≈0.3345,误差为0.002924(置信率为 0.95)。 5
进一步比较两种方法的计算量与M的关系。由于式(6)中的取和数与M2同阶,因此 确定性方法(9)式的计算量与M成正比。至于蒙特卡罗方法的计算量与M的关系,则很明 显,是与M成正比。于是,为了使上述回返概率的计算结果更精确些,比如进一步取M= 2000,即比原来的M大一个数量级,由确定性方法需要计算机CPU时间约260810分钟≈ 半年,而蒙特卡罗方法仅需要744.3分钟≈12.4小时。显然,后者要比前者更实际一些。尤 其是三维以上的随机徘徊问题,情况将更是如此。 3随机误差干扰问题 因素x之值可由试验者使制,对x的响应之指标值为y,由于有随机误差的于扰,y对x 的依赖关系实际上是 y (14) 其中ε为随机误差。对于确定的y=y“和任意的e,用x(e)表示如下方程 h(x, e) (15 中关于x的解。所谓随机误差于扰问題是,扰到这样的x=x,使得随机变量x(e)取x 的可能性最大。 解决上述随机误差干扰问题存在两大困难。第一个困难是,h(x,E)的形式是未知的; 第二个困难是,问题中存在随机误差e的干扰。 1951年, robbins和Mom第一次研究了满足如下条件的随机误差干扰间题9:随机误 差∈只影响y,即h(x,e)具有如下形式: h(x h(x) h(x)为x的递增函数增加速度不快于线性;e服从均值为零的对称分布。所给出的算法 是,对于x的初始近似x1,用下式确定出xM ),m=1,2,…,M-1 其中ym为当x=x时,y的响应使;bn>0,并满足条件 bn=∞ ∑ (18) 当M足够大时,用x作为x“的近似。 用蒙特卡罗方法能否解随机误差干扰问题呢?为了使这件事是可行的,假设存在与随 机误差e无关的两个常数A和B,对于任意的e,方程(15)式在(A,B)上有唯一的根;所要计 算的x为随机变量x(e)关于e的数学期望: x=E(xe)tE (19) 由于只要x(∈)服从正态分布,条件(19)式就一定满足,因此,总的来说,上述条件不仅 比 Robbins-Monr模型的限制条件弱,而且,基本上包括了实际中所可能遇到的大多数问题。 用x表示在(A,B)上服从均匀分布的随机变量,根据对h(x,)的值设,不难确定,对 于任意的ε有 x(e)=A+(B-ae(g(x, e):01x (20) 其中: g(x, e)=(y-h(a, e)(y-h(x, e)) (21) 将式(20)代入到式(19)中,可以进一步得到 6
x"=A+(B-Ae(n(g(x, E)20lx,e) 根据此式,便有解随机误差干扰题的蒙特卡罗方法如下:于(A,B)上随机确定N个试验点 x,x2,…,x,满足均匀分布条件,用 X”=A+(B=)× )≥0 作为x的近似估计,其中y为当x=A时,y的响应值 比较解随机误差干扰问题的上述 Robbins- Monp方法和蒙特卡罗方法,前者曾被荣为随 机通近方面的开创性工作,后者却是一种普通的蒙特卡罗技巧,而且,前者要求满足的条件 较苛刻后者要求满足的条件较一般。这一对比情况同样表明,对于许多其他方法所难以解 决的随机性问题,蒙特卡罗方法可以比较方便地解决 4不公平博弈问题 总共有M个人参加博弈,他们的编号分别为1,2,…,M.博弈所用的工具是,点数为0 1,2,…,M的牌各M张。博弈的方法是,首先将零牌分给每人各一张,剩下的牌机会均等地 分发给每人各M张然后按照如下规则进行抽牌(抽出后放回):由编号为1的人率先在自 己的牌中抽出一张,若为零牌,即为获胜者;若为其他牌,此牌的点数即为下一个在自己牌中 的抽牌者的编号,估此类推,直至抽出零牌为止,谁首先抽出零牌谁就是获胜者,依照上述 博办法,很明显,编号为1的人比其他人有更多的机会获胜,其他人获胜的机会则是均等 的。所谓不公平博弈问题是问,编号为1的人获胜的机会比其他人多多少 用n,表示在某向中编号为讠的人得点数为j的牌的张数,则不难确定,如下线代数方 程组的解x1 11y r2lyn2:…,n2M 0 nvI, ni,', nMLxM 恰为编号为1的人在这局博弈获胜的概率.用P(n)和q(n)依次表示分牌结果为n,j =1,2,…,M的概率和编号为1的人获胜的概率。于是,很明显,编号为1的人获胜的平均 概率为 Q=∑q(n)p(n2) (24) 由于除编号为1的人以外的其他M-1个人,获胜的机会是相同的即获胜的平均概率均为 (1-Q)/(M-1),因此,编号为1的人获胜的机会比其他人多 c=Q 1 M 由以上所述,很明显,解不公平博穿问题的主要内容是,解具有随机系数的线代数方程 组(23)式和确定其解x1=q(n)的数学期望式(24).由此可见,用确定性方法解不公平博弈 问题是非常困难的 用蒙特卡罗方法解不公平博弈问題,情况将如何呢?按照不公平博弈问题中的博弈规 则,若在K,k=1,2,…,M2中以任意的顺序存放自然数1,2,…,M各M个,则有解不公平
博弈问题的蒙特卡罗方法如下 (1)为计算作准备 令n=1;x1=0 (2)开始第n盘的博弈 令i=0 (3)进行分牌 令i=i+1。在牌K,K+1,…,K是机会均等地抽出一张K分给编号为1的人: (M2-i+1)e),t=〔x1 (26) 将牌K与牌K进行对调以保证没有发出的牌为K+1,K2+2,…,KM (4)分牌是否结束 当i<M2时,分牌尚未结束转至步骤式(3);否则,分牌结束。 (5)由编号为1的人率先抽牌 令i=1 (6)确定下一抽牌者或结束此盘博弈 令 〔(M+1)〕 (27) 当j≠0时,下一抽牌者为i=KM(;-1)+,重复步骤(6);当j=0时,此盘博弈结束。 (7)记录贡献 (28) (8)N次博弈是否已全部完成 当n<N时,N次博弈尚未全部完成,令n=n+1,转至步骤(2);否则,N次博弈已全 部完成。 (9)给出计算结果 C 取M=2,4,8,16;N=500,上述解不公平博耷问题的蒙特卡罗方法,在 Cyber0/825机 上进行了计算,计算结果误差和所需要的计算机CP时间,列人在表1中.从所需计算机 时间看,计算上述全部4种情况, 襄1解不公平博弈问题的计算结果 别 Q的误差 所需时间(分) 0.4126 0.004315 0.2l 5638 0.2241 0,003655 0.1133 0833 0.1166 0.002813 0.0577 7.9833 8