数学建模第3讲蒙特卡洛方法初探(2007)

1、蒙特卡洛方法初探实验目的实验目的实验内容实验内容学习计算机模拟的基本过程与方法。学习计算机模拟的基本过程与方法。1 1、模拟的概念。、模拟的概念。4 4、实验作业。、实验作业。3 3、计算机模拟实例。、计算机模拟实例。2 2、产生随机数的计算机命令。、产生随机数的计算机命令。模拟的概念模拟的概念 模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。 模拟的基本思想是建立一个试验模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息模拟的方法模拟的方法1、物理模拟物理模拟： 对实际系统及其过程用功能相似的实物系统去模仿。例

2、如，军事演习、船艇实验、沙盘作业等。 物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难。而且，许多系统无法进行物理模拟，如社会经济系统、生态系统等。 在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。 在一定的假设条件下，运用数学运算模拟系统在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟。计算机上进行的，称为计算机模拟。2、数学模拟数学模拟 计算

3、机模拟可以反复进行，改变系统的结构和系数都比较容易。 蒙特卡洛（蒙特卡洛（Monte CarloMonte Carlo）方法）方法是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样，通过对样本值的统计分析，求得所研究系统的某些参数用蒙特卡洛方法进行计算机模拟的步骤用蒙特卡洛方法进行计算机模拟的步骤：1 设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。2 模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象产生模拟随机数的计算机命令产生模拟随机数的计算机命令 在Matlab软件中，可以直接产生满足各种分布的随机数，命令如下：1产生m*n阶(

4、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(N)R = unidrnd(N,mm,nn)当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。若连续型随机变量X的概率密度函数为 其中 0为常数，则称X服从参数为 的指数指数分布分布。 0001

5、)(/xxexfx 指数分布的期望值为 排队服务系统中顾客到达间隔、质量与可靠性中电子元件的寿命通常服从指数分布。例例 顾客到达某商店的间隔时间服从参数为顾客到达某商店的间隔时间服从参数为10(10(分钟分钟) )的指数分布的指数分布( (指数分布的均值为指数分布的均值为10) - - 指两个顾客到达商店的平均间隔时间指两个顾客到达商店的平均间隔时间是是1010分钟分钟. .即平均即平均1010分钟到达分钟到达1 1个顾客个顾客. . 顾客顾客到达的间隔时间可用到达的间隔时间可用exprnd(10)exprnd(10)模拟。模拟。设离散型随机变量X的所有可能取值为0,1,2,且取各个值的概率为

6、其中 0为常数，则称X服从参数为 的帕松分布帕松分布。, 2 , 1 , 0,!)(kkekXPk帕松分布在排队系统、产品检验、天文、物理等领域有广泛应用。帕松分布的期望值为6产生1个参数为n,p的二项分布的随机数binornd(n,p)，产生mn个参数为n,p的二项分布的随机数binornd(n,p,m,n) 。 掷一枚均匀硬币，正面朝上的次数掷一枚均匀硬币，正面朝上的次数X X服从参数为，服从参数为，p p的二项分布的二项分布,XB(1,p),XB(1,p)频率的稳定性模拟频率的稳定性模拟 1. 1.事件的频率事件的频率 在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，次

7、试验，记记m m是是n n次试验中事件次试验中事件A A发生的次数。发生的次数。 频率频率 f=m/n f=m/n 2.2.频率的稳定性频率的稳定性 掷一枚均匀硬币，记录掷硬币试验中掷一枚均匀硬币，记录掷硬币试验中频率频率P P* *的波动情况。的波动情况。 function liti1(p,mm)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=pronum=1:mm;plot(num,pro)在在MatlabMatlab中编辑中编辑.m.m文件输入以下

8、命令：文件输入以下命令：在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(0.5,1000)在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(0.5,10000)练习练习掷一枚不均匀硬币，正面出现概率掷一枚不均匀硬币，正面出现概率为为0.30.3，记录前，记录前10001000次掷硬币试验中正面次掷硬币试验中正面频率的波动情况，并画图。频率的波动情况，并画图。 在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(0.3,1000)例例2 2 掷两枚不均匀硬币，每枚正面出现概率掷两枚不均匀硬币，每

9、枚正面出现概率为为0.40.4，记录前，记录前10001000次掷硬币试验中两枚都为次掷硬币试验中两枚都为正面频率的波动情况，并画图。正面频率的波动情况，并画图。 在在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,1000

10、0)例例3: 在一袋中有在一袋中有10 10 个相同的球，分别标有号码个相同的球，分别标有号码1,2,1,2,10,10。每次任取一个球，记录其号码后放。每次任取一个球，记录其号码后放回袋中，再任取下一个。这种取法叫做回袋中，再任取下一个。这种取法叫做“有放有放回抽取回抽取”。今有放回抽取。今有放回抽取3 3个球，求这个球，求这3 3个球的个球的号码均为偶数的概率。号码均为偶数的概率。（用频率估计概率）（用频率估计概率） 解：有放回取解：有放回取3个球个球, 所有取法有所有取法有103种种; 有放回有放回取取3个偶数号码的球个偶数号码的球, 所有取法有所有取法有53种种. 所以所以81105)

