《蒙特卡罗方法.docx》由会员分享,可在线阅读,更多相关《蒙特卡罗方法.docx(11页珍藏版)》请在得力文库 - 分享文档赚钱的网站上搜索。
1、蒙特卡罗方法蒙特卡罗(Monte-Carlo,简写为M-C)方法属于计算数学的一个分支,它是在二十世纪四十年月中期 为了适应时原子能事业的进展而进展起来的,但它与一般计算方法有很大区分,一般计算方法对于解决 多维或因素简单的问题特别困难,而蒙特卡罗方法对于解决这方面的问题却比较简洁。因而蒙特卡罗方 法在近十年来进展很快,特殊是随着快速电子计算机的进展,蒙特卡罗方法得到了快速进展与广泛应用。 蒙特卡罗方法也称随机抽样技术(Random Sampling Technique)或统计试验方法(Method of Statistical Test)o蒙特卡罗是欧洲摩纳哥国的一个重要城市,以赌博著称。蒙
2、特卡罗方法是以概率论与数理统计学为 基础的,是通过统计试验达到计算某个量的目的。而赌博时,概率论是一种有力的手段。所以,以蒙特 卡罗作为方法的名字,缘由也许于此。由于蒙特卡罗方法是采用一连串的随机数来求解问题的,因此求解随机过程,放射性衰变和布朗运 动等问题,它是很有效的。它除了在原子能工业广泛应用外,在物理、化学、地质、石油、线性规划、 计算机研制、计算机模拟试验、解决多体问题等领域中都有不同程度上的应用。第一节.蒙持卡罗方法的基本思想、特点及其局限性一、 蒙特卡罗方法的基本思想用下述三个例子,说明蒙特卡罗方法的基本思想。例1.产品合格率的计算某工厂生产一批产品,其合格率表示是:P(合格率)
3、=M(产品中合格严显总数)N(产品的总数)(1.1)为了确定合格率,应检查这批产品的全部,确定其中合格的数目。但是,由于产品数量多,检查全部产 品花费的代价大。因此,通常实行抽取部分产品,在这部分产品中确定其合格的数目。然后用这部分产 品的合格率3(部分产品合格率)=-(1.2)(1.2)N(部分产品的总数)来代替所要计算的合格率P。例如,检查某批产品,当被检查的产品长度介于13.60cm13.90cm内 时,则认为是合格的,否则是次品。分别抽取5件,10件,60件,150件,600件,900件,1200件, 1800件来检查,其状况如下表和图20所示。Z产品抽样统计表图1.1产品抽样统计图1
4、 .对于结定的一个P位数的初值r。(例如取r=0.3415926, 0.27182818, ),对r0进行平方,得到一个2P位数堤,实际上由于计算机字长有限,尾数长度是有限制的,当婿的位数超出计算机尾数 长度的时候,机器自动会舍掉多出的尾数,这样截去尾数避开了任意数的平方的末位数只能消失。,1, 4, 5, 6, 9而不会消失2, 3, 7, 8的系统偏倚性。2 .推断厂是否小于1,假如。21时,则右移4的小数点直至厂并时为止,然后,截去整数部分, 保留小数部分,把保留的小数部分作为(0, 1)区间上的匀称分布的随机数门,这样截去首位避开了小于 1的数的平方有对小数目偏倚的现象。3 .以口为初
5、值,重复1, 2过程产生戏。这一过程始终重复下去,就可产生(。,1)区间上匀称分布 的随机数序列A r2, .。例题例1.取初值YFL=0.31415926所产生的2000个(0, 1)区间上的匀称随机数进行统计特性检验。1)参数检验最大值 max(n,12,,200。)= 0.99907696最小值 min(口,2,r2ooo)= .17911196X 103一阶矩 ooo =0.48966279二阶矩低。=0.31479022二阶中心矩 52()()()=0.845612382)匀称性检验将(0, 1)区间等分成39个小区间,对2000个随机数进行匀称性检验,由计算得: X2(38)=45
6、.316(45.31653.33354)在显著水平a=0.05下,接受匀称性假设。例2.取YFL=0.27182818产生2000个(0, 1)区间上均勾分布的随机数。计算结果最大值 max(n,殳,12000,)=0.999828589最小值 min(n,r2,,r2oooj=0.57747960X IO-3,一阶矩 Gooo =0.49808762.附注1 .设rl为(。,1)区间上的匀称分布的随机数,则Ri=a+(b-a)ri,即为(a, b)区间上匀称分布的随机数。2 .本方法程序简洁,只要一个存贮单元就够了。但在使用上要留意初值的选取,初值选得好可以 产生周期较长、统计特性较好的随机
7、数。产生匀称分布随机数的方法很多,这里就不一一列举了。下面给出用M-C方法计算多重积分的 程序:SUBROUTINE RJ (N, N” YFL, RY)F=0.0DO 5 1=1, Ni5 F=F 十 FX(N)RY = F/FLOAT(Ni)RETURNEND SUBROUTINE RAND(YFL)随机数程序段体部分END程序功能本程序采纳蒙特卡罗(MC)方法求多重积分1=/:何力,,也,的近似值。使用说明1 .子程序语句SUBROUTINE RJ(N,Ni,YFL,RY)2 .哑元说明N:整变量,输入参数,积分重数。N1:整变量,输入参数,一个充分大的正整数。YFL: 同子程序RAND
8、的哑元说明。RY:实变量,输出参数,存放积分结果。3 .所调用的过程SUBROUTINE RAND(YFL), 表示产生(0, 1)区间上匀称分布的随机数的子程序段。FUNCTION FX(N),表示计算被积函数的函数段,由使用者自行编写,其自变量由于程序RAND 产生的随机数求得,N为积分重数。例题取S = 2, 3,,8计算多重积分T:宫灯,的结果。精确结果当产生随机数的初值取YFL=0.31415926时,计算结果如下表:S2 I134516!178i I0.61.01.31.62.0*2.32.6I、N 、1 _1567810000.6440.984(1.2741.6171.9792.
9、3132.64720000.6500.9631.3071.650一一 1.989. 一一一一2.3152.64530000.6500.9661.3221.6531.9862.3162.64640000.6460.979la 3241.6531.9832.3172.647可见,计算结果的误差小于10%。上表可看作八次试验,从结果中看出,随着抽取件数的增多,合格率愈来愈趋于一个稳定值。.9。假如 定义随机变量(b第z次试验为正品,幼=(0,第,次试验为次品。则做了N次试验后,正品个数共为Nr= 2劭3(1这样,(1.2)式可进一步写为P 仪二 2co i v 0 9 N N占人们从阅历中还知道,当
10、N数目越大,r/N作为正品率的估量值就越精确的可能性也越大。类似这种把观看正品消失的频率作为近似概率的例子在生产中是很常见的。例2.射击问题(打靶嬉戏)设r表示射击运动员的弹着点到靶心的距离,g(r)表示击中r处相应的得分数(环数),分布密度函数 f(r)表示该运动员的弹着点分布,它反映运动员射击水平。积分=(1 3)表示这个运动员的射击成果。用概率语言说,就是随机变量g(r)的数学期望值,即Vg=Eg(r)。 现在,假设这个射击运动员射击N次,弹着点依次是口&,,则自然地认为N次射击得分的平均值$8(兀) 1-1(1.4)相当好地代表了这个射击运动员的成果。换句话说,gN是积分(1.3)式的
11、一个估量值(或近似值)。这个 例子通常称为打靶嬉戏,它直观地说明白蒙特卡罗方法计算定积分的基本思想。为进一步阐明这个思想, 我们再举个例子:计算积分1= J: fWdx (OWfG)WD直观上,就是在边长为1的正方形里随机投点,当点落在y=f(x)曲线下面(见图21),对积分值有“贡 献”,否则对积分值无“贡献:为此,设自是0, 1上匀称分布的随机变数,定义1,第2次试验引刀(短)时,0,否则,(1.5)就是积分I的一个估量值。例3.求解三维椭圆型偏微分方程的边值问题考虑三维空间中的一个体积V上的Laplace方程:d2fd2fd2fdx2dyzdz2=0,(1.6)f | S = 0(P)#
12、其中S为边界,要求f在V中的数值。用蒙特卡罗方法求解是:将V用等距dx=Jy=J2的网格剖分,在扩中任一内点Q,要求f(Q)之值。 用一个骰子投掷,其点数1, 2, 3, 4, 5及6设分别表示向X轴正向,负向,y轴正向,负向,2轴正 向,负向移动一步。由Q点动身,按投骰子的方法,在严内进行移动,当此点到达边界时,登记边值 之数值;又回头由Q点按投段子之点数移动,如此连续下去,设得到一系列的边值:外,2,外,则f(Q)的近似值可取I Nf (Q) 2f(Q) =留意,这种作法的来源是依据两个基本假设,其一为方程(1.6)可以用差分格式f 1+ fl- l,j,k + fi,i+ l,k + f
13、bj- lk + fI,j,k+ 1 + fl+j,k - 1) b来近似,其二是假设骰子为匀称的,即消失6个数机会均等。从上述三个例子可以看到,当所要求解的问题是某种大事消失的概率,或者是某个随机变量的期望 值时,它们可以通过某种“试验”的方法,得到这种大事消失的频率,或者这个随机变数的平均值,并 用它们作为问题的解。这就是蒙特卡罗方法的基本思想。由上面例子还可看出,用蒙特卡罗方法求解问题时,首先要建立一个随机模型,然后要制造一系列 的随机数用以模拟这个过程,最终要作统计性的处理。关于建立随机模型,因问题而异。下面来争论随 机数的产生方法。二、随机数和伪随机数的产生1.单位矩形分布最简洁、最
14、重要、最基本的一个概率分布是单位区间(0, 1)上匀称分布(或称单位矩形分布),其分 布密度函数是:(1,当 0 至0,其它点处.2 .随机数由具有单位矩形分布的总体内抽取的简洁子样:Xi, X2,,Xn,简称为随机数序列。其中的每 一个体称为随机数。由于它在蒙特卡罗方法中占有极其重要的地位,我们将用特地的符号家包, 松来表示。连续产生随机数的例子如投骰子,有奖储蓄开奖时的摇号码机等。3 .伪随机数及其产生的方法计算机不会掷骰子,它是采用数论的方法来产生随机数的。由于这种方法属于半阅历性质,因此只 能近似地具备随机数的性质,所以称为伪随机数。最初冯诺伊曼(Von Neumann)建议的“平方取
15、中法” 如下;以十进制数为例,平方取中法是把一个28位的十进制数自乘后,去头裁尾保留出间23个数字, 然后用102s来除。例如=6406, 父2 = 0368, %3 = 1354, % 4 = 8333 . %6 = 4388,讨=41036836, 港= 00135424, %; =01833316, x =694388899相应自6伪随机数序列便是0.6406, 0.0368, 0.1354, 0.8333, 0.4388,。平方取中法的一般形式是(占+1 三110f%|2(modl()2s),/ = _ /+14i+i 一,其中“mod”表示取模运算,X为任意给定的初始值,由2s位十进
16、制数组成。在电子计算机上采纳二 进制的数,此时平方取中法如下叫+i 三2-s%彳(mod22s)Xi为任意给定的2s位二进制数。同平方取中法相类似,乘积取中法的一般形式是x1+2 =C2$X|X|+iD (mod2Xs) 9b 一 4l+2gi+Z -2 2Xi,X2为任意两个初始值,由2s位二进制数组成。乘同余方法的一般形式是 xI+I =ax (modN),Xi为任意初值。乘加同余法的一般形式是科1 三(/+c) (modN) 9其中%。= Seed) (0WSeedW199617) c= 99991,a 24298, N = 199617.关于产生伪随机数的方法,性能以及检查已有很多工作
17、。其实,在大多数计算中,产生匀称分布的 伪随机数已作为内部函数或库文件存在计算机内,随时可以调用。三、蒙特卡罗方法的特点及其局限性对于多种介质,几何上不规章,维数高等类型的问题,蒙特卡罗方法照样可用,不须修改。因此,在肯定意义上说,它是一种很普遍适用的方法。而它的误差 是由概率论中的大数定律limP/oUS 承一n-8 I 2v i=)所掌握的,并由中心极限定理Pr 牌配-P Y芳卜? e a,这表明,不等式X rrP-P JN (.7)以概率1-a成立。通常a称为置信度,1-a称为置信水平。假如必。,则由(L7)式知,蒙特卡罗方法的误 差为U _ Xq。N (1.8)a和Xa的关系可从正态分
18、布N(0, 1)积分表中查到。常用的几组a和Xa如下a0.50.67450.051.960.0113特殊称a=0.5时的误差0.6745a/ J犷为概然误差。再如,取置信水平为95%,则Xa=L96,此时表明误差不等式:以95%的可能性具有精确度为E= 1.960a/07。所以,MC方法对于误差的估量具 有概率性质。即对于这个方法不能断言误差不超过某值,而只能指出误差以某种(如接近1)的概率不超 过某值。还可看出,当给定置信度a后,误差E由。和J犷打算。要减小E,或者是增大N,或者是减 小方差。2。在固定o下,要提高精度一位数字,就要增加100倍工作量,因此,单纯增大N,不是一个 有效的方法。
19、改进的方法之一是削减方差这里有很多技巧,如“重要抽样”、“系统抽样”、“相 关、“对偶变数”、”半解析方法”、“统计估量”、 “分裂与俄国轮盘赌”等等。蒙特卡罗方法主要有如下三个特点:1.收敛速度与问题维数无关从误差表达式(1.8)中,明显看出,在肯定的置信水平下,MC方法的误差,除了与方差d有关 外,只取决于于样容量N,而与于样中的元素所在的集合空间*的组成无关。问题的维数变化,除了引 起抽样时间及计算估量量的时间有变动外,不影响问题的误差。换言之,要达到同一精度,用M-C 方法选取的点数与维数无关;计算时间仅与维数成比例。但一般数值方法,比如在计算多重积分时,达 到同样的误差,点数与维数的
20、事次成比例,即计算量要随维数的基次方增加。这一特性,打算了 MC 方法相宜于多维问题。例如,计算积分:1 of(X =xdx o用mc方法是抽取随机数&,o和i,令co=0,其它.于是P 1N1 = 1(1.9)就是积分Pl的一个估量值。由于均方差b = L, P. =,依据误差不等式,知 22色一小母,于是相对误差|/P1.96,N(1.10)将上述积分推广到s维的情形,取%2, 一, %s) = (%1 + 汇2 + 一 + Xs)/S,积分变成Ps =J 0 Jo13h S xidx1dx2 dx.3 i(1.11)此时,MC方法的估量同样是 NPs =方S小N 1(1.12)其中1 s
21、1,当心您g = K手心0,其它。(1.13)由此可见,用MC方法计算单重积分P和S重积分R,除了确定coi值时与问题的维数有关外(前者是国庐言,后者抵()1 5孔,仅多产生s1个随机数),其它与维数无关,尤其重要的是误差与维 i=数无关。2.受问题的条件限制的影响小比如计算S维中任意一个区域Ds上积分Psd =DsJg(%i,,xdxxdxz -dx(1.14)无论Ds的外形如何特殊,只要能给出描述Ds特性的几何条 件,那么总可以给出类似于(1. 4)式的估量如下Ad = # ,P?)(1.15)其由2是对全部满足(X;,X;X:)Ds的点取的。而其它数值方法受问题的条件限制影响比较大。可见
22、,M-C方法是解决简单的几何模型的一种有效的方法。3 .程序结构简洁,占用内存单元较少最终,在电子计算机上实现MC计算时,程序结构清楚简洁。例如,用式(1.15)计算积分(1.14) 的MC方法程序结构如下面框图所示。其中假设D s居于S维单位超立方体,AN表示每组抽样的数目,那后:为随机数,X。为常数,依据置信水平要求确定,E表示对计算相对误差的要求。4 .蒙特卡罗方法的局限性前面已经指出,蒙特卡罗方法的弱点是收敛速度慢,误差的概率性质,误差与点数(抽样数N)的平 方根成反比,而其它数值方法的误差是准确的,通常状况下,误差与点数成反比。因此,对于维数不高 的问题(一维、二维),数值方法可以结
23、出误差很小的结果,而且计算量小。也就是说,MC方法计算 量相对说来比较大,又给出的是概率误差,STA RTIE E2 = 0IN=n=u您上通常给出二至三位有效数字。除此之外,经证明,只有当系统的大小与粒子的平均自由程可以相比较时, 一般在10个平均自由程左右,MC方法算出的结果较为满足。而对于大系统、深穿透问题(一般指系 统大小为1620个平均自由程以上),算出的结果往往偏低。而对于大系统,数值方法往往很适应,能 算出较好的结果。因此,现在已有人将数值方法与蒙特卡罗方法联合起来使用,相互取长补短,克服这 种局限性,取得了肯定的效果。其次节蒙持卡罗方法计算积分计算多重积分是蒙特卡罗方法的重要应
24、用领域之一。为任一积分,都可看作某个随机变量的数学期 望。因此,以用这个随机变量的算术平均值来近似。一、蒙特卡罗方法求积分例如求积分 =f(x)dxJ (2.1)选取分布密度函数(1,当e1o n, P(x) = I 0,其它.则式(2.1)可写成(2.2)即。是随机变量f(x)的数学期望。现从p(x)抽取随机变量X的N个样本:Xi, i=l, 2,,N,则算术 平均值(2.3)就是。的一个近似估量。对于一般的S重积分5= u G(P)dPJ V(2.4)其中,P=P(Xi, X2,Xs)表示S维空间的点,Vs表示积分区域。用MC方法求解的步骤是:1 .在所求解约区域上构造一个分布密度函数。取
25、Vs上任一概率密度函数f(P),它满足如下条件:f(P)g 当 PWPs, G(P)于0 时.当/(P)妾o,g(P)=g(P)=0,其它.(2.5)则式(2.4)可改写为(2.6)e= jFs g(P)f(P)JP=ECg(P)3即。是随机变量g(P)的数学期望。2 .将g(P)的数学期望值用其算术平均值来近似。现从f(P)抽取随机向量P的N个样本:Pi,i=l, 2,,N,则g、1寺 g(P)IN i= Ig、1寺 g(P)IN i= I(2.7)就是积分。的近似估量。选取f(P)的最简洁方法是取Vs上的匀称分布:1/暝,当 PCPs, f(P) =(2.8)(2.8)I 0,其它.这里,
26、Vs也表示积分区域的体积。这时= ysG(p) (29)明显,存在一个如何选取f(P)的问题。选取最优的估量形式,将是MC方法要争论的问题之一。 二、无偏估量假如随机变量Y,其期望值为0*,即班小明 (2.10)则称Y是0*的一个无偏估量。不难验证,8田)和都是。的无偏估量。三、收敛性及误差估量 由上节可得如下结果:1 .当Nfoo时,抽样(算术)平均以概率为1收敛到e即Proi lim gN = 8 )=1 j /(2.11)2 .假如有不为零的有限方差存在,即f Cg(尸)一尸)dFv+8J 吸(2.12)存在,则在l-a置信水平下,成立不等式:|gN - 8 | VXaOg/也就是说,用
27、MC方法求积分,其误差的阶为。(N2).最终,给出用M-C方法产生匀称分布随机数相计算多重积分的程序。 产生匀称分布随机数的程序:SUBROUTINE RAND(YFL)IF(YFL. EQ. 0.0) YFL=0.31415926 Y = YFL*YFL1 Y = Y*10.0 IF(Y-l.O) 1, 2, 2YFL=Y-FLOAT(IFIX(Y)RETURNEND程序功能 本程序用于产生(0, 1)区间上匀称分布的随机数。 使用说明1 .子程序语句SUBROUTINE RAND(YFL)2 .哑元说明YFL:实变量,输入、输出参数,要求在第一次调用前,对 YFL赋值, 当YFL=0.。时,程序自 带初值:YFL=0.31415926;当使用者自己选定初值时,可将选定的初值赋于YFL。以后每调用一次便 在YFL中产生一个(O, 1)区间上的随机数,同时,YFL又作为产生下一个随机数的初值。方法简介本程序采纳的计算步骤如下:
限制150内