多元统计课程设计momo.doc

用频率法计算二重积分 1 蒙特卡罗模拟思想初步1.1 不规则图形面积的求解 给定一个不规则图形，怎样求解它的面积？这是一个比较棘手的问题。通常我们都会想到“网格法”，即将不规则图形以微小的正方形“网格”化把不规则图形分化成多个小面积的方形，边界部分进行割补并计数：不足半格的计数为0，超过半格的为1。最后将所有的格数加起来，乘以方形面积即可。计算结果的精确度与网格的细密度直接相关。但是，网格的细化程度越高，操作起来便越不方便，如若仅靠眼睛来记录网格数目，实不为一种可取方法。那么，我们可否寻求其他方法求面积呢？先看下面一个例子：古代一位智者求图形面积，他先将图形画在一块方形布上，然后将一袋装有很多大小与重量大致相同的豆子随意倒在布上面，那么落入图形里的豆子个数（重量）与总豆子个数（重量）的比值乘以方布面积就是所求图形面积（见图1）。 图1 从表面上看，这种方法似乎很奏效，方便简洁；从理论上，我们又该怎样解释这种方法的实用性呢？1.2 蒙特卡罗模拟思想 我们先来解释上面方法的合理性。我们可以这样认为：随机抛出一颗豆子，它必定落入方形中，但却未必能落入所画图形中。但是“落入图形内部”这件事必定以一定概率发生。并且图形面积越大，豆子落入其中的可能性越大，“落入图形内部”这件事发生的概率也越大。不难想到概率等于图形面积与方形面积之比。我们也知道，在一定条件下，某一事件在多次试验中出现的频率可以近似看作该事件发生的概率。换句话说，我们可以用频率估计概率。现在我们将一袋豆子随机倒出，豆子总数N可看作试验次数，“落入图形内部”的豆子数n可视作事件发生的次数，很直然地就有频率为n/N，以n/N为豆子“落入图形内部”的概率，则得出上述结论。在这里，我们应该注意几点：a.应保证豆子的“随机性”，如豆子本身性质（形状、重量等）应基本一致且愈小愈好；豆子抛出方向是随机的等。b.豆子的数量应该足够多。用频率估计概率有个前提条件，即试验次数足够多，越多越精确。 在此，我们给出了一种与频率有关的新的计算图形面积的算法。现实中由于种种不确定因素导致计算结果不甚精确，如豆子大小影响随机性。试想，我们可否对这种现实情景进行虚拟的计算机模拟？在计算机里，豆子可由没有大小重量的随机点替代以确保随机性，随机点可由计算机尽可能多的产生以确保足够多数量，这样就可以是计算更为精确。我们称这种模拟思想为蒙特卡罗模拟思想。下面提出蒙特卡罗模拟确定性模型。2 蒙特卡罗模拟确定性模型2.1 曲变梯形面积的近似计算原理对于图示的函数y=f(x) 图像: y y=M y=f(x) .P(x,y) 0 a b x图2怎样求解y=f(x)在区间a,b上的一段与x轴所围部分面积？这实际上就是个一重定积分的求解问题。我们可以这样做：取一个矩形（见图二），使其包含曲边梯形。由计算机产生位于矩形里的随机点，判断该点与曲线的位置关系（曲线下方或上方）；产生足够多的随机点N个，记下落入曲线下方或曲线上的点的数目n，则由大数定律，当N足够大时有下列近似公式： 曲边梯形面积矩形面积nN这样就得到曲边梯形面积的近似公式：曲边梯形面积(b-a)*M*n/N2.2 蒙特卡罗模拟曲边梯形面积算法：下面给出蒙特卡罗模拟曲边梯形面积算法：(1)赋初值:曲线下方点数n=0;输入模拟所需总的试验次数（即总的随机点数）N.(2)对i=1,2,N,执行以下三步: 产生随机数x(i),y(i); 计算f(x(i); 如果y(i)f(x(i),则n=n+1;否则,转.(3)计算曲边梯形的面积:S=(b-a)*M*n/N. 3 基于Monto Carlo模拟思想的二重积分求解3.1 问题的提出 上一小节提到用蒙特卡罗模拟思想求任意区边梯形的面积，实际上也就解决了一重定积分的求解问题。那么，我们也能否根据这种模拟思想，求某个立体图形的体积呢？如果能，我们又该如何模拟呢？进一步来说，如果我们可以求得体积，那么我们就解决了二重定积分的求解问题了。3.2 问题的分析 同样，我们可以依照平面（二维）上的蒙特卡罗模拟，进行空间（三维）上的蒙特卡罗模拟：用一个立方体“罩住”立体图形后，计算机产生位于立方体内部的随机点，点可能落在立体图形的外部，也可能是内部。然后根据落在内部的点的数目、总的随机点数目及立方体的体积得出目标立体图形的体积。由于计算机可以很好很方便地产生随机点，我们有理由相信，这将是种可行的模拟方法。3.3 实例求解现在我们看具体实例。用频率法求解下列二重积分： e-x2+y22dxdy，其中积分区域为圆域 x2+y24。要求给出适合的试验次数，使得误差精度以一定程度保证在0.01内，并说明其原理。 3.3.1 二重积分的估计值 我们先用MATLAB编程绘制出二维函数z=e-x2+y22 在圆域的图像，为了便于清楚观察该曲面，我们同时给出方位角为0度，仰角为90度的视点所观察到的视野图（程序见附录）： 图3 显然，该曲面与五个平面都相切于一点，这样就能让我们更方便地找的一个立方体“罩住”这个立体图形：六个平面x=-2，x=2，y=-2，y=2，z=0所围成的长方体。我们给出这个立方体“罩住”曲面的图（图4）： 图4类似曲边梯形面积算法，我们也给出本例蒙特卡罗模拟立体图形体积算法：(1)赋初值:曲面下方点数n=0;输入模拟所需总的试验次数（即总的随机点数）N.(2)对i=1,2,N,执行以下三步: 产生随机数x(i),y(i),z(i),使得随机点( x(i),y(i),z(i)位于立方体内部或立方体上; 计算z(x(i),y(i); 如果x(i)2+y(i)24且z(i)z(x(i),y(i),则n=n+1;否则,转.(3)计算曲面与xoy平面所围成的立体图的体积:V=16*n/N.在这里给定N=1000,为了减小与精确值的偏差，我们编程重复二十次，取二十次的平均值（程序见附录），得到体积V= 5.4880，也即二重积分的估计值。 3.3.2 用频率法求二重积分的精度 上面只给出了N=1000时的估计值，即使取二十次的平均值，也不能保证精度。但我们知道N越大，精度就越高，结果就越接近精确值，下一步要找到这里的N与精度的关系。譬如说，要以概率保证误差精度在内，产生的随机点个数（即试验次数）N应满足什么条件？这里就要运用概率论与数理统计里的一个重要定理中心极限定理。我们知道，每产生一个随机点，其落在目标立体图形内部的可能性（即概率）Q等于两个体积之比。产生一个随机点，落于内部或外部用随机变量X表示，记“落于内部”为“X=1”，“落于外部”为“X=0”。产生第k个随机点的结果用Xk表示，现在独立地产生N个随机点，Xk（k=1,2,N）与X同分布，且期望和方差分别为Q和Q(1-Q)，k=1NXk 则恰好表示N次试验中“落于内部”的点的个数n，由中心极限定理知：k=1NXk-NQNQ(1-Q) 渐近地服从标准正态分布N（0,1）。保证精度在内等价于k=1NXkN*VQ*V，亦即k=1NXk-NQNQ(1-Q)V*NQ(1-Q)，这样“以概率保证误差精度在内”等价于P( k=1NXk-NQNQ(1-Q)4);z(i)=nan;subplot(1,2,1);mesh(x,y,z);xlabel(x轴);ylabel(y轴);zlabel(z轴);text(0,1,1,曲面z=e-0.5(x2+y2);subplot(1,2,2);mesh(x,y,z);view(0,90);text(0,0,1,改变视点后视野图)立方体“罩住”曲面图像程序：x y=meshgrid(-2:0.01:2);z1=exp(-(x.2+y.2)/2);i=find(x.2+y.24);z1(i)=nan;mesh(x,y,z1);hold on;z2=ones(size(x);mesh(x,y,z2)求二重积分估计值程序：N=1000;x=;y=;z=;n=0;for i=1:N x(i)=unifrnd(-2,2,1,1); y(i)=unifrnd(-2,2,1,1); z(i)=rand; if exp(-(x(i)2+y(i)2)/2)=