Monte Carlo方法及在计算积分和古典概率方面的应用

王坤

【摘 要】介绍了Monte Carlo方法的思想,主要在定积分计算方面介绍了随机投点法和平均值法,并将其推广到二重积分、三重积分和多重积分情形,最后以棋手分奖金问题介绍了Monte Carlo方法在古典概率中的应用.分析了误差,介绍了减少误差的方法.给出这些方法的实例及其Mathematica实现程序. 【期刊名称】《曲靖师范学院学报》 【年(卷),期】2010(029)003 【总页数】6页(P33-38)

【关键词】Monte Carlo方法;积分计算;古典概率;模拟 【作 者】王坤

【作者单位】曲靖师范学院,数学与信息科学学院,云南,曲靖,655011 【正文语种】中 文 【中图分类】O242.2

Monte Carlo方法,源于第二次世界大战美国关于研制原子弹的“曼哈顿计划”.该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城——摩纳哥的Monte Carlo来命名这种方法,为它蒙上了一层神秘色彩.Monte Carlo方法的基本思想很早以前就被人们所发现和利用.早在17世纪,人们就知道用事件发生的“频率”来决定事件的“概率”.19世纪人们用投针试验的方法来确定圆周率.20世纪40年代电子计算机的出现,特别是近年来高速电子计算机的出现,使得数学方法在计算机上

大量、快速地模拟这样的试验成为可能.

Monte Carlo方法研究的问题大致可分为两种类型:一种是问题本身就是随机的,另一种本身属于确定性问题,但可以建立它的解与特定随机变量或随机过程的数字特征或分布函数之间的联系,因而也可用随机模拟方法解决.

文[1～7]介绍了Monte Carlo方法的思想,但没有给出具体的实例及实现过程.本文介绍了Monte Carlo方法的思想,从计算定积分和古典概率两方面的应用进行研究,给出了实例及其Mathematica实现程序. 1.1 Monte Carlo方法思想概述

Monte Carlo方法,有时也称随机模拟(Random Simulation)方法或统计试验(Statistical Testing)方法.它的基本思想是:首先建立一个概率模型或随机过程,使它的参数等于问题的解;然后通过对模型或过程的观察、抽样来计算所求参数的统计特征;最后给出所求解的近似值,而解的精度可用估计值的标准误差来表示. 假设所求的量 x是随机变量ζ的数学期望E(ζ),那么近似确定x的方法是对ζ进行重复抽样,产生相互独立的ζ值的序列ζ1,ζ2,…,ζN并计算其算术平均值: 根据大数定理,当充分大时以概率1成立,即可用作为x的估计值.

Monte Carlo方法以概率统计理论为基础,以随机抽样(随机变量的抽样)为手段,在很多方面有重要的应用.它的优点表现在3个方面:方法和程序的结构简单,易分析、易理解;收敛的概率性和收敛速度与问题的维数无关,很好的避免了维数问题;受问题条件限制的影响较小,很好的提高可行性. 使用Monte Carlo方法的步骤如下:

(1)构造或描述概率过程.对于本身就具有随机性质的问题,主要是正确地描述和模拟这个概率过程.对于本来不是随机性质的确定性问题,就必须事先构造一个人为的概率过程,使它的某些参量正好是所要求解问题的解(将不具有随机性质的问题,转化为随机性质的问题).

(2)实现从已知概率分布中抽样.由于各种概率模型都可以看作是有各种各样的概率分布构成的,因此产生己知概率分布的随机变量,就成为实现Monte Carlo方法模拟试验的基本手段.

(3)建立各种估计量.一般说来,构造了概率模型并能从中抽样后,即能实现模拟试验后,就要确定一个随机变量,作为所要求的问题的解的估计量.在Monte Carlo计算中,使用最多的是无偏估计.建立各种估计量,相当于对模拟试验的结果进行考察和登记,从中得到问题的解.

1.2 Monte Carlo方法的可行性

从Monte Carlo方法的基本思想可以得到它通常的做法,利用数学或物理方法产生[0,1]中均匀分布的随机数,再变换得到任意分布的随机数.随机数个数很大时,可以由大数定理求出事件的概率值.这种做法的可行性主要依据下面的事实:

(1)如果随机变量ε的分布函数是F(x),由于 F(x)非降.对于任意的y(0y}作为 F(x)的反函数.我们考虑随机变量η=F(ζ)的分布,这里假定 F(x)是连续函数,则对于0≤x≤1有: 即η=F(ζ)服从[0,1]上的均匀分布.

(2)反之,如果服从[0,1]上的均匀分布,则对于任意的分布函数 F(x),令ζ=F-1(η),则: 因此ε是服从分布函数F(x)的随机变量.

所以我们只要能够产生[0,1]中均匀分布的随机变量的子样,那么通过(2)式我们就可以得到任意分布函数 F(x)的随机变量的子样.再结合大数定理,就可以运用Monte Carlo方法进行随机模拟,解决一些实际的问题. 2.1 随机投点法

