蒙特卡罗算法.ppt

Monte Carlo Simulation Methods(蒙特卡罗模拟方法)主要内容: 一. M-C方法概述. 二. 随机数的生成. 三. 模拟训练. 四. 实验题目.成信院数学与信息科学系 李胜坤1 蒙特卡洛蒙特卡洛(Monte Carlo)方法，或称计算机方法，或称计算机随机模拟方法，是一种基于随机模拟方法，是一种基于“随机数随机数”的计的计算方法这一方法源于美国在第二次世界大算方法这一方法源于美国在第二次世界大战中研制原子弹的战中研制原子弹的“曼哈顿计划曼哈顿计划”该计划的主持人之一、数学家冯的主持人之一、数学家冯·诺伊曼用驰名世诺伊曼用驰名世界的赌城界的赌城—摩纳哥的摩纳哥的Monte Carlo—来命名来命名这种方法，为它蒙上了一层神秘色彩这种方法，为它蒙上了一层神秘色彩一.M-C方法概述2基本思想很早以前就被人们所发现和利用17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”19世纪人们用投针试验的方法来决定π高速计算机的出现，使得用数学方法在计算机上大量模拟这样的试验成为可能3从Buffon(蒲丰)投针问题谈起 45试验者时间(年)针长投针次数相交次数π的估计值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592926数值积分问题选取(0, 1)中随机数序列x1， x2， x3， …… xn。

则误差约 ，它并不能和一些高级的数值积分算法比拟， 但对多维情况，MC方法却很有吸引力7我们可产生一系列随机数可简单取3个随机数构成一个随机点，即相应地，一般地，8Monte Carlo数值积分的优点与一般的数值积分方法比较，Monte Carlo方法具有以下优点：9M-C的基本思路1.针对实际问题建立一个简单且便于实现的概率统计模型，使所求的量（或解）恰好是该模型某个指标的概率分布或者数字特征2.对模型中的随机变量建立抽样方法，在计算机上进行模拟测试，抽取足够多的随机数，对有关事件进行统计3.对模拟试验结果加以分析，给出所求解的估计及其精度(方差)的估计4.必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率10回顾几种连续型分布1.均匀分布U(a,b)其概率密度函数为有11 均匀性特点:均匀分布随机变量X落在(a,b)内任意子区间的概率只与子区间的长度有关,而与子区间的位置无关. 可假设有这种特性的随机变量服从均匀分布.均匀分布随机变量X的取值具有”均匀性”122. 正态分布正态分布随机变量X的概率密度函数是正态分布由两个参数 和 唯一确定.其中 是X的均值(数学期望): =E(X),它确定了概率曲线的中心位置,而 是X的标准差: ,它确定了概率曲线的”宽窄”程度.13 在许多实际问题中,有一类随机变量可以表示成为许多相互独立的随机变量之和,而其中每个随机变量对总和只起微小的影响,这类随机变量往往服从或近似服从正态分布.在实际应用中,如果我们分析到一个随机变量受到较多独立的微小因素的叠加影响,就可以用正态分布来模拟这个变量.如:工厂产品的测量尺寸,农作物的收获量,某地区成年人的身高,体重等可看成服从正态分布的随机变量.143. 指数分布指数分布随机变量X的概率密度为指数分布常用来描述寿命问题.15二.随机数的生成1.蒙特卡罗模拟的关键是生成优良的随机数。

2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)我们通常称之为伪随机数(pseudo-random numbers)3.在模拟中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布U(0,1)的随机数16U(0,1)随机数的生成乘同余法： 称为种子,a 是乘因子,m是模数17一个简单的例子18一个简单的例子(续)上面的例子中，第一个随机数生成器的周期长度是 10，而后两个生成器的周期长度只有它的一半我们自然希望生成器的周期越长越好，这样我们得到的分布就更接近于真实的均匀分布19线性同余生成器(混合同余法)(Linear Congruential Generator )1.c是非负整数.通过适当选取参数c可以改善2.随机数的统计性质(独立性,均匀性).20常用的线性同余生成器Modulus mMultiplier aReference2^31-1=214748364716807Lewis, Goodman, and Miller39373L’Ecuyer742938285Fishman and Moore950706376Fishman and Moore1226874159Fishman and Moore214748339940692L’Ecuyer214748356340014L’Ecuyer21复杂一些的生成器Multiple recursive generator22算法实现许多程序语言中都自带生成随机数的方法，如 c 中的 random() 函数，Matlab中的rand()函数等。

