蒙特卡罗算法ppt课件

1、Monte Carlo Simulation Methods(蒙特卡罗模拟方法)主要内容: 一. M-C方法概述. 二. 随机数的生成. 三. 模拟训练. 四. 实验标题.成信院数学与信息科学系 李胜坤 蒙特卡洛蒙特卡洛(Monte Carlo)方法，或称计算机方法，或称计算机随机模拟方法，是一种基于随机模拟方法，是一种基于“随机数的计随机数的计算方法。这一方法源于美国在第二次世界大算方法。这一方法源于美国在第二次世界大战中研制原子弹的战中研制原子弹的“曼哈顿方案。该方案曼哈顿方案。该方案的主持人之一、数学家冯的主持人之一、数学家冯诺伊曼用著名世诺伊曼用著名世界的赌城界的赌城摩纳哥的摩纳哥的M

2、onte Carlo来命名来命名这种方法，为它蒙上了一层奥秘颜色。这种方法，为它蒙上了一层奥秘颜色。一.M-C方法概述根本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率来决议事件的“概率。19世纪人们用投针实验的方法来决议。高速计算机的出现，使得用数学方法在计算机上大量模拟这样的实验成为能够。从Buffon(蒲丰)投针问题谈起 220, /2 0, sin,:sin.llaXAX随机投针可以理解成针的中心点与最近的平行线的距离X是均匀地分布在区间 上的r.v.，针与平行线的夹角 是均匀地分布在区间 上的r.v.，且X与 相互独立，于是针与平行线相交的充要条件为 即相交

3、 2sin0022(sin)2lllpP Xdxdaa 于是有：2lap试验者时间(年)针长投针次数相交次数的估计值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.14159292数值积分问题niiixfnxfExxf110)(1)(d)(选取(0, 1)中随机数序列x1， x2， x3， xn。那么误差约 ，它并不能和一些高级的数值积分算法比较， 但对多维情况，MC方法却很有吸引力。n1 niiiizyxfnzyxzyxf110

4、1010),(1ddd),(我们可产生一系列随机数可简单取3个随机数构成一个随机点，即,.,654321),.,(),(654321相应地，niibaxfnxxfab1)(1d)(1 niiidcbayxfnyxyxfabcd1),(1dd),()(1普通地，A)in points randomn over of (average) of measure(d)(fAxfAMonte Carlo数值积分的优点与普通的数值积分方法比较，Monte Carlo方法具有以下优点：1. Monte Carlo一般的数值方法很难推广到高维积分的情形，而方法很容易推广到高维情形2/1/22. ()() dO

5、 nO n一般的数值积分方法的收敛阶为 ，而由中心极限定理可以保证 Monte Carlo 方法的收敛阶为 。此收敛阶与维数无关，且在高维时明显优于一般的数值方法。M-C的根本思绪1.针对实践问题建立一个简单且便于实现的概率统计模型，使所求的量或解恰好是该模型某个目的的概率分布或者数字特征。2.对模型中的随机变量建立抽样方法，在计算机上进展模拟测试，抽取足够多的随机数，对有关事件进展统计3.对模拟实验结果加以分析，给出所求解的估计及其精度(方差)的估计4.必要时，还应改良模型以降低估计方差和减少实验费用，提高模拟计算的效率回想几种延续型分布1.均匀分布U(a,b)其概率密度函数为其他,0,1)

