MonteCarlo算法模拟.ppt

Monte Carlo 模拟,内容提纲,1.引言 2.Monte Carlo模拟基本思想 3.随机数生成函数 4.应用实例举例 5.排队论模拟 6.Monte Carlo模拟求解规划问题,引言(Introduction),Monte Carlo方法：,蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,Monte Carlo是摩纳哥（monaco)的首都，该城以赌博闻名,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,Monte Carlo方法的基本思想,蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。源于美国在第二次世界大战研制原子弹的“曼哈顿计划”，该计划的主持人之一数学家冯·诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。 蒙特卡罗方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用蒲丰投针的方法来计算圆周率，上世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,蒲丰投针实验： 法国科学家蒲丰(Buffon)在1777年提出的蒲丰 投针实验是早期几何概率一个非常著名的例子。蒲丰 投针实验的重要性并非是为了求得比其它方法更精确 的值，而是它开创了使用随机数处理确定性数学问 题的先河，是用偶然性方法去解决确定性计算的前导。由此可以领略到从“概率土壤”上开出的一朵瑰丽的鲜花-蒙特卡罗方法(MC) 蒲丰投针实验可归结为下面的数学问题：平面上画有距离为a的一些平行线，向平面上任意投一根长为l (la)的针，假设针落在任意位置的可能性相同，试求针与平行线相交的概率P(从而求),蒲丰投针实验： 如右图所示，以M表示针落下 后的中点，以x表示M到最近一条平行 线的距离，以表示针与此线的交角： 针落地的所有可能结果满足： 其样本空间视作矩形区域, 面积是: 针与平行线相交的条件： 它是样本空间子集A，面积是： syms l phi; int('l/2*sin(phi)',phi,0,pi) %ans=l 因此，针与平行线相交的概率为： 从而有：,蒲丰投针实验的计算机模拟： format long; %设置15位显示精度 a=1; l=0.6; %两平行线间的宽度和针长 figure; axis(0,pi,0,a/2); %初始化绘图板 set(gca,'nextplot','add'); %初始化绘图方式为叠加 counter=0; n=2010; %初始化计数器和设定投针次数 x=unifrnd(0,a/2,1,n); phi=unifrnd(0,pi,1,n); %样本空间 for i=1:n if x(i)l*sin(phi(i)/2 %满足此条件表示针与线的相交 plot(phi(i),x(i),r.); counter=counter+1; %统计针与线相交的次数 frame(counter)=getframe; %描点并取帧 end end fren=counter/n; pihat=2*l/(a*fren) %用频率近似计算 figure(2) movie(frame,1) %播放帧动画1次,一些人进行了实验，其结果列于下表 ：,蒙特卡罗投点法是蒲丰投针实验的推广： 在一个边长为a的正方形内随机投点， 该点落在此正方形的内切圆中的概率 应为该内切圆与正方形的面积比值， 即 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),基本思想,由上面的例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与之有关的量时，通过某种试验的方法，得出该事件发生的频率，再通过它得到问题的解。这就是蒙特卡罗方法的基本思想。,蒙特卡罗方法的关键步骤在于随机数的产生，计算机产生的随机数都不是真正的随机数(由算法确定的缘故)，如果伪随机数能够通过一系列统计检验，我们也可以将其当作真正的随机数使用。,rand('seed',0.1); rand(1) 每次运行程序产生的值是相同的,rand('state',sum(100*clock)*rand); rand(1) %每次重新启动matlab时，输出的随机数不一样,注意: 产生一个参数为的指数分布的随机数应输入 exprnd(1/),产生m×n阶参数为A1,A2,A3的指定分布'name'的随机数矩阵 random('name',A1,A2,A3,m,n) 举例: 产生2×4阶的均值为0方差为1的正态分布的随机数矩阵 random('Normal',0,1,2,4) 'name'的取值可以是(详情参见help random)： 'norm' or 'Normal' / 'unif' or 'Uniform' 'poiss' or 'Poisson' / 'beta' or 'Beta' 'exp' or 'Exponential' / 'gam' or 'Gamma' 'geo' or 'Geometric' / 'unid' or 'Discrete Uniform' ,非常见分布的随机数的产生,逆变换方法,Acceptance-Rejection 方法 最早由 Von Neumann提出，现在已经广泛应用于各种随机数的生成。,基本思路：,通过一个容易生成的概率分布 g 和一个取舍 准则生成另一个与 g 相近的概率分布 f 。,为要模拟服从给定分布的随机变量,用生成一个易于实现的不可约遍历链 作为随机样本，使其平稳分布为给定分布的方法，称为马氏链蒙特卡罗方法.,马氏链蒙特卡罗方法,1 连续型随机变量(以指数分布为例)：,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 指数分布的随机数 %指数分布随机数产生函数已经提供 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(rFx); x(index(1)-1) %通用离散分布随机数产生程序,例1 在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮,现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。,分析：这是一个概率问题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程.,为了能显示我方20次射击的过程，现采用模拟的方式。,举例,需要模拟出以下两件事：,2 当指示正确时，我方火力单位的射击结果情况,1 观察所对目标的指示正确与否,模拟试验有两种结果，每一种结果出现的概率都是1/2,因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确,模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6),这时可用投掷骰子的方法来确定： 如果出现的是、三个点：则认为没能击中敌人； 如果出现的是、点：则认为毁伤敌人一门火炮； 若出现的是点：则认为毁伤敌人两门火炮,问题分析,i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数； k3：击中敌人两门火炮的射击总数； E：有效射击比率； E1：20次射击平均每次毁伤敌人的火炮数,符号说明,模拟框图,模拟结果,理论计算,结果比较,虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程,用蒙特卡洛方法进行计算机模拟的步骤：,1 设计一个逻辑框图，即模拟模型 2 根据流程图编写程序，模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象 3 分析模拟结果，计算所需要结果.,随机投掷均匀硬币，验证国徽朝上与朝下的 概率是否都是 1/2,n=10000; % 给定试验次数 m=0; for i=1:n x=randperm(2)-1;%randperm(2) returns a 随机排列 of the integers 1:2. y=x(1); if y=0 % 0 表示国徽朝上，1 表示国徽朝下 m=m+1; end end fprintf('国徽朝上的频率为：%fn',m/n);,小实例一：投掷硬币,随机投掷骰子，验证各点出现的概率是否为 1/6,n=10000; m1=0; m2=0; m3=0; m4=0; m5=0;m6=0; for i=1:n x=randperm(6); y=x(1); switch y case 1, m1=m1+1; case 2, m2=m2+1; case 3, m3=m3+1; case 4, m4=m4+1; case 5, m5=m5+1; otherwise, m6=m6+1; end end . % 输出结果,小实例二：投掷骰子,设某班有 m 个学生，则该班至少有两人同一天生日的概率是多少？,小实例五：生日问题,n=1000; p=0; m=50; % 设该班的人数为 50 for t=1:n a=; q=0; for k=1:m b=randperm(365); a=a,b(1); end c=unique(a); % 去掉重复数字从小到大重新排序 if length(a)=length(c) p=p+1; end end fprintf('至少两人同一天生日的概率为：%fn',p/n);,试验五源程序,clear; m = 50; p1= 1:365; p2= 1:365-m, 365*ones(1,m); p = p1./p2; p = 1- prod(p); fprintf('至少两人同一天生日的概率为：%fn',p);,小实例五的理论值计算,排队问题随机模拟,排队论主要研究随机服务系统的工作过程。,在排队系统中，服务对象的到达时间和服务时间都是随机的。排队论通过对随机服务现象的统计研究，找出反映这些随机现象平均特性的规律指标，如排队的长度、等待的时间及服务利用率等。,1 系统的假设： （1） 顾客源是无穷的； （2） 排队的长度没有限制； （3）到达系统的顾客按先后顺序依次进入服务，“先到先服务”。,在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店设： 1顾客到来间隔时间服从参数为0.1的指数分布 对顾客的服务时间服从，上的均匀分布 排队按先到先服务规则，队长无限制,假定一个工作日为8小时，时间以分钟为单位。 1模拟一个工作日内完成服务的个数及顾客平均等待时间t 2模拟100个工作日，求出平均每日完成服务的个数及每日顾客的平 均等待时间。,单服务员的排队模型模拟,w：总等待时间； ci：第i个顾客的到达时刻； bi：第i个顾客开始服务时刻； ei：第i个顾客服务结束时刻； xi：第i-1个顾客与第i个顾客之间到达的间隔时间； yi：对第i个顾客的服务时间。,符号说明,c1,b1,c3,c4,c5,c2,e1,b2,e2,b3,e3,b4,e4,b5,ci=ci-1+ xi ei=bi+yi bi=max(ci,ei-1),t,思路分析,流程框图,用蒙特卡洛法解非线性规划问题,先产生一个随机数作为初始试验点,以后则将上一个试验点的第j个分量随机产生，其它分量不变而产生一新的试验点这样，每产生一个新试验点只需一个新的随机数分量当KMAXK或PMAXP时停止迭代,在Matlab软件包中编程，共需三个文件:randlp.m, mylp.m, lpconst.m.主程序为randlp.m.,% mylp.m function z=mylp(x) %目标函数 z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2); %转化为求最小值问题,% randlp.m function sol,r1,r2=randlp(a,b,n) %随机模拟解非线性规划 debug=1; a=0; %试验点下界 b=10; %试验点上界 n=1000; %试验点个数 r1=unifrnd(a,b,n,1); %a,b均匀分布随机数矩阵 r2=unifrnd(a,b,n,1); sol=r1(1) r2(1); z0=inf; for i=1:n x1=r1(i); x2=r2(i); lpc=lpconst(x1 x2)； if lpc=1 z=mylp(x1 x2)； if zz0 z0=z； sol=x1 x2； end end end,与Monte Carlo方法相似，但理论基础不同的方法“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)近年来也获得迅速发展。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般可比Monte Carlo方法提高数百倍。,