但这些生成器生成的随机数效果很不一样，比如 c 中的函数生成的随机数性质就比较差，如果用 c ，最好自己再编一个程序Matlab 中的 rand() 函数，经过了很多优化可以产生性质很好的随机数，可以直接利用23从U(0,1)到其它概率分布的随机数1.离散型随机数的模拟2.连续型随机数的模拟3.正态随机数的模拟241.离散型随机数的模拟设随机变量 的分布律为令将 作为区间(0,1)的分点.若随机变量 ,有25令则有据此,可得产生 的随机数的具体过程为:每产生一个(0,1)区间上均匀分布随机数 ,若 则令 取值 .26例1:离散型随机变量X有如下分布律: X 0 1 2P(x) 0.3 0.3 0.4设 是(0,1)上均匀分布的随机数,令则 是具有X分布律的随机数.272.连续型随机数的模拟a.逆变换方法(常用) (Inverse Transform Method)b.舍取方法 (Acceptance-Rejection Method)定理: 设随机变量Y的分布函数F(y)是连续函数, 而U是在(0,1)上均匀分布的随机变量, 令, 则Y与X有相同的分布.28293031舍取方法舍取方法(Acceptance-Rejection )方法最早由 Von Neumann提出，现在已经广泛应用于各种随机数的生成。