6、(bxaabxf有),(), c (,badabdcdxcP其中 均匀性特点:均匀分布随机变量X落在(a,b)内恣意子区间的概率只与子区间的长度有关,而与子区间的位置无关. 可假设有这种特性的随机变量服从均匀分布.均匀分布随机变量X的取值具有均匀性2. 正态分布正态分布随机变量X的概率密度函数是Rxxxf,)(21exp21)(2正态分布由两个参数 和 独一确定.其中 是X的均值(数学期望): =E(X),它确定了概率曲线的中心位置,而 是X的规范差: ,它确定了概率曲线的宽窄程度.)(XD 在许多实践问题中,有一类随机变量可以表示成为许多相互独立的随机变量之和,而其中每个随机变量对总和只起微

7、小的影响,这类随机变量往往服从或近似服从正态分布.在实践运用中,假设我们分析到一个随机变量遭到较多独立的微小要素的叠加影响,就可以用正态分布来模拟这个变量.如:工厂产品的丈量尺寸,农作物的收获量,某地域成年人的身高,体重等可看成服从正态分布的随机变量.3. 指数分布指数分布随机变量X的概率密度为0, 00,)(xxexfx指数分布常用来描画寿命问题.二.随机数的生成1.蒙特卡罗模拟的关键是生成优良的随机数。2.在计算机实现中,我们是经过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模拟了随机数的性质(如可以经过 统计检验)。我们通常称之为伪随机数(pseudo-r

8、andom numbers)。3.在模拟中，我们需求产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布U(0,1)的随机数。U(0,1)随机数的生成乘同余法：1101 mod , , /iiiiixaxmuaxxmxm其中 均为整数, 可以任意选取。 称为种子,a 是乘因子,m是模数0 x一个简单的例子1i+116 mod11, u/11 6,11iiixxxam（）0 1 ,x 1,6,3,7,9,10,5,8,4当时 得到序列:,1,6,3.,2.003,1 ,1,3,9.3,2,2,1,3,9,5,42,6,7,10,8,6.axax如果令 得到序列:如果令 得到序列：

9、一个简单的例子(续)上面的例子中，第一个随机数生成器的周期长度是 10，而后两个生成器的周期长度只需它的一半。我们自然希望生成器的周期越长越好，这样我们得到的分布就更接近于真实的均匀分布。0 (max在给定 的情况下，生成器的周期与和初值种子)选择有关。线性同余生成器(混合同余法)(Linear Congruential Generator )111 () mod /iiiixaxcmuxm一般形式：02. 1 mmaxm线性同余器可以达到的最长周期为 ，我们可以通过适当的选择 和，使无论选取怎样的初值 都可以达到最大周期(一般选取 为质数) c是非负整数.经过适中选取参数c可以改善 随机数的

10、统计性质(独立性,均匀性).常用的线性同余生成器Modulus mMultiplier aReference231-1=214748364716807Lewis, Goodman, and Miller39373LEcuyer742938285Fishman and Moore950706376Fishman and Moore1226874159Fishman and Moore214748339940692LEcuyer214748356340014LEcuyer复杂一些的生成器Multiple recursive generator1122102(.) mo,.)d/iiiki kikk

