第六讲 蒙特卡罗方法 与计算机模拟1第六讲 蒙特卡罗与计算机模拟内容：计算机模拟(或称仿真)是一种广义数值计算方法，适合解决一些规模大、难以解析化以及受随机因素影响的不确定数学模型目的：了解蒙特卡罗方法的基本思想，掌握利用Matlab对离散/连续系统进行模拟的方法 要求：掌握Matlab随机数函数，处理应用问题Ø了解蒙特卡罗方法的起源和基本思想Ø了解蒲丰投针实验设计思路和计算机模拟Ø掌握随机数函数 rand unifrnd normrnd exprndØ掌握离散系统和连续系统计算机模拟实例第六讲 蒙特卡罗方法 与计算机模拟2蒙特卡罗方法的起源和基本思想蒙特卡罗方法(Monte Carlo method)，或称计 算机随机模拟方法，是一种基于“随机数”的计算方 法源于美国在第二次世界大战研制原子弹的“曼 哈顿计划”，该计划的主持人之一数学家冯·诺伊曼 用驰名世界的赌城—摩纳哥的Monte Carlo—来命 名这种方法，为它蒙上了一层神秘色彩蒙特卡罗方法的基本思想很早以前就被人们所 发现和利用早在17世纪，人们就知道用事件发生 的“频率”来决定事件的“概率”19世纪人们用蒲丰 投针的方法来计算圆周率π，上世纪40年代电子计 算机的出现，特别是近年来高速电子计算机的出现 ，使得用数学方法在计算机上大量、快速地模拟这 样的试验成为可能。
第六讲 蒙特卡罗方法 与计算机模拟3蒲丰投针实验近似计算圆周率π蒲丰投针实验：法国科学家蒲丰(Buffon)在1777年提 出的蒲丰投针实验是早期几何概率一个 非常著名的例子蒲丰投针实验的重要 性并非是为了求得比其它方法更精确的 π值，而是它开创了使用随机数处理确定性数学问 题的先河，是用偶然性方法去解决确定性计算的前 导，由此可以领略到从“概率土壤”上开出的一朵瑰 丽的鲜花——蒙特卡罗方法(MC)蒲丰投针实验可归结为下面的数学问题：平面 上画有距离为a的一些平行线，向平面上任意投一 根长为l(la)的针，假设针落在任意位置的可能性 相同，试求针与平行线相交的概率P(从而求π)第六讲 蒙特卡罗方法 与计算机模拟4蒲丰投针实验近似计算圆周率π蒲丰投针实验：如右图所示，以M 表示针落下后的中点， 以x表示M到最近一条 平行线的距离，以φ表示针与此线的交角： 针落地的所有可能结果满足： 其样本空间视作矩形区域Ω, 面积是: 针与平行线相交的条件： 它是样本空间Ω子集A，面积是： 积分计算 syms l phi; int('l/2*sin(phi)',phi,0,pi); %ans=l 因此，针与平行线相交的概率为： 从而有： 特别当 时 p为统计频率第六讲 蒙特卡罗方法 与计算机模拟5蒲丰投针实验近似计算圆周率π蒲丰投针实验的计算机模拟： format long; a=1; l=0.6; %显示精度, 线宽和针长 figure; axis([0,pi,0,a/2]); %初始化绘图板 set(gca,'nextplot','add'); %初始化绘图方式为叠加 counter=0; n=2010; %初始化计数器和设定投针次数 x=unifrnd(0,a/2,1,n); phi=unifrnd(0,pi,1,n); %样本空间Ω for i=1:nif x(i)l*sin(phi(i))/2 %满足此条件表示针与线的相交plot(phi(i),x(i),'r.'); frame(i)=getframe; %描点并取帧title(['Current Point ',num2str(i),' Total ',num2str(n)]);counter=counter+1; %统计针与线相交的次数end end fren=counter/n; pihat=2*l/(a*fren) %用频率近似计算π %movie(frame,1) %播放帧动画1次第六讲 蒙特卡罗方法 与计算机模拟6蒲丰投针实验近似计算圆周率π蒲丰投针实验的历史记录： 意大利数学家拉泽里尼得到了准确到6位小数 的π值，不过他的实验因为太准确而受到了质疑第六讲 蒙特卡罗方法 与计算机模拟7蒲丰投针实验计算圆周率π蒙特卡罗投点法是蒲丰投针实验的推广：在一个边长为a的正方形内随机投点，该点落在 此正方形的内切圆中的概率应为该内切圆与正方形 的面积比值，即 n=10000; a=2; m=0; for i=1:nx=rand(1)*a; y=rand(1)*a;if ( (x-a/2)^2+(y-a/2)^2 = (a/2)^2 )m=m+1;end enddisp(['投点法近似计算的π为: ',num2str(4*m/n)]);xyo(a/2,a/2)第六讲 蒙特卡罗方法 与计算机模拟8常见分布的随机数产生语句蒙特卡罗方法的关键步骤在于随机数的产生， 计算机产生的随机数都不是真正的随机数(由算法 确定的缘故)，如果伪随机数能够通过一系列统计 检验，我们也可以将其当作真正的随机数使用： 第六讲 蒙特卡罗方法 与计算机模拟9常见分布的随机数产生语句MATLAB可以直接产生满足各种分布的随机数 具体命令如下: ① 产生m×n阶[0,1]上均匀分布的随机数矩阵rand(m,n) 产生一个[0,1]上均匀分布的随机数rand ② 产生m×n阶[a,b]上均匀分布的随机数矩阵unifrnd (a,b,m, n) 产生一个[a,b]上均匀分布的随机数unifrnd(a,b) ③ 产生一个1:n的随机排列(元素均出现且不重复)p=randperm(n) 注意: randperm(6)与unifrnd (1,6,1, 6)的区别第六讲 蒙特卡罗方法 与计算机模拟10常见分布的随机数产生语句④ 产生m×n阶均值为mu方差为sigma的正态分布 的随机数矩阵normrnd(mu,sigma,m,n) 产生一个均值为mu方差为sigma的正态分布的随机 数 normrnd(mu,sigma) ⑤ 产生m×n阶期望值为mu (mu=1/λ)的指数分布的 随机数矩阵exprnd(mu,m,n) 产生一个期望值为mu的指数分布的随机数exprnd(mu) 注意: 产生一个参数为λ的 指数分布的随机数应输入 exprnd(1/λ)第六讲 蒙特卡罗方法 与计算机模拟11常见分布的随机数产生语句⑥ 产生m×n阶参数为A1,A2,A3的指定分布'name' 的随机数矩阵random('name',A1,A2,A3,m,n) 产生一个参数为为A1,A2,A3的指定分布'name'的随 机数 random('name',A1,A2,A3) 举例: 产生2×4阶的均值为0方差为1的正态分布的 随机数矩阵 random('Normal',0,1,2,4) 'name'的取值可以是(详情参见help random)： 'norm' or 'Normal' / 'unif' or 'Uniform' 'poiss' or 'Poisson' / 'beta' or 'Beta' 'exp' or 'Exponential' / 'gam' or 'Gamma' 'geo' or 'Geometric' / 'unid' or 'Discrete Uniform' ……第六讲 蒙特卡罗方法 与计算机模拟12MATLAB随机数的“重置”问题Matlab的随机数是伪随机数，但在一定的信度 之下可以看作真正的随机数。
问题是rand函数产生 的随机数从一个随机数序列中取出来，而每次启动 Matlab时，rand的状态都会被重置(相当于把序列 的指针移到了随机数序列的开始)，换言之第一次 启动Matlab调用的第n次rand函数与下一次启动调 用的第n个rand函数产生相同的数值如果想打乱这种状态，可以为rand指定一个与 当前时间相关的初始状态，而不用默认状态： rand('state',sum(100*clock));或者 rand('state',sum(100*clock)*rand);第六讲 蒙特卡罗方法 与计算机模拟13非常见分布的随机数的产生对于常见分布随机数，可由相应Matlab函数直 接产生，对于非常见分布随机数可如下处理： [1] 连续型随机变量(以p116指数分布为例)：概率密度函数: 反映随机变量出现在某个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为r Fxinv=inline(Fxinv) x=Fxinv(3,rand) %产生参数 lambda=3 指数分布的随机数 %指数分布随机数产生函数已经提供 exprnd(1/3,1,1)第六讲 蒙特卡罗方法 与计算机模拟15非常见分布的随机数的产生[2] 离散型随机变量(以p117离散分布为例)：x=[2,4,6,8]; px=[0.1,0.4,0.3,0.2]; %以下为程序片段 Fx=0; for n=1:length(px), Fx=[Fx,sum(px(1:n))]; end r=rand; index=find(rFx); x(index(1)-1) %已编写通用离散分布随机数产生程序 scatrnd(x,px,n)。