大数定律积分应用.doc

大数定律的应用
------------蒙特卡罗积分法
一、概述
（一）、大数定律：
设X1，X2,.…是相互独立的，服从同一分布的随机变量序列，且具有数学期望
E（Xk）=u(k=1,2,…).作前n个变量的算术平均值，则对于任意ε>0，有
①

等式①表明，当时这个事件当概率趋于1.即对于任意正数 ε，当n充分大时，不等式 | -µ|<ε成立的概率很大，通俗地说，辛钦大数定理是说，对于独立同分布且具有均值，的随机变量X1 ,…，Xn ,当n很大时他们的算术平均 很可能接近于.
（二）、蒙特卡罗法
蒙特学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法，与它对应的是确定性算法。
定积分的计算是Monte Carlo 方法引入计算数学的开端，在实际中，许多需要计算多重积分的复杂问题，用Monte Carlo 方法一般都能够很有效地予以解决，尽管Monte Carlo 方法给出的计算结果的精确度不是非常高，但它能很快地提供出一个低精度的模拟结果也是很有价值的。在多重积分计算中，由于Monte Carlo 方法的误差与积分重数无关，所以它比常用的均匀网格求积公式要优越。
二、基本原理
根据数学期望的定义，当x在（a,b)上满足均匀分布时，
②

我们在（a,b）之间取n个随机数，n的值非常大，当n的值大到一定程度时我们可以用
近似g(x)在（a,b)上的数学期望E[g(x)],由此我们再利用公式②，即可得到g(x)在区间（a,b)上的积分的近似值。
三、 实例分析
（一）sin函数的积分
首先计算一个比较简单的可以直接进行手算积分的函数的积分

1） 、基本思路
在（0，/2)之间取10000个随机数，然后计算出它们的平均值，由于10000个数据已经足够大，在要求不是非常高的情况下我们可以近似得认为它们的平均值就等于g(x)在（0，/2)之间的数学期望E[g(x)]。最后利用公式②可得
=（/2-0）E[g(x)]
2） 、c程序实现
#include "stdio.h"
#include "time.h"
#include "math.h"
#include "dos.h"
#define PI 3.1415926

void main()
{ double t,x,y,sum=0;
long i;
char end;
printf("wait...");
srand((int)time(0)); /*设置随机数种子，保证每次运行程序的结果不一样 */
for(i=0;i<=100000;i++)
{
x=rand()%10000; /* 取一个在0到10000范围内的整数 */
x=x/10000*PI/2; /*将随机数转换为（0，/2)上的随机数 */
t=sin(x); /*计算sin函数，即求每个随机数对应的函数值 */
sum=sum+t; /* 对随机数的函数值进行求和 */
}

printf("\n%f ",PI/2*sum/--i); /*计算出g(x)的数学期望乘以/2求出积分的*/
/* 似值并输出 */
scanf("%c",&end); /*使程序停在输出完结果的画面， */
}
3) 理论计算
由于sin函数的积分是比较简单的，我们可以直接对其进行积分，求得精确结果
=(-cos)=1
4）蒙特卡罗法估算结果分析
多次运行程序得到运行结果如下
0.983056
0.981828
0.982811
0.982450
0.981573
0.981261
0.981746
0.981057
0.982819
0.082183
从上面10次运行程序输出的结果可以看到用蒙特卡罗积分法得到的估计值基本稳定在
0.982附近，首先这个值已经比较精确，已经基本达到估算的要求。但令一方面，这个值并不是在准确值1附近浮动，而是普遍偏小一点点。
误差的来源应该是多方面的但主要在于以下几点
1、在c程序中产生的随机数并不是真正意义上的完美的随机数，我们利用
x=rand()%10000;
产生令从0到9999之间的整数，然后令
x=x/10000*PI/2;
从而得到（0，/2)区间上的随机数，很显然用这种方法得到的随机数是一些特定的离散的值，着必然会带来一定的误差。
2、在c程序中当数值的长度超过变量的位数是采用的是舍弃最后多余位的方法，这样必然导致每一个函数值都要比真实值略小，虽然这个误差非常小，但是在100000个抽样值中都存在这个问题，所以积少成多，最终的运算结果也会比真实值略小。

（二） 、多项式函数的积分
接下来计算一个积分范围不是从0开始的多项式积分

1）、基本思路
在（1,2)之间取10000个随机数，然后计算出它们的平均值，由于10000个数据已经足够大，在要求不是非常高的情况下我们可以近似得认为它们的平均值就等于g(x)在（1,2)之间的数学期望E[g(x)]。最后利用公式②可得
=（）E[g(x)]
2）、c程序实现
#include "stdio.h"
#include "time.h"
#include "math.h"
#include "dos.h"

void main()
{ double t,x,y,sum=0,result,average;
long i;
char end;
printf("wait...");
srand((int)time(0)); /*设置随机数种子，保证每次运行程序的结果不一样 */

for(i=0;i<=100000;i++)
{
x=rand()%10000; /* 取一个在0到10000范围内的整数 */
x=x/10000+1; /*将随机数转换为（1,2)上的随机数 */
t=x*x*x*x+2*x*x*x+4*x; /*计算该多项式函数，即求每个随机数对应的函数值 */
sum=sum+t; /* 对随机数的函数值进行求和 */
}
average=sum/--i; /* 对随机数的函数值进行求平均值 */
result=1*average; /*利用 =（）E[g(x)] */
/* 求得积分结果*/
printf("\n%f ",result); /*输出结果*/
scanf("%c",&end);
}

