第6讲蒙特卡罗方法与计算机模拟(第1次课)

1、第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟1第六讲 蒙特卡罗与计算机模拟内容：内容：计算机模拟计算机模拟(或称仿真或称仿真)是一种广义数值计算是一种广义数值计算 方法，适合解决一些规模大、难以解析化方法，适合解决一些规模大、难以解析化 以及受随机因素影响的不确定数学模型以及受随机因素影响的不确定数学模型目的：目的：了解蒙特卡罗方法的基本思想，掌握利用了解蒙特卡罗方法的基本思想，掌握利用 Matlab对离散对离散/连续系统进行模拟的方法连续系统进行模拟的方法要求：要求：掌握掌握Matlab随机数函数，处理应用问题随机数函数，处理应用问题了解蒙特卡罗方法的起源和基本思想了解蒙特

2、卡罗方法的起源和基本思想了解蒲丰投针实验设计思路和计算机模拟了解蒲丰投针实验设计思路和计算机模拟掌握掌握随机数随机数函数函数 rand unifrnd normrnd exprnd掌握掌握离散系统离散系统和和连续系统连续系统计算机模拟实例计算机模拟实例第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟2蒙特卡罗方法的起源和基本思想 蒙特卡罗方法蒙特卡罗方法(Monte Carlo method)，或称计，或称计算机随机模拟方法，是一种基于算机随机模拟方法，是一种基于“随机数随机数”的计的计算方法。源于美国在第算方法。源于美国在第二二次世界大战研制原子弹次世界大战研制原子弹的的“曼哈

3、顿计划曼哈顿计划”，该计划的主持人之一数学家该计划的主持人之一数学家冯冯诺伊曼诺伊曼用驰名世界的用驰名世界的赌城赌城摩纳哥的摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色来命名这种方法，为它蒙上了一层神秘色彩。彩。 蒙特卡罗蒙特卡罗方法的方法的基本思想基本思想很早以前就被人们所很早以前就被人们所发现和利用。早在发现和利用。早在17世纪，人们就知道世纪，人们就知道用事件发用事件发生的生的“频率频率”来决定事件的来决定事件的“概率概率”。19世纪人世纪人们用们用蒲丰蒲丰投针投针的方法来的方法来计算计算圆周率圆周率，上，上世纪世纪40年年代电子计算机的出现，特别是近年来代电子计算

4、机的出现，特别是近年来高速电子计高速电子计算机算机的出现，使得用数学方法在计算机上大量、的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。快速地模拟这样的试验成为可能。第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟3蒲丰投针实验近似计算圆周率蒲丰投针实验蒲丰投针实验： 法国科学家蒲丰法国科学家蒲丰(Buffon)在在1777年提年提出的蒲丰投针实验是早期几何概率一个出的蒲丰投针实验是早期几何概率一个非常著名的例子。非常著名的例子。蒲丰投针实验蒲丰投针实验的重要的重要性并非是为了求得比其它方法更精确的性并非是为了求得比其它方法更精确的值，而是它值，而是它开创了开

5、创了使用使用随机数处理确定性数学问随机数处理确定性数学问题题的先河，是用的先河，是用偶然性方法去解决确定性计算偶然性方法去解决确定性计算的的前导，前导，由此由此可以可以领略到从领略到从“概率土壤概率土壤”上开出的上开出的一朵瑰丽的鲜花一朵瑰丽的鲜花蒙特卡罗蒙特卡罗方法方法(MC) 蒲丰投针实验蒲丰投针实验可归结为下面的数学问题可归结为下面的数学问题：平面：平面上画有距离为上画有距离为a的一些平行线，向平面上任意投一的一些平行线，向平面上任意投一根长为根长为l(la)的针，假设针落在任意位置的可能性的针，假设针落在任意位置的可能性相同，试求相同，试求针与平行线相交的概率针与平行线相交的概率P(从

6、而求从而求)第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟4蒲丰投针实验近似计算圆周率蒲丰投针实验蒲丰投针实验： 如右图所示，以如右图所示，以M表示针落下后的中点，表示针落下后的中点，以以x表示表示M到最近一条到最近一条平行线的距离，以平行线的距离，以表示针与此线的交角：表示针与此线的交角：针落地的所有可能结果满足：针落地的所有可能结果满足：其样本空间视作矩形区域其样本空间视作矩形区域, 面积是面积是:针与平行线相交的条件：针与平行线相交的条件：它是样本空间它是样本空间子集子集A，面积是：，面积是：积分计算积分计算 syms l phi; int(l/2*sin(phi),ph

7、i,0,pi); %ans=l因此，针与平行线相交的概率为：因此，针与平行线相交的概率为：从而有：从而有： 特别当特别当 时时 p为统计频为统计频率率02,0 xa( )2Sa 0sin2,0 xl0( )sin2S Al dl( )/ ( ) 2 /p S A Sl a 2l ap2al1 p第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟5蒲丰投针实验近似计算圆周率蒲丰投针实验蒲丰投针实验的计算机模拟：的计算机模拟：format long; a=1; l=0.6; %显示精度显示精度, 线线宽宽和和针长针长figure; axis(0,pi,0,a/2); %初始化绘图板初始

