第3讲蒙特卡洛方法初探(2007)

蒙特卡洛方法初探实验目的实验内容学习计算机模拟的基本过程与方法。1、模拟的概念。4、实验作业。3、计算机模拟实例。2、产生随机数的计算机命令。模拟的概念模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。模拟的基本思想是建立一个试验模型，这个模型包含所研究系统的主要特点．通过对这个实验模型的运行，获得所要研究系统的必要信息模拟的方法1、物理模拟：对实际系统及其过程用功能相似的实物系统去模仿。例如，军事演习、船艇实验、沙盘作业等。物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难。而且，许多系统无法进行物理模拟，如社会经济系统、生态系统等。在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟。2、数学模拟计算机模拟可以反复进行，改变系统的结构和系数都比较容易。蒙特卡洛（MonteCarlo）方法是一种应用随机数来进行计算机模拟的方法．此方法对研究的系统进行随机观察抽样，通过对样本值的统计分析，求得所研究系统的某些参数．用蒙特卡洛方法进行计算机模拟的步骤：[1]设计一个逻辑框图，即模拟模型．这个框图要正确反映系统各部分运行时的逻辑关系。[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(N)R=unidrnd(N,mm,nn)3.产生mn阶均值为，标准差为的正态分布的随机数矩阵：normrnd(,,m,n)产生一个均值为，标准差的正态分布的随机数：normrnd(,)•当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。4．产生mn阶期望值为的指数分布的随机数矩阵：exprnd(,m,n)•若连续型随机变量X的概率密度函数为其中0为常数，则称X服从参数为的指数分布。0001)(/xxexfx•指数分布的期望值为•排队服务系统中顾客到达间隔、质量与可靠性中电子元件的寿命通常服从指数分布。例顾客到达某商店的间隔时间服从参数为10(分钟)的指数分布(指数分布的均值为10)-----指两个顾客到达商店的平均间隔时间是10分钟.即平均10分钟到达1个顾客.顾客到达的间隔时间可用exprnd(10)模拟。•设离散型随机变量X的所有可能取值为0,1,2,…,且取各个值的概率为其中0为常数，则称X服从参数为的帕松分布。,,2,1,0,!)(kkekXPk5．产生mn阶参数为的帕松分布的随机数矩阵：poissrnd(,m,n)•帕松分布在排队系统、产品检验、天文、物理等领域有广泛应用。•帕松分布的期望值为6产生1个参数为n,p的二项分布的随机数binornd(n,p)，产生mn个参数为n,p的二项分布的随机数binornd(n,p,m,n)。掷一枚均匀硬币，正面朝上的次数X服从参数为１，p的二项分布,X~B(1,p)频率的稳定性模拟1.事件的频率在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数。频率f=m/n2.频率的稳定性掷一枚均匀硬币，记录掷硬币试验中频率P*的波动情况。functionliti1(p,mm)pro=zeros(1,mm);randnum=binornd(1,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(0.5,1000)在Matlab命令行中输入以下命令：liti1(0.5,10000)练习掷一枚不均匀硬币，正面出现概率为0.3，记录前1000次掷硬币试验中正面频率的波动情况，并画图。在Matlab命令行中输入以下命令：liti1(0.3,1000)例2掷两枚不均匀硬币，每枚正面出现概率为0.4，记录前1000次掷硬币试验中两枚都为正面频率的波动情况，并画图。在Matlab中编辑.m文件输入以下命令：functionliti2(p,mm)pro=zeros(1,mm);randnum=binornd(1,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)liti2(0.4,100)liti2(0.4,10000)例3:在一袋中有10个相同的球，分别标有号码1,2,…,10。每次任取一个球，记录其号码后放回袋中，再任取下一个。这种取法叫做“有放回抽取”。今有放回抽取3个球，求这3个球的号码均为偶数的概率。（用频率估计概率）解：有放回取3个球,所有取法有103种;有放回取3个偶数号码的球,所有取法有53种.所以81105)(33AP古典概率模拟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根火柴的概率有多大？（用频率估计概率）functionproguji=liti40(mm)%mm是随机实验次数frq=0;randnum=binornd(1,0.5,mm,2*20);proguji=0;fori=1:mma1=0;a2=0;j=1;while(a120)&(a220)ifrandnum(i,j)==1a1=a1+1;elsea2=a2+1;endj=j+1;endifabs(a1-a2)=5frq=frq+1;endend;proguji=frq/mmliti40(100)proguji=0.4800liti40(1000)proguji=0.4970liti40(10000)proguji=0.4910liti40(100000)proguji=0.4984例4两盒火柴，每盒n根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有k根火柴的概率有多大？（用频率估计概率）functionproguji=liti4(n,k,mm)%n是每盒中的火柴数%k是剩余的火柴数%mm是随机实验次数frq=0;randnum=binornd(1,0.5,mm,2*n);proguji=0;fori=1:mma1=0;a2=0;j=1;while(a1n)&(a2n)ifrandnum(i,j)==1a1=a1+1;elsea2=a2+1;endj=j+1;endifabs(a1-a2)=k,frq=frq+1;end%a1=a1,a2=a2,frq%pauseend;proguji=frq/mmliti4(20,5,100)proguji=0.4800liti4(20,5,1000)proguji=0.4970liti4(20,5,10000)proguji=0.4910liti4(20,5,100000)proguji=0.4984几何概率模拟1.定义向任一可度量区域G内投一点，如果所投的点落在G中任意可度量区域g内的可能性与g的度量成正比，而与g的位置和形状无关，则称这个随机试验为几何型随机试验。或简称为几何概型。2.概率计算P(A)=[A的度量]/[S的度量]例5两人约定于12点到1点到某地会面，先到者等20分钟后离去，试求两人能会面的概率？解：设x,y分别为甲、乙到达时刻(分钟)令A={两人能会面}={(x,y)||x-y|≤20，x≤60,y≤60}P(A)=A的面积/S的面积=（602-402）/602=5/9=0.5556functionproguji=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/mmliti5(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.模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜mm?E=(k2+k3)/mm停止硬币正面?YNNY1，2，34，56functionliti6(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.理论计算设：观察所对目标指示正确确观