蒙克卡罗方法

《蒙克卡罗方法》由会员分享，可在线阅读，更多相关《蒙克卡罗方法（114页珍藏版）》请在装配图网上搜索。

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

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

3、以反复进行，改变系统的结构和系数都比较容易。蒙特卡洛（蒙特卡洛（Monte CarloMonte Carlo）方法）方法是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样，通过对样本值的观察统计，求得所研究系统的某些参数 蒙特卡洛方法也称为蒙特卡洛方法也称为随机模拟方法,其起源最早可以追溯到18世纪下半叶的Buffon试验.例 在1777年,法国学者Buffon提出用试验方法求圆周率鸬闹.其原理如下：假设平面上有元数条距离为1的等矩平行线,现向该平面随机地投掷一根长度为KI1的针,则我们可以计算该针与任一平行线相交的概率.此处随机投针可以这样理解z针的中心点与最近的平行线

4、间的距离Z均匀地分布在区间0.1/2上,针与平行线的夹角以不管相交与否)均匀地分布在区间0,而上(见图6。.于是,针与线相交的充要条件是本寸,从而针线相交概率为1用蒙特卡洛方法进行计算机模拟的步骤用蒙特卡洛方法进行计算机模拟的步骤：1 设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。2 模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象产生模拟随机数的计算机命令产生模拟随机数的计算机命令 在Matlab软件中，可以直接产生满足各种分布的随机数，命令如下：2产生m*n阶离散均匀分布的随机数矩阵：R=unidrnd(N)R=unidrnd(N,mm,nn)当只知

5、道一个随机变量取值在（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服从参数为 的指数分指数分布布。0001)(/xx

6、exft指数分布的期望值为 排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。指数分布在排队论、可靠性分析中有广泛应用。注意：注意：Matlab中，产生参数为 的指数分布的命令为exprnd()例例 顾客到达某商店的间隔时间服从参数为顾客到达某商店的间隔时间服从参数为1010的指的指数分布数分布 指数分布的均值为指数分布的均值为1010。指两个顾客到达商店的平均间隔时间是指两个顾客到达商店的平均间隔时间是1010个单个单位时间位时间.即平均即平均1010个单位时间到达个单位时间到达1 1个顾客个顾客.顾客顾客到达的间隔时间可用到达的间隔时间可用exprnd(1

7、0)exprnd(10)模拟。模拟。设离散型随机变量X的所有可能取值为0,1,2,且取各个值的概率为其中 0为常数，则称X服从参数为 的帕松分布帕松分布。,2,1,0,!)(kkekXPk帕松分布在排队系统、产品检验、天文、物理等领域有广泛应用。帕松分布的期望值为1 1 事件的频率事件的频率 在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，次试验，记记m m是是n n次试验中事件次试验中事件A A发生的次数。发生的次数。频率频率 f=m/n f=m/n 2.2.频率的稳定性频率的稳定性 掷一枚均匀硬币，记录掷硬币试验中频率掷一枚均匀硬币，记录掷硬币试验中频率P P*的波动情况

8、。的波动情况。R=binornd(N,P,mm,nn)例例1 1 频率的稳定性频率的稳定性3 3 概率的频率定义概率的频率定义 在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，记次试验，记m m是是n n次次试验中事件试验中事件A A发生的次数。当试验次数发生的次数。当试验次数n n很大时，很大时，如果频率如果频率m/nm/n稳定地在某数值稳定地在某数值p p附近摆动，而且一附近摆动，而且一般地说，随着试验次数的增加，这种摆动的幅度般地说，随着试验次数的增加，这种摆动的幅度越来越小，称数值越来越小，称数值p p为事件为事件A A在这一组不变的条件在这一组不变的条件下发生的概率

9、，记作下发生的概率，记作P(A)=p.P(A)=p.4 频率的基本性质频率的基本性质（1）对任意事件对任意事件A，有，有 1)(0AP（2）1)(SP0)(P（3）若）若A1，A2，An是互不相容的，则是互不相容的，则)()(11nkknkkAPAP 频率定义的意义频率定义的意义:(1)提供了估计概率的方法提供了估计概率的方法;(2)提供了一种检验理论正确与否的准则提供了一种检验理论正确与否的准则.理论依据：理论依据：大数定律大数定律 大量的随机现象中平均结果的稳定性大量的随机现象中平均结果的稳定性 大数定律的客观背景大数定律的客观背景大量抛掷硬币大量抛掷硬币正面出现频率正面出现频率字母使用频

10、率字母使用频率生产过程中的生产过程中的废品率废品率大数定律大数定律贝努里（Bernoulli）大数定律设 nA 是 n 次独立重复试验中事件 A 发生的次数,p 是每次试验中 A 发生的概率，则0有0limpnnPAn或1limpnnPAn在概率的统计定义中，事件 A 发生的频率“稳定于”事件 A 在一次试验中发生的概率是指：nnAnnA频率与 p 有较大偏差pnnA是小概率事件,因而在 n 足够大时,可以用频率近似代替 p.这种稳定称为依概率稳定.贝努里（贝努里（Bernoulli）大数定律的意义大数定律的意义：定义a 是一常数，0limaYPnn(或)1limaYPnn则称随机变量序列,2

11、1nYYY依概率收敛于常数 a,记作aYnPn故pnnnPA,21nYYY是一系列随机变量，设0有若在 Bernoulli 定理的证明过程中，Y n 是相互独立的服从 0-1分布的随机变量序列 Xk 的算术平均值,Y n 依概率收敛于其数学期望 p.结果同样适用于服从其它分布的独立随机变量序列Chebyshev 大数定律,21nXXX相互独立，设随机变量序列(指任意给定 n 1,相互独立)，且具有相同的数学期望和方差nXXX,21,2,1,)(,)(2kXDXEkk则0有01lim1nkknXnP或11lim1nkknXnP定理的意义定理的意义：当 n 足够大时，算术平均值几乎就是一个常数,可

12、以用算术平均值近似地代替数学期望.具有相同数学期望和方差的独立随机变量序列的算术平均值依概率收敛于数学期望.例如要估计某地区的平均亩产量，要例如要估计某地区的平均亩产量，要收割某些有代表性的地块，例如收割某些有代表性的地块，例如n n 块块.计算其平均亩产量，则当计算其平均亩产量，则当n n 较大时，可用较大时，可用它作为整个地区平均亩产量的一个估计它作为整个地区平均亩产量的一个估计.辛钦大数定律辛钦大数定律 设,21nXXX相互独立，服从同一分布，且具有数学期望 E(X k)=,k=1,2,则对任意正数 001lim1nkknXnP,21nXXX相互独立，注3:设随机变量序列,2,1,)(i

13、XEkki则0有01lim1knikinXnP具有相同的分布，且记knikiMXn1111nPA),(21kAAAgnP),(21kg则则22nPAknPkA),(21kxxxg连续，若 大数定律以严格的数学形式表达了随大数定律以严格的数学形式表达了随机现象最根本的性质之一：机现象最根本的性质之一：它是随机现象统计规律的具体表现它是随机现象统计规律的具体表现.大数定律在理论和实际中都有广泛的应用大数定律在理论和实际中都有广泛的应用.平均结果的稳定性平均结果的稳定性例例1 1 频率的稳定性频率的稳定性1 1 事件的频率事件的频率 在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，

14、次试验，记记m m是是n n次试验中事件次试验中事件A A发生的次数。发生的次数。频率频率 f=m/n f=m/n 2.2.频率的稳定性频率的稳定性 掷一枚均匀硬币，记录掷硬币试验中频率掷一枚均匀硬币，记录掷硬币试验中频率P P*的波动情况。的波动情况。R=binornd(N,P,mm,nn)function liti1(n,p,mm)pro=zeros(1,mm);randnum=binornd(n,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

15、中编辑中编辑.m.m文件输入以下命令：文件输入以下命令：在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(1,0.5,1000)在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(1,0.5,10000)练习练习 频率的稳定性频率的稳定性1 1 事件的频率事件的频率 R=binornd(N,P,mm,nn)在一组不变的条件下，重复作在一组不变的条件下，重复作n n次试验，次试验，记记m m是是n n次试验中事件次试验中事件A A发生的次数。发生的次数。频率频率 f=m/n f=m/n 2.2.频率的稳定性频率的稳定性 练习练习掷

16、一枚不均匀硬币，正面出现概率为掷一枚不均匀硬币，正面出现概率为0.30.3，记录前记录前10001000次掷硬币试验中正面频率的波动次掷硬币试验中正面频率的波动情况，并画图。情况，并画图。在在MatlabMatlab命令行中输入以下命令：命令行中输入以下命令：liti1(1,0.3,1000)例例2 2 掷两枚不均匀硬币，每枚正面出现概率掷两枚不均匀硬币，每枚正面出现概率为为0.40.4，记录前，记录前10001000次掷硬币试验中两枚都为次掷硬币试验中两枚都为正面频率的波动情况，并画图。正面频率的波动情况，并画图。在在Matlab中编辑中编辑.m文件输入以下命令：文件输入以下命令：funct

17、ion liti2(n,p,mm)pro=zeros(1,mm);randnum=binornd(n,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)熊宇乐 y=zeros(1,1000);a=binornd(1,0.4,1,1000);b=binornd(1,0.4,1,1000);c=0;d=0;for i=1:1000 c=c+a(1,i).*b(1,i);y(i)=c/i;end y=y;num=1:1000;plot(num,y)孟亚fu

18、nction bino(n,p,m)x=binornd(n,p,1,m);y=binornd(n,p,1,m);for i=1:m if x(i)=1&y(i)=1 s(i)=1;else s(i)=0;endend for i=1:m y(i)=sum(s(1,1:i)/i;endplot(y)liti2(1,0.4,100)liti2(1,0.4,10000)在一袋中有在一袋中有10 10 个相同的球，分别标有号码个相同的球，分别标有号码1,2,1,2,10,10。每次任取一个球，记录其号码。每次任取一个球，记录其号码后放回袋中，再任取下一个。这种取法叫做后放回袋中，再任取下一个。这种取法

19、叫做“有放回抽取有放回抽取”。今有放回抽取。今有放回抽取3 3个球，求个球，求这这3 3个球的号码均为偶数的概率。个球的号码均为偶数的概率。（用频率估计（用频率估计概率）概率）例例3:解：令解：令A=有放回抽取有放回抽取3个球，求这个球，求这3个球的号码个球的号码均为偶数均为偶数=(2,2,2)，(2,2,4)，.，(10,10,10)81)(APfunction 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

20、,3)+1)；if mod(a,2)=1 frq=frq+1 endend;proguji=frq/mm例例4 4 两盒火柴，每盒两盒火柴，每盒2020根。每次随机在任一根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有完而另一盒中至少还有5 5根火柴的概率有多大？根火柴的概率有多大？（用频率估计概率）（用频率估计概率）liti4(20,5,100)proguji=0.4800 liti4(20,5,1000)proguji=0.4970 liti4(20,5,10000)proguji=0.4910 liti4(20,5,100

21、000)proguji=0.4984function proguji=liti4(nn,num,mm)%nn 是每盒中的火柴数%num 是剩余的火柴数%mm 是随机实验次数frq=0;randnum=binornd(1,0.5,mm,2*nn);proguji=0;for i=1:mm a1=0;a2=0;j=1;while(a120)&(a2=5 frq=frq+1;end%a1=a1,a2=a2,frq%pauseend proguji=frq/mm二二.几何概率几何概率1.1.定义定义 向任一可度量区域向任一可度量区域G G内投一点，如果所投内投一点，如果所投的点落在的点落在G G中任意

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

23、60,y y6060P(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在我方某前沿防守地域，敌人

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

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

26、火炮的射击总数；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=0;k3=0;for i=1:mm if randnum1(i)=0 k1=k1+1;

27、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命令行中输入以下命令：命令行中输入以下命令：liti6(0.5,20000)5.5.理论计算理论计算6.6.结果比较结果比较 模拟结

28、果与理论计算近似一致，能更模拟结果与理论计算近似一致，能更加真实地表达实际战斗动态过程加真实地表达实际战斗动态过程 3.某厂生产的灯泡能用某厂生产的灯泡能用1000小时的概率为小时的概率为0.8,能能用用1500小时的概率为小时的概率为0.4,求已用求已用1000小时的灯泡小时的灯泡能用到能用到1500小时的概率小时的概率2.在一袋中有在一袋中有10 个相同的球，分别标有号码个相同的球，分别标有号码1,2,10。今有放回任取两个球，求取得的第一个球号码。今有放回任取两个球，求取得的第一个球号码为奇数，第二个球的号码为偶数的概率。为奇数，第二个球的号码为偶数的概率。1.掷三枚不均匀硬币，每枚正面

29、出现概率为掷三枚不均匀硬币，每枚正面出现概率为0.3，记，记录前录前1000次掷硬币试验中至少两枚都为正面频率的次掷硬币试验中至少两枚都为正面频率的波动情况，并画图。波动情况，并画图。作业作业:某厂生产的灯泡能用1000小时的概率为0.8,能用1500小时的概率为0.4,求已用1000小时的灯泡能用到1500小时的概率解解 令 A 灯泡能用到1000小时 B 灯泡能用到1500小时所求概率为)()(APABPABPAB218.04.0)()(APBP例例 例例:在一袋中有在一袋中有10 个相同的球，分别标有号个相同的球，分别标有号码码1,2,10。今不放回任取两个球，求取。今不放回任取两个球，

30、求取得的第一个球号码为奇数，第二个球的得的第一个球号码为奇数，第二个球的号码为偶数的概率。号码为偶数的概率。解：令解：令A=A=抽取抽取2 2个球，第一个球号码为奇数，第二个球个球，第一个球号码为奇数，第二个球的号码为偶数的号码为偶数=(1,2)=(1,2)，(1,4)(1,4)，.，(9,10)(9,10)185)(AP2022-10-27圆的面积 单位圆的面积等于 计算第一象限内的单位圆的面积，方法为：把它分成n个窄的曲边梯形，计算S大，S小，其中n可以为1000,10000,.实验二实验二:的计算的计算2022-10-27数值积分法4)1(141102102dxxdxx2022-10-2

31、7梯形公式inabaxnabxfbfafdxxfiniiba,)(2/)()()(112022-10-27辛普森公式nabxfxfbfafdxxfniiniiba6)(2)(4)()()(11105.02022-10-27无穷级数法 arctg x=x-x3/3+x5/5-x7/7+x9/9-./4=arctg 1=1-1/3+1/5-1/7+1/9-收敛太慢！|x|应当比 1 小很多，级数收敛才快。/4=arctg 1/2+arctg 1/3 /4=4 arctg 1/5-arctg 1/2392022-10-27 某次随机投点的结果0.20.40.60.810.20.40.60.81 Mo

32、nte-Carlo 方法方法谢文贤谢文贤西北工业大学理学院西北工业大学理学院 高性能计算中心高性能计算中心2004-05-142004-05-14Monte-Carlo 方法方法 一、一、MC 的起源和发展的起源和发展 Buffon试验试验 排队系统模拟：排队系统模拟：M/G/1M/G/1排队系统排队系统 二、二、MC 的原理的原理 三、随机数的产生原理与三、随机数的产生原理与IMSL库库 均匀分布均匀分布U(0,1)的随机数的产生的随机数的产生 其他各种分布的随机数的产生其他各种分布的随机数的产生 随机过程模拟随机过程模拟 四、四、MC的应用举例的应用举例 定积分的定积分的MC计算计算 随机

33、微分方程的数值模拟随机微分方程的数值模拟 五、五、EM算法及其算法及其MCMC方法方法一、一、MC 的起源和发展的起源和发展 随机模拟方法，也称为Monte Carlo方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战进行的研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。冯诺伊曼是公理化方法和计算机体系的领袖人物，Monte Carlo方法也是他的功劳。事实上，Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定

34、事件的“概率”。18世纪下半叶的法国学者Buffon提出用投针试验的方法来确定圆周率的值。这个著名的Buffon试验是Monte Carlo方法的最早的尝试！历史上曾有几位学者相继做过这样的试验。不过呢，他们的试验是费时费力的，同时精度不够高，实施起来也很困难。然而，随着计算机技术的飞速发展，人们不需要具体实施这些试验，而只要在计算机上进行大量的、快速的模拟试验就可以了。在大众的心目中，科学的代言人是心不在焉的牛顿或者爆炸式发型的爱因斯坦，但这只是传统形象，比他们更了解现代计算技术的冯诺伊曼是个衣着考究，风度翩翩的人物，他说：纯粹数学和应用数学的许多分支非常需要计算工具，用以打破目前由于纯粹分

35、析的研究方法不能解决非线性问题而形成的停滞状态。Monte Carlo方法是现代计算技术的最为杰出的成果之一，它在工程领域的作用是不可比拟的。Buffon试验试验 假设平面上有无数条距离为1的等距平行线，现向该平面随机投掷一根长度为 的针()，则我们可计算该针与任一平行线相交的概率。这里，随机投针指的是：针的中心点与最近的平行线间的距离 均匀的分布在区间 。上，针与平行线的夹角 （不管相交与否）均匀的分布在区间 上。因此，针与线相交的充要条件是 l1l21,0 x,021sinxBuffon试验试验 从而针线相交的概率为 根据上式，若我们做大量的投针试验并记录针与线相交的次数，则由大数定理可以

36、估计出针线相交的概率 ，从而得到 的估计值。针与线的位置关系：ldxdlXPpl22sin20sin20 xlp排队系统模拟：排队系统模拟：M/G/1M/G/1排队系统排队系统服务规则：服务规则：先到先服务假设假设：（）顾客到达遵循Poisson分布；（）服务时间服从一般分布；（）到达间隔与服务时间相互独立1排队系统模拟：排队系统模拟：M/G/1M/G/1排队系统排队系统关心的指标：关心的指标：（）时刻t时，系统中的顾客数；即队长的分布；（）顾客的等待时间；（）服务的忙碌程度；（）用最朴素的Monte-Carlo方法可以得到这些指标的估计二、二、MC 的原理的原理 应用Monte Carlo方

37、法求解工程技术问题可以分为两类：确定性问题 随机性问题确定性系统确定性系统随机性系统随机性系统模拟模拟自然界自然界Monte-Carlo模拟，即随机模拟（重复模拟，即随机模拟（重复“试验试验”）重复试验重复试验计算机模拟计算机模拟思路：思路：1、针对实际问题建立一个简单且便于实现的概率统计模型，使问题的解对应于该模型中随机变量的概率分布或其某些数字特征，比如，均值和方差等。所构造的模型在主要特征参量方面要与实际问题或系统相一致的。2、根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，再进行随机模

38、拟试验。思路：思路：3、根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）。4、按照所建立的模型进行仿真试验、计算，求出问题的随机解。5、统计分析模拟试验结果，给出问题的估计以及其精度估计。6、必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。收敛性收敛性:由大数定律,Monte-Carlo模拟的收敛是以概率而言的 误差误差:用频率估计概率时误差的估计，可由中心极限定理，给定置信水平 的条件下，有：模拟次数模拟次数:由误差公式得NU2/1|)(XgVar22/1)(UN三、随机数的产生原

39、理与三、随机数的产生原理与IMSL库库 随机数的产生是实现MC计算的先决条件。而大多数概率分布的随机数的产生都是基于均匀分布U(0,1)的随机数。首先，介绍服从均匀分布U(0,1)的随机数的产生方法。其次，介绍服从其他各种分布的随机数的产生方法。以及服从正态分布的随机数的产生方法。接下来，介绍随机过程模拟。最后，关于随机数的几点注。均匀分布均匀分布U(0,1)的随机数的产生的随机数的产生 产生均匀分布的标准算法在很多高级计算机语言的书都可以看到。算法简单，容易实现。使用者可以自己手动编程实现。IMSL统计库中也提供给我们用于产生均匀分布的各种函数。我们的重点是怎样通过均匀分布产生服从其他分布的

40、随机数。因此，直接使用IMSL提供的可靠安全的标准函数，当然不用费事了。IMSL库中的函数使用库中的函数使用 RNSET:种子的设定 CALL RNSET(ISEED)RNOPT:产生器的类型的设定 CALL RNOPT(IOPT)RNUN/DRNUN:产生均匀分布的随机数 CALL RNUN(NR,R)其他各种分布的随机数的产生其他各种分布的随机数的产生 基本方法有如下三种：逆变换法 合成法 筛选法 逆变换法逆变换法 设随机变量 的分布函数为 ，定义 定理定理 设随机变量 服从 上的均匀分布，则 的分布函数为 。因此，要产生来自 的随机数，只要先产生来自 的随机数，然后计算 即可。其步骤为

41、X xF 10,:inf1yyxFxyFU)1,0(UFX1 xF xF1,0U uF1 uFxuU11,0计算，抽取由合成法合成法 合成法的应用最早见于Butlter 的书中。构思如下：如果 的密度函数 难于抽样，而 关于 的条件密度函数 以及 的密度函数 均易于抽样，则 的随机数可如下产生：可以证明由此得到 的服从 。X xpXYyxpY ygX xyxpyygY抽取由条件分布，抽取的分布由X xp筛选抽样筛选抽样 假设我们要从 抽样，如果可以将 表示成 ，其中 是一个密度函数且易于抽样，而 ，是常数，则 的抽样可如下进行：定理定理 设 的密度函数 ，且 ，其中 ，是一个密度函数。令 和

42、分别服从 和 ，则在 的条件 下，的条件密度为 xp xgxhcxp xph 10 xg1cX 。，回到如果，停止，则如果，抽取，由抽取由1321,01yguyxyguyyhuUX xp xgxhcxph 10 xg1cUY1,0U yh YgU xpYgUxpYY标准正态分布的随机数标准正态分布的随机数 的随机数产生方法很多。简要介绍三种。法法1、变换法（Box 和Muller 1958）设 ，是独立同分布的 变量，令 则 与 独立，均服从标准正态分布。法法2、结合合成法与筛选法。（略）法法3 3、近似方法（利用中心极限定理）即用 个 变量产生一个 变量。其中 是抽自 的随机数，可近似为一个

43、 变量。1,0N1U2U1,0U2122112sinln22cosln22121UUXUUX1X2Xn1,0U1,0N2112unxiu1,0Uniiunu111,0NIMSL库中的函数使用库中的函数使用 RNSET:种子的设定 CALL RNSET(ISEED)RNOPT:产生器的类型的设定 CALL RNOPT(IOPT)RNNOA:产生服从标准正态分布的伪随机数 CALL RNNOA(NR,R)随机过程模拟随机过程模拟 高斯白噪声的产生：一种简单有效的模拟高斯白噪声过程的方法是由独立的单位正态随机序列 ，按照如图所示的方式连接，其中 式中 为白噪声的强度。kktDkD 由于正态随机数的独

44、立性，当 很小时，按(a)、(b)所构成的过程的相关时间很短，从而具有很高的截止频率。当然，过小，样本的计算量过大。因此，选取 适当小即可。ttt关于随机数的几点注关于随机数的几点注 注注1 1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注注2 2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。四、四、MC的

45、应用举例的应用举例 例一、定积分的例一、定积分的MC计算计算 随机投点法随机投点法 样本平均值法样本平均值法 几种降低估计方差的几种降低估计方差的MC方法方法 例二、例二、随机微分方程的数值模拟随机微分方程的数值模拟 系统的可靠性系统的可靠性数值模拟数值模拟计算问题计算问题定积分的定积分的MC计算计算 事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。下面考虑一个简单的定积分 为了说明问题，我们首先介绍两种求 的简单的MC方法，然后给出几种较为复杂而更有效的MC方法。dxxfba 在计算积分上，MC的实用场合是计算重积分 其中 是 维空间的点，当 较大时，用MC方

46、法比一般的数值方法有优点，主要是它的误差与维数 无关。dPPgIkPmmm随机投点法随机投点法 方法简述方法简述：设，有限，并设 是在 上均匀分布的二维随机变量，其联合密度函数为 。则易见 是 中 曲线下方的面积。假设我们向 中进行随机投点，若点落在 下方，(即 称为中的，否则称为不中，则点中的概率为 。若我们进行了 次投点，其中 次中的，则用频率来估计概率 。即 。ab Mxf0Mybxayx0,:,YX,MybxaIabM0,1 dxxfba xfy xfy xfy abMpn0npnnp0 那么我们可以得到 的一个估计 具体试验步骤为nnabM01 用上式来估计的个数统计，和，计算，随机

47、数，个独立地产生0,.,11,02nyxfxfMvyabuaxnivuUniiiiiiiii 注注1 1 随机投点法的思想简单明了，且每次投点结果服从二项分布，故 ，其中 注注2 2 可证 是 的无偏估计。若用估计的标准差来衡量其精度，则估计 的精度的阶为 。注注3 3 这里，定积分的解，就对应我们选定的随机变量的概率值。pnbn,0abMp1121n样本平均值法样本平均值法 基本原理基本原理：对积分 ，设 是 上的一个密度函数，改写 可见，任一积分均可以表示为某个随机变量（函数）的期望。由矩法，若有 个来自 的观测值，则可给出 的一个矩估计。最简单的，若，有限，可取 。dxxfba xgba

48、,XgXfEdxxgxgxfban xgab abxg/1 设 是来自 的随机数，则 的一个估计为 具体步骤为 注注 可证 是 的无偏估计。一般而言，样本均值法要比随机投点法更有效。nxx,.,1baU,niiniiixfnabxgxfn1121 用上式来估计，和计算，随机数个独立地产生nixfuabaxuuUniiin,.,1,.,1,012求解定积分的算例求解定积分的算例 计算定积分 事实上，其精确解为 随机投点法：sjtd_djf.f90，及动态图 样本平均值法：ybpjz_djf.f90 注注 所有计算中的随机数的产生均来自IMSL库。样本数为10000。增加样本数目，可提高计算精度，

49、但计算时间也会提高。0.4sin0.8dxx400.2cos几种降低估计方差的几种降低估计方差的MC方法方法 重要抽样法重要抽样法 特点：相对样本均值法而言，样本均值法是由于假设 是均匀分布的概率密度，故采用的是均匀抽样，各随机数 是均匀分布的随机数，各 对 的贡献是不同，大则贡献大，但在抽样时，这种差别未能体现出来。而重要抽样法，则希望贡献率大的随机数出现的概率大，贡献小的随机数出现概率小，从而提高抽样的效。xgixix2 ixf几种降低估计方差的几种降低估计方差的MC方法方法 重要抽样法重要抽样法 关键因素在于 的选取，使得估计的方差较小。重要抽样法的基本思想，就是通过选取与 形状接近的密

50、度函数 来降低估计的方差。xg xf xg几种降低估计方差的几种降低估计方差的MC方法方法 分层抽样法分层抽样法 同样是利用贡献率大小来降低估计方差的方法。它首先是把样本空间 分成一些小区间 ，且诸 不交，然后在各个小区间内的抽样数由其贡献大小决定。对 贡献大的 抽样多，可提高抽样效率。如果能够提出较好抽样区间的分配和各子区间内抽样次数的分配方案，分层抽样法估计积分可以达到非常令人满意的效果。DmDD,.,1iDDUDiiD几种降低估计方差的几种降低估计方差的MC方法方法 关联抽样法关联抽样法 将需要估计的积分分解成两个积分之差，对 的估计转化为对 ，的估计的差。则相应的，其估计的方差的大小则

51、与 ，的估计的正相关度有关，若两者的相关程度越高，则 的估计方差越小。这便是关联抽样法的基本出发点。2121IIdxxfdxxfdxxfbababa1I2I1I2I随机微分方程的数值模拟随机微分方程的数值模拟 考虑It随机微分方程 一般，可以求精确解析解的It随机微分方程只有两类，一种是线性的It随机微分方程，另一种是可用It微分公式化为线性的非线性It随机微分方程，并且主要是标量方程。而实用中较复杂的It随机微分方程则只能用数值解法。对上式用最简单的Euler-Maruyama近似：tdBXdtXmtdXtttttBXtXmXXttt 求解下面受谐和激励与高斯白噪声激励的Duffing振子：

52、参数取值：将上式化为标准It微分方程，借助 Euler-Maruyama近似，可求得上述Duffing方程的时间历程图，及相轨图。由于需要大量计算时间，我们只给出结果，以作参考。ttxxxxcos322 4.1,1.02 5.50,0.40 xx 时间历程图 相轨图一个元件(或系统)能正常工作的概率称为元件(或系统)的可靠性系统由元件组成,常见的元件连接方式：串联并联1221系统的可靠性计算问题 例例5 5设 两系统都是由 4 个元件组成,每个元件正常工作的概率为 p,每个元件是否正常工作相互独立.两系统的连接方式如下图所示，比较两系统的可靠性.A1A2B2B1S1:)()()()()()(2

53、121212121211BBAAPBBPAAPBBAAPSP)2(22242ppppA1A2B2B1S2:212)()(iiiBAPSP)()(12SPSP22)2(pp)2(22pp222pp五、五、EM算法及其算法及其MCMC方法方法 EM算法：算法：是一种迭代方法，最初由Dempster等提出，并主要应用于较为复杂的后验分布，来计算后验均值或后验众数，即极大似然估计的一种数据添加算法。最大优点是简单和稳定。MCMC算法：算法：当后验分布较为复杂时，对于后验分布的积分计算，像后验均值、后验方差、后验分布的分位数等等，就不得不求助于MCMC算法。他在统计物理学中得到广泛的应用，近年来，迅速发展到Bayes统计、显著性检验、极大似然估计等方面。关于这两方面的详细知识，感兴趣的可自己阅读！