《蒙克卡罗方法》PPT课件.ppt

《《蒙克卡罗方法》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《蒙克卡罗方法》PPT课件.ppt（114页珍藏版）》请在三一办公上搜索。

1、蒙特卡洛方法基本思想,实验目的,实验内容,学习计算机模拟的基本过程与方法。,1、模拟的概念。,4、实验作业。,3、计算机模拟实例。,2、产生随机数的计算机命令。,模拟的概念,模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。,模拟的基本思想是建立一个试验模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息,模拟的方法,1、物理模拟：对实际系统及其过程用功能相似的实物系统去模仿。例如，军事演习、船艇实验、沙盘作业等。,物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难。而且，许多系统无法进行物理模

2、拟，如社会经济系统、生态系统等。,在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。,在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟。现代的数学模拟都是在计算机上进行的，称为计算机模拟。,2、数学模拟,计算机模拟可以反复进行，改变系统的结构和系数都比较容易。,蒙特卡洛（Monte Carlo）方法是一种应用随机数来进行计算机模拟的方法此方法对研究的系统进行随机观察抽样，通过对样本值的观察统计，求得所研究系统的某些参数,蒙特卡洛方法也称为随机模拟方法,其起源最早可

3、以追溯到18世纪下半叶的Buffon试验.,例 在1777年,法国学者Buffon提出用试验方法求圆周率鸬闹.其原理如下：假设平面上有元数条距离为1的等矩平行线,现向该平面随机地投掷一根长度为KI1的针,则我们可以计算该针与任一平行线相交的概率.此处随机投针可以这样理解z针的中心点与最近的平行线间的距离Z均匀地分布在区间0.1/2上,针与平行线的夹角以不管相交与否)均匀地分布在区间0,而上(见图6。.于是,针与线相交的充要条件是本寸,从而针线相交概率为1,用蒙特卡洛方法进行计算机模拟的步骤：,1 设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系。2 模拟随机现象可通过具

4、有各种概率分布的模拟随机数来模拟随机现象,产生模拟随机数的计算机命令,在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),当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时

5、，可以认为该对象服从正态分布。,机械加工得到的零件尺寸的偏差、射击命中点与目标的偏差、各种测量误差、人的身高、体重等，都可近似看成服从正态分布。,若连续型随机变量X的概率密度函数为 其中 0为常数，则称X服从参数为 的指数分布。,指数分布的期望值为,排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。,指数分布在排队论、可靠性分析中有广泛应用。,注意：Matlab中，产生参数为 的指数分布的命令为exprnd(),例 顾客到达某商店的间隔时间服从参数为10的指数分布,指数分布的均值为10。指两个顾客到达商店的平均间隔时间是10个单位时间.即平均10个单位时间到达

6、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很大时，如

7、果频率m/n稳定地在某数值p附近摆动，而且一般地说，随着试验次数的增加，这种摆动的幅度越来越小，称数值p为事件A在这一组不变的条件下发生的概率，记作P(A)=p.,4 频率的基本性质,（1）对任意事件A，有,（2）,（3）若A1，A2，An是互不相容的，则,频率定义的意义:,(1)提供了估计概率的方法;(2)提供了一种检验理论正确与否的准则.,理论依据：大数定律,大量的随机现象中平均结果的稳定性,大数定律的客观背景,大量抛掷硬币正面出现频率,字母使用频率,生产过程中的废品率,大数定律,贝努里（Bernoulli）大数定律,设 nA 是 n 次独立重复试验中事件 A 发生的次数,p 是每次试验中