对于定积分 为使计算机模拟简单起见,设 a,b有限,0≤f(x)≤M,令Ω =

{(x,y)∶a≤x≤b,0≤y≤M},并设(X,Y)是在Ω上均匀分布的二维随机变量,其联合密度函数为是Ω中曲线y=f(x)下方的面积.

假设我们向Ω中进行随机投点.若点落在y =f(x)下方(即 y该方法的具体计算步骤为:独立地产生2n个U(0,1)随机数ui,vi,i=1,…,n;计算xi=a +ui(b-a),yi=Mvi和f(xi);统计f(xi)≥yi的个数 n0;用(3)估计θ.

例1 1777年,法国学者Buffon提出用试验方法求圆周率π的值.原理如下: 假设平面上有无数条距离为1的等距平行线,现向该平面随机地投掷一根长度为l(l≤1)的针.则我们可以计算该针与任一平行线相交的概率.此处随机投针可以这样理解:针中心与最近的平行线间的距离 x均匀地分布在区间[0, 1/2]上,针与平行线间的夹角φ(不管相交与否)均匀地分布在区间[0,π]上(如图1).于是,针与线相交的充要条件是从而针线相交概率为:

而由大数定律可以估计出针线相交的概率其中n为掷针次数,m为针线相交次数,从而圆周率其Mathematica实现语句见附录1. 历史上曾有几位学者做过这样的试验,结果如表1所示.

Buffon投针试验由于当时条件的落后,随机试验中误差比较大,精确度不高(Lazzarini的试验精度很高纯属巧合). 2.2 样本平均值法

对于积分 设 g(x)是(a,b)上的一个密度函数,改写

由矩法,若有 n个来自g(x)的观测值,则可给出θ的一个矩估计,这便是样本平均值法的基本原理.

若a,b有限,可取设x,…,1

该方法的具体计算步骤为:独立地产生 n个U(0,1)随机数 u1,…,un;计算 xi=a+ui(ba)和f(xi),i=1,…,n;用(5)估计θ.

后面将给出一个例子说明此方法的应用. xn是来自U(a,b)的随机数,则θ的一个估

计为 方法1

其中Ω0为 S维单位立方体,{0≤xi≤1,i =1,2,…,S},在Ω0上有:0≤g(x1,x2,…,xi)≤1.很明显.此时积分(5)可以看作为求S+1维空间长方体 V:Ω0×[0,g(x1,x2,…,xs)]的体积.即:

对于这种较为一般形式的多重积分计算问题,采用的还是随机投点. 具体步骤如下:

首先产生S+1个随机数ζi(i=1,2,…,S)及η,构造S+1维随机向量η),然后检验ε是否落后在V中,同理可以推论.检验η≤g(ζ1,ζ2,…,ζs)是否成立,如果在构成的N个随机向量ζ中,有m个随机向量落于V中,那么取作为积分的近似值,即如果积分区域及被积函数不满足上述条件,那么可以通过变换达到所希望的条件. 其中积分区域Ω包含在K维多面体中,此多面体决定于 K个不等式ai≤xi≤bi(i=1,2,…,K).

设函数f(x1,x2,…,xk)在Ω内连续且满足条件:0≤f(x1,x2,…,xk)≤M,N是在K+1维多面体Ωx[0,M]中均匀分布的随机质点的个数,n是在N个随机点之中落入以K维区域V为底以 Z=f(x1,x2,…,xk)为顶之曲顶柱体内的随机点的个数.

这里Ωx[0,M]表示由不等式ai≤xi≤bi(i =1,2,…,K)和0≤Z≤M决定的K+1维多面体.则 K重积分的Monte Carlo近似计算公式为:

例2 在三维空间中,由三个圆柱面:x2+ y2=1,x2+z2=1,y2+z2=1围成一个立体,利用Monte Carlo方法求它的体积.分析:据题意,所求体积中

考虑在空间Ω内随机的产生n个点,落在空间内有m个,则在Mathematica中模拟程序见附录2.

在生产、生活中概率是用于描述不确定性事件发生的可能性的大小,对于一些比较简单的随机事件的概率计算问题,我们可以通过比较常用的概率计算公式进行准确

计算.但是随着研究问题的深入和事件本身的复杂性,直接导致了概率计算的困难,甚至无法计算.但由于计算机技术的发展,近代发展起来的Monte Carlo方法在复杂的事件的概率计算中起了十分重要的作用.

下面的例子说明了Monte Carlo方法在古典概率中的应用.

例3 甲乙两位棋手棋艺相当,现他们在一项奖金为1000元的比赛中相遇,比赛为5局3胜制,已经进行了3局的比赛,结果为甲2胜1负,现因故要停止比赛,问应该如何分配这1000元比赛奖金才算公平?

