蒙特卡洛模拟，准确的说它并不是一种算法，而是一种解决问题的思路，因为算法的实现代码是确定的，而蒙特卡洛模拟的实现代码并不确定，针对不同问题的求解，其实现代码是不一样的，其核心思路是根据要求解的问题建立一个概率事件，并确定该事件的随机输入变量，该然后对这个概率事件有限次地输入随机数，统计输出结果，从而根据输出结果的统计特征得到问题的近似解。下面举一个很简单很容易理解的例子，来讲解蒙特卡洛模拟的思路。

问题：如上图，有一个矩形操场，已知其长100米，宽50米，该操场中间画了一条不规则的红线，求红线左侧和右侧区域的面积分别是多少？

总体求解思路：只要求得两个区域的面积之比，即可以得到各自的面积。

蒙特卡洛模拟思路：

1. 概率事件：每一个学生走进操场时都随机选择一个位置待着，因此他有可能待在操场的任何一个角落（假设同一位置可以待多个人）。

2. 随机输入变量：每一个学生。

3. 有限次输入随机数：有500个学生，依次走进操场，并随机地选择一个位置待着。

4. 统计输出结果：统计待在左侧区域的学生总数为n1，待在右侧区域的学生总数为n2

5. 根据输出结果的统计特征得到问题的近似解：因为每个学生待在左侧还是右侧是随机的，可以认为左侧与右侧的学生人数比n1/n2近似等于左右侧的面积比，所以可以得到左侧面积为100*50*n1/(n1+n2)，右侧面积为100*50*n2/(n1+n2)。

在蒙特卡洛模拟中，模拟次数越多，得到的结果越接近真实解，前提条件是输入的随机变量足够随机，如果变量不够随机，很可能导致离谱的错误。

下面再举一个例子，并使用C++编码进行蒙特卡洛模拟，以求得问题的近似解。

如上图所示，在y=sinx曲线的[0, π]区间中，求解曲线下方阴影部分的面积，也即相当于求函数y=sinx在该区间的定积分。蒙特卡洛模拟的思路与步骤：

1. 概率事件：随机生成x坐标在[0, π]区间、y坐标在[0, 1]区间的点，该点有可能落在阴影区域，也有可能落在阴影之外的区域。

2. 随机变量：点的x坐标与y坐标。

3. 有限次输入随机数：生成多个点，比如100个，1000个，10000个......

4. 统计输出结果：统计落在阴影区域内的点数，以及落在阴影区域之外的点数。

5. 根据落在阴影区域内部的点数与落在外部的点数之比，得到矩形中阴影区域与非阴影区域的近似面积比，从而求得阴影部分的近似面积。

首先，我们来计算一下真实解，阴影部分的面积可按牛顿-莱布尼茨公式计算（这里可能有人会疑惑了，电脑为什么不直接用这个公式计算，我们要知道，电脑是很笨的，它是不会直接使用这个公式来计算定积分的）。所以，问题的真实解是2。

其次，说明一下怎么判断点落在阴影区域内部还是外部。假设通过随机生成x坐标与y坐标，得到点P(x, y)，那么可如此判断：当y≤sinx时，点P落在阴影内部，反之当y>sinx时，点P落在阴影外部。

下面上代码：

   

    

     

      

     
     

      

       //生成[min, max]区间的浮点数
      
     
    

     

      

     
     

      

       float get_random(float min, float max)
      
     
    

     

      

     
     

      

       {
      
     
    

     

      

     
     

      

         
       float x = rand() / 
       double(RAND_MAX);    
       //生成0~1的浮点数
      
     
    

     

      

     
     

      
  
       return (x*(max - min) + min);   
      
     
    

     

      

     
     

      

       }
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

       //蒙特卡洛模拟，cnt为模拟次数
      
     
    

     

      

     
     

      

       void MonteCarlo_test(int cnt)
      
     
    

     

      

     
     

      

       {
      
     
    

     

      

     
     

      
  
       int n1 = 
       0;   
       //阴影区域内部点的计数器
      
     
    

     

      

     
     

      
  
       int n2 = 
       0;   
       //阴影区域外部点的计数器
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
  
       const 
       float PI = 
       3.141592;
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
  
       float s = PI*
       1.0;
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

         srand((unsigned)time(NULL));  
       //设置随机数种子
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
  
       for (
       int i = 
       0; i < cnt; i++)
      
     
    

     

      

     
     

      

         {
      
     
    

     

      

     
     

      
    
       float x = get_random(
       0, PI);  
       //随机生成在[0, π]区间的x坐标
      
     
    

     

      

     
     

      
    
       float y = get_random(
       0, 
       1);   
       //随机生成在[0, 1]区间的y坐标
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
    
       float fx = sin(x);   
       //计算sinx
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

           
       if (y <= fx)  
       //如果y≤sinx则认为点在阴影内部
      
     
    

     

      

     
     

      

             n1++;    
       //阴影区域内部点计数器加1
      
     
    

     

      

     
     

      

           
       else      
       //如果y>sinx则认为点在阴影内部
      
     
    

     

      

     
     

      

             n2++;    
       //阴影区域外部点计数器加1
      
     
    

     

      

     
     

      

         }
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

         
       float s1 = n1 * 
       1.0 / cnt * s;   
       //n1/(n1+n2)*s = n1/cnt*s，即为阴影部分面积
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      
 
      
     
    

     

      

     
     

      

         printf(
       "模拟次数：%d，阴影部分面积：%f\n", cnt, s1);
      
     
    

     

      

     
     

      

       }
      
     
   

运行以上代码，逐渐增加模拟次数cnt，得到以下结果，可以看到，随着模拟次数的增加，得到的近似解会越来越接近真实解。

查看评论