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