蒙特卡罗方法介绍及其建模应用.pptx
2023/2/18 19:25课程说明公用邮箱： key：ahualian2008参考书目:n黄燕、吴平.SAS统计分析及应用，机械工业出版社.n陈杰.Matlab宝典，电子工业出版社.n张文彤等.SPSS11.0统计分析教程，北京希望电子出版社.n薛益、陈立萍.统计建模与R软件，清华大学出版社.第1页/共72页2023/2/18 19:25主要内容蒙特卡洛方法应用实例2排队论模拟介绍3蒙特卡洛方法介绍12009-B 眼科病床安排应用4第2页/共72页2023/2/18 19:25蒙特卡洛方法介绍蒙特卡洛起源与发展1蒙特卡洛模拟误差分析2随机数的产生原理3第3页/共72页2023/2/18 19:25蒙特卡洛起源与发展1第4页/共72页2023/2/18 19:25模拟的概念模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。模拟的基本思想是建立一个试验模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息第5页/共72页2023/2/18 19:25模拟的方法物理模拟对实际系统及其过程用功能相似的实物系统去模仿。例如，军事演习、船艇实验、沙盘作业等 物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难许多系统无法进行物理模拟，如社会经济系统、生态系统等数学模拟在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟计算机模拟可以反复进行，改变系统的结构和系数都比较容易在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。第6页/共72页2023/2/18 19:25Monte Carlo的起源Monte Carlo方法：又称随机模拟方法，对研究的系统进行随机观察抽样,通过对样本值的统计分析，求得所研究系统的某些参数它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来，“曼哈顿计划”主持人之一、数学家：冯诺伊曼用驰名世界的赌城摩纳哥最大的城市Monte Carlo来命名这种方法Monte-Carlo,Monaco第7页/共72页2023/2/18 19:25MC的起源和发展在大众的心目中，科学的代言人是“心不在焉”的牛顿或者“爆炸式“发型的爱因斯坦但这只是传统形象，比他们更了解现代计算技术的冯诺伊曼是个”衣着考究，风度翩翩“的人物，他说：纯粹数学和应用数学的许多分支非常需要计算工具，用以打破目前由于纯粹分析的研究方法不能解决非线性问题而形成的停滞状态Monte Carlo方法是现代计算技术的最为杰出的成果之一，它在工程领域的作用是不可比拟的第8页/共72页2023/2/18 19:25Monte Carlo方法的应用物理：核物理，热力学与统计物理，粒子输运问题等数学：多重积分、解微分方程、非线性方程组求解等工程领域：真空技术，水力学，激光技术等经济学领域：期权定价、项目管理、投资风险决策等其他领域：化学、医学，生物，生产管理、系统科学、公用事业等方面，随着科学技术的发展，其应用范围将更加广泛。第9页/共72页2023/2/18 19:25MC的起源和发展事实上，Monte Carlo方法的基本思想很早以前就被人们所发现和利用：早在17世纪，人们就知道用事件发生的“频率”来近似事件的“概率”18世纪下半叶法国学者Buffon提出用投针试验的方法来确定圆周率的值的Buffon投针试验是Monte Carlo方法的最早的尝试历史上曾有几位学者相继做过这样的试验：试验费时费力精度不够高实施困难随着计算机技术的飞速发展，人们不需要具体实施这些试验，而只要在计算机上进行大量的、快速的模拟试验第10页/共72页2023/2/18 19:25 Buffon 试验假设平面上有无数条距离为1的等距平行线，现向该平面随机投掷一根长度为l的针(l1)，则我们可计算该针与任一平行线相交的概率。这里，随机投针指的是：针的中心点与最近的平行线间的距离X均匀地分布在区间0,1/2上，针与平行线的夹角（不管相交与否）均匀的分布在区间0,上。此时，针与线相交的充要条件是第11页/共72页2023/2/18 19:25从而针线相交的概率为根据上式，若我们做大量的投针试验并记录针与线相交的次数，则由大数定理可以估计出针线相交的概率p，从而得到 的估计值。Buffon 试验第12页/共72页2023/2/18 19:25function 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)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第14页/共72页2023/2/18 19:25建立统计模型，主要特征参量方面要与实际问题或系统相一致，问题的解对应于模型中随机变量的概率分布或其某些数字特征根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数，进而进行随机模拟实验根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）按照所建立模型进行仿真试验、计算，求出问题的随机解统计分析模拟试验结果，给出问题的估计以及其精度估计。必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。用蒙特卡洛方法进行计算机模拟的步骤第15页/共72页2023/2/18 19:25蒙特卡洛模拟误差分析2第16页/共72页2023/2/18 19:25蒙特卡洛模拟的理论基础大数定律-贝努里（Bernoulli）大数定律中心极限定理第17页/共72页2023/2/18 19:25蒙特卡洛模拟的误差分析由中心极限定理可知：这表明，不等式 近似地以概率1成立。上式也表明，收敛到 的阶为O(n-1/2)。通常，蒙特卡罗方法的误差 定义为：第18页/共72页2023/2/18 19:25蒙特卡洛模拟的误差分析 关于蒙特卡罗方法的误差需说明两点：1.蒙特卡罗方法的误差为概率误差，也即蒙特卡罗方法的收敛是概率意义下的收敛，虽然不能断言其误差不超过某个值，但能指出其误差以接近1的概率不超过某个界限 例如：=0.5，误差 此时，误差超过的概率与小于的概率1-相等，都等于0.5。2.误差中的均方差 是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出。第19页/共72页2023/2/18 19:25蒙特卡洛模拟的误差分析 减小方差的各种技巧 显然，当给定置信度后，误差由和n决定。要减小，或者是增大n，或者是减小方差2。在固定的情况下，要把精度提高一个数量级，试验次数n需增加两个数量级。因此，单纯增大n不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于n增大四倍的效果(n=(u/)2)。一般来说，降低方差的技巧，往往会使观察一个子样的时间增加，在固定时间内，使观察的样本数减少。第20页/共72页2023/2/18 19:25蒙特卡洛模拟的效率分析蒙特卡罗方法中效率：用来衡量一种方法的优劣，其由方差和观察一个子样的费用（使用计算机的时间）两者来衡量，它定义为nc:nc=(u/)2 2c 其中c是观察一个子样的平均费用,显然它与 2c成正比。总而言之，作为提高蒙特卡洛方法效率的重要方向，是在减小标准差的同时兼顾考虑费用大小，使 2c尽可能地小。第21页/共72页2023/2/18 19:25蒙特卡洛方法的特点Monte Carlo方法及其程序结构简单产生随机数，通过大量简单重复抽样和简单计算计算相应的值收敛速度与问题维数无关Monte Carlo方法的收敛速度为O(n-1/2)，与一般数值方法相比很慢。因此，用Monte Carlo方法不能解决精确度要求很高的问题Monte Carlo方法误差只与标准差和样本容量n有关，而与样本所在空间无关，即Monte Carlo方法的收敛速度与问题维数无关，而其他数值方法则不然。Monte Carlo方法的适用性强Monte Carlo方法对多维问题的适用性在解题时受问题条件限制的影响较小 例如：要计算s维空间中的任一区域Ds上的积分第22页/共72页2023/2/18 19:25随机数的产生原理3第23页/共72页2023/2/18 19:25常用分布的随机数生成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)来模拟它第24页/共72页2023/2/18 19:25常用分布的随机数生成2.正态分布N(，2)产生m*n阶均值为，标准差为的正态分布的随机数矩阵：R=normrnd(，m，n)适用范围：当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布第25页/共72页2023/2/18 19:25常用分布的随机数生成3.指数分布E()产生m*n阶均值为的指数分布的随机数矩阵：R=exprnd(，m，n)适用范围：排队服务系统中顾客到达间隔、质量与可靠性中电子元件的寿命通常服从指数分布。例：顾客到达某商店的间隔时间服从参数为10(分钟)的指数分布(指数分布的均值为10)指两个顾客到达商店的平均间隔时间是10分钟.即平均10分钟到达1个顾客.顾客到达的间隔时间可用exprnd(10)模拟。第26页/共72页2023/2/18 19:25常用分布的随机数生成4.泊松分布()产生m*n阶均值为的泊松分布的随机数矩阵R=poissrnd(，m，n)适用范围：泊松分布在排队系统、产品检验、天文、物理等领域有广泛应用5.二项分布 B(n,p)产生mn个参数为n,p的二项分布的随机数R=binornd(n,p,m,n)第27页/共72页2023/2/18 19:25一般分布随机数产生方法基本方法有如下三种：逆变换法复合抽样方法筛选法 第28页/共72页2023/2/18 19:25逆变换法(直接抽样方法)设随机变量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)即可。其步骤为第29页/共72页2023/2/18 19:251离散型分布即=infx:F(x)u其中令=1时为了实现由任意离散型分布的随机抽样，直接抽样方法是非常理想的!第30页/共72页2023/2/18 19:251离散型分布例1.掷骰子点数的抽样按照离散分布的直接抽样:(1)由U(0,1)抽取u即：等价于：也可使用如下更简单的方法第31页/共72页2023/2/18 19:25function 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离散型分布第32页/共72页2023/2/18 19:251离散型分布例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第33页/共72页2023/2/18 19:25cdfplot(liti12(1000)cdfplot(liti12(1000)+20)cdfplot(liti12(1000)+500)第34页/共72页2023/2/18 19:252连续分布对于连续型分布，如果分布函数F(x)的反函数F1(x)能够解析表示，则直接抽样方法是:第35页/共72页2023/2/18 19:25 在a，b上均匀分布的分布函数为：则 (1)由U(0,1)抽取u例3.在a，b上均匀分布的抽样第36页/共72页2023/2/18 19:25指数分布为连续型分布，其一般形式如下：其分布函数为：则 (1)由U(0,1)抽取u 因为1-u 也是(0,1)上均匀随机数，可将上式简化为 例4.指数分布第37页/共72页2023/2/18 19:25Randnum=（-2）*log(unifrnd(0,1,1,1000);cdfplot(Randnum)例5.产生指数分布 的随机数第38页/共72页2023/2/18 19:25 设X分布函数为F(x)，X1,Xn独立且与同分布，试设X FX(x),Y FY(y),且相互独立,M=maxX,Y，N=minX,Y,求 M,N 的分布函数.例6.极值分布第39页/共72页2023/2/18 19:25推广至相互独立的 n 个随机变量的情形：相互独立，且设则当X1,Xn相互独立相同分布函数F(x)时，有 FM(z)=F(z)n FN(z)=1-1-F(z)n第40页/共72页2023/2/18 19:25需要指出的是，当X1,Xn相互独立且具有相同分布函数F(x)时,常称M=max(X1,Xn)，N=min(X1,Xn)为极值.适用范围：由于一些灾害性的自然现象，如地震、洪水等等都是极值，研究极值分布具有重要的意义和实用价值.第41页/共72页2023/2/18 19:25例7 设系统 L 由相互独立的 n 个元件组成，连 接方式为：(1)串联；并联；如果 n 个元件的寿命分别为且求在以上 2 种组成方式下，系统 L 的寿命 X 的密度函数.第42页/共72页2023/2/18 19:25解:(1)(2)第43页/共72页2023/2/18 19:25例例8 设Xi 分布函数为生成n=20的1行10000列随机数，并画经验分布函数曲线。n=20Randnum=1-(1-unifrnd(0,1,1,10000).(1/n);cdfplot(Randnum)第44页/共72页2023/2/18 19:25逆变换法(直接抽样方法)连续性分布函数的直接抽样方法对于分布函数的反函数容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的：分布函数无法用解析形式表达，因而无法给出反函数的解析形式分布函数有解析形式，但是反函数的解析形式给不出来反函数有解析形式，但运算量很大下面叙述的抽样方法是能够克服这些困难的比较好的方法。第45页/共72页2023/2/18 19:25复合抽样方法复合抽样方法的基本思想是由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)第46页/共72页2023/2/18 19:25特殊复合分布在实际问题中，经常有这样的随机变量，它服从的分布与某随机变量的取值有关，而该随机变量服从一确定分布，例如，概率密度这是一个复合分布，其中fn(x)为与n有关的概率密度，pn0，n1，且第47页/共72页2023/2/18 19:25设Y为一离散型随机变量，它可能的取值为1,2,n,，取这些值的概率分别为p1,p2,pn,，Y的分布函数为：fn(x)为给定Y=n时X的条件密度。该复合分布f(x)的抽样方法为:首先从离散分布F(y)中抽样N然后再从密度函数fN(x)中抽样确定XfN 第48页/共72页2023/2/18 19:25 设有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 混合分布抽样第49页/共72页2023/2/18 19:25分布密度函数：的抽样方法为：1.首先从Y的离散分布中抽样N，N=1或2。根据:得第50页/共72页2023/2/18 19:25即Y为离散型随机变量，取1，2两个值，取1的概率为0.3，取2的概率为0.7。第51页/共72页2023/2/18 19:25分布密度函数：的抽样方法为：2.所以从 f(x)的抽样如下进行从 U(0,1)中抽取u第52页/共72页2023/2/18 19:253.画 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第53页/共72页2023/2/18 19:25Liti19(100)Liti19(1000)Liti19(10000)第54页/共72页2023/2/18 19:25指数函数分布的一般形式为例10 指数函数分布的抽样则使用复合抽样方法，抽取服从该分布的样本，生成10000个随机数，画经验分布函数，n=5.引入如下两个密度函数：第55页/共72页2023/2/18 19:25对应分布函数为使用复合抽样方法，首先从f1(y)中抽取y从U(0,1)中抽取u,令第56页/共72页2023/2/18 19:25指数分布,均值为1/y再由f2(x|yf1)中抽取xf 第57页/共72页2023/2/18 19:25function 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)第58页/共72页2023/2/18 19:25 筛 选 抽 样定理:设 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的抽样可如下进行：第59页/共72页2023/2/18 19:25例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)第60页/共72页2023/2/18 19:252.筛选抽样方法：取：则抽样框图为 显然，没有必要舍弃u1u2的情况，此时，只需取：亦即第61页/共72页2023/2/18 19:25function 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)第62页/共72页2023/2/18 19:25随机向量的抽样：直接抽样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的条件密度函数，即第63页/共72页2023/2/18 19:25根据上述边缘密度函数和条件密度函数，二维分布f(x,y)的抽样方法为:首先由 fl(x)中抽取 xf1，再由 f2(y|xf1)中抽样yf2(xf1,yf2)就是该二维分布的一个抽样例12 对下面二维分布进行抽样将 f(x,y)写为：其中第64页/共72页2023/2/18 19:25直接抽样方法：1.首先由 fl(x)中抽取 xf1:2.再由 f2(y|xf1)中抽样yf2:第65页/共72页2023/2/18 19:25function 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)第66页/共72页2023/2/18 19:25随机向量的抽样：筛选抽样定 理：设 (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)可如下进行：第67页/共72页2023/2/18 19:25例13 生成单位圆内均匀分布的10000个随机数，并画散点图。相当于-1，1-1，1上均匀分布。第68页/共72页2023/2/18 19:25mm=10000；xRandnum=zeros(1,mm);yRandnum=zeros(1,mm);ii=1;while iimm Randnum1=unifrnd(-1,1);Randnum2=unifrnd(-1,1);s=Randnum12+Randnum22；if s=1;xRandnum(1,ii)=Randnum1;yRandnum(1,ii)=Randnum2;ii=ii+1;endendplot(xRandnum,yRandnum,.)第69页/共72页2023/2/18 19:25实习一1.生成单位球内均匀分布的1行10000列随机数，并画散点图；2.设密度函数为 生成 (n=2，期望=5)的随机数10000个，并绘制经验分布图。为常数第70页/共72页2023/2/18 19:25实习一3.设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为生成10000个随机数，并画散点图。第71页/共72页2023/2/18 19:25南京信息工程大学感谢您的观看！第72页/共72页
