《第讲蒙特卡洛方法初探.ppt》由会员分享,可在线阅读,更多相关《第讲蒙特卡洛方法初探.ppt(63页珍藏版)》请在三一文库上搜索。
1、蒙特卡洛方法初探,实验目的,实验内容,学习计算机模拟的基本过程与方法。,1、模拟的概念。,4、实验作业。,3、计算机模拟实例。,2、产生随机数的计算机命令。,模拟的概念,模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程,以寻求过程规律的一种方法。,模拟的基本思想是建立一个试验模型,这个模型包含所研究系统的主要特点通过对这个实验模型的运行,获得所要研究系统的必要信息,模拟的方法,1、物理模拟: 对实际系统及其过程用功能相似的实物系统去模仿。例如,军事演习、船艇实验、沙盘作业等。,物理模拟通常花费较大、周期较长,且在物理模型上改变系统结构和系数都较困难。而且,许多系统无法进行物理模拟
2、,如社会经济系统、生态系统等。,在实际问题中,面对一些带随机因素的复杂系统,用分析方法建模常常需要作许多简化假设,与面临的实际问题可能相差甚远,以致解答根本无法应用。这时,计算机模拟几乎成为唯一的选择。,在一定的假设条件下,运用数学运算模拟系统的运行,称为数学模拟。现代的数学模拟都是在计算机上进行的,称为计算机模拟。,2、数学模拟,计算机模拟可以反复进行,改变系统的结构和系数都比较容易。,蒙特卡洛(Monte Carlo)方法是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样,通过对样本值的统计分析,求得所研究系统的某些参数,用蒙特卡洛方法进行计算机模拟的步骤:,1 设计
3、一个逻辑框图,即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。 2 模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象,产生模拟随机数的计算机命令,在Matlab软件中,可以直接产生满足各种分布的随机数,命令如下:,1产生m*n阶(a,b)均匀分布U(a,b)的随机数矩阵: unifrnd (a,b,m, n); 产生一个a,b均匀分布的随机数:unifrnd (a,b),当只知道一个随机变量取值在(a,b)内,但不知道(也没理由假设)它在何处取值的概率大,在何处取值的概率小,就只好用U(a,b)来模拟它。,2产生mm*nn阶离散均匀分布的随机数矩阵: R = unidrnd
4、(N) R = unidrnd(N,mm,nn),当研究对象视为大量相互独立的随机变量之和,且其中每一种变量对总和的影响都很小时,可以认为该对象服从正态分布。,若连续型随机变量X的概率密度函数为 其中 0为常数,则称X服从参数为 的指数分布。,指数分布的期望值为,排队服务系统中顾客到达间隔、质量与可靠性中电子元件的寿命通常服从指数分布。,例 顾客到达某商店的间隔时间服从参数为10(分钟)的指数分布(指数分布的均值为10),- 指两个顾客到达商店的平均间隔时间是10分钟.即平均10分钟到达1个顾客. 顾客到达的间隔时间可用exprnd(10)模拟。,设离散型随机变量X的所有可能取值为0,1,2,
5、且取各个值的概率为 其中 0为常数,则称X服从参数为 的帕松分布。,帕松分布在排队系统、产品检验、天文、物理等领域有广泛应用。,帕松分布的期望值为,6 产生1个参数为n,p的二项分布的随机数binornd(n,p),产生mn个参数为n,p的二项分布的随机数binornd(n,p,m,n) 。,掷一枚均匀硬币,正面朝上的次数X服从参数为,p的二项分布,XB(1,p),频率的稳定性模拟,1.事件的频率 在一组不变的条件下,重复作n次试验,记m是n次试验中事件A发生的次数。 频率 f=m/n,2.频率的稳定性,掷一枚均匀硬币,记录掷硬币试验中频率P*的波动情况。,function liti1(p,m
6、m) pro=zeros(1,mm); randnum = binornd(1,p,1,mm) a=0; for i=1:mm a=a+randnum(1,i); pro(i)=a/i; end pro=pro num=1:mm; plot(num,pro),在Matlab中编辑.m文件输入以下命令:,在Matlab命令行中输入以下命令:,liti1(0.5,1000),在Matlab命令行中输入以下命令:,liti1(0.5,10000),练习掷一枚不均匀硬币,正面出现概率为0.3,记录前1000次掷硬币试验中正面频率的波动情况,并画图。,在Matlab命令行中输入以下命令:,liti1(0
7、.3,1000),例2 掷两枚不均匀硬币,每枚正面出现概率为0.4,记录前1000次掷硬币试验中两枚都为正面频率的波动情况,并画图。,在Matlab中编辑.m文件输入以下命令:,function liti2(p,mm) pro=zeros(1,mm); randnum = binornd(1,p,2,mm);a=0; for i=1:mm a=a+randnum(1,i)*randnum(2,i); pro(i)=a/i; end pro=pro,num=1:mm;plot(num,pro),liti2(0.4,100),liti2(0.4,10000),例3: 在一袋中有10 个相同的球,分
8、别标有号码1,2,10。每次任取一个球,记录其号码后放回袋中,再任取下一个。这种取法叫做“有放回抽取”。今有放回抽取3个球,求这3个球的号码均为偶数的概率。(用频率估计概率),解:有放回取3个球, 所有取法有103种; 有放回取3个偶数号码的球, 所有取法有53种. 所以,古典概率模拟,function proguji=liti3(n,mm) frq=0; randnum=unidrnd(n,mm,3);proguji=0; for i=1:mm a=(randnum(i,1)+1)*(randnum(i,2)+1)*(randnum(i,3)+1); if mod(a,2)=1 frq=fr
9、q+1 end end; proguji=frq/mm,例4 两盒火柴,每盒20根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有5根火柴的概率有多大?(用频率估计概率),function proguji=liti40(mm) %mm 是随机实验次数 frq=0; randnum=binornd(1,0.5,mm,2*20);proguji=0; for i=1:mm a1=0;a2=0;j=1; while (a1=5 frq=frq+1; end end; proguji=frq/mm, liti40(100) proguji = 0.4800 liti40(10
10、00) proguji = 0.4970 liti40(10000) proguji = 0.4910 liti40(100000) proguji = 0.4984,例4 两盒火柴,每盒n根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有k根火柴的概率有多大?(用频率估计概率),function proguji=liti4(n,k,mm) %n是每盒中的火柴数 %k 是剩余的火柴数 %mm 是随机实验次数 frq=0; randnum=binornd(1,0.5,mm,2*n);proguji=0; for i=1:mm a1=0;a2=0;j=1; while (
11、a1=k , frq=frq+1; end % a1=a1,a2=a2,frq % pause end; proguji=frq/mm, liti4(20,5,100) proguji = 0.4800 liti4(20,5,1000) proguji = 0.4970 liti4(20,5,10000) proguji = 0.4910 liti4(20,5,100000) proguji = 0.4984,几何概率模拟,1.定义,向任一可度量区域G内投一点,如果所投的点落在G中任意可度量区域g内的可能性与g的度量成正比,而与g的位置和形状无关,则称这个随机试验为几何型随机试验。或简称为几何
12、概型。,2. 概率计算,P(A)=A的度量/S的度量,例5 两人约定于12点到1点到某地会面,先到者等20分钟后离去,试求两人能会面的概率?,解:设x, y分别为甲、乙到达时刻(分钟),令A=两人能会面=(x,y)|x-y|20,x60,y60,P(A)=A的面积/S的面积 =(602-402)/602=5/9=0.5556,function proguji=liti5(mm) %mm 是随机实验次数 frq=0; randnum1=unifrnd(0,60,mm,1); randnum2=unifrnd(0,60,mm,1); randnum=randnum1-randnum2; progu
13、ji=0; for ii=1:mm if abs(randnum(ii,1)=20 frq=frq+1; end end proguji=frq/mm,liti5(10000) proguji = 0.5557,例6 在我方某前沿防守地域,敌人以一个炮排(含两门火炮)为单位对我方进行干扰和破坏为躲避我方打击,敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现,我方指挥所对敌方目标的指示有50是准确的,而我方火力单位,在指示正确时,有1/3的概率能毁伤敌人一门火炮,有1/6的概率能全部消灭敌人,现在希望能用某种方式把我方将要对敌人实施的1次打击结果显现出来,利用频率稳定性,确定有效射击(
14、毁伤一门炮或全部消灭)的概率.,复杂概率模拟,分析: 这是一个复杂概率问题,可以通过理论计算得到相应的概率.,为了直观地显示我方射击的过程,现采用模拟的方式。,需要模拟出以下两件事:,1. 问题分析,1 观察所对目标的指示正确与否,模拟试验有两种结果,每一种结果出现的概率都是1/2,因此,可用投掷一枚硬币的方式予以确定,当硬币出现正面时为指示正确,反之为不正确,2 当指示正确时,我方火力单位的射击结果情况,模拟试验有三种结果:毁伤一门火炮的可能性为1/3(即2/6),毁伤两门的可能性为1/6,没能毁伤敌火炮的可能性为1/2(即3/6),这时可用投掷骰子的方法来确定: 如果出现的是、三个点:则认
15、为没能击中敌人; 如果出现的是、点:则认为毁伤敌人一门火炮; 若出现的是点:则认为毁伤敌人两门火炮,2. 符号假设,i:要模拟的打击次数; k1:没击中敌人火炮的射击总数; k2:击中敌人一门火炮的射击总数; k3:击中敌人两门火炮的射击总数 E:有效射击(毁伤一门炮或两门炮)的概率,3. 模拟框图,function liti6(p,mm) efreq=zeros(1,mm); randnum1 = binornd(1,p,1,mm); randnum2 = unidrnd(6,1,mm);k1=0;k2=0;k3=0; for i=1:mm if randnum1(i)=0 k1=k1+1;
16、 else if randnum2(i)=3 k1=k1+1; elseif randnum2(i)=6 k3=k3+1; else k2=k2+1; end end efreq(i)=(k2+k3)/i; end num=1:mm;plot(num,efreq),在Matlab中编辑.m文件输入以下命令:,在Matlab命令行中输入以下命令:,liti6(0.5,2000),在Matlab命令行中输入以下命令:,liti6(0.5,20000),5. 理论计算,6. 结果比较,模拟结果与理论计算近似一致,能更加真实地表达实际战斗动态过程,三.蒙特卡洛模拟的理论基础与模拟结果的误差,大数定律,
17、中心极限定理,大数定律,贝努里(Bernoulli) 大数定律,设 nA 是 n 次独立重复试验中事件 A 发生的 次数, p 是每次试验中 A 发生的概率,则,有,或,贝努里(Bernoulli) 大数定律的意义:,设 X1,X2,XN是来自总体X (EX) 的简单随机样本,即X1,X2,XN 独立同分 布,则,辛钦 大数定律,即,中心极限定理,则对于任意实数 x ,若令,则,等价于,于是,这表明,不等式,近似地以概率1成立。上式也表明,,收敛到,的阶为O(n -1/2)。通常,蒙特卡罗方法的误差 定义为,关于蒙特卡罗方法的误差需说明两点: (1),蒙特卡罗方法的误差为概率误差, 也即蒙特卡
18、罗方法的收敛是概率意义下的 收敛,虽然不能断言其误差不超过某个 值,但能指出其误差以接近1的概率不超过 某个界限。,如=0.5,误差,此时,误差,超过的概率与小于的概率1-相等,都 等于0.5。,来代替,在计算所求量的同时,可计算出 。,关于蒙特卡罗方法的误差需说明两点: (2)误差中的均方差 是未知的,必须 使用其估计值,显然,当给定置信度后,误差由 和 n决定。要减小 ,或者是增大n ,或者是减 小方差 2。在 固定的情况下,要把精度提 高一个数量级,试验次数n需增加两个数量 级。因此,单纯增大n不是一个有效的办法。,四、减小方差的各种技巧,另一方面,如能减小估计的均方差 , 比如降低一半
19、,那误差就减小一半,这相 当于n增大四倍的效果(n=(u /)2)。因此 降低方差的各种技巧,引起了人们的普遍 注意。,一般来说,降低方差的技巧,往往会使观 察一个子样的时间增加。在固定时间内,使 观察的样本数减少。所以,一种方法的优 劣,需要由方差和观察一个子样的费用(使 用计算机的时间)两者来衡量。这就是蒙特 卡罗方法中效率的概念。它定义为nc ,其中c 是观察一个子样的平均费用。显然 nc= (u /)2 2c 它与 2c成正比。,四、效率,总而言之,作为提高蒙特卡洛方法效率的 重要方向,是在减小标准差的同时兼顾考虑费 用大小,使 2c尽可能地小。,nc= (u /)2 2c),五、蒙特
20、卡洛方法的特点,Monte Carlo方法及其程序结构简单, 产生随机数,计算相应的值。即通过大 量的简单重复抽样和简单计算实现该方法。,2)收敛速度与问题维数无关, Monte Carlo方法的收敛速度为O(n -1/2), 与一般数值方法相比很慢。因此,用Monte Carlo方法不能解决精确度要求很高的问题, 从,可见,Monte Carlo方法的误差只与标准差 和样本容量n有关,而与样本所在空间无关, 即Monte Carlo方法的收敛速度与问题维数无 关。而其他数值方法则不然。因此,这就决定 了Monte Carlo方法对多维问题的适用性。, 在解题时受问题条线限制的影响较小,例 如
21、要计算s维空间中的任一区域Ds上的积分,3) Monte Carlo方法的适用性强,时,无论区域Ds的形状多么特殊,只要能给出 描述Ds的几何特征的条件,就可以从Ds中均匀 产生n个点 , 得到积分的近似值,其中Ds为区域Ds的体积。这是数值方法难以作 到的。,3. 某厂生产的灯泡能用1000小时的概率为0.8, 能用1500小时的概率为0.4 , 求已用1000小时的灯泡能用到1500小时的概率(频率估计概率)。,2. 在一袋中有10 个相同的球,分别标有号码1,2, 10。今有放回任取两个球,求取得的第一个球号码 为奇数,第二个球的号码为偶数的概率(频率估计概率),掷三枚不均匀硬币,每枚正面出现概率为0.3,记 录前1000次掷硬币试验中至少两枚都为正面频率的 波动情况,并画图。,作业:,4 : 两船欲停靠同一个码头, 设两船到达码头的时间各不相干,而且到达码头的时间在一昼夜内是等可能的. 如果两船到达码头后需在码头停留的时间分别是1 小时与2 小 时,试求在一昼夜内,任一船到达时,需要等待空出码头的概率. (频率估计概率),5 :在0,1,2,3,9中不重复地任取4个数,求它们能排成首位非零的四位偶数的概率. (频率估计概率),6:从1,2, ,10十个数字中有放回地任取5个数字, 求取出的5个数字中按由小到大排列, 中间的那个数等于 4 的概率. (频率估计概率),
链接地址:https://www.31doc.com/p-2574011.html