50蒙特卡罗算法

蒙特卡罗（Montecarlo）方法【实验目的】z了解Montecarlo随机摸拟方法的具体应用z学习使用MATLAB软件中有关随机数生成函数【实验准备与内容】1．Montecarlo方法蒙特卡罗(MonteCarlo)方法，又称计算机随机模拟方法，是一种基于随机数的计算方法。这一方法源于美国在第二次世界大战研制原子弹的曼哈顿计划。该计划的主持人之一数学家冯诺伊曼用驰名世界的赌城-摩纳哥的MonteCarlo来命名这种方法。MonteCarlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的频率来决定事件的概率。1777年，蒲丰（Buffon）提出著名的Buffon投针试验永来近似计算圆周率π。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得人们在计算机上利用数学方法大量、快速地模拟这样的试验成为可能。目前这一方法已经广泛地运用到数学、物理、管理、生物遗传、社会科学等领域，并显示出特殊的优越性。2．伪随机数实际应用中的随机数通常都是某些数学公式计算而产生的伪随机数，这样的伪随机数从数学意义上讲不是严格的随机数，但是，只要伪随机数能够通过随机数的一系列统计检验，我们就可以将它当作真随机数而放心使用。这样我们就可以方便、经济、重复地产生随机数。理论上要求伪随机数产生器具备以下特征：良好的统计分布特性，高效率的伪随机数产生，伪随机数产生的循环周期长，伪随机数可以重复产生等。到目前为止，已经提出了各种分布的伪随机数产生方法，在这里，我们不准备从数学上介绍随机数的产生原理，而只给出伪随机数的MATLAB产生函数。z生成各类分布的随机数统一生成函数random('name',A1,A2,A3,m,n)——生成以A1,A2,A3为参数分布为name的nm×阶随机数矩阵，其中，‘name’为包含特定分布名称的字符串(如表*.*)，A1,A2,A3是分布参数矩阵，根据分布不同，各参数的含义也不相同，且其中一些参数也不是必须的。z针对不同分布的随机数产生函数MATLAB还提供了不同分布的随机数产生函数，这些函数的具体使用格式及参数含义较random函数更明确（如表*.*）。表*.*特定分布字符串‘name’分布名称随机数产生函数函数格式'beta‘β分布BetarndBetarnd(A,B,m,n)'bino‘二项分布BinorndBinornd(N,P,m,n)'chi2'2χ分布chi2rndchi2rnd(V,m,n)'exp‘指数分ExprndExprnd(MU,m,n)布'f'F分布frndFrnd(V1,V2,m,n)'gam'γ分布GamrndGamrnd(A,B,m,n)'geo'几何分布GeorndGeornd(P,m,n)'hyge'超几何分布HygerndHygernd(M,K,N,mm,nn)'logn'对数正态分布LognrndLognrnd(Mu,Sigma,m,n)'nbin'负二项分布NbinrndNbinrnd((R,P,m,n)'ncf‘非中心F分布NcfrndNcfrnd(Nu1,Nu2,Delta,m,n)'nct‘非中心t分布NctrndNctrnd(V,Delta,m,n)'ncx2'非中心分布2χncx2rndncx2rnd(V,Delta,m,n)'norm'正态分布NormrndNormrnd(Mu,Sigma,m,n)'poiss'泊松分布PoissrndPoissrnd(Lambda,m,n)'rayl'Rayleigh分布raylrndRaylrnd(B,m,n)'t't分布TrndTrnd(V,m,n)'unif'连续均匀分布unifrndUnifrnd(A,B,m,n)'unid'离散均匀分布UnidrndUnidrnd(N,m,n)'weib‘Weibull分布weibrndWeibrnd(A,B,m,n)为了说明函数的使用方法，我们举两个例子。例：用random函数生成均值为0，标准差为2的正态分布随机数。rn=random('norm',0,2,10000,1);％生成正态分布随机数hist(rn，30)％绘出随机数的直方图-10-50510020040060080010001200例：生成区间[1,5]的服从均匀分布的随机数rn=unifrnd(1,5,10000,1);％生成均匀分布随机数hist(rn，30)％绘出随机数的直方图123450501001502002503003504003.Buffon投针试验在平滑的桌面上画一组相距2a的平行线束，向此桌面上投一枚长为2b的细针，为避免针与两平行线同时相交的复杂情况，假定。现在，我们来求针与平行线相交的概率。0ba设M为针的中点，y为M与昀近平行线的距离，x为与针与平行线的交角（π≤≤x0），于是，针与平行线相交的充要条件是xbysin0≤≤,故相交的概率为图中曲线xbysin=与x轴所夹图形的面积占图中长方形面积的百分比。xbysin=oaπxy即abxdxbapπππ2sin10==∫从而pab2=π用N表示投针的次数，n表示其中针与平行线相交的次数，由贝努里（Bernoulli）定理知，当N充分大时，频率接近于概率，即Nnp≈，从而nabN2≈π。这就是Buffon用随机实验近似求圆周率π的基本公式。我们不妨用产生随机数的办法来模拟投针实验：对每一次投针实验，针与平行线的交角x（π≤≤x0）和针的中点M与昀近平行线的距离y（ay≤≤0）是一随机变量，其分别服从],0[π和上的均匀分布。因此，N次投针实验可以通过以上产生随机数的函数生成，然后根据针与直线相交的充要条件，即可判断针与直线是否相交，从而可以计算出n的值，达到近似计算圆周率],0[axbysin0≤≤π的目的。具体程序如下function[pai,number]=buffon(a,b,N)%2a，2b分别为平行线间距和针长，N为投针次数x=unifrnd(0,pi,N,1);y=unifrnd(0,a,N,1);number=0;%相交计数器fori=1:Nify(i)=b*sin(x(i))number=number+1;endendpai=2*b*N/(a*number)其计算结果如下：abN相交次数π的估计值45361000051393.11344536100000509923.137745365000002547223.14074．积分计算定积分的计算是MonteCarlo方法引入计算数学的开端，在实际问题中，许多需要计算多重积分的复杂问题，用MonteCarlo方法一般都能够很有效地予以解决，尽管MonteCarlo方法计算结果的精度不很高，但它能很快地提供出一个低精度的模拟结果也是很有价值的。而且，在多重积分中，由于MonteCarlo方法的计算误差与积分重数无关，因此它比常用的均匀网格求积公式要优越。4．1随机投点法（1）设计算的定积分为，其中为有限数，被积函数是连续随机变量dxxfIba∫=)(ba,)(xfς的概率密度函数，因此满足如下条件：)(xf)(xf非负，且（1）1)(=∫+∞∞−dxxf显然I是一个概率积分，其积分值等于概率)(baP≤ς。下面按给定分布随机投点的办法，给出如下MonteCarlo近似求积算法：)(xfstep1:产生服从给定分布的随机变量值Nixi,...,2,1,=；step2:检查是否落入积分区间。如果条件ixbxai≤满足，则记录落入积分区间一次。ix假设在N次实验以后，落入积分区间的总次数为n，那么用ixNnI=_作为概率积分的近似值，即NnI≈（2）如果要计算的定积分为，其中为有限数，但被积函数不是某随机变量的的概率密度函数。在这种情况下，当然我们可以对积分作变换，将积分变换成满足（1）的条件。但在变换以后，需要产生新的分布随机变量，因此常会遇到很多困难和比较复杂的计算。dxxfIba∫=)(ba,)(xf上述求积方法需要产生给定分布的随机变量，它适合于解决特殊类型的概率积分。如果只用随机数完成随机投点，那么，下面的方法可以解决较为广泛的一类积分问题。（3）设是[0,1]上的连续函数，且)(xf1)(0≤≤xg，需计算积分，如下图所示阴影部分的面积。在图中单位正方形内均匀地投点dxxfIba∫=)(),(ης，则该随机点落入曲线)(xfy=下面的概率为IdxxfxfyP==≤∫10)()}({因此，给出如下MonteCarlo近似求积算法：step1:产生两组[0，1]区间内的随机数Niyxii,...,2,1,,=，并把作为随机点),(iiyx),(ης的可取坐标；step2:检查（，是否落入ix)iyΩ内，如果条件)(iixfy≤满足，则记录（，落入积分区间一次。ix)iy假设在N次实验以中，落入Ω内总次数为n，那么量NnI=_近似等于随机点落入内的概率，即ΩNnI≈假如所需计算积分，其中为有限数，被积函数有界，并用dxxfIba∫=)(ba,)(xfM和分别表示其昀大值和昀小值。作变换，m*)(xabax−+=])(([1)(***mxabafmMxf−−+−=，此时)()())((10*abmdxxfabmMI−+−−=∫且有，即转化为上面讨论过的情况。]1,0[,1)(0**∈≤≤xxf4．2平均值法任取一组相互独立、同分布的随机变量}{iξ，iξ在[a,b]内服从分布率，令)(xp)()()(*xpxfxp=，则也是一组相互独立、同分布的随机变量，而且)}({*ipξ∫∫===babaiIdxxfdxxpxppE)()()()}({**ξ由强大数定理1)(1lim1*=⎟⎠⎞⎜⎝⎛=∑=∞→NiiNIpNPξ若记∑=−=NiiIpN1*)(1ξ，则−I依概率1收敛到I，平均值法就是用−I作为I的近似值。假如所需计算积分为，其中被积函数在[a,b]内可积，任意选择一个有简单办法可以进行抽样的概率密度函数，使其满足条件：dxxfIba∫=)()(xpz,0)(≠xp当时（0)(≠xfbxa≤≤）z1)(=∫dxxpba记⎪⎩⎪⎨⎧=≠=.0p(x),0,0)(,)()()(*xpxpxfxp则所求积分为dxxpxpIba∫=)()(*因而MonteCarlo近似求积算法为：step1:产生服从分布律的随机数)(xpNixi,...,2,1,=；step2:计算均值∑=−=NiipNI1*)(1ξ，即有II≈−。如果为有限数，那么可以取为均匀分布：ba,)(xp⎪⎩⎪⎨⎧≤≤−=其它,0,1)(bxaabxp此时dxabxfabIba∫−−=1)()(，具体的算法为：step1:产生[a,b]上的均匀随机数Nixi,...,2,1,=；step2:计算均值∑=−−=NiixfNabI1)()(，且有−≅II。4．3多重积分的计算根据上面定积分的Montecarlo计算方法，我们可以很容易地将它推广到重积分的近似计算。为了避免重复，在此只给出具体的二重和三重积分计算公式。)2(≥nnz二重积分的计算∑∫∫=≅NiiiDyxfNDdxdyyxf1),(||),(其中，是平面区域D中的均匀分布的N个独立选取的随机点列，表示区域D的面积。),(iiyx||Dz三重积分的计算∑∫∫∫=ΩΩ≅NiiiizyxfNdzdxdyzyxf1),,(||),,(其中，是空间区域),,(iiizyxΩ中的均匀分布的N个独立选取的随机点列，||Ω表示区域的Ω体积。很显然，该方法特别适合于多重积分，特别是积分区域比较复杂的重积分计算。4．4误差估计MonteCarlo求积方法常常以某个随机变量θ的简单子样Nθθθ,,,21L的算术平均值∑=−=NiiN11θθ作为积分I的近似值，由强大数定理知道，如果随机变量序列Nθθθ,,,21L相互独立、同分布、期望值存在，那么当时，以概率1收敛到∞→N−θI。按照中心极限定理，只要随机变量序列Nθθθ,,,21L相互独立、同分布、数学期望存在，且具有有限标准差0≠σ，那么，当当时，随机变量∞→NNIY/σθ