分析:平均分对甲欠公平,全归甲则对乙欠公平.合理的分法是按一定的比例分配. 现在我们用计算机模拟两位棋手后面的比赛,就可以知道奖金分配方案.由于两位棋手的棋艺相当,可以假定他们在以下每局的比赛胜负的机会各半.Mathematica中函数产生随机数0或1,0与1出现的机会各占一半,可以用随机数1表示甲棋手胜,而随机数0表示乙胜.(也可以用[0,1]中的随机实数来模拟两人的胜负,随机数大于0.5表示甲胜,否则乙胜)连续模拟1000次(或更多次数)每次模拟到甲乙两方有一方胜了三局为止.按所说方案分配奖金,1000次模拟结束后,计算两棋手每次的平均奖金,就是该棋手应得的奖金. 模拟结果:

甲:750,乙:250(程序见附录1)

最终以甲分到;乙分到即甲750元,乙250元.

实际上,因为比赛只需进行两局.则可分出胜负.结果无非是以下4种情况之一: 甲甲、甲乙、乙甲、乙乙.

上面4种情况可看出,甲获胜的概率为乙获胜的概率为在Mathematica中模拟程序见附录3. 5.1 收敛性

蒙特卡罗方法是由随机变量 X的简单子样

x1,x2,…,xn的算术平均值 作为所求解的近似值.由大数定律可知,如 x1,x2,…,xn独立同分布,且具有有限期望值(E(X)< ∞),则即随机变量X的简单子样的算术平均值,当子样数N充分大时,以概率1收敛于它的期望值 E(X). 5.2 误差

蒙特卡罗方法的近似值与真值的误差问题,概率论的中心极限定理给出了答案.该定理指出,如果随机变量序列 X1,X2,…,XN独立同分布,且具有有限非零的方差σ2,即 f(X)是 X的分布密度函数.则 当N充分大时,有如下的近似式

其中α称为置信度,1-α称为置信水平.这表明,不等式近似地以概率1-α成立,且误差收敛速度的阶为 O(N-1/2).

通常,Monte Carlo方法的误差ε定义为

上式中λα与置信度α是一一对应的,根据问题的要求确定出置信水平后,查标准正态分布表,就可以确定出λα.

下面给出几个常用的α与λα的数值:

关于蒙特卡罗方法的误差需说明两点:第一,蒙特卡罗方法的误差为概率误差,这与其他数值计算方法是有区别的.第二,误差中的均方差σ是未知的,必须使用其估计值来代替,在计算所求量的同时,可计算出

例4 用平均值法估计圆周率π,并考虑置信度为5%,精度要求为0.01的情况下所需的试验次数.

解:易知故考虑令 X～ U(0,1),令其期望值为 因此其中xi是[0,1]区间上均匀分布的随机数. 此时,α=0.05,Zα/2=1.96,ε=0.01/4,所以 5.3 减小方差的各种技巧

显然,当给定置信度α后,误差ε由σ和N决定.要减小ε,或者是增大 N,或者是减小

方差σ2.在σ固定的情况下,要把精度提高一个数量级,试验次数N需增加两个数量级.因此,单纯增大N不是一个有效的办法.

另一方面,如能减小估计的均方差σ,比如降低一半,那误差就减小一半,这相当于 N增大4倍的效果.因此降低方差的各种技巧,引起了人们的普遍注意.

一般来说,降低方差的技巧,往往会使观察一个子样的时间增加.在固定时间内,使观察的样本数减少.所以,一种方法的优劣,需要由方差和观察一个子样的费用(使用计算机的时间)两者来衡量.这就是蒙特卡罗方法中效率的概念,它定义为σ2c,其中c是观察一个子样的平均费用.显然σ2c越小,方法越有效.

总的来说,增大样本n的值对计算机要求较高;减小方差的技巧都只具有指导思想上的意义.对于实际的计算问题,往往要求对涉及的随机变量有先验的了解,或者对发生的物理过程的性态有一定的认识.通过利用这些预知的信息采取相应的手段减小误差,提高精度. 附录

【相关文献】

[1]徐钟济.蒙特卡罗方法[M].上海科学技术出版社, 1985:171～188.

[2]茆诗松,王静龙,濮晓龙.高等数理统计[M].北京:高等教育出版社,2006:415～454. [3]周铁,徐树方,张平文,等.计算方法[M].吉林:清华大学出版社,2006:299～353. [4]李尚志,陈发来,张韵华,等.数学实验[M].北京:高等教育出版社,2004:23～30.

[5]王岩.Monte Carlo方法应用研究[J].云南大学学报(自然科学版),2006,28(S1):23～26. [6]薛毅,陈立萍.统计建模与R软件[M].北京:清华大学出版社,2008:476～485. [7]杨自强.你也需要蒙特卡罗方法——提高应用水平的若干技巧[J].数理统计与管理,2007,27(2):355～376.