MonteCarlo算法模拟

1、 内容提纲 1.引言 2.Monte Carlo模拟基本思想 3.随机数生成函数 4.应用实例举例 5.排队论模拟 6.Monte Carlo模拟求解规划问题 Monte Carlo方法：方法： 蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四 十年代中期为了适应当时原子能事业的发展而发展起来的。十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法，亦称统计模拟方法， statistical simulation method 利用随机数进行数值模拟的方法利用随机数进行数值模拟的方法 M

2、onte Carlo名字的由来：名字的由来： Monte Carlo是摩纳哥（是摩纳哥（monaco)的首都，该城以赌博闻名的首都，该城以赌博闻名 Nicholas Metropolis (1915-1999) Monte-Carlo, Monaco Monte CarloMonte Carlo方法的基本思想方法的基本思想 蒙特卡罗方法，或称计算机随机模拟方法，是一种基蒙特卡罗方法，或称计算机随机模拟方法，是一种基 于于“随机数随机数”的计算方法。源于美国在第的计算方法。源于美国在第二二次世界大战研次世界大战研 制原子弹的制原子弹的“曼哈顿计划曼哈顿计划”，该计划的主持人之一数学家该计划的主持

3、人之一数学家 冯冯诺伊曼用驰名世界的赌城诺伊曼用驰名世界的赌城摩纳哥的摩纳哥的Monte CarloMonte Carlo来来 命名这种方法，为它蒙上了一层神秘色彩。命名这种方法，为它蒙上了一层神秘色彩。 蒙特卡罗蒙特卡罗方法的基本思想很早以前就被人们所发现和方法的基本思想很早以前就被人们所发现和 利用。早在利用。早在1717世纪，人们就知道用事件发生的世纪，人们就知道用事件发生的“频率频率”来来 决定事件的决定事件的“概率概率”。1 19 9世纪人们用世纪人们用蒲丰蒲丰投针的方法来投针的方法来计计 算算圆周率圆周率，上，上世纪世纪4040年代电子计算机的出现，特别是近年代电子计算机的出现，特

4、别是近 年来高速电子计算机的出现，使得用数学方法在计算机上年来高速电子计算机的出现，使得用数学方法在计算机上 大量、快速地模拟这样的试验成为可能。大量、快速地模拟这样的试验成为可能。 蒲丰投针实验蒲丰投针实验： 法国科学家蒲丰法国科学家蒲丰(Buffon)(Buffon)在在17771777年提出的蒲丰年提出的蒲丰 投针实验是早期几何概率一个非常著名的例子。投针实验是早期几何概率一个非常著名的例子。蒲丰蒲丰 投针实验投针实验的重要性并非是为了求得比其它方法更精确的重要性并非是为了求得比其它方法更精确 的的值，而是它开创了使用随机数处理确定性数学问值，而是它开创了使用随机数处理确定性数学问 题的

5、先河，是用偶然性方法去解决确定性计算的前导。题的先河，是用偶然性方法去解决确定性计算的前导。 由此由此可以可以领略到从领略到从“概率土壤概率土壤”上开出的一朵瑰丽的鲜上开出的一朵瑰丽的鲜 花花- -蒙特卡罗蒙特卡罗方法方法( (MC)MC) 蒲丰投针实验可归结为下面的数学问题：平面上画蒲丰投针实验可归结为下面的数学问题：平面上画 有距离为有距离为a a的一些平行线，向平面上任意投一根长为的一些平行线，向平面上任意投一根长为l l (la)(la)的针，假设针落在任意位置的可能性相同，试求针的针，假设针落在任意位置的可能性相同，试求针 与平行线相交的概率与平行线相交的概率P(P(从而求从而求)

6、) 蒲丰投针实验蒲丰投针实验： 如右图所示，以如右图所示，以M M表示针落下表示针落下 后的中点，以后的中点，以x x表示表示M M到最近一条平行到最近一条平行 线的距离，以线的距离，以表示针与此线的交角：表示针与此线的交角： 针落地的所有可能结果满足：针落地的所有可能结果满足： 其样本空间视作矩形区域其样本空间视作矩形区域, , 面积是面积是: : 针与平行线相交的条件：针与平行线相交的条件： 它是样本空间它是样本空间子集子集A A，面积是：，面积是： syms l phi; int(l/2*sin(phi),phi,0,pi) %ans=l 因此，针与平行线相交的概率为：因此，针与平行线相

