资源描述
Monte Carlo Simulation Methods(蒙特卡罗模拟方法蒙特卡罗模拟方法)一.概述与思想 二.随机数的生成.三.实例-港口模型 四.作业引例：引例：电梯系统电梯系统 然而然而这种做法可能是难以接受的这种做法可能是难以接受的，因为在收集统计，因为在收集统计数据时要再三惊扰乘客，并且电梯运行模式的不断变数据时要再三惊扰乘客，并且电梯运行模式的不断变化也会使乘客感到迷惑化也会使乘客感到迷惑 我们可以我们可以提出若干供选择的电梯运行模式提出若干供选择的电梯运行模式，如设定，如设定停偶数层、奇数层的电梯或直达电梯停偶数层、奇数层的电梯或直达电梯理论上，对每种理论上，对每种供选择的模式都能够做若干次试验供选择的模式都能够做若干次试验，以确定哪一种模式，以确定哪一种模式能为那些要到达特定楼层的乘客提供最好的服务。能为那些要到达特定楼层的乘客提供最好的服务。一 概述与思想 所以在某些情况下，所以在某些情况下，对对象的行为进行直接观测对对象的行为进行直接观测或重复试验可能是不可行的或重复试验可能是不可行的。与此有关的与此有关的另一个问题是大城市交通控制系统可另一个问题是大城市交通控制系统可供选择的运行模式的检验供选择的运行模式的检验，为了做试验而不停地改变，为了做试验而不停地改变单行道的交通方向和配置交通信号将是不现实的单行道的交通方向和配置交通信号将是不现实的 在对象的行为不能做分析性的解释，或数据无法直在对象的行为不能做分析性的解释，或数据无法直接收集的情况下，建模者接收集的情况下，建模者可以用某种方式间接地模拟可以用某种方式间接地模拟其行为其行为，试验所研究的供选择的各种方案，以估计它，试验所研究的供选择的各种方案，以估计它们怎样影响对象的行为，然后收集数据来确定哪种方们怎样影响对象的行为，然后收集数据来确定哪种方案是最好的案是最好的 例如，为了得到一艘拟建造的潜艇受到的阻力，例如，为了得到一艘拟建造的潜艇受到的阻力，造一个原型是不可行的，我们可以按比例建一个模型，造一个原型是不可行的，我们可以按比例建一个模型，去模拟实际的潜艇的行为去模拟实际的潜艇的行为 这里将研究另外一种形式的模拟这里将研究另外一种形式的模拟蒙特卡罗蒙特卡罗(MonteCarlo)模拟，模拟，一般是借助于计算机完成的一般是借助于计算机完成的 蒙特卡洛蒙特卡洛(Monte Carlo)方法，或称计算机随机方法，或称计算机随机模拟方法，是一种基于模拟方法，是一种基于“随机数随机数”的计算方法。的计算方法。这一方法源于美国在第二次世界大战中研制原这一方法源于美国在第二次世界大战中研制原子弹的子弹的“曼哈顿计划曼哈顿计划”。该计划的主持人之一、数。该计划的主持人之一、数学家冯学家冯诺伊曼用驰名世界的赌城诺伊曼用驰名世界的赌城摩纳哥的摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层来命名这种方法，为它蒙上了一层神秘色彩。神秘色彩。从Buffon(蒲丰)投针问题谈起 试验者时间(年)针长投针次数相交次数的估计值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592921 1 确定行为的模拟确定行为的模拟例：曲线下的面积例：曲线下的面积 本节以曲线下的面本节以曲线下的面积为例说明蒙特卡罗积为例说明蒙特卡罗模拟在确定行为建模模拟在确定行为建模中的应用中的应用 下面的算法给出了用蒙特卡罗方法求曲线下面积下面的算法给出了用蒙特卡罗方法求曲线下面积的计算机模拟的计算格式的计算机模拟的计算格式 在给定区间上曲线在给定区间上曲线y=cosy=cosx下面积的真值是下面积的真值是2 2注意到即使对于注意到即使对于产生的相当多的点数，产生的相当多的点数，误差也是可观的误差也是可观的对单变量函数，一般说对单变量函数，一般说来，来，蒙特卡罗方法无法与在数值分析中学到的积分方法相比蒙特卡罗方法无法与在数值分析中学到的积分方法相比，没没有误差界以及难以求出函数的上界有误差界以及难以求出函数的上界M M也是它的缺点也是它的缺点然而，蒙特卡然而，蒙特卡罗方法可以推广到多变量函数，在那里它变得更加实用罗方法可以推广到多变量函数，在那里它变得更加实用做任何的蒙特卡罗模拟，都要用到做任何的蒙特卡罗模拟，都要用到随机数随机数2.2.同样抛同样抛100100次的结果也不近相同。次的结果也不近相同。蒙特卡罗模拟是一种随机模型蒙特卡罗模拟是一种随机模型1.1.抛抛100100次硬币得到次硬币得到5151个正面，并且接下来的个正面，并且接下来的1010次次(即使不太可能刚巧即使不太可能刚巧1010次次)全为正面的情况是可全为正面的情况是可能出现的，这样，用能出现的，这样，用110110次的结果进行估计实际次的结果进行估计实际上比用上比用100100次要差次要差 要记住，要记住，对于根据模拟结果的预测寄予太多的对于根据模拟结果的预测寄予太多的信任是有危险的，信任是有危险的，特别是在模拟中包含的假设没有特别是在模拟中包含的假设没有清楚表明的时候清楚表明的时候还有，由于用了大量的数据和庞还有，由于用了大量的数据和庞大的计算，再加上非专业人员理解模拟模型和计算大的计算，再加上非专业人员理解模拟模型和计算机输出相对容易，所以常会导致对模拟结果的过分机输出相对容易，所以常会导致对模拟结果的过分相信相信二.随机数的生成1.蒙特卡罗模拟的关键是生成优良的随机数。关键是生成优良的随机数。2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。我们通常称之为伪随机数(pseudo-random numbers)。3.在模拟中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布均匀分布U(0,1)U(0,1)的随机数。的随机数。U(0,1)随机数的生成乘同余法：称为种子,a 是乘因子,m是模数一个简单的例子一个简单的例子(续)上面的例子中，第一个随机数生成器的周期长度是 10，而后两个生成器的周期长度只有它的一半。我们自然希望生成器的周期越长越好，这样我们得到的分布就更接近于真实的均匀分布。线性同余生成器(混合同余法)(Linear Congruential Generator)1.c是非负整数.通过适当选取参数c可以改善随机数的统计性质(独立性,均匀性).常用的线性同余生成器Modulus mMultiplier aReference231-1=214748364716807Lewis,Goodman,and Miller39373LEcuyer742938285Fishman and Moore950706376Fishman and Moore1226874159Fishman and Moore214748339940692LEcuyer214748356340014LEcuyer复杂一些的生成器Multiple recursive generator算法实现许多程序语言中都自带生成随机数的方法，如 c 中的 random()函数，Matlab中的rand()函数等。但这些生成器生成的随机数效果很不一样，比如 c 中的函数生成的随机数性质就比较差，如果用 c，最好自己再编一个程序。Matlab 中的 rand()函数，经过了很多优化。可以产生性质很好的随机数，可以直接利用。从U(0,1)到其它概率分布的随机数1.离散型随机数的模拟2.连续型随机数的模拟3.正态随机数的模拟1.离散型随机数的模拟设随机变量 的分布律为令将 作为区间(0,1)的分点.若随机变量 ,有令则有据此,可得产生 的随机数的具体过程为:每产生一个(0,1)区间上均匀分布随机数 ,若 则令 取值 .例1:离散型随机变量X有如下分布律:X 0 1 2P(x)0.3 0.3 0.4设 是(0,1)上均匀分布的随机数,令则 是具有X分布律的随机数.2.连续型随机数的模拟a.逆变换方法(常用)(Inverse Transform Method)b.舍取方法 (Acceptance-Rejection Method)定理:设随机变量Y的分布函数F(y)是连续函数,而U是在(0,1)上均匀分布的随机变量,令,则Y与X有相同的分布.舍取方法舍取方法(Acceptance-Rejection)方法最早由 Von Neumann提出，现在已经广泛应用于各种随机数的生成。基本思路：通过一个容易生成的概率分布 g 和一个取舍准则生成另一个与 g 相近的概率分布 f。具体步骤：下面我们验证由上述步骤生成的随机数 Y 确实具有密度函数 f(x)所以为了提高舍取法的效率，我们应该使 c 的取值尽可能的小，也就是使 f 和 g 的分布更为相近。3.标准正态分布随机数的生成正态分布是概率统计中最重要的分布，在此我们着重讨论如何生成标准正态分布随机数。引理：Box-Muller 算法利用中心极限定理设 是n个相互独立的在(0,1)上均匀分布的随机变量,由中心极限定理知 渐近服从正态分布N(0,1).一般取n=10即可,若取n=12,则上式简化为 再由公式 即可得到正态分布 的随机数.Matlab程序Function r=rnd-u(a,b)%产生在a,b间均匀分布的随机数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实例一 例例（港口系统）（港口系统）考察一个带有船只卸货设备的小港考察一个带有船只卸货设备的小港口，任何时间仅能为一艘船只卸货船只进港是为了口，任何时间仅能为一艘船只卸货船只进港是为了卸货，卸货，相邻两艘船到达的时间间隔在相邻两艘船到达的时间间隔在1515分钟到分钟到145145分钟分钟之间变化之间变化一艘船只卸货的时间由所卸货物的类型决一艘船只卸货的时间由所卸货物的类型决定，在定，在4545分钟到分钟到9090分钟之间变化分钟之间变化 需要回答以下问题：需要回答以下问题：1 1每艘船只在港口的平均时间和最长时间是多少每艘船只在港口的平均时间和最长时间是多少?2 2若一艘船只的等待时间是从到达到开始卸货的若一艘船只的等待时间是从到达到开始卸货的时间，每艘船只的平均等待时间和最长等待时间是多时间，每艘船只的平均等待时间和最长等待时间是多少少?3 3卸货设备空闲时间的百分比是多少卸货设备空闲时间的百分比是多少?4 4船只排队最长的长度是多少船只排队最长的长度是多少?3 3 实例二实例二 为了得到一些合理的答案，为了得到一些合理的答案，利用计算机或可编程利用计算机或可编程计算器来模拟港口的活动计算器来模拟港口的活动假定相邻两艘船到达的假定相邻两艘船到达的时间间隔和每艘船只卸货的时问在它们各自的时间时间间隔和每艘船只卸货的时问在它们各自的时间区问内均匀分布区问内均匀分布，例如两艘船到达的时间间隔可以，例如两艘船到达的时间间隔可以是是1515到到145145之间的任何整数，且这个区间内的任何整之间的任何整数，且这个区间内的任何整数等可能地出现数等可能地出现 在给出模拟这个港口系统的一般算法之前，在给出模拟这个港口系统的一般算法之前，考虑考虑有有5 5艘船只的假想情况对每艘船只有以下数据：艘船只的假想情况对每艘船只有以下数据：设想码头设备的拥有者关心他们设想码头设备的拥有者关心他们提供服务的质量提供服务的质量，并且要并且要评价各种管理模式以确定为了改善服务是否值评价各种管理模式以确定为了改善服务是否值得增加费用得增加费用做一些统计可以帮助对服务质量的评价做一些统计可以帮助对服务质量的评价。例如，呆在港口时间最长的船只是船例如，呆在港口时间最长的船只是船5 5，呆了，呆了130130分分钟，而平均是钟，而平均是8989分钟分钟(表表5-14)5-14)通常顾客对等待时间通常顾客对等待时间的长短非常在乎，例中最长的等待时间为的长短非常在乎，例中最长的等待时间为5555分钟，而分钟，而平均是平均是2626分钟如果排队太长有些顾客会改到别处去分钟如果排队太长有些顾客会改到别处去做生意，例中最长的队是做生意，例中最长的队是2 2 用下面的蒙特卡罗模拟算法可以做这些统计，对各用下面的蒙特卡罗模拟算法可以做这些统计，对各种管理模式进行估价种管理模式进行估价 现在假定你是码头设备拥有者的顾问，如果能现在假定你是码头设备拥有者的顾问，如果能够雇用更多的劳动力，或者得到更好的卸货设备，够雇用更多的劳动力，或者得到更好的卸货设备，使卸货时间减少到每艘船使卸货时间减少到每艘船35357575分钟分钟，会有什么影，会有什么影响响?表表5-165-16给出了基于模拟算法的结果。给出了基于模拟算法的结果。从表从表5-165-16可以看到，每艘船的卸货时间减少可以看到，每艘船的卸货时间减少15152020分钟，使得分钟，使得船只船只呆在港口的时间特别是等待时间缩短呆在港口的时间特别是等待时间缩短了，而了，而设备空闲时问设备空闲时问的百分比却增加了近一倍的百分比却增加了近一倍船主对此是满意的，因为这提高了船主对此是满意的，因为这提高了长期行驶时每艘船运送货物的效率长期行驶时每艘船运送货物的效率 这样，人港贸易好像会增加如果这样，人港贸易好像会增加如果贸易量增加使贸易量增加使得相邻两艘船到达的时间间隔缩短到得相邻两艘船到达的时间间隔缩短到1010到到120120分钟之分钟之间间，模拟结果如表，模拟结果如表5-175-17从这个表可以看到，随着从这个表可以看到，随着贸易量的增加，船只又要在港口呆更长的时间，但贸易量的增加，船只又要在港口呆更长的时间，但设备空闲时间少多了，于是船主和设备拥有者都随设备空闲时间少多了，于是船主和设备拥有者都随着贸易量的增加而受益着贸易量的增加而受益 假定我们假定我们对两艘船到达的时间间隔和每艘船的卸对两艘船到达的时间间隔和每艘船的卸货时间分别在货时间分别在15between15betweeni145145和和45unload45unloadi9090内均匀分布不满意内均匀分布不满意，决定收集港口系统的经验数据，决定收集港口系统的经验数据，并将结果并人我们的模型设想观测了利用港口卸并将结果并人我们的模型设想观测了利用港口卸货的货的12001200艘船，收集的数据见表艘船，收集的数据见表5-185-18 最后，按照表最后，按照表5-195-19和表和表5-205-20给出的规则生成给出的规则生成betweenbetweeni和和unloadunloadi。(i=1,2,3,n)，将线性样条子模型并入，将线性样条子模型并入港口系统的模拟模型用这些子模型，表港口系统的模拟模型用这些子模型，表5-215-21给出了每给出了每次次100100艘船共艘船共6 6次独立模拟的结果次独立模拟的结果.作业一例 模拟求近似圆周率在边长为1的正方形内有一半径为0.5的内切圆.现在模拟产生在正方形内均匀分布的点n个.如果有m个在圆内,则圆面积与正方形的面积比可近似为m/n.即/4m/n 4m/nn=10000m=0;For i=1:n if rand2+rand2=0,则可以通过模拟估算.构造一个矩形包含曲边梯形,dmax f(x).产生n(足够大)个在矩形区域内的点,如果落在由函数f(x)构成的曲边梯形内的点为m个,则所求定积分为 .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采用前面讲的方法:作业 三 赶上火车的概率赶上火车的概率问题描述:如图,一列火车从A站开往B站,某人每天赶往B站上这趟火车.他已了解到:(1)火车从A站到B站的运行时间是均值为30分钟,标准差为2分钟的随机变量;(2)火车在下午大约1点离开A站,离开时刻的频率分布如下:出发时刻午后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问他能赶上火车的概率是多少?变量说明 :火车从A站出发的时刻;:火车从A站到B站的运行时间,单位:分钟;:他到达B站的时刻问题分析与假设:此题包含多个随机因数.这里假设 ,都是随机变量,其中 服从正态分布.模型建立很显然,他能及时赶上火车的条件是:+.为了简化计算,将下午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.1如果r为在(0,1)均匀分布的随机数,为了模拟随机变量 ,可以通过如下方法:其中,和 分别用来模拟随机变量 和 .function pro_huoche%模拟赶上火车的概率n=input(输入模拟次数:)m=0;for i=1:n r1=rand;if r1=0.7 t1=0;elseif r1=0.9 t1=5;else t1=10;end r3=rand;if r3=0.3 t3=28;elseif r3=0.7 t3=30;elseif r3=0.9 t3=32;else t3=34;end t2=30+randn*2;if t3(t1+t2)m=m+1;endendp=m/nend
展开阅读全文