8、 A 发生的概率，则,有,或,贝努里（Bernoulli）大数定律的意义：,定义,a 是一常数，,(或,故,在 Bernoulli 定理的证明过程中，Y n 是相互独立的服从 0-1分布的随机变量序列 Xk 的算术平均值,Y n 依概率收敛于其数学期望 p.,结果同样适用于服从其它分布的独立随机变量序列,Chebyshev 大数定律,(指任意给定 n 1,相互独立)，且具有相同的数学期望和方差,或,定理的意义：,当 n 足够大时，算术平均值几乎就是一个常数,可以用算术平均值近似地代替数学期望.,具有相同数学期望和方差的独立随机变量序列的算术平均值依概率收敛于数学期望.,例如要估计某地区的平均亩

9、产量，要收割某些有代表性的地块，例如n 块.计算其平均亩产量，则当n 较大时，可用它作为整个地区平均亩产量的一个估计.,辛钦大数定律,具有相同的分布，且,记,则,则,大数定律以严格的数学形式表达了随机现象最根本的性质之一：,它是随机现象统计规律的具体表现.,大数定律在理论和实际中都有广泛的应用.,平均结果的稳定性,例1 频率的稳定性,1 事件的频率 在一组不变的条件下，重复作n次试验，记m是n次试验中事件A发生的次数。频率 f=m/n,2.频率的稳定性,掷一枚均匀硬币，记录掷硬币试验中频率P*的波动情况。R=binornd(N,P,mm,nn),function liti1(n,p,mm)pr

10、o=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),在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.频率的

11、稳定性,练习掷一枚不均匀硬币，正面出现概率为0.3，记录前1000次掷硬币试验中正面频率的波动情况，并画图。,在Matlab命令行中输入以下命令：,liti1(1,0.3,1000),例2 掷两枚不均匀硬币，每枚正面出现概率为0.4，记录前1000次掷硬币试验中两枚都为正面频率的波动情况，并画图。,在Matlab中编辑.m文件输入以下命令：,function 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=

12、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;endy=y;num=1:1000;plot(num,y),孟亚function bino(n,p,m)x=binornd(n,p,1,m);y=binornd(n,p,1,m);for i=1:m if x(i)=1endplot(y),liti2(1,0.4,100),liti2(1,0.4,10000),在

13、一袋中有10 个相同的球，分别标有号码1,2,10。每次任取一个球，记录其号码后放回袋中，再任取下一个。这种取法叫做“有放回抽取”。今有放回抽取3个球，求这3个球的号码均为偶数的概率。（用频率估计概率）,例3:,解：令A=有放回抽取3个球，求这3个球的号码均为偶数=(2,2,2)，(2,2,4)，.，(10,10,10),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

14、(a,2)=1 frq=frq+1 endend;proguji=frq/mm,例4 两盒火柴，每盒20根。每次随机在任一盒中取出一根火柴。问其中一盒中火柴被取完而另一盒中至少还有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,100000)proguji=0.4984,function proguji=liti4(nn,num,mm)%nn 是每盒中的火柴数%num 是剩余的火柴数%mm 是随

