【第四讲】随机模拟讲义稿

1、第四讲 随机模拟如果认为计算机在数学建模中的作用仅仅是根据已经建立的模型，代入数据算出最后结果的话，那就错了。其实计算机可以在建模阶段帮助我们选捧合理的模型，对模型进行定量的测试从而帮助我们评价和改进。计算机模拟中的蒙特卡罗方法又称随机抽样技巧或统计试验方法。半个多世纪以来，由于科学技术的发展和电子计算机的发明 ，这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域

2、日趋广泛。4.1相关原理及计算实验1.模拟法分类（1）运筹对策法：主要用于军事对策和企业管理对策。如现代化战争的军事演习、新式武器的试验等。最早于 40 年代末美国纽曼等人首先用运筹模拟法解决了核屏蔽实验问题。（2）蒙特卡罗法：蒙特卡罗方法又称随机抽样技巧或统计试验方法，与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。例1：设总计投了M根针，落入阴影部分N根，则不规则图形的面积为11针在平行线间的位置 图（4.1） 图（4.2）例2. （蒲丰投针实验） 为了求得圆周率值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（ la

3、）的平行线相交的频率代替概率P，再利用准确的关系式：求出值其中为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。一些人进行了实验，其结果列于下表 ：表4-1实验者年份投计次数的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929解：设针投到地面上的位置可以用一组参数（x,）来描述，x为针中心的坐标，为针与平行线的夹角，如图所示。任意投针，就是意味着x与都是任意取的，但x的范围限于0，a，夹角的范围限于0，。在此情况下，

4、针与平行线相交的数学条件是如何产生任意的（x,）？x在0，a上任意取值，表示x在0，a上是均匀分布的，其分布密度函数为：类似地，的分布密度函数为：因此，产生任意的（x,）的过程就变成了由f1(x)抽样x及由f2()抽样的过程了。由此得到：其中1，2均为（0,1）上均匀分布的随机变量。每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,），然后定义描述针与平行线相交状况的随机变量s(x,)，为 如果投针次，则 是针与平行线相交概率的估计值。事实上， 于是有 由以上例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通

5、过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。 （3）系统模拟法：是用数字对含有随机变量的系统进行模拟，可看作是蒙特卡洛法的应用。一般说来，蒙特卡洛法用于静态计算，而系统模拟法用于动态模型计算。我们主要讨论此法。 2.随机数的产生（1）0,1区间上均匀分布随机数的产生在连续型随机变量的分布中，最简单且最基本的分布是单位均匀分布。定义:设为0,1上服从均匀分布的随机变量，即的分布密度函数与分布函数分别为：， 则的样本值，即以即以等概率取自0,1的一串数称为0,1上均匀分布的随机数。 随机数的产生方法主要有以下

6、几种：物理方法：一是放射性物质随机蜕变；二是电子管回路的热噪声。（如可将热噪声源装于计算机外部，按其噪声电压的大小表示不同的随机数。此法产生的随机性最好，但产生过程复杂。）查随机数表（Rand Table）（1955年由美国兰德公司编制，有随机数100万个。）随机数表中的数字具有均匀的随机性，没有周期性。使用时，可根据需要任取一段（横或竖）。如需20个，便可从中取（顺次）20个，要几位取几位，随机数表无所谓位数，不能四舍五入。由递推公式（如同余数公式）在计算机内产生伪随机数：由于第i+1个随机数是由第i个按一定公式推算出来的，故并非真正的随机数。（2）任意概率分布随机数的产生以上介绍了均匀分布

7、R 的随机数的产生方法，那么任意分布X的随机数如何产生？ 我们说，X的随机数可以利用得到。那么X与R间的关系又是什么？ 定理：设R是服从0，1区间上均匀分布的随机变量，X的分布函数为，则。例3：利用0，1区间均匀分布的随机数表示服从负指数分布的随机数。 解：设X 服从负指数分布，则，由，知 所以即为所求。例4：求a,b区间上的均匀分布的随机数。解：设X 服从a,b区间上的均匀分布，则，令，知即为所求。4.2 模拟实例例1：某生产电子产品的企业，要对某型号的产品平均无故障运行时间做出估计。该产品由A、B、C三个部件串联而成。因此，当这三个部件中任何一个部件发生故障而失效时，则该电子产品也即告失效

