一蒙特卡洛随机模拟

系列一蒙特卡洛随机模拟实验目的：学会用计算机随机模拟方法来解决随机性问题蒙特卡洛模拟法简介蒙特卡洛(MonteCarlo)方法是一种应用随机数来进行计算机摸拟的方法。此方法对研究对象进行随机抽样，通过对样本值的观察统计，求得所研究系统的某些参数。作为随机模拟方法，起源可追溯到18世纪下半叶蒲峰实验。蒙特卡洛模拟法的应用领域蒙特卡洛模拟法的应用领域主要有：直接应用蒙特卡洛模拟：应用大规模的随机数列来模拟复杂系统，得到某些参数或重要指标。蒙特卡洛积分：利用随机数列计算积分，维数越高，积分效率越高。蒙特卡洛模拟法求解步骤应用此方法求解工程技术问题可以分为两类：确定性问题和随机性问题。解题步骤如下：1.根据提出的问题构造一个简单、适用的概率模型或随机模型，使问题的解对应于该模型中随机变量的某些特征(如概率、均值和方差等)，所构造的模型在主要特征参量方面要与实际问题或系统相一致2.根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，方可进行随机模拟试验。根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样(包括直接抽样、分层抽样、相关抽样、重要抽样等)。按照所建立的模型进行仿真试验、计算，求出问题的随机解。统计分析模拟试验结果，给出问题的概率解以及解的精度估计。在可靠性分析和设计中，用蒙特卡洛模拟法可以确定复杂随机变量的概率分布和数字特征，可以通过随机模拟估算系统和零件的可靠度，也可以模拟随机过程、寻求系统最优参数等。一.预备知识：随机数的产生提示：均匀分布U(0,1)的随机数可由C语言或Matlab自动产生，在此基础上可产生其他分布的随机数.逆变换法：设随机变量U服从(0,1)上的均匀分布，则X=F-'(U)的分布函数为F(x)步骤：(1)产生U(0J)的随机数U；②计算X=F-1(U),则X服从F(x)分布.问题：练习用此方法产生常见分布随机数例如“指数分布，均匀分布U(a,b)”.还有其它哪种常见分布的随机数可用此方法方便产生？产生离散分布随机数己知离散随机变量X的概率分布:P(X=xk)=I\,(K=1,2…)，产生随机变量X的随机数可采用如下算法：将区间［0.1］依次分为长度为Pi,p?,・••的小区间L，L，・••；产生［0,1］均匀分布随机数R,若Rclk则令X=xk,重复(b),即得离散随机变量X的随机数序列.问题：(1)下表给出了离散分布X的概率分布表，试产生100个随机数X的概率分布表Xk13579Pk0.50.050.20.050.2用此方法给出100个二项分布B(20,0.1)的随机数及10个泊松分布P(l)的随机数.正态分布的抽样提示：设U］,U2是独立同分布的U(0Q变量，令X］=(-21nU］)”2cos(2^u2)X2=(-21nU1)1/2sin(2MJ2)则X.与X,独立，均服从标准正态分布.步骤：(1)由U(0J)独立抽取Ui=g=U2(2)用(*)式计算^,X2.用此方法可同时产生两个标准正忐分布的随机数问题：有关随机数产生方法很多，查阅相关材料进行系统总结.二.随机决策问题某小贩每天以一元的价格购进一种鲜花，卖出价为b元/束，当天卖不出去的花全部损失，顾客一天内对花的需求量是随机变量，服从泊松分布，P(X=k)=e-4—,k=0,1,2,...,,其中常数;I由多口销传量的平均值来估计，问小贩每天应购进多少束鲜花?(准则:期望收入，(①最局)问题：(1)在给定b=1.25,2=50的值后，画出目标函数S(u)连线散点图，观察单调性，给出最优决策U*：。选取其他的b,4,再观察S(u)的单调性；用计算机模拟方法来求出最优决策对固定的U,例如，u=40,对随机变量X模拟100次，每次模拟得到一个收入，求出100个收入的平均值，即得到在决策u=40情况下的可能收入；对所有的可能的u,重复(3),从中找最大的，并与⑴的结果相比较.蒙特卡罗计算兀值思路：图1表示一个内接于正方形的圆的半径R.圆的面积是刀R，,正方形的面积是(2R)2.圆和正方形的面枳的比值就是将上述比值乘以4,就能获得〃值.(2R)-4