15、机实验次数frq=0;randnum=binornd(1,0.5,mm,2*nn);proguji=0;for i=1:mm a1=0;a2=0;j=1;while(a1=5 frq=frq+1;end%a1=a1,a2=a2,frq%pauseend proguji=frq/mm,二.几何概率,1.定义,向任一可度量区域G内投一点，如果所投的点落在G中任意可度量区域g内的可能性与g的度量成正比，而与g的位置和形状无关，则称这个随机试验为几何型随机试验。或简称为几何概型。,2.概率计算,1.P(A)=A的度量/S的度量,两人约定于12点到1点到某地会面，先到者等20分钟后离去，试求两人能会面的

16、概率？,例5,解：设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;proguji=0;for ii=1:mm if abs(randnum(ii,1)=20 frq=frq+1;endendproguji=frq/mm,

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

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

19、3：击中敌人两门火炮的射击总数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;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)

20、,在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，记录前10

21、00次掷硬币试验中至少两枚都为正面频率的波动情况，并画图。,作业:,某厂生产的灯泡能用1000小时的概率为0.8,能用1500小时的概率为0.4,求已用1000小时的灯泡能用到1500小时的概率,解 令 A 灯泡能用到1000小时 B 灯泡能用到1500小时,所求概率为,例:,在一袋中有10 个相同的球，分别标有号码1,2,10。今不放回任取两个球，求取得的第一个球号码为奇数，第二个球的号码为偶数的概率。,解：令A=抽取2个球，第一个球号码为奇数，第二个球的号码为偶数=(1,2)，(1,4)，.，(9,10),圆的面积,单位圆的面积等于 计算第一象限内的单位圆的面积，方法为：把它分成n个窄的曲

22、边梯形，计算S大，S小，其中n可以为1000,10000,.,实验二:的计算,数值积分法,梯形公式,辛普森公式,无穷级数法,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/239,某次随机投点的结果,Monte-Carlo 方法,谢文贤西北工业大学理学院 高性能计算中心2004-05-14,Monte-Carlo 方法,一、MC 的起源和发展 Buffon试验 排队系统模拟：M/G/1

23、排队系统二、MC 的原理三、随机数的产生原理与IMSL库 均匀分布U(0,1)的随机数的产生 其他各种分布的随机数的产生 随机过程模拟四、MC的应用举例 定积分的MC计算 随机微分方程的数值模拟五、EM算法及其MCMC方法,一、MC 的起源和发展,随机模拟方法，也称为Monte Carlo方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战进行的研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。冯诺伊曼是公理化方法和计算机体系的领袖人物，Monte Carlo方法也是他的功劳。,

24、事实上，Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。18世纪下半叶的法国学者Buffon提出用投针试验的方法来确定圆周率的值。这个著名的Buffon试验是Monte Carlo方法的最早的尝试！,历史上曾有几位学者相继做过这样的试验。不过呢，他们的试验是费时费力的，同时精度不够高，实施起来也很困难。然而，随着计算机技术的飞速发展，人们不需要具体实施这些试验，而只要在计算机上进行大量的、快速的模拟试验就可以了。,在大众的心目中，科学的代言人是心不在焉的牛顿或者爆炸式发型的爱因斯坦，但这只是传统形象，比他们更了

25、解现代计算技术的冯诺伊曼是个衣着考究，风度翩翩的人物，他说：纯粹数学和应用数学的许多分支非常需要计算工具，用以打破目前由于纯粹分析的研究方法不能解决非线性问题而形成的停滞状态。Monte Carlo方法是现代计算技术的最为杰出的成果之一，它在工程领域的作用是不可比拟的。,Buffon试验,假设平面上有无数条距离为1的等距平行线，现向该平面随机投掷一根长度为 的针()，则我们可计算该针与任一平行线相交的概率。这里，随机投针指的是：针的中心点与最近的平行线间的距离 均匀的分布在区间。上，针与平行线的夹角（不管相交与否）均匀的分布在区间 上。因此，针与线相交的充要条件是,Buffon试验,从而针线相

26、交的概率为 根据上式，若我们做大量的投针试验并记录针与线相交的次数，则由大数定理可以估计出针线相交的概率，从而得到 的估计值。针与线的位置关系：,排队系统模拟：M/G/1排队系统,服务规则：先到先服务假设：（）顾客到达遵循Poisson分布；（）服务时间服从一般分布；（）到达间隔与服务时间相互独立,1,排队系统模拟：M/G/1排队系统,关心的指标：（）时刻t时，系统中的顾客数；即队长的分布；（）顾客的等待时间；（）服务的忙碌程度；（）用最朴素的Monte-Carlo方法可以得到这些指标的估计,二、MC 的原理,应用Monte Carlo方法求解工程技术问题可以分为两类：确定性问题 随机性问题,

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

28、直接抽样、分层抽样、相关抽样、重要抽样等）。4、按照所建立的模型进行仿真试验、计算，求出问题的随机解。5、统计分析模拟试验结果，给出问题的估计以及其精度估计。6、必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。,收敛性:由大数定律,Monte-Carlo模拟的收敛是以概率而言的误差:用频率估计概率时误差的估计，可由中心极限定理，给定置信水平 的条件下，有：模拟次数:由误差公式得,三、随机数的产生原理与IMSL库,随机数的产生是实现MC计算的先决条件。而大多数概率分布的随机数的产生都是基于均匀分布U(0,1)的随机数。首先，介绍服从均匀分布U(0,1)的随机数的产生方法。其

29、次，介绍服从其他各种分布的随机数的产生方法。以及服从正态分布的随机数的产生方法。接下来，介绍随机过程模拟。最后，关于随机数的几点注。,均匀分布U(0,1)的随机数的产生,产生均匀分布的标准算法在很多高级计算机语言的书都可以看到。算法简单，容易实现。使用者可以自己手动编程实现。IMSL统计库中也提供给我们用于产生均匀分布的各种函数。我们的重点是怎样通过均匀分布产生服从其他分布的随机数。因此，直接使用IMSL提供的可靠安全的标准函数，当然不用费事了。,IMSL库中的函数使用,RNSET:种子的设定 CALL RNSET(ISEED)RNOPT:产生器的类型的设定 CALL RNOPT(IOPT)R

30、NUN/DRNUN:产生均匀分布的随机数 CALL RNUN(NR,R),其他各种分布的随机数的产生,基本方法有如下三种：逆变换法 合成法 筛选法,逆变换法,设随机变量 的分布函数为，定义 定理 设随机变量 服从 上的均匀分布，则 的分布函数为。因此，要产生来自 的随机数，只要先产生来自 的随机数，然后计算 即可。其步骤为,合成法,合成法的应用最早见于Butlter 的书中。构思如下：如果 的密度函数 难于抽样，而 关于 的条件密度函数 以及 的密度函数 均易于抽样，则 的随机数可如下产生：可以证明由此得到 的服从。,筛选抽样,假设我们要从 抽样，如果可以将 表示成，其中 是一个密度函数且易于

31、抽样，而，是常数，则 的抽样可如下进行：定理 设 的密度函数，且，其中，是一个密度函数。令 和 分别服从 和，则在 的条件下，的条件密度为,标准正态分布的随机数,的随机数产生方法很多。简要介绍三种。法1、变换法（Box 和Muller 1958）设，是独立同分布的 变量，令 则 与 独立，均服从标准正态分布。法2、结合合成法与筛选法。（略）法3、近似方法（利用中心极限定理）即用 个 变量产生一个 变量。其中 是抽自 的随机数，可近似为一个 变量。,IMSL库中的函数使用,RNSET:种子的设定 CALL RNSET(ISEED)RNOPT:产生器的类型的设定 CALL RNOPT(IOPT)R

32、NNOA:产生服从标准正态分布的伪随机数 CALL RNNOA(NR,R),随机过程模拟,高斯白噪声的产生：一种简单有效的模拟高斯白噪声过程的方法是由独立的单位正态随机序列，按照如图所示的方式连接，其中式中 为白噪声的强度。,由于正态随机数的独立性，当 很小时，按(a)、(b)所构成的过程的相关时间很短，从而具有很高的截止频率。当然，过小，样本的计算量过大。因此，选取 适当小即可。,关于随机数的几点注,注1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪

33、随机数。但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。,四、MC的应用举例,例一、定积分的MC计算 随机投点法 样本平均值法 几种降低估计方差的MC方法 例二、随机微分方程的数值模拟 系统的可靠性数值模拟计算问题,定积分的MC计算,事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。下面考虑一个简单的定积分为了说明问题，我们首先介绍两种求 的简单的MC方法，然后给出几种较为复杂而更有效的MC方法。,在计算积分上，MC的实用场合

34、是计算重积分其中 是 维空间的点，当 较大时，用MC方法比一般的数值方法有优点，主要是它的误差与维数 无关。,随机投点法,方法简述：设，有限，并设 是在 上均匀分布的二维随机变量，其联合密度函数为。则易见 是 中 曲线下方的面积。假设我们向 中进行随机投点，若点落在 下方，(即 称为中的，否则称为不中，则点中的概率为。若我们进行了 次投点，其中 次中的，则用频率来估计概率。即。,那么我们可以得到 的一个估计 具体试验步骤为,注1 随机投点法的思想简单明了，且每次投点结果服从二项分布，故，其中注2 可证 是 的无偏估计。若用估计的标准差来衡量其精度，则估计 的精度的阶为。注3 这里，定积分的解，

35、就对应我们选定的随机变量的概率值。,样本平均值法,基本原理：对积分，设 是 上的一个密度函数，改写可见，任一积分均可以表示为某个随机变量（函数）的期望。由矩法，若有 个来自 的观测值，则可给出 的一个矩估计。最简单的，若，有限，可取。,设 是来自 的随机数，则 的一个估计为 具体步骤为注 可证 是 的无偏估计。一般而言，样本均值法要比随机投点法更有效。,求解定积分的算例,计算定积分事实上，其精确解为 随机投点法：sjtd_djf.f90，及动态图样本平均值法：ybpjz_djf.f90注 所有计算中的随机数的产生均来自IMSL库。样本数为10000。增加样本数目，可提高计算精度，但计算时间也会

36、提高。,几种降低估计方差的MC方法,重要抽样法 特点：相对样本均值法而言，样本均值法是由于假设 是均匀分布的概率密度，故采用的是均匀抽样，各随机数 是均匀分布的随机数，各 对 的贡献是不同，大则贡献大，但在抽样时，这种差别未能体现出来。而重要抽样法，则希望贡献率大的随机数出现的概率大，贡献小的随机数出现概率小，从而提高抽样的效。,几种降低估计方差的MC方法,重要抽样法关键因素在于 的选取，使得估计的方差较小。重要抽样法的基本思想，就是通过选取与 形状接近的密度函数 来降低估计的方差。,几种降低估计方差的MC方法,分层抽样法 同样是利用贡献率大小来降低估计方差的方法。它首先是把样本空间 分成一些

37、小区间，且诸 不交，然后在各个小区间内的抽样数由其贡献大小决定。对 贡献大的 抽样多，可提高抽样效率。如果能够提出较好抽样区间的分配和各子区间内抽样次数的分配方案，分层抽样法估计积分可以达到非常令人满意的效果。,几种降低估计方差的MC方法,关联抽样法 将需要估计的积分分解成两个积分之差，对 的估计转化为对，的估计的差。则相应的，其估计的方差的大小则与，的估计的正相关度有关，若两者的相关程度越高，则 的估计方差越小。这便是关联抽样法的基本出发点。,随机微分方程的数值模拟,考虑It随机微分方程 一般，可以求精确解析解的It随机微分方程只有两类，一种是线性的It随机微分方程，另一种是可用It微分公式

38、化为线性的非线性It随机微分方程，并且主要是标量方程。而实用中较复杂的It随机微分方程则只能用数值解法。对上式用最简单的Euler-Maruyama近似：,求解下面受谐和激励与高斯白噪声激励的Duffing振子：参数取值：将上式化为标准It微分方程，借助 Euler-Maruyama近似，可求得上述Duffing方程的时间历程图，及相轨图。由于需要大量计算时间，我们只给出结果，以作参考。,时间历程图 相轨图,一个元件(或系统)能正常工作的概率称为元件(或系统)的可靠性,系统由元件组成,常见的元件连接方式：,串联,并联,设 两系统都是由 4 个元件组成,每个元件正常工作的概率为 p,每个元件是否正常工作相互独立.两系统的连接方式如下图所示，比较两系统的可靠性.,S1:,S2:,五、EM算法及其MCMC方法,EM算法：是一种迭代方法，最初由Dempster等提出，并主要应用于较为复杂的后验分布，来计算后验均值或后验众数，即极大似然估计的一种数据添加算法。最大优点是简单和稳定。MCMC算法：当后验分布较为复杂时，对于后验分布的积分计算，像后验均值、后验方差、后验分布的分位数等等，就不得不求助于MCMC算法。他在统计物理学中得到广泛的应用，近年来，迅速发展到Bayes统计、显著性检验、极大似然估计等方面。关于这两方面的详细知识，感兴趣的可自己阅读！,