计算机技术-系统仿真monte carlo模拟

MonteCarlo模拟

内容提纲1.引言2.MonteCarlo模拟基本思想3.随机数生成函数4.应用实例举例5.排队论模拟6.MonteCarlo模拟求解规划问题引言(Introduction)MonteCarlo方法：蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法，statisticalsimulationmethod利用随机数进行数值模拟的方法MonteCarlo名字的由来：MonteCarlo是摩纳哥（monaco)的首都，该城以赌博闻名NicholasMetropolis(1915-1999)Monte-Carlo,MonacoMonteCarlo方法的基本思想

蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。源于美国在第二次世界大战研制原子弹的“曼哈顿计划”，该计划的主持人之一数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法，为它蒙上了一层神秘色彩。蒙特卡罗方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用蒲丰投针的方法来计算圆周率π，上世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。蒲丰投针实验：法国科学家蒲丰(Buffon)在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子。蒲丰投针实验的重要性并非是为了求得比其它方法更精确的π值，而是它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导。由此可以领略到从“概率土壤”上开出的一朵瑰丽的鲜花-蒙特卡罗方法(MC)蒲丰投针实验可归结为下面的数学问题：平面上画有距离为a的一些平行线，向平面上任意投一根长为l(l<a)的针，假设针落在任意位置的可能性相同，试求针与平行线相交的概率P(从而求π)蒲丰投针实验：如右图所示，以M表示针落下后的中点，以x表示M到最近一条平行线的距离，以φ表示针与此线的交角：针落地的所有可能结果满足：其样本空间视作矩形区域Ω,面积是:针与平行线相交的条件：它是样本空间Ω子集A，面积是：symslphi;int('l/2*sin(phi)',phi,0,pi)%ans=l因此，针与平行线相交的概率为：从而有：蒲丰投针实验的计算机模拟：formatlong;%设置15位显示精度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);

%样本空间Ωfori=1:nifx(i)<l*sin(phi(i))/2

%满足此条件表示针与线的相交

plot(phi(i),x(i),‘r.’);counter=counter+1;%统计针与线相交的次数frame(counter)=getframe;%描点并取帧endendfren=counter/n;

pihat=2*l/(a*fren)

%用频率近似计算πfigure(2)movie(frame,1)%播放帧动画1次一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929蒙特卡罗投点法是蒲丰投针实验的推广：

在一个边长为a的正方形内随机投点，该点落在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值，即n=10000;a=2;m=0;fori=1:nx=rand(1)*a;y=rand(1)*a;if((x-a/2)^2+(y-a/2)^2<=(a/2)^2)m=m+1;endenddisp(['投点法近似计算的π为:',num2str(4*m/n)]);xyo(a/2,a/2)基本思想

由上面的例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与之有关的量时，通过某种试验的方法，得出该事件发生的频率，再通过它得到问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法的关键步骤在于随机数的产生，计算机产生的随机数都不是真正的随机数(由算法确定的缘故)，如果伪随机数能够通过一系列统计检验，我们也可以将其当作真正的随机数使用。rand('seed',0.1);rand(1)％每次运行程序产生的值是相同的rand('state',sum(100*clock)*rand);rand(1)%每次重新启动matlab时，输出的随机数不一样注意:产生一个参数为λ的指数分布的随机数应输入exprnd(1/λ)产生m×n阶参数为A1,A2,A3的指定分布'name'的随机数矩阵random('name',A1,A2,A3,m,n)举例:产生2×4阶的均值为0方差为1的正态分布的随机数矩阵random('Normal',0,1,2,4)'name'的取值可以是(详情参见helprandom)：'norm'or'Normal'/'unif'or'Uniform''poiss'or'Poisson'/'beta'or'Beta''exp'or'Exponential'/'gam'or'Gamma''geo'or'Geometric'/'unid'or'DiscreteUniform'……非常见分布的随机数的产生逆变换方法Acceptance-Rejection方法

最早由VonNeumann提出，现在已经广泛应用于各种随机数的生成。

基本思路：通过一个容易生成的概率分布g和一个取舍准则生成另一个与g相近的概率分布f。

为要模拟服从给定分布的随机变量,用生成一个易于实现的不可约遍历链作为随机样本，使其平稳分布为给定分布的方法，称为马氏链蒙特卡罗方法.马氏链蒙特卡罗方法[1]连续型随机变量(以指数分布为例)：symstxlambda;Fx=int('lambda*exp(-lambda*t)',t,0,x)%分布函数symsr;Fxinv=finverse(Fx,x);%求反函数Fxinv=subs(Fxinv,x,r)%替换反函数变量x为rFxinv=inline(Fxinv)x=Fxinv(3,rand)%产生参数lambda=3指数分布的随机数%指数分布随机数产生函数已经提供exprnd(1/3,1,1)[2]离散型随机变量(以离散分布为例)：x=[2,4,6,8];px=[0.1,0.4,0.3,0.2];%以下为程序片段Fx=0;forn=1:length(px),Fx=[Fx,sum(px(1:n))];endr=rand;index=find(r<Fx);x(index(1)-1)%通用离散分布随机数产生程序例1

在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点．

经过长期观察发现，我方指挥所对敌方目标的指示有50％是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮．

现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。分析：这是一个概率问题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程.

为了能显示我方20次射击的过程，现采用模拟的方式。举例

需要模拟出以下两件事：