7、交的概率为： 从而有：从而有： 0, 2/0ax 2/aS 0, 2/sin0lx 0 2/sinldlAS alSASp/2)(/ )( apl /2 蒲丰投针实验蒲丰投针实验的计算机模拟：的计算机模拟： format long; %format long; %设置设置1515位显示精度位显示精度 a=1;a=1; l=0.6;l=0.6; % %两平行线间的宽度两平行线间的宽度和和针长针长 figure; axis(0,pi,0,a/2); %figure; axis(0,pi,0,a/2); %初始化绘图板初始化绘图板 set(gca,nextplot,set(gca,nextplot,

8、a adddd);); % %初始化绘图方式为叠加初始化绘图方式为叠加 counter=0;counter=0; n=n=20102010; ; % %初始化计数器和设定初始化计数器和设定投针次数投针次数 x=unifrnd(0,a/2,1,n);x=unifrnd(0,a/2,1,n); p phihi=unifrnd(0,pi,1,n);=unifrnd(0,pi,1,n); % %样本空间样本空间 for i=1:nfor i=1:n if x(i)l if x(i)l* *sin(sin(p phihi(i)/2(i)/2 % %满足此条件表示针与线的相交满足此条件表示针与线的相交 p

9、lot(phi(i),x(i),plot(phi(i),x(i),r.r.) ); ; counter=counter+1; counter=counter+1; % %统计针与线相交的次数统计针与线相交的次数 frame(frame(countercounter)=getframe; %)=getframe; %描点并取帧描点并取帧 end end endend fren=counter/n;fren=counter/n; pihat=2pihat=2* *l/(al/(a* *fren)fren) % %用用频率频率近似计算近似计算 figure(2)figure(2) movie(fra

10、me,1) %movie(frame,1) %播放帧动画播放帧动画1 1次次 一些人进行了实验，其结果列于下表一些人进行了实验，其结果列于下表 ： 实验者年份投计次数的实验值 沃尔弗(Wolf)185050003.1596 斯密思(Smith)185532043.1553 福克斯(Fox)189411203.1419 拉查里尼 (Lazzarini) 190134083.1415929 蒙特卡罗投点法是蒙特卡罗投点法是蒲丰投针实验蒲丰投针实验的推广：的推广： 在一个边长为在一个边长为a的正方形内随机投点，的正方形内随机投点， 该点落在此正方形的内切圆中的概率该点落在此正方形的内切圆中的概率 应

11、为该内切圆与正方形的面积比值，应为该内切圆与正方形的面积比值， 即即 n=10000; a=2; m=0; for i=1:n x=rand(1)*a; y=rand(1)*a; if ( (x-a/2)2+(y-a/2)2 = (a/2)2 ) m=m+1; end end disp(投点法近似计算的投点法近似计算的为: ,num2str(4*m/n); x y o (a/2,a/2) /4a:a/2 2 2 基本思想基本思想 由上面的例子可以看出，当所求问题的解由上面的例子可以看出，当所求问题的解 是某个事件的概率，或者是某个随机变量的数是某个事件的概率，或者是某个随机变量的数 学期望，或

12、者是与之有关的量时，通过某种试学期望，或者是与之有关的量时，通过某种试 验的方法，得出该事件发生的频率，再通过它验的方法，得出该事件发生的频率，再通过它 得到问题的解。这就是蒙特卡罗方法的基本思得到问题的解。这就是蒙特卡罗方法的基本思 想。想。 蒙特卡罗方法的关键步骤在于随机数的产蒙特卡罗方法的关键步骤在于随机数的产 生，计算机产生的随机数都不是真正的随机数生，计算机产生的随机数都不是真正的随机数 ( (由算法确定的缘故由算法确定的缘故) )，如果伪随机数能够通过，如果伪随机数能够通过 一系列统计检验，我们也可以将其当作真正的一系列统计检验，我们也可以将其当作真正的 随机数使用。随机数使用。

13、rand(seed,0.1); rand(1) 每次运行程序产生的值是相同的每次运行程序产生的值是相同的 rand(state,sum(100*clock)*rand); rand(1) %每次重新启动每次重新启动matlab时，输出的随机数不一样时，输出的随机数不一样 注意注意: : 产生一个参数为产生一个参数为的指数分布的随机数应输入的指数分布的随机数应输入 exprnd(1/ exprnd(1/) ) e0 ( ) 00 t x f x x 产生产生m mn n阶阶参数为参数为A1,A2,A3A1,A2,A3的指定分布的指定分布namename的的随机数矩阵随机数矩阵 random(na

14、me,A1,A2,A3,m,n) random(name,A1,A2,A3,m,n) 举例举例: : 产生产生2 24 4阶的均值为阶的均值为0 0方差为方差为1 1的正态分布的随机数矩阵的正态分布的随机数矩阵 random(Normal,0,1,2,4)random(Normal,0,1,2,4) namename 的取值可以是的取值可以是( (详情参见详情参见help random)help random)： norm or Normal / unif or Uniformnorm or Normal / unif or Uniform poiss or Poisson / beta or

15、 Betapoiss or Poisson / beta or Beta exp or Exponential / gam or Gammaexp or Exponential / gam or Gamma geo or Geometric / unid or Discrete Uniformgeo or Geometric / unid or Discrete Uniform 非常见分布的随机数的产生非常见分布的随机数的产生 逆变换方法逆变换方法 1 1 ( ) (0,1) ( ) F x UuFu 由定理，要产生来自的随机数，只要先 产生来自随机数，然后计算 即 可。具体步骤如下： (1)

16、 (0,1)U上生成均匀分布的随机数 。 -1 (2) , ( ) ( ) XXUFFx计算则为来自 分布的随机数. Acceptance-Rejection 方法方法 最早由 Von Neumann提出，现在已经广 泛应用于各种随机数的生成。 基本思路： 通过一个容易生成的概率分布 g 和一个取舍 准则生成另一个与 g 相近的概率分布 f 。 为要模拟服从给定分布的随机变量,用生成 一个易于实现的不可约遍历链 作为随机样本，使其平稳分布为给定分布 的方法，称为马氏链蒙特卡罗方法马氏链蒙特卡罗方法. 0,nXX n 马氏链蒙特卡罗方法马氏链蒙特卡罗方法 1 连续型随机变量连续型随机变量(以指数

17、分布为例以指数分布为例)： e01-e0 ( )( )( ) 0000 tx xtx f tf t dt F x tx syms t x lambda; Fx=int(lambda*exp(-lambda*t),t,0,x) %分布函数分布函数 syms r; Fxinv=finverse(Fx,x); %求反函数求反函数 Fxinv=subs(Fxinv,x,r) %替换反函数变量替换反函数变量x为为r Fxinv=inline(Fxinv) x=Fxinv(3,rand) %产生参数产生参数 lambda=3 指数分布的随机数指数分布的随机数 %指数分布随机数产生函数已经提供指数分布随机数

18、产生函数已经提供 exprnd(1/3,1,1) 2 离散型随机变量离散型随机变量(以离散分布为例以离散分布为例)： x=2,4,6,8; px=0.1,0.4,0.3,0.2; %以下为程序片段以下为程序片段 Fx=0; for n=1:length(px),Fx=Fx,sum(px(1:n);end r=rand; index=find(rMAXK或PMAXP时停止迭代 框框 图图 初始化:给定MAXK,MAXP;k=0,p=0,Q:大整数 xj=aj+R(bj-aj) j=1,2,n j=0 j=j+1,p=p+1 PMAXP? YN xj=aj+R(bj-aj) gi(X)0? i=1,2n YN jMAXK? YN 输出X,Q,停止 Y N 例例 max 2121 2 2 2 1 382xxxxxxz s.t 103 21 xx 0 1 x 0 2 x 在Matlab软件包中编程，共需三个文件:randlp.m, mylp.m, lpconst.m.主程序为randlp.m. % mylp.m% mylp.m function z=mylp(x) %目标函数 z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-