8、。如果根据该产品投入运行后再对其无故障运行时间做出估计，则费用较高。现在企业已经得到每一种部件的有关运行试验记录资料，其中包括用来确定部件失效时间的概率分布。表4-2 部件失效概率分布和随机数取值表A部件B部件C部件失效时间概率随机数失效时间概率随机数失效时间概率随机数40.1011020.05010560.2012050.2113030.1061570.3215060.3316040.2163580.25517570.2618050.3366580.15769080.15819560.256690100.1910090.05960070.19100表4-3产品失效时间仿真表序号ABC产品失效

9、时间随机数失效时间随机数失效时间随机数失效时间1336244528462506726859313519479944828204869455969177286630588620657245957126580243852174915541599105103865155885111250830463128582343674139285550165147972748494155968061366161152640664179795451565183964757385经计算产品失效时间为4.6小时。例2：（存贮系统的模拟）有某种货物的存贮系统，市场对这种货物的需求量和订货提前期都是随机的，它们的概率分布

10、如下：表4-4 表4-5需求量概率累积概率00.020.0210.080.1020.220.3230.340.6640.180.8450.090.9360.071.00提前期概率累积概率10.230.2320.450.6830.170.8540.090.9450.061.00现在考虑订货、存贮、缺货损失三项费用：订货费用每次25元，订货量每次20单位，订货点为15单位。（即存货低于15单位时订货，但已订货未到前不再订）存贮费每件每周10元，缺货损失费每件每周500元。对于缺货，货到后不补，设开始时存货为20单位。试利用所给随机数R1（在下表内）模拟需求量，R2（50，86，15）模拟订货提前期

11、。模拟14周的运行情况：并求订货费用、存贮费用、缺货费用以及周平均费用。表4-6周需求到货量存储量订货缺货量随机数R1需求量是否订货提前期R2提前期020168416252313502390584593265812467242074431789561186498147109461112820112895151363312151140131可求得：订货费用25×3=75，存贮费用10×200=2000，缺货费用5001=500× 1周平均费用 例3.（排队系统模拟）有一银行营业点打算添置一台自动存取款机（12小时服务），顾客按一定的间隔时间到来，排队接受服务，先来者

12、先用，后来者后用，顾客不愿在队列中等待太久，否则会离去。管理人员想了解等待时间超过3分钟的顾客的比例为多少，若该比例太大，则考虑再增设一台机器。 表4-7 顾客到达间隔时间统计表间隔（分）人数频率累积频率随机数对应范围1180.180.180.000.172170.170.350.180.343150.150.500.350.494120.120.620.500.615100.100.720.620.71690.090.810.720.80780.080.890.810.88850.050.940.890.93920.020.960.940.951010.010.970.961110.010.

13、980.971210.010.990.981310.011.000.99总计1001.00表4-8 顾客用机时间统计表用机时间人数频率累积频率随机数对应范围1480.480.480.000.472200.200.680.480.673160.160.840.680.834120.120.960.840.95520.020.980.960.97620.021.000.980.99总计1001.00产生01（0.00-0.99）间隔两组均匀分布的随机数。 一组用于模拟顾客到达间隔时间，另一组模拟顾客用机时间。由第一组产生的一个随机数代表当前到达存取款机的一位顾客，若此随机数的值为0.70，可以确定

14、所模拟的该顾客到达的时间与前一位顾客到达时的间隔时间为5分钟。由第二组产生的一个随机数代表正在使用存取款机的一位顾客，若此随机数的值为0.90，可以确定所模拟的该顾客使用存取款机的时间为4分钟。表4-9 手工模拟步骤与结果（假设模拟开始时间为0） 顾客编号到达间隔随机数间隔时间用机时间随机数用机时间到达时间开始使用时间离开时间等待时间10.5740.502446020.0310.8945610130.9590.311141415040.3830.8031717200用手工方法模拟10位顾客，模拟开始时间设为0. 手工模拟的过程前4位顾客如表所示。后6位顾客的模拟留给读者作为练习来完成。如果排队

15、等待时间大于或等于3分钟的人数为M，模拟的总人数为T， 当T足够大时，M/T即可作为结果。 例4.单服务员的排队模型：在某商店有一个售货员，顾客陆续到来，售货员逐个接待顾客，顾客排队等待，业务完成后顾客离开商店。设：1. 顾客到来间隔时间服从参数为0.1的指数分布；2. 顾客服务时间服从4,15上的均匀分布；3. 排队按先来后到原则，队长无限制。假定时间以min为单位，一个工作日为8小时，（1） 模拟一个工作日内完成服务的个数及顾客平均等待时间。（2） 模拟100个工作日，并求出每日内完成服务的个数及顾客平均等待时间。解：设：总等待时间；：第个顾客的到达时刻；：第个顾客开始服务时刻；：第个顾客

