蒙特卡罗随机模拟投点法.doc

1、1蒙 特 卡 罗 随 机 模 拟 投 点 法 在 数 字 积 分 中 的应 用数 学 与 应 用 数 学 0901 班 ： 张 瑞 宸指 导 老 师 ： 任 明 慧摘 要 ：本 文 首 先 介 绍 了 蒙 特 卡 罗 方 法 的 产 生 和 发 展 ，然 后 分 析 了 蒙 特 卡 罗 方法 计 算 数 值 积 分 的 理 论 原 理 ，最 后 给 出 了 蒙 特 卡 罗 方 法 计 算 数 值 积 分 的MATLAB 编 程 实 现 ，全 文 主 要 是 讨 论 了 蒙 特 卡 罗 方 法 在 定 积 分 计 算 的 应 用 。而 蒙 特 卡 罗 的 优 点 ：可以计算被积函数非常复杂的定积

2、分、重积分，并且维数没有限制，这是别的数值积分方法还未达到的。 蒙 特 卡 罗 的 缺 点 ：收 敛 速 度 慢 ，误差 一 般 较 大 ，且 是 概 率 的 误 差 ，不 是 真 正 的 误 差 。关 键 词 ：蒙 特 卡 罗 方 法 ，均 值 估 计 法 ，数 值 积 分 ，Matlab 编 程 Abstract：This paper first introduces the emergence and development of the Monte Carlo method, and then analyze the theoretical principles of Monte Ca

3、rlo numerical integration method, Full-text mainly discussed the application of the Monte Carlo method in the definite integral. The advantages of Monte Carlo: can be calculated the integrable functions very complex definite integral, Multiple integrals, and dimension no limit, other numerical integ

4、ration methods have not yet reached. Monte Carlo Disadvantages: slow convergence speed, error generally higher, and the probability of error, not a real error.Keywords: Monte Carlo method，Mean estimation method，numerical integral，Matlab programming 0 引言历史上有记载的蒙特卡罗试验始于十八世纪末期（约 1777 年） ，当时布丰（Buffon）为了

5、计算圆周率，设计了一个“投针试验” ，后文会给出。虽然方法已经存在了 200 多年，此方法命名为蒙特卡罗则是在二十世纪四十年，美国原子弹计划的一个子项目需要使用蒙特卡罗方法模拟中子对某种特殊材料的穿透作用。出于保密缘故，每个项目都要一个代号，传闻命名代号时，项目负责人之一 von Neumann 灵犀一点选择摩洛哥著名赌城 蒙特卡罗（Monte 蒙特卡罗随机模拟投点法在数值积分中的应用2Carlo） 作为该项目名称，自此这种方法也就被命名为 Monte Carlo 方法广为流传。蒙特卡罗方法，又名随机模拟法或统计实验法它是以概率统计理论为基础，依据大数定律（样本均值替代总体均值）利用电子计算机

6、数字模拟技术，解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。通常蒙特卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗方法是一种有效的求出数值解的方法。一般蒙特卡罗法在数学中最常见的应用就是蒙特卡罗随机投点和蒙特卡罗数值积分。1 蒙特卡罗方法的产生与发展蒙特卡罗方法是在二战期间产生和发展起来的他的奠基者是美籍匈牙利人数学家冯诺伊曼（J.Von Neuman

7、n 1903-1957）由于通常计算量相当大而电子计算机在当时还没有出现，所有运算只能用手工进行，故而相当长的时间里蒙特卡罗方法难以推广。1.1 蒙特卡罗方法的产生作为蒙特卡罗方法的最初应用，是解决蒲丰氏问题1777 年，法国数学家Buffon 提出利用投针实验求解的问题：设平面上有无数多条距离为 1 的等距平行线，现向该平面随机的投掷一根长度为 的针。随机投针是指针的中心点于最近的平行线间的距离 均匀l x分布在0,1/2上，针与平行线的夹角 （不管相交与否）均匀分布在 上，0,如图一所示，故而，故其 ， 概率密度函数分别为1(0,)2xU(0,)。1(),()2px故我们得到针与线相交的充

