Monte Carlo(蒙特卡洛方法)

1、Monte Carlo Simulation Methods(蒙特卡罗模拟方法蒙特卡罗模拟方法) 一. 概述与思想 二. 随机数的生成. 三. 实例-港口模型 四. 作业引例：引例：电梯系统电梯系统 然而然而这种做法可能是难以接受的这种做法可能是难以接受的，因为在收集统计，因为在收集统计数据时要再三惊扰乘客，并且电梯运行模式的不断变数据时要再三惊扰乘客，并且电梯运行模式的不断变化也会使乘客感到迷惑化也会使乘客感到迷惑 我们可以我们可以提出若干供选择的电梯运行模式提出若干供选择的电梯运行模式，如设定，如设定停偶数层、奇数层的电梯或直达电梯停偶数层、奇数层的电梯或直达电梯理论上，对每种理论上，对每

2、种供选择的模式都能够做若干次试验供选择的模式都能够做若干次试验，以确定哪一种模式，以确定哪一种模式能为那些要到达特定楼层的乘客提供最好的服务。能为那些要到达特定楼层的乘客提供最好的服务。一 概述与思想 所以在某些情况下，所以在某些情况下，对对象的行为进行直接观测对对象的行为进行直接观测或重复试验可能是不可行的或重复试验可能是不可行的。 与此有关的与此有关的另一个问题是大城市交通控制系统可另一个问题是大城市交通控制系统可供选择的运行模式的检验供选择的运行模式的检验，为了做试验而不停地改变，为了做试验而不停地改变单行道的交通方向和配置交通信号将是不现实的单行道的交通方向和配置交通信号将是不现实的

3、在对象的行为不能做分析性的解释，或数据无法直在对象的行为不能做分析性的解释，或数据无法直接收集的情况下，建模者接收集的情况下，建模者可以用某种方式间接地模拟可以用某种方式间接地模拟其行为其行为，试验所研究的供选择的各种方案，以估计它，试验所研究的供选择的各种方案，以估计它们怎样影响对象的行为，然后收集数据来确定哪种方们怎样影响对象的行为，然后收集数据来确定哪种方案是最好的案是最好的 例如，为了得到一艘拟建造的潜艇受到的阻力，例如，为了得到一艘拟建造的潜艇受到的阻力，造一个原型是不可行的，我们可以按比例建一个模型，造一个原型是不可行的，我们可以按比例建一个模型，去模拟实际的潜艇的行为去模拟实际的

4、潜艇的行为 这里将研究另外一种形式的模拟这里将研究另外一种形式的模拟蒙特卡罗蒙特卡罗(MonteCarlo)模拟，模拟，一般是借助于计算机完成的一般是借助于计算机完成的 蒙特卡洛蒙特卡洛(Monte Carlo)方法，或称计算机随机方法，或称计算机随机模拟方法，是一种基于模拟方法，是一种基于“随机数随机数”的计算方法。的计算方法。 这一方法源于美国在第二次世界大战中研制原这一方法源于美国在第二次世界大战中研制原子弹的子弹的“曼哈顿计划曼哈顿计划”。该计划的主持人之一、。该计划的主持人之一、数学家冯数学家冯诺伊曼用驰名世界的赌城诺伊曼用驰名世界的赌城摩纳哥的摩纳哥的Monte Carlo来命名这

5、种方法，为它蒙上了一层来命名这种方法，为它蒙上了一层神秘色彩。神秘色彩。从Buffon(蒲丰)投针问题谈起 220, /2 0, sin,:sin.llaXAX随机投针可以理解成针的中心点与最近的平行线的距离X是均匀地分布在区间 上的r.v.，针与平行线的夹角 是均匀地分布在区间 上的r.v.，且X与 相互独立，于是针与平行线相交的充要条件为 即相交 2sin0022(sin)2lllpP Xdxdaa 于是有：2lap试验者时间(年)针长投针次数相交次数的估计值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.75

6、10304893.15951Lazzarini19250.83340818083.141592921 1 确定行为的模拟确定行为的模拟例：曲线下的面积例：曲线下的面积 本节以曲线下的面本节以曲线下的面积为例说明蒙特卡罗积为例说明蒙特卡罗模拟在确定行为建模模拟在确定行为建模中的应用中的应用 下面的算法给出了用蒙特卡罗方法求曲线下面积下面的算法给出了用蒙特卡罗方法求曲线下面积的计算机模拟的计算格式的计算机模拟的计算格式 在给定区间上曲线在给定区间上曲线y=cosy=cosx下面积的真值是下面积的真值是2 2注意到即使对于注意到即使对于产生的相当多的点数，产生的相当多的点数，误差也是可观的误差也是可