16、服务结束时刻；则我们首先来模拟整个流程：初始化：令i=1,ei-1=0,w=0产生间隔时间随机数xiExp(0.1)，ci=xi, bi=xi产生服务时间随机数yiU4,15，ei= bi+yi累计等待时间wi= w+bi-ci准备下一次服务i=i+1产生间隔时间随机数xiExp(0.1)，ci=ci-1+xi确定开始服务时间bi= max(ci,ei-1)bi>480i=i-1,t=w/i结束服务，完成个数m=i，平均等待时间为w结束NoYes首先来模拟整个流程：用Matlab编程如下，见排队.m，经模拟计算得到个工作日内完成服务的个数m=45人，顾客平均等待时间t=33min；100

17、个工作日每日内完成服务的个数m=43人，顾客平均等待时间t=27min。Matlab程序为：Simu1.mcleari=2;w=0;e(i-1)=0;x(i)=exprnd(10);c(i)=x(i);b(i)=x(i);while b(i)<=480y(i)=unifrnd(4,15);e(i)=b(i)+y(i);w=w+b(i)-c(i);i=i+1;x(i)=exprnd(10);c(i)=c(i-1)+x(i);b(i)=max(c(i),e(i-1);endi=i-2;t=w/im=iSimu2.mclearcs=100;for j=1:cs j; w(j)=0; i=2;x

18、(i)=exprnd(10);c(i)=x(i);b(i)=x(i);while b(i)<=480y(i)=unifrnd(4,15);e(i)=b(i)+y(i);w(j)=w(j)+b(i)-c(i);i=i+1;x(i)=exprnd(10);c(i)=c(i-1)+x(i);b(i)=max(c(i),e(i-1);endi=i-2;t(j)=w(j)/i;m(j)=i;endpt=0;pm=0;for j=1:cs pt=pt+t(j); pm=pm+m(j);endpt=pt/cspm=pm/cs4.3 蒙特卡洛模拟求积分我们知道在近似计算中Simpson方法是一种常用的数

19、值积分方法，但对于不规则的被积函数及多维情形这种方法往往无能为力。而蒙特卡洛方法计算这种积分就较为方便了，虽然其精度较差，但在某些场合往往有其出现的必要性，当然过后我们会对这一点进行讨论。使用蒙特卡洛模拟计算积分有两种方法，以下我们分别予以介绍。1.随机投点法随机投点法有一个更为有趣的称呼Hit or Miss Monte Carlo积分法。考虑计算积分，其中为计算上述积分，我们考虑取矩形区域，设为在上均匀分布的随机变量，它有概率密度则落入S的概率满足即S（HitMiss那么我们只要能估计出即可给出的近似值。为估计，产生个相互独立的在上均匀分布的随机数 。用表示满足的数目，称为“Hit”的数目

20、，根据大数定律，可以用估计，从而得到的近似值性质：1）2）3）对于任意给定的，为使只须4）对于充分大的，积分有以下近似的的置信区间其中，为标准正态分布的分布函数。（证明略）例5：试用蒙特卡洛模拟中的随机投点法计算积分Mat lab程序如下：p=0; nt=0; %the total number of the test nofin=0; %the number of point in the area figure; for x=0:.01:5; for tmp=0:.1:50 y=500*rand(); y0=4*x3; nt=nt+1; if(y<=y0) nofin=nofin+1

21、; plot(x,y,'r'); hold on end end end sofsq=x*y0; %计算矩形的面积p=nofin/nt; s=p*sofsq2. 样本均值积分法计算积分的另一种方法是将其表示成某个随机变量的数学期望，我们只要将积分改写为这里假定为任意概率密度函数且满足，这样就可以将积分表示为如下的数学期望一般情况下，我们可以取X为上均匀分布的随机变量，即于是有为了估计积分，只要产生N个上均匀分布的随机数，并用样本均值来估计积分即可。性质：1）2）3），即平均值法得到的积分估计比随机投点法得到的更有效。例6. 试用蒙特卡洛模拟中的平均值法计算积分解：对于积分不妨改

22、写为则积分为服从上的均匀分布的随机变量的函数的数学期望，即接下来我们产生N个上的均匀分布的随机数，并求其函数的平均值，并令即可求得积分的近似值。Mat lab程序为：sum=0;for tmp=1:1:50000; x=5*rand(); y=20*x3; sum=sum+y;endms=sum/500004.4 蒙特卡罗方法的收敛性，误差 1. 收敛性我们来讨论平均值法的收敛性，由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，XN的算术平均值： 作为所求解的近似值。由大数定律可知，如X1，X2，XN独立同分布，且具有有限期望值（E(X)<），则 即随机变量X的简单子样的

23、算术平均值，当子样数充分大时，以概率1收敛于它的期望值E(X)。2 误差蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2 ，即f (X)是X的分布密度函数。则当N充分大时，有如下的近似式其中称为置信度，1称为置信水平。这表明，不等式近似地以概率1成立，且误差收敛速度的阶为。通常，蒙特卡罗方法的误差定义为上式中与置信度是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出。下面给出几个常用的与的数值：0.50.050.0030.67451.963关于蒙特卡罗方法的误差需说明两

24、点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出 。 3 减小方差的各种技巧 显然，当给定置信度后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。在固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。 4. 效率一般来说，降低方差的技巧，往往会使观察一个子样的时

25、间增加。在固定时间内使观察的样本数减少。所以，一种方法的优劣需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为，其中c是观察一个子样的平均费用。显然越小，方法越有效。 5 方法特点优点1) 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。2) 受几何条件限制小。在计算s维空间中的任一区域Ds上的积分时，无论区域Ds的形状多么特殊，只要能给出描述