11、ixa xa xa xmuxxxxm需(要选取种子12-1-1 (,.) 1ijiii kkkxmxxxmm每个有种选择，所以向量 可以取个不同的值，所以这样的随机数生成器的最大周期可以达到 ，大大提高了简单同余生成器的周期。算法实现许多程序文语中都自带生成随机数的方法，如 c 中的 random() 函数，Matlab中的rand()函数等。但这些生成器生成的随机数效果很不一样，比如 c 中的函数生成的随机数性质就比较差，假设用 c ，最好本人再编一个程序。Matlab 中的 rand() 函数，经过了很多优化。可以产生性质很好的随机数，可以直接利用。从U(0,1)到其它概率分布的随机数1.

12、离散型随机数的模拟2.延续型随机数的模拟3.正态随机数的模拟1.离散型随机数的模拟设随机变量 的分布律为令将 作为区间(0,1)的分点.假设随机变量 ,有X), 2 , 1(ipxXPii),2, 1(,)(, 0)0(1npnPPnii) 1 , 0( UU),2, 1( ,)1()()()1(npnPnPnPUnPPn)(nP令那么有据此,可得产生 的随机数的详细过程为:每产生一个(0,1)区间上均匀分布随机数 ,假设 那么令 取值 .nx)() 1(nxXnPUnPXnnpxXPX)() 1(nPUnPU例1:离散型随机变量X有如下分布律: X 0 1 2P(x) 0.3 0.3 0.4

13、设 是(0,1)上均匀分布的随机数,令那么 是具有X分布律的随机数.iiiiUUUx6 . 0, 26 . 03 . 0, 13 . 00, 0NUUU,21Nxxx,212.延续型随机数的模拟a.逆变换方法(常用) (Inverse Transform Method)b.舍取方法 (Acceptance-Rejection Method)定理: 设随机变量Y的分布函数F(y)是延续函数, 而U是在(0,1)上均匀分布的随机变量, 令, 那么Y与X有一样的分布.)(1UFX1 1 ( ) (0,1) ( ) F xUuFu由 定 理， 要 产 生 来 自的 随 机 数 ， 只 要 先产 生 来

14、 自随 机 数， 然 后 计 算 即可 。 具 体 步 骤 如 下 ：(1) (0,1)U上生 成均 匀 分 布 的 随 机 数。-1(2) , ( ) ( ) XXUFFx计算则为来自 分布的随机数.11() ()()( )( )FUP XxP FUxP UF xF x 由 的 定 义 和 均 匀 分 布 的 分 布 函明数 可证：得 ：-11 ( , ) 0 ( ) 1, ( )() 01 (0,1) ( - ) ( , ) XU a bxaxaF xaxbbaxbFyaba yyUUab a UU a b例 ： 设，则其分布函数为 ，生成随机数，则是来自的随机数。/12 exp( ) (

15、 )1, 0 ( )log(1) log(1) ( ) 1 xXXF xexFyyXUUUU 例 ：设服从指数分布，则的分布函数为：通过计算得，则：服从指数分布 其中服从均匀分布又因为 和有着同样的分布，所以也可以取： log( )XU 舍取方法舍取方法(Acceptance-Rejection )方法最早由 Von Neumann提出，如今曾经广泛运用于各种随机数的生成。根本思绪：经过一个容易生成的概率分布 g 和一个取舍准那么生成另一个与 g 相近的概率分布 f 。详细步骤： ( ) ( ) ( )( ) 1, f xg xf xcg xcx 假设和均为集合 上的概率密度函数。且满足： 。

16、易于从g(x)中生成样本X.用如下步骤生成有Y:1. ( ) ;2. (0,1),U3. ()/() g xXUXUf Xcg X生成 的样本生成U且 与 独立；如果，则取Y=X,返回步骤1， 否则舍去X,返回步骤1.下面我们验证由上述步骤生成的随机数 Y 确实具有密度函数 f(x) ()(|()/()( , ()/() ()/()( )( )1/( )( )()( )( )( )AAAP YAP XA Uf Xcg XP XAUf Xcg XP Uf Xcg Xf xg x dxccg xf xP YAcg x dxf x dxcg xx对于任何的 ，我们有：）分母即Y的概率密度函数为f(

17、).()/()1/ ( ) ( )(0,1) P Uf Xcg Xcf xcg xU由上面的证明过程我们可以看出每次接受的概率为 ）。也就是说为了生成一个 的样本，我们需要平均生成个和 的样本。所以为了提高舍取法的效率，我们应该使 c 的取值尽能够的小，也就是使 f 和 g 的分布更为相近。3.规范正态分布随机数的生成正态分布是概率统计中最重要的分布，在此我们着重讨论如何生成规范正态分布随机数。引理：121/21121/221212, (0,1) ( 2ln)cos(2) ( 2ln)sin(2) U UUZUUZUUZZ 设 是独立同分布的变量，令：则证明：与 独由随机立，且均服从向量的变换

18、公标准式易正态分布。得，略。Box-Muller 算法121. ,U U生成两个独立的均匀分布随机数 1/21121/221212 ( 2ln)cos(2) ( 2ln)sin(2) ZUUZUUZZ 2.令 则，为独立的正态随机数利用中心极限定理设 是n个相互独立的在(0,1)上均匀分布的随机变量,由中心极限定理知 渐近服从正态分布N(0,1).普通取n=10即可,假设取n=12,那么上式简化为 再由公式 即可得到正态分布 的随机数.nUUU,2112/ )2(1nnUXniiXY),(2N1216iiUXMatlab程序Function r=rnd-u(a,b)%产生在a,b间均匀分布的随

19、机数r=a+(b-a)*rand;returnMatlab程序Function r=rnd-beta(lamda)%模拟指数分布%lamda表示指数分布的参数r=-log(rand)/lamda;returnMatlab程序Function y= rnd(mean, segema)%模拟均值为mean,方差为segma的正态分布r=rand(1,12);x=sum(r)-6;y=segma*x+mean;return三. 模拟训练例1. 模拟求近似圆周率在边长为1的正方形内有一半径为0.5的内切圆.如今模拟产生在正方形内均匀分布的点n个.如果有m个在圆内,那么圆面积与正方形的面积比可近似为m/

20、n.即/4m/n 4m/nn=10000m=0;For i=1:n if rand2+rand2 =0,那么可以经过模拟估算.构造一个矩形包含曲边梯形,dmax f(x).产生n(足够大)个在矩形区域内的点,假设落在由函数f(x)构成的曲边梯形内的点为m个,那么所求定积分为 .dabnmxfba)()(n=106;a=0;b=1;d=1;m=0;for i=1:n x=a+rand*(b-a); y=d*rand; if y=x2 m=m+1; endends=m/n*(b-a)*dn=106;x=rand(1,n);y=x.2;s=sum(y)/n采用前面讲的方法:niiixfnxfExxf

21、110)(1)(d)(例3. 渡口模型问题描画:一个渡口的渡船营运者拥有一只甲板长32米,可以并排停放两列车辆的渡船.他在思索怎样在甲板上安排过河车辆的位置,才干平安地运过最多数量的车辆.分析:怎样安排过河车辆,关怀一次可以运多少辆各类车.预备任务:察看数日,发现每次情况不尽一样,得到以下数据和情况:(1)车辆随机到达,构成一个等待上船的车列;(2)来到车辆,轿车约占40%,卡车约占55%,摩托车约占5%;(3)轿车车身长3.55.5米,卡车车身长为810米.问题分析:这是一个机理较复杂的随机问题,是遵先到先效力的随机排队问题.处理方法:采用模拟模型方法.因此需思索以下问题: (1)应该怎样安

22、排摩托车? (2)下一辆到达的车是什么类型? (3)怎样描画一辆车的车身长度? (4)如何安排到达车辆参与甲板上两列车队中的哪一列中去?本实验主要模拟装载车辆的情况,暂时不思索渡船的平安.模型建立:设到达的卡车,轿车长度分别为随机变量 , .结合实践,这里无妨设卡车,轿车的车身长度 , 均服从正态分布.由于卡车车身长810米,所以卡车车长 的均值为(8+10)/2=9米,由概率知识中的3原那么,其规范差为(9-8)/3=1/3,所以得到 .同理可得 .1L2L1L2L1L)91, 9(1NL)91, 5 . 4(2NL模拟程序设计:由以上的分析,程序设计时应划分的主要模块(函数)如下:(1)确

23、定下一辆到达车辆的类型:(2)根据车的类型确定到达车辆的长度;(3)根据一定的停放规那么,确定放在哪一列.function r =makeid%模拟下一辆到达车的类型t=rand;if t=0.55 r=1;%到达卡车elseif tL(2)v if L(2)+newlen32v pos=2;v elsev full=1;v endvelsev if L(1)+newlen32v pos=1;v elsev full=1;v endvendvfunction sim_dukou v%渡口模型的模拟vn=input(输入模拟次数:);vif isempty(n)|(n500)v n=500;ve

24、ndvN=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/n四.实验标题赶上火车的概率问题描画:如图,一列火车从A站开往B站,某人每天赶往B站上这趟火车.他已了解到:(1)火车从A站到B站的运转时间是均值为30分钟,规范差为2分钟的随机变量;(2)火车在下午大约1点分开A站