11、(33AP古典概率模拟古典概率模拟 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=frq+1 endend; proguji=frq/mm例例4 4 两盒火柴，每盒两盒火柴，每盒2020根。每次随机在任根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有被取完而另一盒中至少还有5 5

12、根火柴的概根火柴的概率有多大？（用频率估计概率）率有多大？（用频率估计概率） 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 (a120)&(a2=5 frq=frq+1; endend; proguji=frq/mm liti40(100)proguji = 0.4800 liti40(1000)proguji = 0.4970 liti40(10000)proguji = 0.4910 liti4

13、0(100000)proguji = 0.4984例例4 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; whi

14、le (a1n)&(a2=k , frq=frq+1; end% a1=a1,a2=a2,frq % pauseend; 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.1.定义定义 向任一可度量区域向任一可度量区域G G内投一点，如果所投的点内投一点，如果所投的点落在落在G G中任意可度量区域中任意可度量区域g

15、 g内的可能性与内的可能性与g g的度的度量成正比，而与量成正比，而与g g的位置和形状无关，则称这的位置和形状无关，则称这个随机试验为几何型随机试验。或简称为几何个随机试验为几何型随机试验。或简称为几何概型。概型。2. 概率计算概率计算 P(A)=A的度量的度量/S的度量的度量例例5 两人约定于两人约定于12点到点到1点到某地会面，先到者点到某地会面，先到者等等20分钟后离去，试求两人能会面的概率？分钟后离去，试求两人能会面的概率？ 解：设解：设x, y分别为甲、乙到达时刻分别为甲、乙到达时刻(分钟分钟)令令A=两人能会面两人能会面=(x,y)|x-y|20，x60,60,y y6060P(

16、A)=A的面积的面积/S的面积的面积=（602-402）/602=5/9=0.5556function proguji=liti5(mm)%mm 是随机实验次数frq=0;randnum1=unifrnd(0,60,mm,1);randnum2=unifrnd(0,60,mm,1);randnum=randnum1-randnum2;proguji=0;for ii=1:mm if abs(randnum(ii,1)=20 frq=frq+1; endendproguji=frq/mmliti5(10000)proguji = 0.5557例例6 6在我方某前沿防守地域，敌人以一个炮排（含两门

17、火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点 经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的概率能毁伤敌人一门火炮，有1/6的概率能全部消灭敌人 现在希望能用某种方式把我方将要对敌人实施的1次打击结果显现出来，利用频率稳定性，确定有效射击(毁伤一门炮或全部消灭)的概率.复杂概率模拟复杂概率模拟分析分析： 这是一个复杂概率问题，可以通过理论计算得到相应的概率. 为了直观地显示我方射击的过程，现采用模拟的方式。 需要模拟出以下两件事： 1. 1. 问题分析问题分析1 1 观察所对目标的指示正确与否观察

18、所对目标的指示正确与否 模拟试验有两种结果，每一种结果出现的概率都是1/2 因此，可用投掷一枚硬币的方式予以确定可用投掷一枚硬币的方式予以确定,当硬币出现正面时为指示正确，反之为不正确 2 2 当指示正确时，我方火力单位的射击结果当指示正确时，我方火力单位的射击结果情况情况 模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6) 这时可用投掷骰子的方法来确定可用投掷骰子的方法来确定：如果出现的是、三个点：则认为没能击中敌人；如果出现的是、点：则认为毁伤敌人一门火炮；若出现的是点：则认为毁伤敌人两门火炮2. 2. 符号假

19、设符号假设i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数E：有效射击(毁伤一门炮或两门炮)的概率3. 3. 模拟框图模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1imm?E=(k2+k3)/mm 停止硬币正面?YNNY1，2，34，56function liti6(p,mm)efreq=zeros(1,mm);randnum1 = binornd(1,p,1,mm);randnum2 = unidrnd(6,1,mm);k1=0;k2=

20、0;k3=0;for i=1:mm if randnum1(i)=0 k1=k1+1; 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)在在MatlabMatlab中编辑中编辑.m.m文件输入以下命令：文件输入以下命令：在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti6(0.5,2000)在在MatlabMatlab命令行中输入以下命令：命令行中输入以

21、下命令：liti6(0.5,20000)5. 5. 理论计算理论计算6. 6. 结果比较结果比较 模拟结果与理论计算近似一致，能更加真模拟结果与理论计算近似一致，能更加真实地表达实际战斗动态过程实地表达实际战斗动态过程 三三. .蒙特卡洛模拟的理论基础与模拟结果的误差蒙特卡洛模拟的理论基础与模拟结果的误差大数定律大数定律中心极限定理中心极限定理大数定律大数定律贝努里（Bernoulli） 大数定律设 nA 是 n 次独立重复试验中事件 A 发生的次数, p 是每次试验中 A 发生的概率，则0有0limpnnPAn或1limpnnPAn在概率的统计定义中，事件 A 发生的频率“ 稳定于”事件 A

22、 在一次试验中发生的概率是指：nnAnnA频率与 p 有较大偏差pnnA是小概率事件, 因而在 n 足够大时, 可以用频率近似代替 p . 这种稳定称为依概率稳定.贝努里（贝努里（Bernoulli） 大数定律的意义大数定律的意义： 设 X1，X2，XN是来自总体X (EX)的简单随机样本，即X1，X2，XN 独立同分布，则 EXXNXPNiiN111)(limXEXPNN辛钦 大数定律即中心极限定理中心极限定理 设随机变量序列,21nXXX相互独立，服从同一分布，且有期望和方差：, 2 , 1,0)(,)(2kXDXEkk则对于任意实数 x ,xtnkkndtexnnXP21221lim若令

23、) 1 , 0(1NnnXnkknkknXnX11则等价于) 1 , 0(/NnXn于是1)(unuXPN这表明，不等式近似地以概率1成立。上式也表明，nuXnnX收敛到的阶为O(n -1/2)。通常，蒙特卡罗方法的误差定义为nu 关于蒙特卡罗方法的误差需说明两点：（1），蒙特卡罗方法的误差为概率误差，也即蒙特卡罗方法的收敛是概率意义下的收敛，虽然不能断言其误差不超过某个值，但能指出其误差以接近1的概率不超过某个界限。如=0.5，误差./6745. 0n此时，误差超过的概率与小于的概率1-相等，都等于0.5。 来代替，在计算所求量的同时，可计算出 。2112)1(1niiniiXnXn 关于蒙

24、特卡罗方法的误差需说明两点：（2）误差中的均方差 是未知的，必须使用其估计值 显然，当给定置信度后，误差由 和n决定。要减小 ，或者是增大n ，或者是减小方差 2。在 固定的情况下，要把精度提高一个数量级，试验次数n需增加两个数量级。因此，单纯增大n不是一个有效的办法。四、减小方差的各种技巧四、减小方差的各种技巧 nu 另一方面，如能减小估计的均方差 ，比如降低一半，那误差就减小一半，这相当于n增大四倍的效果(n=(u /)2)。因此降低方差的各种技巧，引起了人们的普遍注意。nu 一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，

25、需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为nc ，其中c是观察一个子样的平均费用。显然 nc= (u /)2 2c它与 2c成正比。四、效率四、效率 总而言之，作为提高蒙特卡洛方法效率的重要方向，是在减小标准差的同时兼顾考虑费用大小，使 2c尽可能地小。nc= (u /)2 2c)五、蒙特卡洛方法的特点五、蒙特卡洛方法的特点1) Monte Carlo方法及其程序结构简单 产生随机数，计算相应的值。即通过大量的简单重复抽样和简单计算实现该方法。2）收敛速度与问题维数无关）收敛速度与问题维数无关 Monte Carlo方法的收敛速度为O

26、(n -1/2)，与一般数值方法相比很慢。因此，用Monte Carlo方法不能解决精确度要求很高的问题 从nu可见，Monte Carlo方法的误差只与标准差和样本容量n有关，而与样本所在空间无关，即Monte Carlo方法的收敛速度与问题维数无关。而其他数值方法则不然。因此，这就决定了Monte Carlo方法对多维问题的适用性。 在解题时受问题条线限制的影响较小，例如要计算s维空间中的任一区域Ds上的积分ssDdxdxdxxxxggs2121),( 3） Monte Carlo方法的适用性强时，无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生n个点 , 得到积分的近似值),()()(2)(1isii