26、Ds的几何特征的条件，就可以从Ds中均匀产生N个点，得到积分的近似值。其中Ds为区域Ds的体积。这是数值方法难以作到的。另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。 3) 收敛速度与问题的维数无关。由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比

27、如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。4) 具有同时计算多个方案与多个未知量的能力。对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不

28、同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。 5) 误差容易确定。对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。 6) 程序结构简单，易于实现。 在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。缺点1) 收敛速度慢。如前所述，蒙特卡罗方法的收敛速度为，一般不容易得到精确度较高的近似结果。对于维数少（三

29、维以下）的问题，不如其他方法好。 2) 误差具有概率性。由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。 6蒙特卡罗方法的主要应用范围 蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。习题41. 假设某个银行只有一个营业窗口开放，顾客随机性的进入银行，然后按照先后次序排队等候服务。营业员服务顾客的时间也是一个随机变量，假定顾客到达的时间和每位顾客接受服务的时间都已给出，并且已填写在下面的表格中，采

30、用“手工”的方式模拟这个系统。直到第10位顾客服务完成后离开系统。表4-10顾客序号到达时间服务时间服务开始时间服务结束时间等待时间空闲时间滞留时间12326439641035132617572248275931410343根据仿真结果回答下列问题： （1） 顾客在系统中的平均滞留时间（从进入到离开的时间）（2） 顾客的平均排队等待时间（3） 营业员处于空闲状态的百分比2. 有一单服务台的排队系统，根据经验资料知道到达的间隔时间和服务时间的概率分布如下表，其他条件符合标准情形。表4-11到达间隔概率累积概率对应随机数(a)20.40.40.0-0.460.30.70.4-0.7100.20.9

31、0.7-0.9140.11.00.9-1.0表4-12服务间隔概率累积概率对应随机数(a)10.40.40.0-0.430.40.80.4-0.850.21.00.8-1.0(1)今由随机数表任选两组随机数 RNa：902,321,211,021,198,383,107,799,439 RNb：612,484,048,605,583,773,054,853,313,200 试根据这两组随机数分别产生表示开始十位顾客的到达间隔时间的随机数AT和表示服务时间的随机数ST。(2)模拟这个排队系统的运行情况，这一阶段共运行多少分钟？(3)求系统空闲的概率P0。表4-13顾客RNaAT进入系统时刻RNb

32、ST服务时间系统空闲时间开始结束1061230311290214144843141733212160481171842112186053182150212205833212461982227733242773832240541272881072268535283339799103631313637510439642200142432. 一个大型超市每日都从农村采购新鲜农产品出售，正常情况下每公斤可获例1元。如果采购数量过多，次日只能减价出售，每公斤将亏损0.4元，现在该市采用以下采购策略：以前一天的市场需求量作为当天的采购量。据统计分析，每天平均需求量为100 kg,标准差为30 kg。在这种情况下，该超市经营一个月能获多少利润？ 3. 某企业计划投资建立某产品生产线，为此必须对该产品未来能实现多少利润进行预测和分析。建立该生产线需投资5万元。产品能实现多少利润主要受以下三个不确定因素的影响：售价、成本与年销售量。经过有关生产、计划、销售人员分析，考虑到原材料供应、市场竞争和价格浮动等因素的作用，初步估计售价、成本与年销售量可能出现的情况及其发生概率如下表。试用蒙特卡罗模拟法分析建立此产品生产线的未来盈亏状况。售价（元）发生概率成本（元）发生概率年销售量（万件）发生概率50.320.13.50.260.5