蒙克卡罗方法

蒙特卡洛方法基本思想,实验目的,实验内容,学习计算机模拟的基本过程与方法。,1、模拟的概念。,4、实验作业。,3、计算机模拟实例。,2、产生随机数的计算机命令。,模拟的概念,模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。,模拟的基本思想是建立一个试验模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息,模拟的方法,1、物理模拟：对实际系统及其过程用功能相似的实物系统去模仿。例如，军事演习、船艇实验、沙盘作业等。,物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难。而且，许多系统无法进行物理模拟，如社会经济系统、生态系统等。,在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。,在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟。,2、数学模拟,计算机模拟可以反复进行，改变系统的结构和系数都比较容易。,蒙特卡洛（MonteCarlo）方法是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样，通过对样本值的观察统计，求得所研究系统的某些参数,蒙特卡洛方法也称为随机模拟方法,其起源最早可以追溯到18世纪下半叶的Buffon试验.,例在1777年,法国学者Buffon提出用试验方法求圆周率鸬闹.其原理如下：假设平面上有元数条距离为1的等矩平行线,现向该平面随机地投掷一根长度为KI1的针,则我们可以计算该针与任一平行线相交的概率.此处随机投针可以这样理解z针的中心点与最近的平行线间的距离Z均匀地分布在区间0.1/2上,针与平行线的夹角以不管相交与否)均匀地分布在区间0,而上(见图6。.于是,针与线相交的充要条件是本寸,从而针线相交概率为1,用蒙特卡洛方法进行计算机模拟的步骤：,1设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。2模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象,产生模拟随机数的计算机命令,在Matlab软件中，可以直接产生满足各种分布的随机数，命令如下：,2产生m*n阶离散均匀分布的随机数矩阵：R=unidrnd(N)R=unidrnd(N,mm,nn),当只知道一个随机变量取值在（a，b）内，但不知道（也没理由假设）它在何处取值的概率大，在何处取值的概率小，就只好用U（a，b）来模拟它。,1产生m*n阶a，b均匀分布U（a，b）的随机数矩阵：unifrnd(a,b,m,n)产生一个a，b均匀分布的随机数：unifrnd(a,b),当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。,机械加工得到的零件尺寸的偏差、射击命中点与目标的偏差、各种测量误差、人的身高、体重等，都可近似看成服从正态分布。,若连续型随机变量X的概率密度函数为其中0为常数，则称X服从参数为的指数分布。,指数分布的期望值为,排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。,指数分布在排队论、可靠性分析中有广泛应用。,注意：Matlab中，产生参数为的指数分布的命令为exprnd(),例顾客到达某商店的间隔时间服从参数为10的指数分布,指数分布的均值为10。指两个顾客到达商店的平均间隔时间是10个单位时间.即平均10个单位时间到达1个顾客.顾客到达的间隔时间可用exprnd(10)模拟。,设离散型随机变量X的所有可能取值为0,1,2,且取各个值的概率为其中0为常数，则称X服从参数为的帕松分布。,帕松分布在排队系统、产品检验、天文、物理等领域有广泛应用。,帕松分布的期望值为,1事件的频率在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数。频率f=m/n,2.频率的稳定性,掷一枚均匀硬币，记录掷硬币试验中频率P*的波动情况。R=binornd(N,P,mm,nn),例1频率的稳定性,3概率的频率定义,在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数。当试验次数n很大时，如果频率m/n稳定地在某数值p附近摆动，而且一般地说，随着试验次数的增加，这种摆动的幅度越来越小，称数值p为事件A在这一组不变的条件下发生的概率，记作P(A)=p.,4频率的基本性质,（1）对任意事件A，有,（2）,（3）若A1，A2，An是互不相容的，则,频率定义的意义:,(1)提供了估计概率的方法;(2)提供了一种检验理论正确与否的准则.,理论依据：大数定律,大量的随机现象中平均结果的稳定性,大数定律的客观背景,大量抛掷硬币正面出现频率,字母使用频率,生产过程中的废品率,大数定律,贝努里（Bernoulli）大数定律,设nA是n次独立重复试验中事件A发生的次数,p是每次试验中A发生的概率，则,有,或,贝努里（Bernoulli）大数定律的意义：,定义,a是一常数，,(或,故,在Bernoulli定理的证明过程中，Yn是相互独立的服从0-1分布的随机变量序列Xk的算术平均值,Yn依概率收敛于其数学期望p.,结果同样适用于服从其它分布的独立随机变量序列,Chebyshev大数定律,(指任意给定n1,相互独立)，且具有相同的数学期望和方差,或,定理的意义：,当n足够大时，算术平均值几乎就是一个常数,可以用算术平均值近似地代替数学期望.,具有相同数学期望和方差的独立随机变量序列的算术平均值依概率收敛于数学期望.,例如要估计某地区的平均亩产量，要收割某些有代表性的地块，例如n块.计算其平均亩产量，则当n较大时，可用它作为整个地区平均亩产量的一个估计.,辛钦大数定律,具有相同的分布，且,记,则,则,大数定律以严格的数学形式表达了随机现象最根本的性质之一：,它是随机现象统计规律的具体表现.,大数定律在理论和实际中都有广泛的应用.,平均结果的稳定性,例1频率的稳定性,1事件的频率在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数。频率f=m/n,2.频率的稳定性,掷一枚均匀硬币，记录掷硬币试验中频率P*的波动情况。R=binornd(N,P,mm,nn),functionliti1(n,p,mm)pro=zeros(1,mm);randnum=binornd(n,p,1,mm)a=0;fori=1:mma=a+randnum(1,i);pro(i)=a/i;endpro=pronum=1:mm;plot(num,pro),在Matlab中编辑.m文件输入以下命令：,在Matlab命令行中输入以下命令：,liti1(1,0.5,1000),在Matlab命令行中输入以下命令：,liti1(1,0.5,10000),练习频率的稳定性,1事件的频率R=binornd(N,P,mm,nn)在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数。频率f=m/n,2.频率的稳定性,练习掷一枚不均匀硬币，正面出现概率为0.3，记录前1000次掷硬币试验中正面频率的波动情况，并画图。,在Matlab命令行中输入以下命令：,liti1(1,0.3,1000),例2掷两枚不均匀硬币，每枚正面出现概率为0.4，记录前1000次掷硬币试验中两枚都为正面频率的波动情况，并画图。,在Matlab中编辑.m文件输入以下命令：,functionliti2(n,p,mm)pro=zeros(1,mm);randnum=binornd(n,p,2,mm)；a=0;fori=1:mma=a+randnum(1,i)*randnum(2,i);pro(i)=a/i;endpro=pro，num=1:mm;plot(num,pro),熊宇乐y=zeros(1,1000);a=binornd(1,0.4,1,1000);b=binornd(1,0.4,1,1000);c=0;d=0;fori=1:1000c=c+a(1,i).*b(1,i);y(i)=c/i;endy=y;num=1:1000;plot(num,y),孟亚functionbino(n,p,m)x=binornd(n,p,1,m);y=binornd(n,p,1,m);fori=1:mifx(i)=1endplot(y),liti2(1,0.4,100),liti2(1,0.4,10000),在一袋中有10个相同的球，分别标有号码1,2,10。每次任取一个球，记录其号码后放回袋中，再任取下一个。这种取法叫做“有放回抽取”。今有放回抽取3个球，求这3个球的号码均为偶数的概率。（用频率估计概率）,例3:,解：令A=有放回抽取3个球，求这3个球的号码均为偶数=(2,2,2)，(2,2,4)，.，(10,10,10),functionproguji=liti3(n,mm)frq=0;randnum=unidrnd(n,mm,3);proguji=0;fori=1:mma=(randnum(i,1)+1)*(randnum(i,2)+1)*(randnum(i,3)+1)；ifmod(a,2)=1frq=frq+1endend;proguji=frq/mm,例4两盒火柴，每盒20根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有5根火柴的概率有多大？（用频率估计概率）,liti4(20,5,100)proguji=0.4800liti4(20,5,1000)proguji=0.4970liti4(20,5,10000)proguji=0.4910liti4(20,5,100000)proguji=0.4984,functionproguji=liti4(nn,num,mm)%nn是每盒中的火柴数%num是剩余的火柴数%mm是随机实验次数frq=0;randnum=binornd(1,0.5,mm,2*nn);proguji=0;fori=1:mma1=0;a2=0;j=1;while(a1=5frq=frq+1;end%a1=a1,a2=a2,frq%pauseendproguji=frq/mm,二.几何概率,1.定义,向任一可度量区域G内投一点，如果所投的点落在G中任意可度量区域g内的可能性与g的度量成正比，而与g的位置和形状无关，则称这个随机试验为几何型随机试验。或简称为几何概型。,2.概率计算,1.P(A)=A的度量/S的度量,两人约定于12点到1点到某地会面，先到者等20分钟后离去，试求两人能会面的概率？,例5,解：设x,y分别为甲、乙到达时刻(分钟),令A=两人能会面=(x,y)|x-y|20，x60,y60,P(A)=A的面积/S的面积=（602-402）/602=5/9=0.5556,functionproguji=liti5(mm)%mm是随机实验次数frq=0;randnum1=unifrnd(0,60,mm,1);randnum2=unifrnd(0,60,mm,1);randnum=randnum1-randnum2;proguji=0;forii=1:mmifabs(randnum(ii,1)=20frq=frq+1;endendproguji=frq/mm,liti5(10000)proguji=0.5557,例6在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人,现在希望能用某种方式把我方将要对敌人实施的1打击结果显现出来，利用频率稳定性，确定有效射击的概率,分析：这是一个概率问题，可以通过理论计算得到相应的概率和期望值.,为了能显示我方射击的过程，现采用模拟的方式。,需要模拟出以下两件事：,1.问题分析,1观察所对目标的指示正确与否,模拟试验有两种结果，每一种结果出现的概率都是1/2,因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确,2当指示正确时，我方火力单位的射击结果情况,模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6),这时可用投掷骰子的方法来确定：如果出现的是、三个点：则认为没能击中敌人；如果出现的是、点：则认为毁伤敌人一门火炮；若出现的是点：则认为毁伤敌人两门火炮,2.符号假设,i：要模拟的打击次数；k1：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数E：有效射击比率；,3.模拟框图,functionliti6(p,mm)efreq=zeros(1,mm);randnum1=binornd(1,p,1,mm);randnum2=unidrnd(6,1,mm);k1=0;k2=0;k3=0;fori=1:mmifrandnum1(i)=0k1=k1+1;elseifrandnum2(i)=3k1=k1+1;elseifrandnum2(i)=6k3=k3+1;elsek2=k2+1;endendefreq(i)=(k2+k3)/i;endnum=1:mm;plot(num,efreq),在Matlab中编辑.m文件输入以下命令：,在Matlab命令行中输入以下命令：,liti6(0.5,2000),在Matlab命令行中输入以下命令：,liti6(0.5,20000),5.理论计算,6.结果比较,模拟结果与理论计算近似一致，能更加真实地表达实际战斗动态过程,3.某厂生产的灯泡能用1000小时的概率为0.8,能用1500小时的概率为0.4,求已用1000小时的灯泡能用到1500小时的概率,2.在一袋中有10个相同的球，分别标有号码1,2,10。今有放回任取两个球，求取得的第一个球号码为奇数，第二个球的号码为偶数的概率。,掷三枚不均匀硬币，每枚正面出现概率为0.3，记录前1000次掷硬币试验中至少两枚都为正面频率的波动情况，并画图。,作业:,某厂生产的灯泡能用1000小时的概率为0.8,能用1500小时的概率为0.4,求已用1000小时的灯泡能用到1500小时的概率,解令A灯泡能用到1000小时B灯泡能用到1500小时,所求概率为,例:,在一袋中有10个相同的球，分别标有号码1,2,10。今不放回任取两个球，求取得的第一个球号码为奇数，第二个球的号码为偶数的概率。,解：令A=抽取2个球，第一个球号码为奇数，第二个球的号码为偶数=(1,2)，(1,4)，.，(9,10),2020/4/26,圆的面积,单位圆的面积等于计算第一象限内的单位圆的面积，方法为：把它分成n个窄的曲边梯形，计算S大，S小，其中n可以为1000,10000,.,实验二:的计算,2020/4/26,数值积分法,2020/4/26,梯形公式,2020/4/26,辛普森公式,2020/4/26,无穷级数法,arctgx=x-x3/3+x5/5-x7/7+x9/9-./4=arctg1=1-1/3+1/5-1/7+1/9-收敛太慢！|x|应当比1小很多，级数收敛才快。/4=arctg1/2+arctg1/3/4=4arctg1/5-arctg1/239,2020/4/26,某次随机投点的结果,Monte-Carlo方法,谢文贤西北工业大学理学院高性能计算中心2004-05-14,Monte-Carlo方法,一、MC的起源和发展Buffon试验排队系统模拟：M/G/1排队系统二、MC的原理三、随机数的产生原理与IMSL库均匀分布U(0,1)的随机数的产生其他各种分布的随机数的产生随机过程模拟四、MC的应用举例定积分的MC计算随机微分方程的数值模拟五、EM算法及其MCMC方法,一、MC的起源和发展,随机模拟方法，也称为MonteCarlo方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战进行的研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯诺伊曼用驰名世界的赌城摩纳哥的MonteCarlo来命名这种方法，为它蒙上了一层神秘色彩。冯诺伊曼是公理化方法和计算机体系的领袖人物，MonteCarlo方法也是他的功劳。,事实上，MonteCarlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。18世纪下半叶的法国学者Buffon提出用投针试验的方法来确定圆周率的值。这个著名的Buffon试验是MonteCarlo方法的最早的尝试！,历史上曾有几位学者相继做过这样的试验。不过呢，他们的试验是费时费力的，同时精度不够高，实施起来也很困难。然而，随着计算机技术的飞速发展，人们不需要具体实施这些试验，而只要在计算机上进行大量的、快速的模拟试验就可以了。,在大众的心目中，科学的代言人是心不在焉的牛顿或者爆炸式发型的爱因斯坦，但这只是传统形象，比他们更了解现代计算技术的冯诺伊曼是个衣着考究，风度翩翩的人物，他说：纯粹数学和应用数学的许多分支非常需要计算工具，用以打破目前由于纯粹分析的研究方法不能解决非线性问题而形成的停滞状态。MonteCarlo方法是现代计算技术的最为杰出的成果之一，它在工程领域的作用是不可比拟的。,Buffon试验,假设平面上有无数条距离为1的等距平行线，现向该平面随机投掷一根长度为的针()，则我们可计算该针与任一平行线相交的概率。这里，随机投针指的是：针的中心点与最近的平行线间的距离均匀的分布在区间。上，针与平行线的夹角（不管相交与否）均匀的分布在区间上。因此，针与线相交的充要条件是,Buffon试验,从而针线相交的概率为根据上式，若我们做大量的投针试验并记录针与线相交的次数，则由大数定理可以估计出针线相交的概率，从而得到的估计值。针与线的位置关系：,排队系统模拟：M/G/1排队系统,服务规则：先到先服务假设：（）顾客到达遵循Poisson分布；（）服务时间服从一般分布；（）到达间隔与服务时间相互独立,1,排队系统模拟：M/G/1排队系统,关心的指标：（）时刻t时，系统中的顾客数；即队长的分布；（）顾客的等待时间；（）服务的忙碌程度；（）用最朴素的Monte-Carlo方法可以得到这些指标的估计,二、MC的原理,应用MonteCarlo方法求解工程技术问题可以分为两类：确定性问题随机性问题,确定性系统,随机性系统,模拟,自然界,Monte-Carlo模拟，即随机模拟（重复“试验”）,重复试验,计算机模拟,思路：,1、针对实际问题建立一个简单且便于实现的概率统计模型，使问题的解对应于该模型中随机变量的概率分布或其某些数字特征，比如，均值和方差等。所构造的模型在主要特征参量方面要与实际问题或系统相一致的。2、根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，再进行随机模拟试验。,思路：,3、根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）。4、按照所建立的模型进行仿真试验、计算，求出问题的随机解。5、统计分析模拟试验结果，给出问题的估计以及其精度估计。6、必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。,收敛性:由大数定律,Monte-Carlo模拟的收敛是以概率而言的误差:用频率估计概率时误差的估计，可由中心极限定理，给定置信水平的条件下，有：模拟次数:由误差公式得,三、随机数的产生原理与IMSL库,随机数的产生是实现MC计算的先决条件。而大多数概率分布的随机数的产生都是基于均匀分布U(0,1)的随机数。首先，介绍服从均匀分布U(0,1)的随机数的产生方法。其次，介绍服从其他各种分布的随机数的产生方法。以及服从正态分布的随机数的产生方法。接下来，介绍随机过程模拟。最后，关于随机数的几点注。,均匀分布U(0,1)的随机数的产生,产生均匀分布的标准算法在很多高级计算机语言的书都可以看到。算法简单，容易实现。使用者可以自己手动编程实现。IMSL统计库中也提供给我们用于产生均匀分布的各种函数。我们的重点是怎样通过均匀分布产生服从其他分布的随机数。因此，直接使用IMSL提供的可靠安全的标准函数，当然不用费事了。,IMSL库中的函数使用,RNSET:种子的设定CALLRNSET(ISEED)RNOPT:产生器的类型的设定CALLRNOPT(IOPT)RNUN/DRNUN:产生均匀分布的随机数CALLRNUN(NR,R),其他各种分布的随机数的产生,基本方法有如下三种：逆变换法合成法筛选法,逆变换法,设随机变量的分布函数为，定义定理设随机变量服从上的均匀分布，则的分布函数为。因此，要产生来自的随机数，只要先产生来自的随机数，然后计算即可。其步骤为,合成法,合成法的应用最早见于Butlter的书中。构思如下：如果的密度函数难于抽样，而关于的条件密度函数以及的密度函数均易于抽样，则的随机数可如下产生：可以证明由此得到的服从。,筛选抽样,假设我们要从抽样，如果可以将表示成，其中是一个密度函数且易于抽样，而，是常数，则的抽样可如下进行：定理设的密度函数，且，其中，是一个密度函数。令和分别服从和，则在的条件下，的条件密度为,标准正态分布的随机数,的随机数产生方法很多。简要介绍三种。法1、变换法（Box和Muller1958）设，是独立同分布的变量，令则与独立，均服从标准正态分布。法2、结合合成法与筛选法。（略）法3、近似方法（利用中心极限定理）即用个变量产生一个变量。其中是抽自的随机数，可近似为一个变量。,IMSL库中的函数使用,RNSET:种子的设定CALLRNSET(ISEED)RNOPT:产生器的类型的设定CALLRNOPT(IOPT)RNNOA:产生服从标准正态分布的伪随机数CALLRNNOA(NR,R),随机过程模拟,高斯白噪声的产生：一种简单有效的模拟高斯白噪声过程的方法是由独立的单位正态随机序列，按照如图所示的方式连接，其中式中为白噪声的强度。,由于正态随机数的独立性，当很小时，按(a)、(b)所构成的过程的相关时间很短，从而具有很高的截止频率。当然，过小，样本的计算量过大。因此，选取适当小即可。,关于随机数的几点注,注1由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注2应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。,四、MC的应用举例,例一、定积分的MC计算随机投点法样本平均值法几种降低估计方差的MC方法例二、随机微分方程的数值模拟系统的可靠性数值模拟计算问题,定积分的MC计算,事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。下面考虑一个简单的定积分为了说明问题，我们首先介绍两种求的简单的MC方法，然后给出几种较为复杂而更有效的MC方法。,在计算积分上，MC的实用场合是计算重积分其中是维空间的点，当较大时，用MC方法比一般的数值方法有优点，主要是它的误差与维数无关。,随机投点法,方法简述：设，有限，并设是在上均匀分布的二维随机变量，其联合密度函数为。则易见是中曲线下方的面积。假设我们向中进行随机投点，若点落在下方，(即称为中的，否则称为不中，则点中的概率为。若我们进行了次投点，其中次中的，则用频率来估计概率。即。,那么我们可以得到的一个估计具体试验步骤为,注1随机投点法的思想简单明了，且每次投点结果服从二项分布，故，其中注2可证是的无偏估计。若用估计的标准差来衡量其精度，则估计的精度的阶为。注3这里，定积分的解，就对应我们选定的随机变量的概率值。,样本平均值法,基本原理：对积分，设是上的一个密度函数，改写可见，任一积分均可以表示为某个随机变量（函数）的期望。由矩法，若有个来自的观测值，则可给出的一个矩估计。最简单的，若，有限，可取。,设是来自的随机数，则的一个估计为具体步骤为注可证是的无偏估计。一般而言，样本均值法要比随机投点法更有效。,求解定积分的算例,计算定积分事实上，其精确解为随机投点法：sjtd_djf.f90，及动态图样本平均值法：ybpjz_djf.f90注所有计算中的随机数的产生均来自IMSL库。样本数为10000。增加样本数目，可提高计算精度，但计算时间也会提高。,几种降低估计方