7、观的对单变量函数，一般说对单变量函数，一般说来，来，蒙特卡罗方法无法与在数值分析中学到的积分方法相比蒙特卡罗方法无法与在数值分析中学到的积分方法相比，没没有误差界以及难以求出函数的上界有误差界以及难以求出函数的上界M M也是它的缺点也是它的缺点然而，蒙特卡然而，蒙特卡罗方法可以推广到多变量函数，在那里它变得更加实用罗方法可以推广到多变量函数，在那里它变得更加实用Monte Carlo数值积分的优点与一般的数值积分方法比较，Monte Carlo方法具有以下优点：1. Monte Carlo一般的数值方法很难推广到高维积分的情形，而方法很容易推广到高维情形2/1/22. ()() dO nO n

8、一般的数值积分方法的收敛阶为 ，而由中心极限定理可以保证 Monte Carlo 方法的收敛阶为 。此收敛阶与维数无关，且在高维时明显优于一般的数值方法。 蒙特卡罗方法是随机的，为使预测值与真值之差蒙特卡罗方法是随机的，为使预测值与真值之差变小变小需要做大量的试验需要做大量的试验在最终的估计中，在最终的估计中，要讨论要讨论保证预先给定的置信水平所要求的试验次数，保证预先给定的置信水平所要求的试验次数，需要需要统计学的背景知识，统计学的背景知识，然而作为一般准则，然而作为一般准则，结果的精结果的精度提高一倍度提高一倍( (即误差减少一半即误差减少一半) )，试验次数大约需要，试验次数大约需要增至

9、增至4 4倍倍做任何的蒙特卡罗模拟，都要用到做任何的蒙特卡罗模拟，都要用到随机数随机数2.2.同样抛同样抛100100次的结果也不近相同。次的结果也不近相同。蒙特卡罗模拟是一种随机模型蒙特卡罗模拟是一种随机模型1.1.抛抛100100次硬币得到次硬币得到5151个正面，并且接下来的个正面，并且接下来的1010次次( (即使不太可能刚巧即使不太可能刚巧1010次次) )全为正面的情况是可全为正面的情况是可能出现的，这样，用能出现的，这样，用110110次的结果进行估计实际次的结果进行估计实际上比用上比用100100次要差次要差 要记住，要记住，对于根据模拟结果的预测寄予太多的对于根据模拟结果的预

10、测寄予太多的信任是有危险的，信任是有危险的，特别是在模拟中包含的假设没有特别是在模拟中包含的假设没有清楚表明的时候清楚表明的时候还有，由于用了大量的数据和庞还有，由于用了大量的数据和庞大的计算，再加上非专业人员理解模拟模型和计算大的计算，再加上非专业人员理解模拟模型和计算机输出相对容易，所以常会导致对模拟结果的过分机输出相对容易，所以常会导致对模拟结果的过分相信相信二.随机数的生成1.蒙特卡罗模拟的关键是生成优良的随机数。关键是生成优良的随机数。2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。

11、我们通常称之为伪随机数(pseudo-random numbers)。3.在模拟中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布均匀分布U(0,1)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,4

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

13、一般选取 为质数) c是非负整数.通过适当选取参数c可以改善 随机数的统计性质(独立性,均匀性).常用的线性同余生成器Modulus mMultiplier aReference231-1=214748364716807Lewis, Goodman, and Miller39373LEcuyer742938285Fishman and Moore950706376Fishman and Moore1226874159Fishman and Moore214748339940692LEcuyer214748356340014LEcuyer复杂一些的生成器Multiple recursive ge

