蒙特卡罗算法（或蒙特卡洛方法）- Monte Carlo method

是以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。

以概率和统计的理论、方法为基础的一种计算方法，将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解，故又称统计模拟法或统计试验法。 

蒙特卡洛方法的基本思想

 

当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。 有一个例子可以使你比较直观地了解蒙特卡洛方法：假设我们要计算一个不规则图形的面积，那么图形的不规则程度和分析性计算（比如，积分）的复杂程度是成正比的。蒙特卡洛方法是怎么计算的呢？假想你有一袋豆子，把豆子均匀地朝这个图形上撒，然后数这个图形之中有多少颗豆子，这个豆子的数目就是图形的面积。当你的豆子越小，撒的越多的时候，结果就越精确。在这里我们要假定豆子都在一个平面上，相互之间没有重叠。

 

蒙特卡洛方法的工作过程

 

在解决实际问题的时候应用蒙特卡洛方法主要有两部分工作：

1. 用蒙特卡洛方法模拟某一过程时，需要产生各种概率分布的随机变量。
2. 用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。

蒙特卡洛方法分子模拟计算的步骤

使用蒙特卡洛方法进行分子模拟计算是按照以下步骤进行的：

1. 使用随机数发生器产生一个随机的分子构型。
2. 对此分子构型的其中粒子坐标做无规则的改变，产生一个新的分子构型。
3. 计算新的分子构型的能量。
4. 比较新的分子构型于改变前的分子构型的能量变化，判断是否接受该构型。
   • 若新的分子构型能量低于原分子构型的能量，则接受新的构型，使用这个构型重复再做下一次迭代。
   • 若新的分子构型能量高于原分子构型的能量，则计算玻尔兹曼因子，并产生一个随机数。
     • 若这个随机数大于所计算出的玻尔兹曼因子，则放弃这个构型，重新计算。
     • 若这个随机数小于所计算出的玻尔兹曼因子，则接受这个构型，使用这个构型重复再做下一次迭代。
5. 如此进行迭代计算，直至最后搜索出低于所给能量条件的分子构型结束。

蒙特卡洛方法应用

圆周率

蒙特卡洛方法可用于近似计算圆周率：让计算机每次随机生成两个0到1之间的数，看以这两个实数为横纵坐标的点是否在单位圆内。生成一系列随机点，统计单位圆内的点数与总点数，（圆面积和正方形面积之比为PI:4，PI为圆周率），当随机点取得越多（但即使取10的9次方个随机点时，其结果也仅在前4位与圆周率吻合）时，其结果越接近于圆周率。实际上，计算机产生的随机数只能精确到某位数，并不能产生任意实数（例如无理数等等）；上述做法将平面分割成一个个网格，由此计算出来的面积当然与圆或多或少有差距。

蒙特卡罗算法 - 应用程序---Monte Carlo法( 蒙特卡罗)积分

C源代码：

#include<iostream.h>
#include<time.h>

/* srand(),rand() */
#include<math.h>

/* sin() */
#include<stdlib.h>

/* RAND_MAX */

using namespace std ;

double FuncOne(double x);
/*被积一元函数*/
double MonteCarloOne(double(*p)(double),double a,double b,double L,double H);
/*一元函数积分调用接口*/
double FuncTwo(double x,double y);
/*被积二元函数*/
double Func1(double x);
/* 积分下限y=Func1(x) */
double Func2(double x);
/* 积分上限y=Func2(x) */
double MonteCarloTwo(double(*func1)(double),double(*func2)(double),double(*p)(double,double),double a,double b,double l,double h,double L,double H);
/*二元函数积分调用接口.在a
double FuncThree(double x,double y,double z);
/*被积分三元函数*/
double Func3(double x,double y);
/* 积分下限z=Func3(x,y) */
double Func4(double x,double y);
/* 积分上限z=Func4(x,y) */
double MonteCarloThree(double(*func1)(double),double(*func2)(double),double(*func3)(double,double),double(*func4)(double,double),double(*p)(double,double,double),double a,double b,double c,double d,double l,double h,double L,double H);
/*三元函数积分调用接口.在a

double FuncOne(double x)
{
    return sin(x);
}

double MonteCarloOne(double(*p)(double),double a,double b,double L,double H)
{
    unsigned N=100000000 ;
    /* 1e8 */
    double x,y ;
    /*随机点*/
    double t ;
    /*换元*/
    unsigned M=0 ;
    /*处于阴影部分内的点的个数*/
    srand((unsigned)time(NULL));
    /*产生随机数种子*/
    for(unsigned int i=0;i
    {
        t=double(rand())/RAND_MAX ;
        /* 产生0~1的随机浮点数 */
        y=double(rand())/RAND_MAX ;
        x=a+(b-a)*t ;
        /* 得到a~b的随机浮点数 */
        /*随机点处于阴影部分内*/
        if(y*(H-L)<(p(x)-L))M++;
    }
    double result=(H-L)*(b-a)*M/N+(b-a)*L ;
    /* double型乘以unsigned型为double型 */
    return result ;
    /*积分结果*/
}

