对蒙特卡洛模拟的理解

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

问题：如上图，有一个矩形操场，已知其长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，得到以下结果，可以看到，随着模拟次数的增加，得到的近似解会越来越接近真实解。

模拟次数	得到结果
10	1.256637
50	2.261946
100	2.042035
500	2.035752
1000	1.976061
10000	2.016588
20000	2.007320
50000	2.000440