8、要条件是 2sinlx蒙特卡罗随机模拟投点法在数值积分中的应用3图1.1 针与线相交的几种情况则针与线相交的概率是 (sin)2lpx= si sin200(l lpdx= 2.sinlld所以得到圆周率 。2(si)llpxp假如我们能做大量的投针实验并记录下针与线的相交次数，则可以根据大数定律估计出针线相交的概率 P。投针实验 N 次可能有 n 次使针与任意平行线相交, 那么 ，显然，试验次数 N 越多，P 的近似程度越好。历史上曾有npN几位学者做过这样的投针试验，并用手工计算出 值，结果参见表1。表 1实验者 时间（年） 针长 l投针次数 相交次数 估计值Wolf 1850 0.80

9、5000 2532 3.15956Smith 1855 0.60 3024 1218 3.15665Fox 1884 0.75 1030 489 3.15951Lazzarini 1925 0.83 3408 1808 3.14159292应该 指出，上述试验的精度 一般不会很高。譬如，假 设21,pl则。则由中心极限定理， 如果试验次数为 N，则 p 的估计值 渐进服从 ，近似为 。因(1),pN0.231(.6,)N蒙特卡罗随机模拟投点法在数值积分中的应用4此，若要以95%的概率保证 p 的精确到三位有效数字，即 与 p 的差距小于0.001，则 N 必须满足 N 。2251.960.3/

10、.18.90重复进行上千次的投针实验和手工计算，要消耗大量的人力、财力，蒙特卡罗方法虽然能解决此类问题，但得不到推广应用。1.2 蒙特卡罗方法的发展20世纪40年代以后，随着电子计算机的出现和发展，人们有可能用计算机来模拟这类实验和计算。计算机具有计算速度高和存储容量大的特点，采用数字模拟技术可以代替许多实际上非常庞大而复杂的实验，并迅速将实验结果进行运算处理，于是 Monte Carlo 方法重新被提起，引起世人重视，应用日渐广泛。实际上，采用 Monte Carlo 方法在计算机上建立模型来解决 Buffon 问题是非常简单的。我们在计算机上进行模拟试验，给定 ，我们可以在计算机上随机产生

11、 和l x，然后判断 是否成立。若成立，则针线相交，否则不交。假如我们2sinlx在计算机上独立的产生 N 对这样的 和 ，并记录下 成立的次数,记为x2sinlxn。则 估计值可取为 ，这就是随机模拟计算的结果。2l借助 Matlab 在计算机上进行投针实验， （相应程序见附录1）所取得的计算结果参见表2。表 2N =0.4l=0.6l=0.8l=0.9l=1.0l50000 3.1252442 3.1583934 3.1392246 3.14642309 3.1464351100000 3.1241457 3.1211978 3.1269543 3.1324504 3.1286664200

12、000 3.1414435 3.1396353 3.1373472 3.1366262 3.1371811500000 3.1519392 3.1470787 3.1455443 3.1422826 3.1425492800000 3.1483358 3.1465609 3.1440825 3.1423004 3.14302941000000 3.1435296 3.1459072 3.1412646 3.1419370 3.1423368蒙特卡罗随机模拟投点法在数值积分中的应用52000000 3.1497197 3.1462907 3.1418722 3.1410570 3.1407306

13、5000000 3.1393330 3.1398735 3.1410784 3.1416836 3.1412106其中 = 3.1416。不难看出，利用计算机运算不仅结果精确，而且迅速。随着电子计算机的普及，蒙特卡罗方法作为一种独到的方法得到开发，并首先应用到核武器的试验与研制中。尤其是各种可是编程方法的不断涌现，更显示出蒙特卡罗方法的最独到的优点，即形象直观地用数学方法在电子计算机上实现数字模拟实验。2 蒙特卡罗方法计算数值积分的理论原理随机模拟是一种随机试验的方法，也称为蒙特卡罗方法（Monte Carlo）。这种方法利用随机试验，根据频率与概率、平均值与期望值等之间的关系，1推断出预期的

14、结果。2.1 蒙特卡罗算法的理论原理随机模拟法的基本原理非常简单，下面给以直观的说明。首先构造一个定积分 ，其中积分 （即 1/41024dx102dxS单位圆的面积）用随机模拟方法求得。图 2.1 中粗线是 1/4 单位圆，如果向图2.1 中边长为 1 的正方形里随机投 块小石头，n当 很大时小石头会大致均匀地分布在正方形n中，数一下落在 1/4 单位圆内的小石头，假定有 个，那么 / 就可以看作 1/4 单位圆面积kkn/4 的近似值。小石头的位置坐标可以用产生均匀分布随机数的程序得到，记为是0,1区间均匀分布随机数（ 相互独立） ，记录满足),.21(,niyxi iyx,的数量 ，即得

15、 =4 / 。2iikkn我们用概率论中的大数定律来说明这个直观认识的原理。图2.1 xy0 11蒙特卡罗随机模拟投点法在数值积分中的应用6大数定律（伯努利（Bernoull）定理） 设 是 次独立重复试验中事2kn件 A 发生的次数， 是事件 A 每次试验中发生的概率，则对任意的正数 ，有p 1limpnkPn（1）若规定“向图中正方形随机投一块小石头落在四分之一单位圆里”为事件A 发生，则 A 发生的概率 应该等于四分之一单位圆面积，随机投 块石头就p n是独立重复做 次试验，事件 A 发生 次，由（1）式， 无限变大时 / 与nkk之差小于任意一个数 的概率趋于 1。p用计算机算一下就会

16、发现，即使 很大结果也不好，并且很不稳定，远不n如常规数值积分的几种方法。实际上，随机模拟法很少用来做定积分，它的特点是能够方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。2.2 蒙特卡罗方法计算定积分的理论原理通常有两类办法计算定积分。（1）随机投点法在“投石算面积”的例子中，事件 A 在每次试验中发生的概率 是四分p之一单位圆面积，即（2）dxdyxp102102次试验由计算机完成，采用0,1区间上的均匀分布产生相互独立的随机n数。记这样产生的 个点的坐标为 。事件 A 发生等价于n),.21),(niyxi满足 ， A 发生的次数是满足 点的个),(iyx21i

17、ix 2iix),.1(n数 。由伯努利定理， 可以用 / 近似替代。kpkn这种方法可以推广如下。对任意区间 内的连续函数 ，满足,ba)(f，为计算定积分 ，先由计算机产生 个点的坐标dxfc)(,c)0adxf)(n，其中 分别为 和 区间上的均匀分布随机数，,iy.21niyx,c蒙特卡罗随机模拟投点法在数值积分中的应用7然后记录 个点中满足 的数目 ，则n)(iixfykcabncdabba )()(（3）（2）均值估计法这种方法依据概率论的一下两个定理：大数定理（辛钦定理） 若随机变量 相互独立，服从同一2 nY,.21个分布，且具有数学期望 ,则对任意的正数 ，有iEY).,1(

18、n（4）lm1iinP随机变量函数的期望 若随机变量 的概率分布密度是2 X，则随机变量的函数 的期望为)(bxap)(XfYdxpfEba)(（5）于是，当 为 区间均匀分布的随机变量时，X,ba )/(1)ab， （ 5）式给出)(bxa)()(xfEabdxfba（6）只要产生 区间相互独立、均匀分布的随机数 ， 就,ba i),.21(n)(iixfy相互独立。根据辛钦定理，当 很大时期望 就可以用 的平nXfEYi均值近似，所以由（6）式得到niiba xfadxf1)()()(（7）与随机投点法相比，均值估计法没有 的限制，只需dxfc)(,c)0蒙特卡罗随机模拟投点法在数值积分中