double FuncTwo(double x,double y)
{
    return x+y ;
}

double Func1(double x)
{
    return x ;
}

double Func2(double x)
{
    return x*x ;
}

double MonteCarloTwo(double(*func1)(double),double(*func2)(double),double(*p)(double,double),double a,double b,double l,double h,double L,double H)
{
    unsigned N=100000000 ;
    /* 1e8 */
    double x,y,z ;
    /*随机点*/
    double t,s ;
    /*换元*/
    unsigned m=0,M=0 ;
    /*处于阴影部分内的点的个数*/
    srand((unsigned)time(NULL));
    /*产生随机数种子*/
    for(unsigned int i=0;i
    {
        t=double(rand())/RAND_MAX ;
        /* 产生0~1的随机浮点数 */
        s=double(rand())/RAND_MAX ;
        x=a+(b-a)*t ;
        /* 得到a~b的随机浮点数 */
        y=l+(h-l)*s ;
        /* 得到l~h的随机浮点数 */
        if(y>func1(x)&&y
        {
            m++;
            /* 注意,m是与一元积分相对应的 */
            z=double(rand())/RAND_MAX ;
            /*随机点处于阴影部分内*/
            if(z*(H-L)<(p(x,y)-L))M++;
        }
    }
    double result=(H-L)*(b-a)*(h-l)*M/N+(b-a)*(h-l)*L*m/N ;
    /* double型乘以unsigned型为double型 */
    return result ;
}

double FuncThree(double x,double y,double z)
{
    return x+y+z ;
}

double Func3(double x,double y)
{
    return 0 ;
}

double Func4(double x,double y)
{
    return FuncTwo(x,y);
}

double MonteCarloThree(double(*func1)(double),double(*func2)(double),double(*func3)(double,double),double(*func4)(double,double),double(*p)(double,double,double),double a,double b,double c,double d,double l,double h,double L,double H)
{
    unsigned N=100000000 ;
    /* 1e8 */
    double x,y,z,u ;
    /*随机点*/
    double t,s,k ;
    /*换元*/
    unsigned m=0,M=0 ;
    /*处于阴影部分内的点的个数*/
    srand((unsigned)time(NULL));
    /*产生随机数种子*/
    for(unsigned int i=0;i
    {
        t=double(rand())/RAND_MAX ;
        /* 产生0~1的随机浮点数 */
        s=double(rand())/RAND_MAX ;
        k=double(rand())/RAND_MAX ;
        x=a+(b-a)*t ;
        /* 得到a~b的随机浮点数 */
        y=c+(d-c)*s ;
        /* 得到c~d的随机浮点数 */
        if(y>func1(x)&&y
        {
            k=double(rand())/RAND_MAX ;
            z=l+(h-l)*k ;
            /* 得到l~h的随机浮点数 */
            /*随机点处于阴影部分内*/
            if(z>func3(x,y)&&z
            {
                m++;
                /* 注意,m是与二元积分对应的 */
                u=double(rand())/RAND_MAX ;
                if(u*(H-L)<(p(x,y,z)-L))
                M++;
            }
        }
    }
    double result=(H-L)*(b-a)*(d-c)*(h-l)*M/N+(b-a)*(d-c)*(h-l)*L*m/N ;
    /* double型乘以unsigned型为double型 */
    return result ;
}

int main()
{
    double a,b ;
    /* x积分限 */
    double l,h ;
    /* y的值域或y积分限 */
    double result ;
    /*积分结果*/
    
    a=0,b=3.1415926 ;
    l=0,h=1 ;
    result=MonteCarloOne(FuncOne,a,b,l,h);
    cout<<<"\N" ;
    
    a=1.0,b=2.0 ;
    l=1.0,h=4.0 ;
    double L=2.0,H=6.0 ;
    /* z的值域 */
    result=MonteCarloTwo(Func1,Func2,FuncTwo,a,b,l,h,L,H);
    cout<<<"\N" ;
    
    a=1.0,b=2.0 ;
    double c=1.0,d=4.0 ;
    l=0.0,h=6.0 ;
    L=2.0,H=12.0 ;
    result=MonteCarloThree(Func1,Func2,Func3,Func4,FuncThree,a,b,c,d,l,h,L,H);
    cout<<<"\N" ;
    
    return 0 ;
}

 

  运行结果:
  2.00011
  3.35045
  20.994
  准确值分别为:

  2
  3.35
  20.9964286