蒙特卡罗方法介绍及其建模应用PartI.ppt

Monte-Carlo 方法介绍及其建模应用,朱连华南京信息工程大学数学与统计学院,2023/10/30*,南京信息工程大学,课程说明,公用邮箱：key：ahualian2008参考书目:黄燕、吴平.SAS统计分析及应用，机械工业出版社.陈杰.Matlab宝典，电子工业出版社.张文彤等.SPSS11.0统计分析教程，北京希望电子出版社.薛益、陈立萍.统计建模与R软件，清华大学出版社.,2023/10/30*,南京信息工程大学,主要内容,蒙特卡洛方法应用实例,2,排队论模拟介绍,3,蒙特卡洛方法介绍,1,2009-B 眼科病床安排应用,4,2023/10/30*,南京信息工程大学,蒙特卡洛方法介绍,蒙特卡洛起源与发展,1,蒙特卡洛模拟误差分析,2,随机数的产生原理,3,2023/10/30*,南京信息工程大学,蒙特卡洛起源与发展,1,2023/10/30*,南京信息工程大学,模拟的概念,模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。模拟的基本思想是建立一个试验模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息,2023/10/30*,南京信息工程大学,模拟的方法,物理模拟对实际系统及其过程用功能相似的实物系统去模仿。例如，军事演习、船艇实验、沙盘作业等 物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难许多系统无法进行物理模拟，如社会经济系统、生态系统等数学模拟在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟计算机模拟可以反复进行，改变系统的结构和系数都比较容易在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。,2023/10/30*,南京信息工程大学,蒙特卡洛（Monte Carlo）方法,蒙特卡洛方法（Monte Carlo，简写为MC）是一种应用随机数来进行计算机模拟的方法 对研究的系统进行随机观察抽样通过对样本值的统计分析，求得所研究系统的某些参数,2023/10/30*,南京信息工程大学,MC的起源和发展,Monte Carlo方法：一种基于“随机数”的随机模拟方法，源于美国在第一次世界大战进行的研制原子弹的“曼哈顿计划”该计划主持人之一、数学家：冯诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩冯诺伊曼是公理化方法和计算机体系的领袖人物,Monte Carlo方法也是他的重要贡献,2023/10/30*,南京信息工程大学,MC的起源和发展,事实上，Monte Carlo方法的基本思想很早以前就被人们所发现和利用：早在17世纪，人们就知道用事件发生的“频率”来近似事件的“概率”18世纪下半叶法国学者Buffon提出用投针试验的方法来确定圆周率的值的Buffon投针试验是Monte Carlo方法的最早的尝试历史上曾有几位学者相继做过这样的试验：试验费时费力精度不够高实施困难随着计算机技术的飞速发展，人们不需要具体实施这些试验，而只要在计算机上进行大量的、快速的模拟试验,2023/10/30*,南京信息工程大学,MC的起源和发展,在大众的心目中，科学的代言人是“心不在焉”的牛顿或者“爆炸式“发型的爱因斯坦但这只是传统形象，比他们更了解现代计算技术的冯诺伊曼是个”衣着考究，风度翩翩“的人物，他说：纯粹数学和应用数学的许多分支非常需要计算工具，用以打破目前由于纯粹分析的研究方法不能解决非线性问题而形成的停滞状态Monte Carlo方法是现代计算技术的最为杰出的成果之一，它在工程领域的作用是不可比拟的,2023/10/30*,南京信息工程大学,Buffon 试验,假设平面上有无数条距离为1的等距平行线，现向该平面随机投掷一根长度为l的针(l1)，则我们可计算该针与任一平行线相交的概率。这里，随机投针指的是：针的中心点与最近的平行线间的距离X均匀地分布在区间0,1/2上，针与平行线的夹角（不管相交与否）均匀的分布在区间0,上。此时，针与线相交的充要条件是,2023/10/30*,南京信息工程大学,从而针线相交的概率为根据上式，若我们做大量的投针试验并记录针与线相交的次数，则由大数定理可以估计出针线相交的概率p，从而得到 的估计值。,Buffon 试验,2023/10/30*,南京信息工程大学,function piguji=buffon(llength,mm)%llength 是针的长度%mm 是随机实验次数frq=0;xrandnum=unifrnd(0,0.5,1,mm);phi=unifrnd(0,pi,1,mm);for ii=1:mm if(xrandnum(1,ii)=(llength*sin(phi(1,ii)/2)frq=frq+1;end end piguji=2*llength/(frq/mm),2023/10/30*,南京信息工程大学,buffon(.6,1000)piguji=3.1662 buffon(.6,10000)piguji=3.1072 buffon(.6,100000)piguji=3.1522 buffon(.6,1000000)piguji=3.1386 buffon(.6,1000000)piguji=3.1451 buffon(.6,1000000)piguji=3.1418 buffon(.6,1000000)piguji=3.1448 buffon(.6,1000000)piguji=3.1405 buffon(.6,1000000)piguji=3.1394,2023/10/30*,南京信息工程大学,建立统计模型，主要特征参量方面要与实际问题或系统相一致，问题的解对应于模型中随机变量的概率分布或其某些数字特征根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数，进而进行随机模拟实验根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）按照所建立模型进行仿真试验、计算，求出问题的随机解统计分析模拟试验结果，给出问题的估计以及其精度估计。必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。,用蒙特卡洛方法进行计算机模拟的步骤,2023/10/30*,南京信息工程大学,蒙特卡洛模拟误差分析,2,2023/10/30*,南京信息工程大学,蒙特卡洛模拟的理论基础,大数定律-贝努里（Bernoulli）大数定律中心极限定理,2023/10/30*,南京信息工程大学,蒙特卡洛模拟的误差分析,由中心极限定理可知：这表明，不等式 近似地以概率1成立。上式也表明，收敛到 的阶为O(n-1/2)。通常，蒙特卡罗方法的误差 定义为：,2023/10/30*,南京信息工程大学,蒙特卡洛模拟的误差分析,关于蒙特卡罗方法的误差需说明两点：1.蒙特卡罗方法的误差为概率误差，也即蒙特卡罗方法的收敛是概率意义下的收敛，虽然不能断言其误差不超过某个值，但能指出其误差以接近1的概率不超过某个界限 例如：=0.5，误差 此时，误差超过的概率与小于的概率1-相等，都等于0.5。2.误差中的均方差 是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出。,2023/10/30*,南京信息工程大学,蒙特卡洛模拟的误差分析,减小方差的各种技巧 显然，当给定置信度后，误差由 和n决定。要减小，或者是增大n，或者是减小方差2。在 固定的情况下，要把精度提高一个数量级，试验次数n需增加两个数量级。因此，单纯增大n不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于n增大四倍的效果(n=(u/)2)。一般来说，降低方差的技巧，往往会使观察一个子样的时间增加，在固定时间内，使观察的样本数减少。,2023/10/30*,南京信息工程大学,蒙特卡洛模拟的效率分析,蒙特卡罗方法中效率：用来衡量一种方法的优劣，其由方差和观察一个子样的费用（使用计算机的时间）两者来衡量，它定义为nc:nc=(u/)2 2c 其中c是观察一个子样的平均费用,显然它与 2c成正比。总而言之，作为提高蒙特卡洛方法效率的重要方向，是在减小标准差的同时兼顾考虑费用大小，使 2c尽可能地小。,2023/10/30*,南京信息工程大学,蒙特卡洛方法的特点,Monte Carlo方法及其程序结构简单产生随机数，通过大量简单重复抽样和简单计算计算相应的值收敛速度与问题维数无关Monte Carlo方法的收敛速度为O(n-1/2)，与一般数值方法相比很慢。因此，用Monte Carlo方法不能解决精确度要求很高的问题Monte Carlo方法误差只与标准差和样本容量n有关，而与样本所在空间无关，即Monte Carlo方法的收敛速度与问题维数无关，而其他数值方法则不然。Monte Carlo方法的适用性强Monte Carlo方法对多维问题的适用性在解题时受问题条件限制的影响较小 例如：要计算s维空间中的任一区域Ds上的积分,2023/10/30*,南京信息工程大学,随机数的产生原理,3,2023/10/30*,南京信息工程大学,常用分布的随机数生成,1.均匀分布U(a,b)产生m*n阶(a,b)均匀分布的随机数矩阵R=unifrnd(a,b,m,n)产生m*n阶(0,N)离散均匀分布的随机数矩阵R=unidrnd(N)R=unidrnd(N,mm,nn)适用范围：当只知道一个随机变量取值在(a,b)内，但不知道（也没理由假设）它在何处取值的概率大，在何处取值的概率小，就只好用U(a,b)来模拟它,2023/10/30*,南京信息工程大学,常用分布的随机数生成,2.正态分布N(，2)产生m*n阶均值为，标准差为的正态分布的随机数矩阵：R=normrnd(，m，n)适用范围：当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布,2023/10/30*,南京信息工程大学,常用分布的随机数生成,3.指数分布E()产生m*n阶均值为的指数分布的随机数矩阵：R=exprnd(，m，n)适用范围：排队服务系统中顾客到达间隔、质量与可靠性中电子元件的寿命通常服从指数分布。例：顾客到达某商店的间隔时间服从参数为10(分钟)的指数分布(指数分布的均值为10)指两个顾客到达商店的平均间隔时间是10分钟.即平均10分钟到达1个顾客.顾客到达的间隔时间可用exprnd(10)模拟。,2023/10/30*,南京信息工程大学,常用分布的随机数生成,4.泊松分布()产生m*n阶均值为的泊松分布的随机数矩阵R=poissrnd(，m，n)适用范围：泊松分布在排队系统、产品检验、天文、物理等领域有广泛应用5.二项分布 B(n,p)产生mn个参数为n,p的二项分布的随机数R=binornd(n,p,m,n),2023/10/30*,南京信息工程大学,一般分布随机数产生方法,基本方法有如下三种：逆变换法复合抽样方法筛选法,2023/10/30*,南京信息工程大学,逆变换法(直接抽样方法),设随机变量X 的分布函数为F(x)，定义 F-1(y)=infx:F(x)y,0y 1定理 设随机变量U服从(0,1)上的均匀分布，则X=F-1(U)的分布函数为F(x)。因此，要产生来自F(x)的随机数，只要先产生来自U(0,1)的随机数，然后计算F-1(u)即可。其步骤为,2023/10/30*,南京信息工程大学,1离散型分布,即,=infx:F(x)u,其中令=1时,为了实现由任意离散型分布的随机抽样，直接抽样方法是非常理想的!,2023/10/30*,南京信息工程大学,1离散型分布,例1.掷骰子点数的抽样按照离散分布的直接抽样:(1)由U(0,1)抽取u即：等价于：也可使用如下更简单的方法,2023/10/30*,南京信息工程大学,function discreterandom=liti11(mm)Random=unifrnd(0,1,1,mm);for i=1:mm if(floor(6*Random(1,i)=6*Random(1,i)Random(1,i)=6*Random(1,i);elseRandom(1,i)=floor(6*Random(1,i)+1;endendcdfplot(Random),1离散型分布,2023/10/30*,南京信息工程大学,1离散型分布,例2:生成1行1000列的110上离散均匀分布的随机数；并画经验分布函数曲线。生成1行1000列2130上离散均匀分布的随机数；并画经验分布函数曲线。生成1行1000列501510上离散均匀分布的随机数;并画经验分布函数曲线。function Random=liti12(mm)Random=unifrnd(0,1,1,mm);for i=1:mm if(floor(10*Random(1,i)=10*Random(1,i)Random(1,i)=10*Random(1,i);else Random(1,i)=floor(10*Random(1,i)+1;endend,2023/10/30*,南京信息工程大学,cdfplot(liti12(1000)cdfplot(liti12(1000)+20)cdfplot(liti12(1000)+500),2023/10/30*,南京信息工程大学,2连续分布,对于连续型分布，如果分布函数F(x)的反函数F1(x)能够解析表示，则直接抽样方法是:,2023/10/30*,南京信息工程大学,在a，b上均匀分布的分布函数为：则(1)由U(0,1)抽取u,例3.在a，b上均匀分布的抽样,2023/10/30*,南京信息工程大学,指数分布为连续型分布，其一般形式如下：其分布函数为：则(1)由U(0,1)抽取u 因为1-u 也是(0,1)上均匀随机数，可将上式简化为,例4.指数分布,2023/10/30*,南京信息工程大学,Randnum=（-2）*log(unifrnd(0,1,1,1000);cdfplot(Randnum),例5.产生指数分布 的随机数,2023/10/30*,南京信息工程大学,设X分布函数为F(x)，X1,Xn独立且与同分布，试,设X FX(x),Y FY(y),且相互独立,M=maxX,Y，N=minX,Y,求 M,N 的分布函数.,例6.极值分布,2023/10/30*,南京信息工程大学,推广至相互独立的 n 个随机变量的情形：,则,当X1,Xn相互独立相同分布函数F(x)时，有,FM(z)=F(z)n FN(z)=1-1-F(z)n,2023/10/30*,南京信息工程大学,需要指出的是，当X1,Xn相互独立且具有相同分布函数F(x)时,常称,M=max(X1,Xn)，N=min(X1,Xn),为极值.,适用范围：由于一些灾害性的自然现象，如地震、洪水等等都是极值，研究极值分布具有重要的意义和实用价值.,2023/10/30*,南京信息工程大学,例7 设系统 L 由相互独立的 n 个元件组成，连 接方式为：,串联；并联；,如果 n 个元件的寿命分别为,且,求在以上 2 种组成方式下，系统 L 的寿命 X 的密度函数.,2023/10/30*,南京信息工程大学,解:(1),(2),2023/10/30*,南京信息工程大学,例8 设Xi 分布函数为,生成n=20的1行10000列随机数，并画经验分布函数曲线。,n=20Randnum=1-(1-unifrnd(0,1,1,10000).(1/n);cdfplot(Randnum),2023/10/30*,南京信息工程大学,逆变换法(直接抽样方法),连续性分布函数的直接抽样方法对于分布函数的反函数容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的：分布函数无法用解析形式表达，因而无法给出反函数的解析形式分布函数有解析形式，但是反函数的解析形式给不出来反函数有解析形式，但运算量很大下面叙述的抽样方法是能够克服这些困难的比较好的方法。,2023/10/30*,南京信息工程大学,逆变换法(复合抽样方法),复合抽样方法的基本思想是由kahn提出的。考虑如下复合分布：其中f2(x|y)为给定Y=y时X的条件密度，F1(y)为Y的分布函数如果X密度函数f(x)难于抽样，而X关于Y的条件密度函数f2(x|y)以及Y的分布F1(y)均易于抽样，则X的随机数抽样：首先从分布F1(y)中抽样YF1，然后再从密度函数f2(x|YF1)中抽样确定Xf2(x|YF),2023/10/30*,南京信息工程大学,逆变换法(特殊复合分布),在实际问题中，经常有这样的随机变量，它服从的分布与某随机变量的取值有关，而该随机变量服从一确定分布，例如，概率密度这是一个复合分布，其中fn(x)为与n有关的概率密度，pn0，n1，且,2023/10/30*,南京信息工程大学,设Y为一离散型随机变量，它可能的取值为1,2,n,，取这些值的概率分别为p1,p2,pn,，Y的分布函数为：fn(x)为给定Y=n时X的条件密度。该复合分布f(x)的抽样方法为:首先从离散分布F(y)中抽样N然后再从密度函数fN(x)中抽样确定XfN,2023/10/30*,南京信息工程大学,设有X的密度函数为：生成10000个随机数，画经验分布函数，并画分布函数曲线。分析：它相当于设Y为离散型随机变量，取1，2两个值，取1的概率为0.3，取2的概率为0.7。当Y取1时，X的条件密度函数为:f(x|Y=1)=2e-2x,x0;当Y取2时，X的条件密度函数为:f(x|Y=2)=e-x,x0.,例9 混合分布抽样,2023/10/30*,南京信息工程大学,分布密度函数：的抽样方法为：,1.首先从Y的离散分布中抽样N，N=1或2。根据:,得,2023/10/30*,南京信息工程大学,即,Y为离散型随机变量，取1，2两个值，取1的概率为0.3，取2的概率为0.7。,2023/10/30*,南京信息工程大学,分布密度函数：的抽样方法为：,2.所以从 f(x)的抽样如下进行,从 U(0,1)中抽取u,2023/10/30*,南京信息工程大学,3.画 f(x)对应的分布函数图,ezplot(F(x),a,b)表示在axb绘制显函数F(x)的函数图，其中,function liti19(mm)R=unifrnd(0,1,mm,1);R1=exprnd(0.5,mm,1);R2=exprnd(1,mm,1);xR=zeros(mm,1);for ii=1:mm if R(ii,1)=0.3 xR(ii,1)=R1(ii,1);else xR(ii,1)=R2(ii,1);endend cdfplot(xR);hold onezplot(0.7*(1-exp(-x)+0.3*(1-exp(-2*x),0,10)hold off,2023/10/30*,南京信息工程大学,Liti19(100),Liti19(1000),Liti19(10000),2023/10/30*,南京信息工程大学,指数函数分布的一般形式为,例10 指数函数分布的抽样,则,使用复合抽样方法，抽取服从该分布的样本，生成10000个随机数，画经验分布函数，n=5.,引入如下两个密度函数：,2023/10/30*,南京信息工程大学,对应分布函数为,使用复合抽样方法，首先从f1(y)中抽取y,从U(0,1)中抽取u,令,2023/10/30*,南京信息工程大学,指数分布,均值为1/y,再由f2(x|yf1)中抽取xf,2023/10/30*,南京信息工程大学,function liti110(n,mm)R1=unifrnd(0,1,mm,1);R2=unifrnd(0,1,mm,1);x=zeros(mm,1);y=1./R1.(1/n)x=-log(R2)./ycdfplot(x),function liti110(n,mm)R1=unifrnd(0,1,mm,1);x=zeros(mm,1);y=1./R1.(1/n);x=exprnd(1./y);cdfplot(x),使用复合抽样方法，首先从f1(y)中抽取y:再由f2(x|yf1)中抽取X:,Liti110(5,10000),2023/10/30*,南京信息工程大学,筛 选 抽 样,定理:设 X 的密度函数f(x)，且可将其表示成f(x)=ch(x)g(x)，其中0g(x)1，c1是常数，h(.)是一个密度函数,令U和Y 分别服从U(0,1)和h(y)，则在U g(Y)的条件下，Y的条件密度为:依据上述定理，若h(.)易于抽样，则X的抽样可如下进行：,2023/10/30*,南京信息工程大学,例11 令圆半径为R0，该圆上的点到圆心的距离为r，r的密度函数和分布函数分别为：生成10000个随机数，画经验分布函数。1.直接抽样方法：缺点：开方运算在计算机 上很费时间,function liti111_1(R0,mm)R=unifrnd(0,1,mm,1);x=R0*sqrt(R);cdfplot(x),Liti111_1(3,10000),2023/10/30*,南京信息工程大学,2.筛选抽样方法：取：则抽样框图为,显然，没有必要舍弃u1u2的情况，此时，只需取：亦即,2023/10/30*,南京信息工程大学,function liti111_2(R0,mm)R1=unifrnd(0,1,mm,1);R2=unifrnd(0,1,mm,1);x=zeros(mm,1);for ii=1:mm if R1(ii,1)=R2(ii,1)x(ii,1)=R0*R2(ii,1);else x(ii,1)=R0*R1(ii,1);endend cdfplot(x),Liti111_2(3,10000),2023/10/30*,南京信息工程大学,随机向量的抽样：直接抽样,1.分量X与Y相互独立，随机向量(X,Y)的抽样 若X,Y相互独立,分布函数分别为FX(x),FY(y),则从FX(x)中抽样x,从FY(y)中抽样y,得到二维随机变量(X,Y)的抽样(x,y).2.按照条件分布，抽取随机向量(X,Y)的样本 二维随机向量(X,Y)的密度函数f(x,y)，将其表示如下：其中 fl(x)，f2(y|x)分别为X 的边缘密度函数和给定 X=x 条件下Y的条件密度函数，即,2023/10/30*,南京信息工程大学,根据上述边缘密度函数和条件密度函数，二维分布f(x,y)的抽样方法为:首先由 fl(x)中抽取 xf1，再由 f2(y|xf1)中抽样yf2(xf1,yf2)就是该二维分布的一个抽样,例12 对下面二维分布进行抽样,将 f(x,y)写为：其中,2023/10/30*,南京信息工程大学,直接抽样方法：1.首先由 fl(x)中抽取 xf1:2.再由 f2(y|xf1)中抽样yf2:,2023/10/30*,南京信息工程大学,function liti112(mm)R1=unifrnd(0,1,mm,1);R2=unifrnd(0,1,mm,1);x=1./R1;y=-R1.*log(R2);plot(x,y,.);axis(0.0 10 0 10)hold onhold off,liti112(1000),2023/10/30*,南京信息工程大学,随机向量的抽样：筛选抽样,定理：设(X,Y)的密度函数为 f(x,y)，且f(x,y)=ch(x,y)g(x,y)，其中0g(x,y)1，c1，h(x,y)是一个密度函数。令U和Z=(Xh,Yh)分别服从U(0,1)和h(x,y)，则在Ug(Z)的条件下，Z的条件密度为：若h(.)易于抽样，则从f(x,y)中的抽样(X,Y)可如下进行：,2023/10/30*,南京信息工程大学,例13 生成单位圆内均匀分布的10000个随机数，并画散点图。,相当于-1，1-1，1上均匀分布。,2023/10/30*,南京信息工程大学,实习一,1.生成单位球内均匀分布的1行10000列随机数，并画散点图；2.设密度函数为 生成(n=2，期望=5)的随机数10000个，并绘制经验分布图。,为常数,2023/10/30*,南京信息工程大学,实习一,3.设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为生成10000个随机数，并画散点图。,