19、的应用8计算 ，不需要产生随机数 ，也不需要作 的比较，显然大为方)(ixf iyiy)(ixf便。2.3 蒙特卡罗方法计算重积分的理论原理用随机模拟作数值积分的优点不仅在于计算简单，尤其是它可以方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。而通过上文对蒙特卡罗计算定积分的两种方法的比较，我们可知均值估计法比投点法更好操作，经过计算还可知均值估计法精确度（若选择的密度函数适合）会更高。重积分本文讨论均值估计法。考虑 维积分：k kkdxxfI .),.(. 2121其中 为 维积分区域， ( ， )。kViiiba:,.选取 上的一个概率密度函数 ， ( ，),.2

20、1kxpViiibxa:)，且 满足 ，ki,.21),.(21kxp 1.,(. 21 kkd，故：,.(21kx kkk dxxpxpfI .),.(),.(. 212121其中, ),.(21kxfE ),.(21k其中 , mi kxpf12),.( ),.(21kx（8）这里 是落在 区域中的随机投点数。m对于二重积分： dxyfI),(其中 ， 。记 ，设 ，bxa: gyxc)()(21 ,dcbaVS分别是 ， 的测度，若选取 上的均匀分布，则VS是,蒙特卡罗随机模拟投点法在数值积分中的应用9。由几何概率得 ，即 ，于是：y)(x, 1),(Syxg nmSVnSVdxygyx

21、fI),(,mif1),(iyxfS1,miVyxfnS1),(综上可得， （9）miVfnI1),(3 蒙特卡罗方法 MATLAB 编程实例用随机模拟作数值积分主要利用均匀分布产生的随机数及相应的判断、计算和简单运算。MATLAB 提供的命令unifrnd(a,b,m,n)产生 行 列 区间上均匀分布随机数。当 ， 时，可用 rand(m,n)。mn,ba0a1b例 1 用蒙特卡罗方法计算 。解(1) 随机投点法。 MATLAB 程序如下：n=10000;x=rand(2,n);k=0;for i=1:nif x(1,i)2+x(2,i)2=feval(fai1,x(k),y(k)=feval(fai2,x(k)s=s+feval(f,x(k),y(k);endends=s/n*(b-a)*(d-c);return function s=mtcx(f,fai1,fai2,a,b,c,d,n)if nargin8 n=10000;endx=unifrnd(a,b,1,n);y=unifrnd(c,d,1,n);