14、nerator1122102(.) mo,.)d/iiiki kikkixa xa xa xmuxxxxm需(要选取种子12-1-1 (,.) 1ijiii kkkxmxxxmm每个有种选择，所以向量 可以取个不同的值，所以这样的随机数生成器的最大周期可以达到 ，大大提高了简单同余生成器的周期。算法实现许多程序语言中都自带生成随机数的方法，如 c 中的 random() 函数，Matlab中的rand()函数等。但这些生成器生成的随机数效果很不一样，比如 c 中的函数生成的随机数性质就比较差，如果用 c ，最好自己再编一个程序。Matlab 中的 rand() 函数，经过了很多优化。可以产生性

15、质很好的随机数，可以直接利用。从U(0,1)到其它概率分布的随机数1.离散型随机数的模拟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

16、有如下分布律: X 0 1 2P(x) 0.3 0.3 0.4设 是(0,1)上均匀分布的随机数,令则 是具有X分布律的随机数.iiiiUUUx6 . 0, 26 . 03 . 0, 13 . 00, 0NUUU,21Nxxx,212 2 随机行为的模拟随机行为的模拟 本节的目的说明如何本节的目的说明如何对简单的随机行为建模对简单的随机行为建模，以形成直觉和理解考察以形成直觉和理解考察3 3个简单的随机模型：个简单的随机模型： 1 1抛一枚正规的硬币抛一枚正规的硬币 2 2掷一个或一对正规的骰子掷一个或一对正规的骰子 3 3掷一个或一对不正规的骰子掷一个或一对不正规的骰子 抛一枚硬币时得到正面

17、或反面的机会是抛一枚硬币时得到正面或反面的机会是1/2，于是，于是抛很多次时出正面次数的比例接近抛很多次时出正面次数的比例接近0.50.5设设x为为0,10,1内的内的随机数随机数，f(x)定义如下：当定义如下：当00 x0.5, 0.5, f(x)= =正面，正面，否则否则 f(x)=反面，反面， f(x)将结果是正面或反面赋值到将结果是正面或反面赋值到00，11内的一个数，随机赋值时我们可利用这个函数的累内的一个数，随机赋值时我们可利用这个函数的累积性质抛很多次时能够得到如下的出现百分比：积性质抛很多次时能够得到如下的出现百分比：2.1 2.1 一枚正规的硬币一枚正规的硬币2.3 2.3

18、掷一个不正规的骰子掷一个不正规的骰子2.连续型随机数的模拟a.逆变换方法(常用) (Inverse Transform Method)b.舍取方法 (Acceptance-Rejection Method)定理: 设随机变量Y的分布函数F(y)是连续函数, 而U是在(0,1)上均匀分布的随机变量, 令, 则Y与X有相同的分布.)(1UFX1 1 ( ) (0,1) ( ) F xUuFu由 定 理， 要 产 生 来 自的 随 机 数 ， 只 要 先产 生 来 自随 机 数， 然 后 计 算 即可 。 具 体 步 骤 如 下 ：(1) (0,1)U上生 成均 匀 分 布 的 随 机 数。-1(2

19、) , ( ) () 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( ) ( )1, 0 ( )log(1) log(1) ( ) 1 xXXF xexFyyXUUUU 例 ：设服从指数分布，则的分布函数为：通过计算

20、得，则：服从指数分布 其中服从均匀分布又因为 和有着同样的分布，所以也可以取： log( )XU 舍取方法舍取方法(Acceptance-Rejection )方法最早由 Von Neumann提出，现在已经广泛应用于各种随机数的生成。基本思路：通过一个容易生成的概率分布 g 和一个取舍准则生成另一个与 g 相近的概率分布 f 。具体步骤： ( ) ( ) ( )( ) 1, f xg xf xcg xcx 假设和均为集合 上的概率密度函数。且满足： 。易于从g(x)中生成样本X.用如下步骤生成有Y:1. ( ) ;2. (0,1),U3. ()/() g xXUXUf Xcg X生成 的样本

21、生成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( ).()/()1/ ( ) ( )(0,1) P Uf Xcg Xcf xcg xU由上面的证明过程我们可以看出每次接受的概率为 ）。也就是说

22、为了生成一个 的样本，我们需要平均生成个和 的样本。所以为了提高舍取法的效率，我们应该使 c 的取值尽可能的小，也就是使 f 和 g 的分布更为相近。3.标准正态分布随机数的生成正态分布是概率统计中最重要的分布，在此我们着重讨论如何生成标准正态分布随机数。引理：121/21121/221212, (0,1) ( 2ln)cos(2) ( 2ln)sin(2) U UUZUUZUUZZ 设 是独立同分布的变量，令：则证明：与 独由随机立，且均服从向量的变换公标准式易正态分布。得，略。Box-Muller 算法121. ,U U生成两个独立的均匀分布随机数 1/21121/221212 ( 2ln

23、)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间均匀分布的随机数r=a+(b-a)*rand;returnMatlab程序Function r=rnd-beta(lamda)%模拟指数分布%lamda表示指

24、数分布的参数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分钟分钟之间变化之间变化一艘船只卸货的时间由所

25、卸货物的类型决一艘船只卸货的时间由所卸货物的类型决定，在定，在4545分钟到分钟到9090分钟之间变化分钟之间变化 需要回答以下问题：需要回答以下问题： 1 1每艘船只在港口的平均时间和最长时间是多少每艘船只在港口的平均时间和最长时间是多少? ? 2 2若一艘船只的等待时间是从到达到开始卸货的若一艘船只的等待时间是从到达到开始卸货的时间，每艘船只的平均等待时间和最长等待时间是多时间，每艘船只的平均等待时间和最长等待时间是多少少? ? 3 3卸货设备空闲时间的百分比是多少卸货设备空闲时间的百分比是多少? ? 4 4船只排队最长的长度是多少船只排队最长的长度是多少? ?3 3 实例二实例二 为了得

26、到一些合理的答案，为了得到一些合理的答案，利用计算机或可编程利用计算机或可编程计算器来模拟港口的活动计算器来模拟港口的活动假定相邻两艘船到达的假定相邻两艘船到达的时间间隔和每艘船只卸货的时问在它们各自的时间时间间隔和每艘船只卸货的时问在它们各自的时间区问内均匀分布区问内均匀分布，例如两艘船到达的时间间隔可以，例如两艘船到达的时间间隔可以是是1515到到145145之间的任何整数，且这个区间内的任何整之间的任何整数，且这个区间内的任何整数等可能地出现数等可能地出现 在给出模拟这个港口系统的一般算法之前，在给出模拟这个港口系统的一般算法之前，考虑考虑有有5 5艘船只的假想情况对每艘船只有以下数据：

27、艘船只的假想情况对每艘船只有以下数据： 设想码头设备的拥有者关心他们设想码头设备的拥有者关心他们提供服务的质量提供服务的质量，并且要并且要评价各种管理模式以确定为了改善服务是否值评价各种管理模式以确定为了改善服务是否值得增加费用得增加费用做一些统计可以帮助对服务质量的评价做一些统计可以帮助对服务质量的评价。 例如，呆在港口时间最长的船只是船例如，呆在港口时间最长的船只是船5 5，呆了，呆了130130分分钟，而平均是钟，而平均是8989分钟分钟( (表表5-14)5-14)通常顾客对等待时间通常顾客对等待时间的长短非常在乎，例中最长的等待时间为的长短非常在乎，例中最长的等待时间为5555分钟，

28、而分钟，而平均是平均是2626分钟如果排队太长有些顾客会改到别处去分钟如果排队太长有些顾客会改到别处去做生意，例中最长的队是做生意，例中最长的队是2 2 用下面的蒙特卡罗模拟算法可以做这些统计，对各用下面的蒙特卡罗模拟算法可以做这些统计，对各种管理模式进行估价种管理模式进行估价 现在假定你是码头设备拥有者的顾问，如果能现在假定你是码头设备拥有者的顾问，如果能够雇用更多的劳动力，或者得到更好的卸货设备，够雇用更多的劳动力，或者得到更好的卸货设备，使卸货时间减少到每艘船使卸货时间减少到每艘船35357575分钟分钟，会有什么影，会有什么影响响? ?表表5-165-16给出了基于模拟算法的结果。给出

29、了基于模拟算法的结果。从表从表5-165-16可以看到，每艘船的卸货时间减少可以看到，每艘船的卸货时间减少15152020分钟，使得分钟，使得船只船只呆在港口的时间特别是等待时间缩短呆在港口的时间特别是等待时间缩短了，而了，而设备空闲时问设备空闲时问的百分比却增加了近一倍的百分比却增加了近一倍船主对此是满意的，因为这提高了船主对此是满意的，因为这提高了长期行驶时每艘船运送货物的效率长期行驶时每艘船运送货物的效率 这样，人港贸易好像会增加如果这样，人港贸易好像会增加如果贸易量增加使贸易量增加使得相邻两艘船到达的时间间隔缩短到得相邻两艘船到达的时间间隔缩短到1010到到120120分钟分钟之间之间

30、，模拟结果如表，模拟结果如表5-175-17从这个表可以看到，随从这个表可以看到，随着贸易量的增加，船只又要在港口呆更长的时间，着贸易量的增加，船只又要在港口呆更长的时间，但设备空闲时间少多了，于是船主和设备拥有者都但设备空闲时间少多了，于是船主和设备拥有者都随着贸易量的增加而受益随着贸易量的增加而受益 假定我们假定我们对两艘船到达的时间间隔和每艘船的卸对两艘船到达的时间间隔和每艘船的卸货时间分别在货时间分别在15between15betweeni145145和和45unload45unloadi9090内均匀分布不满意内均匀分布不满意，决定收集港口系统的经验数据，决定收集港口系统的经验数据，

31、并将结果并人我们的模型设想观测了利用港口卸并将结果并人我们的模型设想观测了利用港口卸货的货的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个,则所求定积分为 .dab