蒙特卡罗方法初探——附C语言实例程序详解
2007-04-30 17:04
蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在本世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。
传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。这也是我们采用该方法的原因。
蒙特卡罗方法的基本原理及思想如下：
当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。
蒙特卡罗解题三个主要步骤：
构造或描述概率过程：
对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。
实现从已知概率分布抽样：
构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。
建立各种估计量：
一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。
蒙特卡罗方法与一般计算方法有很大区别，一般计算方法对于解决多维或因素复杂的问题非常困难，而蒙特卡罗方法对于解决这方面的问题却比较简单。其特点如下：
1、直接追踪粒子，物理思路清晰，易于理解。
2、采用随机抽样的方法，较真切的模拟粒子输运的过程，反映了统计涨落的规律。
3、不受系统多维、多因素等复杂性的限制，是解决复杂系统粒子输运问题的好方法。
4、MC程序结构清晰简单。
5、研究人员采用MC方法编写程序来解决粒子输运问题，比较容易得到自己想得到的任意中间结果，应用灵活性强。
6、MC方法主要弱点是收敛速度较慢和误差的概率性质，其概率误差正比于，如果单纯以增大抽
样粒子个数N来减小误差，就要增加很大的计算量。
近十年来，蒙特卡罗方法发展很快，从1983年到1988年期刊论文数量增长了五倍，有几本好书是关于电子? 光子蒙特卡罗问题的[注1]，蒙特卡罗方法的代码被认为是黑匣子，它已成为计算数学中不可缺少的组成部分，这主要是因为以下原因：
1、传统的分析方法受到了问题复杂性的限制。
2、MC方法直观，对实验者很有吸引力。
3、计算机变得更快更便宜。
4、量子理论的发展为我们提供了辐射与物质相互作用的截面数据。
（以上内容来自网络）
简而言之，蒙特卡罗方法就是用频率来代替概率，当实验样本足够大的时候，就可以得到比较精确的结果。
下面就利用蒙特卡罗方法作一个简单的任务，从算法和编程中来理解蒙特卡罗方法的基本原理。本人从接触到蒙特卡罗方法到撰写这篇日志刚好一天的时间，其中有很多误解和不当之处，敬请各路高手不吝赐教。
题目：利用蒙特卡罗方法模拟计算圆周率PI。
分析：我们根据计算圆的面积公式来求解圆周率PI，当圆的半径为1时，圆的面积刚好等于圆周率PI。
根据蒙特卡罗方法的思想，我们在以坐标原点为圆心作一个半径为1的单位圆，再作一个正方形与此圆相切（其实是否相切不重要，此处只是为了计算方便）。在这个正方形内任取M个点，并比较判断是否落在圆内，将落在圆内的点计作N，那么M与N的比值就可以看成正方形和圆的面积之比。即圆周率PI=4*N/M。
下面是本题目的C语言源程序：
/*------------------------------------*/
#include"stdio.h"
#include"time.h"
#include"stdlib.h"
#include"graphics.h"
main()
{int gd=DETECT,gm=0;
long m=0,n=0,i;
double xi,yi,y,y1;
initgraph(&gd,&gm,""); /*设置图形显示模式*/
setbkcolor(BLUE);
/*-----------画坐标、正方形和圆----------------*/
line(200,50,200,400);line(200,50,205,60); line(200,50,195,60);
line(50,200,400,200);line(400,200,390,195); line(400,200,390,205);
rectangle(100,100,300,300);
circle(200,200,100);
/*--------------------------------------------*/
/*----------用时间作为产生随机数的种子--------*/
srand((unsigned)time(NULL));
for(i=0;i
{xi=(rand()%(1000-0+1)+0)/500.0-1.0;
yi=(rand()%(1000-0+1)+0)/500.0-1.0;
putpixel((int)(xi*100+200),(int)(yi*100+200),(int)RED); /*在正方形区域内画点*/
/*printf("xi=%f,yi=%f\n",xi,yi);*/
if(yi*yi>=(-(1-xi*xi))&&yi*yi
m++,n++;
else m++;
}
y=4.0*n/m; /*计算PI值*/
/*outtextxg(100,100,"PI=");*/
printf("\n\tPI=%f\n",y);
getch();
closegraph();
}
/*-------------------------------------------------*/
程序十分简单，即产生随机数、画点、判别区域、统计，最后根据统计结果模拟计算PI值。 我调试的结果如下：
循环1000次：3.232 3.060 3.072 3.080 3.224
10000次：3.114 3.134 3.137 3.162 3.136
10万次：3.145 3.150 3.151 3.147 3.146
100万次：3.143 3.146 3.145 3.147 3.147
1000万次：3.146 （时间长，只测了一次，呵呵）
结果分析：
1、循环的次数越多，结果越精确。
2、C语言rand()函数产生的并不是真正懂得随机数，而是伪随机数，会影响结果的精度。有兴趣的朋友不妨使用其它能产生真正意义的语言或算法试试。
3、算法过于简单，计算机在处理数据时会产生很多的数据丢失，影响结果精度。
此文作为学习交流之用，程序在TC2.0环境下调试通过。时间仓促，语焉不详，望各路高手不吝赐教。
#include
#include
#include
#define L 40
double funiseed (double iseed) //同余法产生随机数 {
double mult=1277,modulo,rmodulo;
modulo=rmodulo=pow(2,17);
iseed=iseed*mult;
if(iseed
else iseed=iseed-modulo;
iseed=iseed/rmodulo;
return iseed;
}
main()
{
int i,j,la[L+1][L+1],ip[L+1],im[L+1],count=0;
int ici,ien,mcs,mcsmax,n0,mcstep;
float w[9],m,a;
double jkt,iseed;
//初始化数组 for(i=0;i
for(j=0;j
la[i][j] = -1;
a=1./(L*L);
}
}
//************************************
for(i=1;i
ip[i]=i+1;
im[i]=i-1;
}
ip[L]=1;
im[1]=L;
//***********************************
printf("输入样本步数：");
scanf("%d",&mcstep);
printf("输入n0的值：");
scanf("%d",&n0);
printf("输入jkt的值：");
scanf("%lf",&jkt);
蒙特卡罗方法初探——附C语言实例程序详解
2007-04-30 17:04
蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在本世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。
传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。这也是我们采用该方法的原因。
蒙特卡罗方法的基本原理及思想如下：
当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。
蒙特卡罗解题三个主要步骤：
构造或描述概率过程：
对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。
实现从已知概率分布抽样：
构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。
建立各种估计量：
一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。
蒙特卡罗方法与一般计算方法有很大区别，一般计算方法对于解决多维或因素复杂的问题非常困难，而蒙特卡罗方法对于解决这方面的问题却比较简单。其特点如下：
1、直接追踪粒子，物理思路清晰，易于理解。
2、采用随机抽样的方法，较真切的模拟粒子输运的过程，反映了统计涨落的规律。
3、不受系统多维、多因素等复杂性的限制，是解决复杂系统粒子输运问题的好方法。
4、MC程序结构清晰简单。
5、研究人员采用MC方法编写程序来解决粒子输运问题，比较容易得到自己想得到的任意中间结果，应用灵活性强。
6、MC方法主要弱点是收敛速度较慢和误差的概率性质，其概率误差正比于，如果单纯以增大抽
样粒子个数N来减小误差，就要增加很大的计算量。
近十年来，蒙特卡罗方法发展很快，从1983年到1988年期刊论文数量增长了五倍，有几本好书是关于电子? 光子蒙特卡罗问题的[注1]，蒙特卡罗方法的代码被认为是黑匣子，它已成为计算数学中不可缺少的组成部分，这主要是因为以下原因：
1、传统的分析方法受到了问题复杂性的限制。
2、MC方法直观，对实验者很有吸引力。
3、计算机变得更快更便宜。
4、量子理论的发展为我们提供了辐射与物质相互作用的截面数据。
（以上内容来自网络）
简而言之，蒙特卡罗方法就是用频率来代替概率，当实验样本足够大的时候，就可以得到比较精确的结果。
下面就利用蒙特卡罗方法作一个简单的任务，从算法和编程中来理解蒙特卡罗方法的基本原理。本人从接触到蒙特卡罗方法到撰写这篇日志刚好一天的时间，其中有很多误解和不当之处，敬请各路高手不吝赐教。
题目：利用蒙特卡罗方法模拟计算圆周率PI。
分析：我们根据计算圆的面积公式来求解圆周率PI，当圆的半径为1时，圆的面积刚好等于圆周率PI。
根据蒙特卡罗方法的思想，我们在以坐标原点为圆心作一个半径为1的单位圆，再作一个正方形与此圆相切（其实是否相切不重要，此处只是为了计算方便）。在这个正方形内任取M个点，并比较判断是否落在圆内，将落在圆内的点计作N，那么M与N的比值就可以看成正方形和圆的面积之比。即圆周率PI=4*N/M。
下面是本题目的C语言源程序：
/*------------------------------------*/
#include"stdio.h"
#include"time.h"
#include"stdlib.h"
#include"graphics.h"
main()
{int gd=DETECT,gm=0;
long m=0,n=0,i;
double xi,yi,y,y1;
initgraph(&gd,&gm,""); /*设置图形显示模式*/
setbkcolor(BLUE);
/*-----------画坐标、正方形和圆----------------*/
line(200,50,200,400);line(200,50,205,60); line(200,50,195,60);
line(50,200,400,200);line(400,200,390,195); line(400,200,390,205);
rectangle(100,100,300,300);
circle(200,200,100);
/*--------------------------------------------*/
/*----------用时间作为产生随机数的种子--------*/
srand((unsigned)time(NULL));
for(i=0;i
{xi=(rand()%(1000-0+1)+0)/500.0-1.0;
yi=(rand()%(1000-0+1)+0)/500.0-1.0;
putpixel((int)(xi*100+200),(int)(yi*100+200),(int)RED); /*在正方形区域内画点*/
/*printf("xi=%f,yi=%f\n",xi,yi);*/
if(yi*yi>=(-(1-xi*xi))&&yi*yi
m++,n++;
else m++;
}
y=4.0*n/m; /*计算PI值*/
/*outtextxg(100,100,"PI=");*/
printf("\n\tPI=%f\n",y);
getch();
closegraph();
}
/*-------------------------------------------------*/
程序十分简单，即产生随机数、画点、判别区域、统计，最后根据统计结果模拟计算PI值。 我调试的结果如下：
循环1000次：3.232 3.060 3.072 3.080 3.224
10000次：3.114 3.134 3.137 3.162 3.136
10万次：3.145 3.150 3.151 3.147 3.146
100万次：3.143 3.146 3.145 3.147 3.147
1000万次：3.146 （时间长，只测了一次，呵呵）
结果分析：
1、循环的次数越多，结果越精确。
2、C语言rand()函数产生的并不是真正懂得随机数，而是伪随机数，会影响结果的精度。有兴趣的朋友不妨使用其它能产生真正意义的语言或算法试试。
3、算法过于简单，计算机在处理数据时会产生很多的数据丢失，影响结果精度。
此文作为学习交流之用，程序在TC2.0环境下调试通过。时间仓促，语焉不详，望各路高手不吝赐教。
#include
#include
#include
#define L 40
double funiseed (double iseed) //同余法产生随机数 {
double mult=1277,modulo,rmodulo;
modulo=rmodulo=pow(2,17);
iseed=iseed*mult;
if(iseed
else iseed=iseed-modulo;
iseed=iseed/rmodulo;
return iseed;
}
main()
{
int i,j,la[L+1][L+1],ip[L+1],im[L+1],count=0;
int ici,ien,mcs,mcsmax,n0,mcstep;
float w[9],m,a;
double jkt,iseed;
//初始化数组 for(i=0;i
for(j=0;j
la[i][j] = -1;
a=1./(L*L);
}
}
//************************************
for(i=1;i
ip[i]=i+1;
im[i]=i-1;
}
ip[L]=1;
im[1]=L;
//***********************************
printf("输入样本步数：");
scanf("%d",&mcstep);
printf("输入n0的值：");
scanf("%d",&n0);
printf("输入jkt的值：");
scanf("%lf",&jkt);