基本思路：通过一个容易生成的概率分布 g 和一个取舍准则生成另一个与 g 相近的概率分布 f 32具体步骤：33下面我们验证由上述步骤生成的随机数 Y 确实具有密度函数 f(x) 34所以为了提高舍取法的效率，我们应该使 c 的取值尽可能的小，也就是使 f 和 g 的分布更为相近353.标准正态分布随机数的生成正态分布是概率统计中最重要的分布，在此我们着重讨论如何生成标准正态分布随机数引理：36Box-Muller 算法37利用中心极限定理设 是n个相互独立的在(0,1)上均匀分布的随机变量,由中心极限定理知 渐近服从正态分布N(0,1).一般取n=10即可,若取n=12,则上式简化为 再由公式 即可得到正态分布 的随机数.38Matlab程序Function r=rnd-u(a,b)%产生在[a,b]间均匀分布的随机数r=a+(b-a)*rand;return39Matlab程序Function r=rnd-beta(lamda)%模拟指数分布%lamda表示指数分布的参数r=-log(rand)/lamda;return40Matlab程序Function y= rnd(mean, segema)%模拟均值为mean,方差为segma的正态分布r=rand(1,12);x=sum(r)-6;y=segma*x+mean;return41三. 模拟训练例1. 模拟求近似圆周率在边长为1的正方形内有一半径为0.5的内切圆.现在模拟产生在正方形内均匀分布的点n个.如果有m个在圆内,则圆面积与正方形的面积比可近似为m/n.即л/4≈m/n л≈4m/n42n=10000m=0;For i=1:n if rand^2+rand^2 <=1 m=m+1; endendMypi=4*m/n43例2. 用M-C法估算定积分. 求定积分 .分析:对于 ,如果f(x)>=0,则可以通过模拟估算.构造一个矩形包含曲边梯形,d>max f(x).产生n(足够大)个在矩形区域内的点,如果落在由函数f(x)构成的曲边梯形内的点为m个,则所求定积分为 .44n=10^6;a=0;b=1;d=1;m=0;for i=1:n x=a+rand*(b-a); y=d*rand; if y<=x^2 m=m+1; endends=m/n*(b-a)*d45n=10^6;x=rand(1,n);y=x.^2;s=sum(y)/n采用前面讲的方法:46例3. 渡口模型问题描述:一个渡口的渡船营运者拥有一只甲板长32米,可以并排停放两列车辆的渡船.他在考虑怎样在甲板上安排过河车辆的位置,才能安全地运过最多数量的车辆.分析:怎样安排过河车辆,关心一次可以运多少辆各类车.准备工作:观察数日,发现每次情况不尽相同,得到下列数据和情况:(1)车辆随机到达,形成一个等待上船的车列;(2)来到车辆,轿车约占40%,卡车约占55%,摩托车约占5%;(3)轿车车身长3.5~5.5米,卡车车身长为8~10米.47问题分析:这是一个机理较复杂的随机问题,是遵”先到先服务”的随机排队问题.解决方法:采用模拟模型方法.因此需考虑以下问题: (1)应该怎样安排摩托车? (2)下一辆到达的车是什么类型? (3)怎样描述一辆车的车身长度? (4)如何安排到达车辆加入甲板上两列车队中的哪一列中去?本实验主要模拟装载车辆的情况,暂时不考虑渡船的安全.48模型建立:设到达的卡车,轿车长度分别为随机变量 , .结合实际,这里不妨设卡车,轿车的车身长度 , 均服从正态分布.由于卡车车身长8~10米,所以卡车车长 的均值为(8+10)/2=9米,由概率知识中的”3σ”原则,其标准差为(9-8)/3=1/3,所以得到 .同理可得 .49模拟程序设计:由以上的分析,程序设计时应划分的主要模块(函数)如下:(1)确定下一辆到达车辆的类型:(2)根据车的类型确定到达车辆的长度;(3)根据一定的停放规则,确定放在哪一列.50function r =makeid%模拟下一辆到达车的类型t=rand;if t<=0.55 r=1;%到达卡车elseif t<=0.95 r=2;%到达轿车else r=3;%到达摩托车end51function len=getlength(id) %根据车的类型,产生车长随机数switch id case 1, len=min([9+randn*(1/3),10]); case 2, len=min([4.5+randn*(1/3),5.5]); case 3; len=0;end52vfunction [full,pos]=getiffull(L,newlen)v%增加车长为len后是否可行(是否满)v%pos表示加到那一列去vfull=0;vpos=0;vif L(1)>L(2)v if L(2)+newlen<32v pos=2;v elsev full=1;v endvelsev if L(1)+newlen<32v pos=1;v elsev full=1;v endvend53vfunction sim_dukou v%渡口模型的模拟vn=input(‘输入模拟次数:');vif isempty(n)|(n<500)v n=500;vendvN=zeros(1,3);vfor i=1:nv isfull=0;v L=[0,0];v while ~isfullv id=makeid;v N(id)=N(id)+1;v newlen=getlength(id);v [isfull,pos]=getiffull(L,newlen);v if ~isfullv L(pos)=L(pos)+newlen;v endv endvendvdisp(‘平均每次渡船上的车数')vmean_n=N/n54四.实验题目赶上火车的概率问题描述:如图,一列火车从A站开往B站,某人每天赶往B站上这趟火车.他已了解到:(1)火车从A站到B站的运行时间是均值为30分钟,标准差为2分钟的随机变量;(2)火车在下午大约1点离开A站,离开时刻的频率分布如下:55出发时刻午后1:00午后1:05午后1:10频率 0.7 0.2 0.1此人到达B站的时刻频率分布为:时刻午后1:28午后1:30午后1:32午后1:34频率0.30.40.20.1问他能赶上火车的概率是多少?56变量说明 :火车从A站出发的时刻; :火车从A站到B站的运行时间,单位:分钟; :他到达B站的时刻问题分析与假设:此题包含多个随机因数.这里假设 , , 都是随机变量,其中 服从正态分布.57模型建立很显然,他能及时赶上火车的条件是: < + .为了简化计算,将下午1点记为0时刻. 和 的分布律如下: /min 0 5 10 P(t) 0.7 0.2 0.1 /min 28 30 32 34 P(t) 0.3 0.4 0.2 0.158如果r为在(0,1)均匀分布的随机数,为了模拟随机变量 , ,可以通过如下方法: 其中, 和 分别用来模拟随机变量 和 . 59vfunction pro_huochev%模拟赶上火车的概率vn=input(‘输入模拟次数:')vm=0;vfor i=1:nv r1=rand;v if r1<=0.7v t1=0;v elseif r1<=0.9v t1=5;v else v t1=10;v endv r3=rand;v if r3<=0.3v t3=28;v elseif r3<=0.7v t3=30;v elseif r3<=0.9v t3=32;v elsev t3=34;v endv t2=30+randn*2;v if t3<(t1+t2)v m=m+1;v endvendvp=m/nvend60。