3) 理论计算
由于多项式函数的积分是比较简单的，我们可以直接对其进行积分，求得精确结果
=
=-
=19.7

4）蒙特卡罗法估算结果分析
多次运行程序得到运行结果如下
18.846181
18.832073
18.817216
18.806864
18.775687
18.790586
18,822360
18.874818
18.825712
18.836064
从以上的实验数据可以看出使用蒙特卡罗法估算的积分值和准确值是比较接近的，但还是存在一些误差，估算的结果普遍偏小一点，这可能和c程序中随机数的生成方式有关。另一方面在c程序中当数值的长度超过变量的位数是采用的是舍弃最后多余位的方法，这样必然导致每一个函数值都要比真实值略小，虽然这个误差非常小，但是在100000个抽样值中都存在这个问题，所以积少成多，最终的运算结果也会比真实值略小。
（三） 、复杂函数的积分

1）、基本思路
在（3，5)之间取10000个随机数，然后计算出它们的平均值，由于10000个数据已经足够大，在要求不是非常高的情况下我们可以近似得认为它们的平均值就等于g(x)在（1,2)之间的数学期望E[g(x)]。最后利用公式②可得

2) 、程序实现

#include "stdio.h"
#include "time.h"
#include "math.h"
#include "dos.h"

void main()
{ double t,x,y,sum=0,result,average;
long i;
char end;
printf("wait...");
srand((int)time(0));
for(i=0;i<=100000;i++)
{
x=rand()%10000;
x=x/10000+1;
t=x*x*x*x+2*x*x*x+4*x;
sum=sum+t;

}
average=sum/--i;
result=2*average;
printf("\n%f ",result);
scanf("%c",&end);
}
程序原理和上一个基本类似，不再对代码进行注释。
3）、蒙特卡罗法估算结果分析
多次运行程序得到运行结果如下
77.542380
77.875053
77.511215
77.511215
77.534110
77.616808
77.761919
77.680132
77.863036
77.800019
从以上的估算数据中可以发现，真实值基本上稳定在77.7附近,可以认为该复杂积分的积分结果就是77.7
（四） 、幂指数积分

1）、基本思路
在（1,3)之间取10000个随机数，然后计算出它们的平均值，由于10000个数据已经足够大，在要求不是非常高的情况下我们可以近似得认为它们的平均值就等于g(x)在（1,2)之间的数学期望E[g(x)]。最后利用公式②可得

2)、程序实现

#include "stdio.h"
#include "time.h"
#include "math.h"
#include "dos.h"

void main()
{ double t,x,y,sum=0,result,average;
long i;
char end;
printf("wait...");
srand((int)time(0));
for(i=0;i<=100000;i++)
{
x=rand()%10000;
x=x/10000+1;
t=pow(x,x);
sum=sum+t;

}
average=sum/--i;
result=2*average;
printf("\n%f ",result);
scanf("%c",&end);
}
程序原理和上一个基本类似，不再对代码进行注释。
3）、蒙特卡罗法估算结果分析
多次运行程序得到运行结果如下
12.757260
12.863323
12.804183
12.727121
12.893228
12.785942
12.793447
12.756400
12.874630
12.799330

从以上的估算数据中可以发现，真实值基本上稳定在12,8附近,可以认为该幂指数积分的积分结果就是77.7

总结：以上三个例子都是一元函数积分的情况，一元函数的积分程序基本类似，只需修改一下具体的函数即可，下面我们来研究比较复杂一些的多元函数积分的情况。

（四） 、二元积分

1） 、基本思路
在（1,4）之间取10000个x的随机数，在（2,5）之间取10000个y的随机数然后计算出每一对（x,y)对应的函数值g(x,y),接着求出它们的平均值，由于10000个数据已经足够大，在要求不是非常高的情况下我们可以近似得认为它们的平均值就等于g(xy)在之间的数学期望E[g(x)]。最后利用公式②可得

2) 程序实现

#include "stdio.h"
#include "time.h"
#include "math.h"
#include "dos.h"

void main()
{ double t,x,y,sum=0,result,average;
long i;
char end;
printf("wait...");
srand((int)time(0));
for(i=0;i<=100000;i++)
{
x=rand()%10000; /* 产生一个x的随机数 */
x=3*x/10000+1; /* 将x转换到（1,4）上 */
y=rand()%10000; /*产生一个y的随机数 */
y=3*x/10000+2; /*将y转换到（2,5）上 */
t=exp(sin(x))*pow(y,x); /* 计算每一对（x,y)对应的函数值 */
sum=sum+t; /* 求和 */

}
average=sum/--i; /* 求平均数 */
result=3*3*average; /* 求得最终结果 */
printf("\n%lf ",result);
scanf("%c",&end);
}

3） 、 蒙特卡罗法估算结果分析
76.023007
76.067155
76.075529
76.069914
75.982527
76.095394
75.946521
76.029258
76.058545
76.017462
从以上的估算数据中可以发现，真实值基本上稳定在76.0附近,可以认为该幂指数积分的积分结果就是76.0

四、 总结
通过以上的一些实例，我们可以发现用蒙特卡罗积分法对一些无法直接进行积分的函数计算积分值是一种简单可行的方法，虽然需要大量的抽样，但是对计算机来说这是非常简单的。另一方面虽然这种方法计算得到的积分值和实际存在一定的误差，但是只要加大抽样的数量，还是可以吧误差控制在一定范围里面的，所以蒙特卡罗积分法是一种简单可行的用计算机进行积分运算的方法。