[2]当指示正确时，我方火力单位的射击结果情况[1]观察所对目标的指示正确与否模拟试验有两种结果，每一种结果出现的概率都是1/2．因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确．模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)．这时可用投掷骰子的方法来确定：如果出现的是１、２、３三个点：则认为没能击中敌人；如果出现的是４、５点：则认为毁伤敌人一门火炮；若出现的是６点：则认为毁伤敌人两门火炮．问题分析i：要模拟的打击次数；k1：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数；E：有效射击比率；E1：20次射击平均每次毁伤敌人的火炮数．符号说明模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜20?E=(k2+k3)/20E1=0*k1/20+1*k2/20+2*k3/20停止硬币正面?YNNY1，2，34，56模拟结果理论计算结果比较虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程．用蒙特卡洛方法进行计算机模拟的步骤：[1]设计一个逻辑框图，即模拟模型．[2]根据流程图编写程序，模拟随机现象．可通过具有各种概率分布的模拟随机数来模拟随机现象．[3]分析模拟结果，计算所需要结果.

随机投掷均匀硬币，验证国徽朝上与朝下的概率是否都是1/2

n=10000;%

给定试验次数m=0;fori=1:nx=randperm(2)-1;%randperm(2)returnsa随机排列oftheintegers1:2.

y=x(1);ify==0%0表示国徽朝上，1表示国徽朝下m=m+1;endendfprintf('国徽朝上的频率为：%f\n',m/n);小实例一：投掷硬币随机投掷骰子，验证各点出现的概率是否为1/6

n=10000;m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;fori=1:nx=randperm(6);y=x(1);switchycase1,m1=m1+1;case2,m2=m2+1;case3,m3=m3+1;case4,m4=m4+1;case5,m5=m5+1;otherwise,m6=m6+1;endend...%

输出结果小实例二：投掷骰子

设某班有m

个学生，则该班至少有两人同一天生日的概率是多少？小实例五：生日问题解：设一年为365天，且某一个学生的生日出现在一年中的每一天都是等可能的，则班上任意两个学生的生日都不相同的概率为：所以，至少有两个学生同一天生日的概率为：n=1000;p=0;m=50;%

设该班的人数为50fort=1:na=[];q=0;fork=1:mb=randperm(365);a=[a,b(1)];endc=unique(a);%去掉重复数字从小到大重新排序iflength(a)~=length(c)p=p+1;endendfprintf('至少两人同一天生日的概率为：%f\n',p/n);试验五源程序clear;m=50;p1=1:365;p2=[1:365-m,365*ones(1,m)];p=p1./p2;p=1-prod(p);fprintf('至少两人同一天生日的概率为：%f\n',p);小实例五的理论值计算排队问题随机模拟排队论主要研究随机服务系统的工作过程。在排队系统中，服务对象的到达时间和服务时间都是随机的。排队论通过对随机服务现象的统计研究，找出反映这些随机现象平均特性的规律指标，如排队的长度、等待的时间及服务利用率等。[1]系统的假设：（1）顾客源是无穷的；（2）排队的长度没有限制；（3）到达系统的顾客按先后顺序依次进入服务，“先到先服务”。在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客．当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店．设：1．顾客到来间隔时间服从参数为0.1的指数分布．２．对顾客的服务时间服从［４，１５］上的均匀分布．３．排队按先到先服务规则，队长无限制．

假定一个工作日为8小时，时间以分钟为单位。[1]模拟一个工作日内完成服务的个数及顾客平均等待时间t．[2]模拟100个工作日，求出平均每日完成服务的个数及每日顾客的平均等待时间。单服务员的排队模型模拟w：总等待时间；ci：第i个顾客的到达时刻；bi：第i个顾客开始服务时刻；ei：第i个顾客服务结束时刻；xi：第i-1个顾客与第i个顾客之间到达的间隔时间；yi：对第i个顾客的服务时间。符号说明c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci=ci-1+xiei=bi+yibi=max(ci,ei-1)t思路分析初始化：令i=1，ei-1=0，w=0产生间隔时间随机数xi~参数为0.1的指数分布ci=xi,bi=xi

产生服务时间随机数yi~[4,15]的均匀分布ei=bi+yi累计等待时间：w=w+bi-ci准备下一次服务：i=i+1产生间隔时间随机数xi~参数为0.1的指数分布ci=ci-1+xi

确定开始服务时间：bi=max(ci,ei-1)bi＞480?YNi=i-1,t=w/i输出结果：完成服务个数：m=i

平均等待时间：t停止流程框图用蒙特卡洛法解非线性规划问题基本假设试验点的第j个分量xj服从[aj,bj]内的均匀分布．符号假设求解过程先产生一个随机数作为初始试验点,以后则将上一个试验点的第j个分量随机产生，其它分量不变而产生一新的试验点．这样，每产生一个新试验点只需一个新的随机数分量．当K>MAXK或P>MAXP时停止迭代．框图初始化:给定MAXK,MAXP;k=0,p=0,Q:大整数xj=aj+R(bj-aj)j=1,2,…,nj=0j=j+1,p=p+1P>MAXP?YNxj=aj+R(bj-aj)gi(X)≥0?i=1,2…nYNj<n?f(X)≥Q?YNX*=X,Q=f(X)k=k+1k>MAXK?YN输出X,Q,停止YN

在Matlab软件包中编程，共需三个Ｍ－文件:r