蒙特卡罗（Monte carlo）方法

蒙特卡罗（Monte carlo）方法 【实验目的】 z 了解 Monte carlo 随机摸拟方法的具体应用 z 学习使用 MATLAB 软件中有关随机数生成函数 【实验准备与内容】 1．Monte carlo 方法 蒙特卡罗(Monte Carlo)方法，又称计算机随机模拟方法，是一种基于"随机数"的计算方法。这 一方法源于美国在第二次世界大战研制原子弹的"曼哈顿计划"。该计划的主持人之一数学家冯诺伊 曼用驰名世界的赌城-摩纳哥的 Monte Carlo 来命名这种方法。 Monte Carlo 方法的基本思想很早以前就被人们所发现和利用。早在 17 世纪，人们就知道用事 件发生的"频率"来决定事件的"概率"。1777 年，蒲丰（Buffon）提出著名的 Buffon 投针试验永来近 Yi Xuming 矩形 Yi Xuming 铅笔 Yi Xuming 矩形 似计算圆周率 π。本世纪 40 年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得人 们在计算机上利用数学方法大量、快速地模拟这样的试验成为可能。目前这一方法已经广泛地运用 到数学、物理、管理、生物遗传、社会科学等领域，并显示出特殊的优越性。 2．伪随机数 实际应用中的随机数通常都是某些数学公式计算而产生的伪随机数，这样的伪随机数从数学意 义上讲不是严格的随机数，但是，只要伪随机数能够通过随机数的一系列统计检验，我们就可以将 它当作真随机数而放心使用。这样我们就可以方便、经济、重复地产生随机数。理论上要求伪随机 数产生器具备以下特征：良好的统计分布特性，高效率的伪随机数产生，伪随机数产生的循环周期 长，伪随机数可以重复产生等。到目前为止，已经提出了各种分布的伪随机数产生方法，在这里， 我们不准备从数学上介绍随机数的产生原理，而只给出伪随机数的 MATLAB 产生函数。 z 生成各类分布的随机数统一生成函数 random('name',A1,A2,A3,m,n)——生成以 A1,A2,A3 为参数分布为 name 的 nm× 阶随机数矩阵， 其中，‘name’为包含特定分布名称的字符串(如表*.*)，A1,A2,A3 是分布参数矩阵，根据分布不同， Yi Xuming 矩形 Yi Xuming 椭圆形 各参数的含义也不相同，且其中一些参数也不是必须的。 z 针对不同分布的随机数产生函数 MATLAB还提供了不同分布的随机数产生函数，这些函数的具体使用格式及参数含义较 random 函数更明确（如表*.*）。 表*.* 特定分布 字符串 ‘name’ 分布名 称 随机数 产生函数 函数格式 'beta‘ β分布 Betarnd Betarnd(A,B,m,n) 'bino‘ 二项分 布 Binornd Binornd(N,P,m,n) 'chi2' 2χ 分布 chi2rnd chi2rnd(V,m,n) 'exp‘ 指数分 Exprnd Exprnd(MU,m,n) 布 'f' F 分布 frnd Frnd(V1,V2,m,n) 'gam' γ 分布 Gamrnd Gamrnd(A,B,m,n) 'geo' 几何分 布 Geornd Geornd(P,m,n) 'hyge' 超几何 分布 Hygernd Hygernd(M,K,N,mm,nn) 'logn' 对数正 态分布 Lognrnd Lognrnd(Mu,Sigma,m,n) 'nbin' 负二项 分布 Nbinrnd Nbinrnd((R,P,m,n) 'ncf‘ 非中心F 分布 Ncfrnd Ncfrnd(Nu1,Nu2,Delta,m,n) 'nct‘ 非中心 t 分布 Nctrnd Nctrnd(V,Delta,m,n) 'ncx2' 非中心 分布 2χ ncx2rnd ncx2rnd(V,Delta,m,n) 'norm' 正态分 布 Normrnd Normrnd(Mu,Sigma,m,n) 'poiss' 泊松分 布 Poissrnd Poissrnd(Lambda,m,n) 'rayl' Rayleigh 分布 raylrnd Raylrnd(B,m,n) 't' t 分布 Trnd Trnd(V,m,n) 'unif' 连续均 匀分布 unifrnd Unifrnd(A,B,m,n) 'unid' 离散均 匀分布 Unidrnd Unidrnd(N,m,n) 'weib‘ Weibull 分布 weibrnd Weibrnd(A,B,m,n) 为了说明函数的使用方法，我们举两个例子。 例：用 random 函数生成均值为 0， 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 差为 2 的正态分布随机数。 rn=random('norm',0,2,10000,1); ％生成正态分布随机数 hist(rn，30) ％绘出随机数的直方图 -10 -5 0 5 10 0 200 400 600 800 1000 1200 例：生成区间[1,5]的服从均匀分布的随机数 rn=unifrnd(1,5,10000,1); ％生成均匀分布随机数 hist(rn，30) ％绘出随机数的直方图 1 2 3 4 5 0 50 100 150 200 250 300 350 400 3. Buffon 投针试验 在平滑的桌面上画一组相距 2a 的平行线束，向此桌面上投一枚长为 2b 的细针，为避免针与两 平行线同时相交的复杂情况，假定 。现在，我们来求针与平行线相交的概率。 0>> ba 设 M 为针的中点，y 为 M 与最近平行线的距离，x 为与针与平行线的交角（ π≤≤ x0 ），于是，针 与平行线相交的充要条件是 xby sin0 ≤≤ ,故相交的概率为图中曲线 xby sin= 与 x轴所夹图形的面积占图 Yi Xuming 矩形 中长方形面积的百分比。 xby sin= o a π x y 即 a bxdxb a p ππ π 2sin1 0 == ∫ 从而 pa b2=π 用 N 表示投针的次数，n 表示其中针与平行线相交的次数，由贝努里（Bernoulli）定理知，当 N 充 分大时，频率接近于概率，即 N np ≈ ，从而 na bN2≈π 。这就是 Buffon 用随机实验近似求圆周率π的基 本公式。 我们不妨用产生随机数的办法来模拟投针实验：对每一次投针实验，针与平行线的交角 x （ π≤≤ x0 ）和针的中点 M 与最近平行线的距离 y（ ay ≤≤0 ）是一随机变量，其分别服从 ],0[ π 和 上的均匀分布。因此，N 次投针实验可以通过以上产生随机数的函数生成，然后根据针与直线相 交的充要条件 ，即可判断针与直线是否相交，从而可以计算出 n 的值，达到近似计算圆 周率 ],0[ a xby sin0 ≤≤ π的目的。具体程序如下 function [pai,number]=buffon(a,b,N) % 2a，2b 分别为平行线间距和针长，N 为投针次数 x=unifrnd(0,pi,N,1); y=unifrnd(0,a,N,1); number=0; % 相交计数器 for i=1:N if y(i)<=b*sin(x(i)) number=number+1; end end pai=2*b*N/(a*number) 其计算结果如下： a b N 相交次 数 π的估计 值 45 36 10000 5139 3.1134 45 36 100000 50992 3.1377 45 36 500000 254722 3.1407 4．积分计算 定积分的计算是 Monte Carlo 方法引入计算数学的开端，在实际问题中，许多需要计算多重积 分的复杂问题，用 Monte Carlo 方法一般都能够很有效地予以解决，尽管 Monte Carlo 方法计算结果 的精度不很高，但它能很快地提供出一个低精度的模拟结果也是很有价值的。而且，在多重积分中， 由于 Monte Carlo 方法的计算误差与积分重数无关，因此它比常用的均匀网格求积公式要优越。 4．1 随机投点法 （1） 设计算的定积分为 ，其中 为有限数，被积函数 是连续随机变量dxxfI b a ∫= )( ba, )(xf ς 的概率密度 函数，因此 满足如下条件： )(xf )(xf 非负，且 （1） 1)( =∫+∞ ∞− dxxf 显然 I 是一个概率积分，其积分值等于概率 )( baP <≤ ς 。 下面按给定分布 随机投点的办法，给出如下 Monte Carlo 近似求积算法： )(xf step1: 产生服从给定分布的随机变量值 Nixi ,...,2,1, = ； Yi Xuming 矩形 Yi Xuming 椭圆形 step2: 检查 是否落入积分区间。如果条件ix bxa i <≤ 满足，则 记录 混凝土 养护记录下载土方回填监理旁站记录免费下载集备记录下载集备记录下载集备记录下载 落入积分区间一次。 ix 假设在 N 次实验以后， 落入积分区间的总次数为 n，那么用 ix N nI =_ 作为概率积分的近似值，即 N nI ≈ （2） 如果要计算的定积分为 ，其中 为有限数，但被积函数 不是某随机变量的的 概率密度函数。在这种情况下，当然我们可以对积分作变换，将积分变换成满足（1）的条件。但 在变换以后，需要产生新的分布随机变量，因此常会遇到很多困难和比较复杂的计算。 dxxfI b a ∫= )( ba, )(xf 上述求积方法需要产生给定分布的随机变量，它适合于解决特殊类型的概率积分。如果只用随 机数完成随机投点，那么，下面的方法可以解决较为广泛的一类积分问题。 （3） 设 是[0,1]上的连续函数，且)(xf 1)(0 ≤≤ xg ，需计算积分 ，如下图所示阴影部分的面积。 在图中单位正方形内均匀地投点 dxxfI b a ∫= )( ),( ης ，则该随机点落入曲线 )(xfy = 下面的概率为 IdxxfxfyP ==≤ ∫1 0 )()}({ 因此，给出如下 Monte Carlo 近似求积算法： step1: 产生两组[0，1]区间内的随机数 Niyx ii ,...,2,1,, = ，并把 作为随机 点 ),( ii yx ),( ης 的可取坐标； step2: 检查（ ， 是否落入ix )iy Ω内，如果条件 )( ii xfy ≤ 满足，则记录（ ， 落入积分区间一次。 ix )iy 假设在 N 次实验以中，落入Ω内总次数为 n，那么量 N nI =_ 近似等于随机点落入 内的概率，即 Ω N nI ≈ 假如所需计算积分 ，其中 为有限数，被积函数 有界，并用dxxfI b a ∫= )( ba, )(xf M 和 分别表示其最 大值和最小值。作变换 ， m *)( xabax −+= ])(([1)( *** mxabaf mM xf −−+−= ，此时 )()())(( 1 0 * abmdxxfabmMI −+−−= ∫ 且有 ，即转化为上面讨论过的情况。 ]1,0[,1)(0 ** ∈≤≤ xxf 4．2 平均值法 任取一组相互独立、同分布的随机变量 }{ iξ ， iξ 在[a,b]内服从分布率 ，令)(xp )( )()(* xp xfxp = ，则 也是一组相互独立、同分布的随机变量，而且 )}({ * ip ξ ∫ ∫ === b a b a i IdxxfdxxpxppE )()()()}({ ** ξ Yi Xuming 椭圆形 由强大数定理 1)(1lim 1 * =⎟⎠ ⎞⎜⎝ ⎛ =∑ =∞→ N i iN Ip N P ξ 若记 ∑ = −= N i i IpN 1 * )(1 ξ ，则 −I 依概率 1 收敛到 I ，平均值法就是用 −I 作为 I 的近似值。 假如所需计算积分为 ，其中被积函数在[a,b]内可积，任意选择一个有简单办法可以进行 抽样的概率密度函数 ，使其满足条件： dxxfI b a ∫= )( )(xp z ,0)( ≠xp 当 时（0)( ≠xf bxa ≤≤ ） z 1)( =∫ dxxp b a 记 ⎪⎩ ⎪⎨ ⎧ = ≠= . 0p(x) ,0 ,0)(, )( )( )(* xpxp xf xp 则所求积分为 dxxpxpI b a ∫= )()(* 因而 Monte Carlo 近似求积算法为： step1: 产生服从分布律 的随机数)(xp Nixi ,...,2,1, = ； step2: 计算均值 ∑ = − = N i ipN I 1 * )(1 ξ ，即有 II ≈− 。 如果 为有限数，那么 可以取为均匀分布： ba, )(xp ⎪⎩ ⎪⎨ ⎧ ≤≤−= 其它,0 ,1)( bxaabxp 此时 dxabxfabI b a ∫ −−= 1)()( ，具体的算法为： step1: 产生[a,b]上的均匀随机数 Nixi ,...,2,1, = ； step2: 计算均值 ∑ = − −= N i ixfN abI 1 )()( ，且有 −≅ II 。 4．3 多重积分的计算 根据上面定积分的 Monte carlo 计算方法，我们可以很容易地将它推广到 重积分的近似计 算。为了避免重复，在此只给出具体的二重和三重积分计算公式。 )2( ≥nn z 二重积分的计算 ∑∫∫ = ≅ N i ii D yxf N Ddxdyyxf 1 ),(||),( 其中， 是平面区域 D 中的均匀分布的 N 个独立选取的随机点列， 表示区域 D 的面积。 ),( ii yx || D z 三重积分的计算 ∑∫∫∫ =Ω Ω≅ N i iii zyxfN dzdxdyzyxf 1 ),,(||),,( 其中， 是空间区域),,( iii zyx Ω中的均匀分布的 N 个独立选取的随机点列， || Ω 表示区域的Ω体积。 很显然，该方法特别适合于多重积分，特别是积分区域比较复杂的重积分计算。 4．4 误差估计 Monte Carlo 求积方法常常以某个随机变量θ的简单子样 Nθθθ ,,, 21 L 的算术平均值 ∑ = −= N i iN 1 1 θθ 作为积分 I 的近似值，由强大数定理知道，如果随机变量序列 Nθθθ ,,, 21 L 相互独立、同分布、期望值存在，那么当 时， 以概率 1 收敛到∞→N −θ I 。按照中心极限定理，只要随机变量序列 Nθθθ ,,, 21 L 相互独立、同分布、 数学期望存在，且具有有限标准差 0≠σ ，那么，当当 时，随机变量 ∞→N N IY /σ θ−= − 渐进地服从标准正态分布 ，即 )1,0(N dxetYP t x ∫ ∞− −≅< α πα 2 2 2 1)( 因此，对任何 ， 0>αt απ σθ α α α −=≅⎟⎟⎠ ⎞⎜⎜⎝ ⎛ <−=< ∫ −− 12 2||)|(| 0 2 2 dxe N tIPtYP t x 这表示，不等式 N tI σθ α<−− || 成立的概率近似地等于 α−1 。其中当α 很小时，α 称为显著水平， α−1 称为 置信水平。因此 Monte Carlo 方法的误差ε可以写为 N t σε α= ε，N 的值应取为（在概率意义下的统计估计值） 由此可以确定，对于给定的 2 22 ε σαtN = 但是，在一般意义下，我们无法知道理论值σ ，所以在实际计算中，只能用标准差的近似估计值 来 代替 −σ σ 。具体地说，对于随机投点法， )1( −−− −= IIσ ；对于平均值法， ∑∑ == − −= N i i N i i xpN xp N 1 2* 1 2* ))(1()(1σ 。 5．库存管理问题 某种商品进货价格为 a 元，出售价格为 b 元，假设对该商品每天早晨进货配齐 n 个，每日顾客 相互独立地到来，平均每日 m 人，且服从普阿松（Poisson）分布 !)( k emkP mk − = ，其中 k 为每日到来 的人数。问当 a=2、b＝3、m=10 时，每天早晨该商品备齐多少个可得到最大利润。 假设购买此商品的顾客每人只购一个，而且备齐的商品如果当天卖不出去则不可在第二天出售， 若用 表示顾客为 人时的日利润，则 ),( kna k ⎩⎨ ⎧ ≥− <−= nknanb nknakb kna , , ),( 利用 MATLAB 给定的 Poisson 随机数产生函数 poissrnd，生成含有 N 个元素随机数序列 ，这样 平均利润为 N iik 1}{ = N kna N i i∑ ),( 。下面给出了 MATLAB 程序： Yi Xuming 矩形 function y=lirun(a,b,m,n,N) % N 为产生随机数的个数 totol=0; ％totol 为总利润 rand('state',sum(100*clock)); randnumber=poissrnd(m,N,1); ％产生 Poisson 随机数 for i=1:N if randnumber(i) 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 中选出一个较优的方案。 )4.22 2 方案 1：按前一天的销售量作为当天的货存量； 方案 2：按前两天的平均销售量作为当天的货存量。 试用 Monte Carlo 方法确定该选择以上哪一种方案？