8、化绘图板set(gca,nextplot,add); %初始化绘图方式为叠加初始化绘图方式为叠加counter=0; n=2010; %初始化计数器和设定初始化计数器和设定投针次数投针次数x=unifrnd(0,a/2,1,n); phi=unifrnd(0,pi,1,n); %样本空间样本空间for i=1:n if x(i)l*sin(phi(i)/2 %满足此条件表示针与线的相交满足此条件表示针与线的相交 plot(phi(i),x(i),r.); frame(i)=getframe; %描点并取帧描点并取帧 title(Current Point ,num2str(i), Total

9、,num2str(n); counter=counter+1; %统计针与线相交的次数统计针与线相交的次数 endendfren=counter/n; pihat=2*l/(a*fren) %用用频率频率近似计算近似计算%movie(frame,1) %播放帧动画播放帧动画1次次第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟6蒲丰投针实验近似计算圆周率蒲丰投针实验蒲丰投针实验的历史记录：的历史记录： 意大利数学家拉泽里尼得到了准确到意大利数学家拉泽里尼得到了准确到6位小数位小数的的值，不过他的实验因为值，不过他的实验因为太准确太准确而受到了质疑而受到了质疑第六讲第六讲 蒙特卡罗

10、方法蒙特卡罗方法与计算机模拟与计算机模拟7蒲丰投针实验计算圆周率蒙特卡罗投点法蒙特卡罗投点法是是蒲丰投针实验蒲丰投针实验的推广：的推广： 在一个边长为在一个边长为a的正方形内随机投点，该点落的正方形内随机投点，该点落在此正方形的内切圆中的概率应为该内切圆与正在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值，即方形的面积比值，即n=10000; a=2; m=0; for i=1:n x=rand(1)*a; y=rand(1)*a; if ( (x-a/2)2+(y-a/2)2 = (a/2)2 ) m=m+1; endenddisp(投点法近似计算的投点法近似计算的为: ,num2s

11、tr(4*m/n);22/2:/4aaxyo(a/2,a/2)第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟8常见分布的随机数产生语句 蒙特卡罗方法的蒙特卡罗方法的关键步骤关键步骤在于在于随机数随机数的产生，的产生，计算机产生的随机数都不是计算机产生的随机数都不是真正的随机数真正的随机数(由算法由算法确定的缘故确定的缘故)，如果，如果伪随机数伪随机数能够通过一系列统计能够通过一系列统计检验，我们也可以将其检验，我们也可以将其当作真正的随机数当作真正的随机数使用：使用： 第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟9常见分布的随机数产生语句MATLAB可以可以直

12、接产生直接产生满足各种分布的随机数满足各种分布的随机数具体具体命令如下命令如下: 产生产生mn阶阶0,1上上均匀分布均匀分布的随机数矩阵的随机数矩阵 rand(m,n)产生一个产生一个0,1上上均匀分布均匀分布的随机数的随机数 rand 产生产生mn阶阶a,b上上均匀分布均匀分布的随机数矩阵的随机数矩阵 unifrnd (a,b,m, n) 产生一个产生一个a,b上上均匀分布均匀分布的随机数的随机数 unifrnd(a,b) 产生一个产生一个1:n的随机排列的随机排列(元素均出现且不重复元素均出现且不重复) p=randperm(n)注意注意: randperm(6)与与unifrnd (1,

13、6,1, 6)的区别的区别第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟10常见分布的随机数产生语句 产生产生mn阶阶均值为均值为mu方差为方差为sigma的的正态正态分布分布的随机数矩阵的随机数矩阵 normrnd(mu,sigma,m,n)产生一个均值为产生一个均值为mu方差为方差为sigma的的正态分布正态分布的随机的随机数数 normrnd(mu,sigma) 产生产生mn阶期望值为阶期望值为mu (mu=1/)的的指数分布指数分布的随机数矩阵的随机数矩阵 exprnd(mu,m,n)产生一个期望值为产生一个期望值为mu的指数分布的随机数的指数分布的随机数 exprnd

14、(mu)注意注意: 产生一个参数为产生一个参数为的的指数分布的随机数应输入指数分布的随机数应输入 exprnd(1/)e0( )00txf xx第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟11常见分布的随机数产生语句 产生产生mn阶阶参数为参数为A1,A2,A3的指定分布的指定分布name的的随机数矩阵随机数矩阵 random(name,A1,A2,A3,m,n)产生一个参数为为产生一个参数为为A1,A2,A3的指定分布的指定分布name的的随机数随机数 random(name,A1,A2,A3)举例举例: 产生产生24阶的均值为阶的均值为0方差为方差为1的正态分布的的正态分

15、布的随机数矩阵随机数矩阵 random(Normal,0,1,2,4)name的取值可以是的取值可以是(详情参见详情参见help random)：norm or Normal / unif or Uniformpoiss or Poisson / beta or Betaexp or Exponential / gam or Gammageo or Geometric / unid or Discrete Uniform第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟12MATLAB随机数的“重置”问题 Matlab的随机数是的随机数是伪随机数伪随机数，但在一定的信度，但在一定的

16、信度之下可以看作真正的随机数。问题是之下可以看作真正的随机数。问题是rand函数函数产产生的随机数从一个随机数序列中取出来生的随机数从一个随机数序列中取出来，而每次，而每次启动启动Matlab时，时，rand的状态都会被重置的状态都会被重置(相当于把相当于把序列的指针移到了随机数序列的开始序列的指针移到了随机数序列的开始)，换言之第，换言之第一次启动一次启动Matlab调用的第调用的第n次次rand函数与下一次启函数与下一次启动调用的第动调用的第n个个rand函数函数产生相同的数值产生相同的数值。 如果想如果想打乱这种状态打乱这种状态，可以为，可以为rand指定一个指定一个与与当前时间相关当前

17、时间相关的初始状态，而不用默认状态：的初始状态，而不用默认状态：rand(state,sum(100*clock);或者或者rand(state,sum(100*clock)*rand);第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟13非常见分布的随机数的产生 对于常见分布随机数，可由相应对于常见分布随机数，可由相应Matlab函数直函数直接产生，对于非常见分布随机数可如下处理：接产生，对于非常见分布随机数可如下处理：1 连续型随机变量连续型随机变量(以以p116指数分布为例指数分布为例)：e01-e0( )( )( )0000txxtxf tf t dt F xtx 概率密

18、度函数概率密度函数: :反映随机变量出现在某个反映随机变量出现在某个t*位置的可能性位置的可能性f(t*)( (概概率率) ),直观上体现为稠密程度直观上体现为稠密程度( (概率密度概率密度); ); 累积分布函数累积分布函数: :对概率密度函数对概率密度函数积分积分得到得到累积分布函数累积分布函数, ,反映出现反映出现小于等于某个随机变量小于等于某个随机变量x*的所有可能性之和的所有可能性之和F(x*).常用累积分布常用累积分布反函数反函数生成服从该分布的随机变量生成服从该分布的随机变量. .第六讲第六讲 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟14非常见分布的随机数的产生syms t x lambda;Fx=int(lambda*exp(-lambda*t),t,0,x) %累积分布函数累积分布函数syms r;Fxinv=finverse(Fx,x); %求累积分布函数反函数求累积分布函数反函数Fxinv=subs(Fxinv,x,r) %替换反函数变量替换反函数变量x为为rFxinv=inline(Fxinv)x=Fxinv(3,rand) %产生参数产生参数 lambda=3 指数分布的随机数指数分布的随机数%指数分布随机数产生函数已经