2R2R图1例用蒙特卡罗方法计算力值一重定积分的蒙特卡罗算法问题描述：假设函数f(x)在［a,b］内有界连续，且f(x)ZO,求解定积分I=「f(x)dxJa为计算出其值，可构造概率模型如下：取一个边长分别为b-a和c的矩形D,使曲边梯形在矩形域之内，如图2,并在矩形内随机投点，假设随机点均匀地落在整个矩形之内,则落在图中灰色区域内的随机点数k与投点总数N之比k/N就近似地等于曲线下方面枳(即阴影面枳)与矩形面积之比，从而得出近似积分a)c京京e1.——求例由于厂子是非初等函数，我们很难求出其原函数，所以用牛顿■莱布尼茨公式无法求解，但可以运用蒙特卡罗方法求出其近似值.将上述方法推广到一般情况：假设函数f(x)在［a,b］内有界连续，对于定枳分I=［f(x)dx，为计算出其值，可构造如下•概率模型:取一个边长分别为b-a和c-d的矩形D,使曲线［a,b］段的值在矩形域之内，如图3,并在矩形内随机投点，假设随机点均匀地落在整个矩形之内，则落在图中x轴上下灰色区域内的随机点数m与n的差与投点总数p之比(m-n)/P就近似地等于曲线上下方面枳之差(即阴影面积之差)与矩形面枳之比，从而得出近似积分I/%2(b-a)(c-d).图3二重积分的蒙特卡罗算法问题描述：实际计算中常常要遇到如j］f(x,y)dxdy的二重积分，发现被积函数的原函数往*D往很难求出，或者原函数根本就不是初等函数，对于这样的重积分，蒙特卡罗方法也有成熟的计算方法.方法1:步骤：取一个包含D的矩形区域Q：a<x<b,c<y<d,面积A=(b-a)(d-c)；(无,乂)，i=1,2,・••,n,为Q上的均匀分布随机数列，不妨设(&乂),(i= 为落Ak在D中的ii个随机数，则11充分大时，有f(x,y)dxdyb—习f(*,y)d n1=1方法2：对二重积分1= 假设f(&y)为区域A上的有界函数，且f(x,y)20,几何意义对应的是以f(x,y)为曲面顶，A为底的曲顶柱体C的体枳.因此，用均匀随机数计算二重枳分的蒙特卡罗方法基本思路为：假设曲顶柱体C包含在己知体积为Vd的几何体D的内部，在D内产生N个均匀随机点，统计出在C内部的随机点数目Nc,则I=^Nc..

例：计算JJ(1+小一x^-y2一jL+y^dxdy,其中A={(x,y)|x2+y2<1}A分析：该二重积分可以看作以l+小一妒—y2-jV+y?为顶的曲顶主体的体积,此曲顶柱体在一个边长为2的立方体内，用数学分析方法可计算出其精确值为兀.5用蒙特卡洛方法计算体积z>y/x2+y2,z<1+1—x2—y2冰激凌锥含于体积=8的六面体Q={(x,y,z)-1<x<l,-l<y<l,0<z<2)蒙特卡罗求一元函数的最值数学分析中求一元函数f(x)在闭区间［a,b］的最值的方法为：找出函数f(x)的所有驻点和不可导点，然后计算出其相应函数值，并将其与区间端点函数值比较得出函数的最大值或最小值，但实际问题中方程f(x)=0很难求出其根，所以驻点往往很难求出.若利用蒙特卡罗方法，在［a,b］上产生n个随机数，计算出这些数的函数值并作比较,则可得出最值的近似值.例求f(x)=(l-x3)sin3x在［-盆,2勿］上的最大值-6-4 -2 0 2 4 6x(1-x°)sin(3x)-6-4 -2 0 2 4 6x(1-x°)sin(3x)150100<-v=(1-x3<-v=(1-x3)*sin(3*x)0-50-100图4蒙特卡罗方法求一元函数的根根据数学分析相关知识，求一元函数f(x)=0的根有很多方法，如一般迭代法，牛顿切线法等等，这些方法讨论根的收敛性，与初始迭代值相密切相关通常要给出一个合适的初始迭代值知是比较困难的，利用蒙特卡罗方法可以摆脱根的收敛性对初始值Xg的依赖性.具体方法如下：要解方程f(x)=0,XG［a,b］,其中函数f(x)为连续函数，为指定精度令Xd=m=a,X*=b,k=l,F_(1)=FO,步骤如下：(1)^k=k+l,在［X顷，担］内产生n个随机数*(i=l,2,...n)，计算并比较出这n个随机数的函数绝对值的最小值f(&J=得事|1'(而)|｝，业为i的某个取值，令Ff(k)=min｛Fxi林(k-l),f(&)｝⑵若写_«)<言成立，则终止计算，令。硕=跖，根就是捻ot；若Fx^(k)>^fiF_(k)=F_(k-l),则令k=k—1,转至(1)；若E_(k)>S且％颇«)<岌心(卜1)，则令。皿=自，转至(3).(3)令d=cIq/k,Xdami=fgjt—d，纨？=。硕+d，转至(1)•说明：a)F°是人为给定的一个很大的正数，do«b-afido>Ob)k=k+1表示重新赋值给k,使k的值增加1,对k=k-1同理.C)区间［Xm'XJ一定在定义域［a,b］之内.此迭代步骤能使函数值序列Fxmi(l>Fx >.Fxmk,最终使岌心时成立，得出函数f(x)=0达到精度要求的根板心例求f(x)=e-°—tanx+800=0在(0,乙)上的实根，(假定|f(x)|<10'\则认为x为根)2分析：对上面的方程，如取初值05,0.7,0.9,1.1,1.3,1.5,用牛顿迭代公式在区间内找不到根，若用蒙特卡罗方法，则不需要给定初值.蒙特卡罗方法解线性规划用蒙特卡罗方法可解约束规划问题：minf(X),XeEngi(X)Z0,i=l,2,...,ms.t.( ,・■■a,Sbj,j=l,2,・..,n基本思想：在估计的区域｛(改，也，…，％Xj£［aj,bj］,j=1,2,...n)｝内随机取若干个试验点，然后从试验点中找出可行点，再从可行点中选择出使函数值最小的点.符号假设和说明：设试验点的第j个分量Xj服从［apbj内的均匀分布；P：试验点总数；MAX{P)：最大试验点总数；K：可行点总数：MAX{K}：最大可行点数；X*：迭代产生的最优点；R：在［0,1］上的均匀随机数：Q：迭代产生的最小值f(X*),其初始