蒙特卡洛算法简单理解与demo

所谓蒙特卡洛算法，实际上就是用频率估计概率。

首先我们知道一个边长为2的正方形面积是2*2=4，他的内接圆的面积是π*1，那么我们在这样一个正方形内随机生成10000个点，落在圆里的点的个数/10000就应该是π/4，所以我们可以估计出π是落在圆里的点的个数*4/10000.

直接看程序:

import numpy as npradius=1test_num=100000#random generate 10000 points in a sqare width is 2,and the center is in originsqare_x=2*np.random.random_sample(test_num)-1sqare_y=2*np.random.random_sample(test_num)-1incircle_point_num=0#count how many points in the incirclefor point_count in range(len(sqare_x)):    if(sqare_x[point_count]*sqare_x[point_count]+               sqare_y[point_count]*sqare_y[point_count]<               radius*radius):#judge if the distance lower than radius        incircle_point_num+=1print 'PI: '+str(4.0*incircle_point_num/test_num)

这里随机生成了10000个在边长为2，中心在原点的正方形中的点，根据坐标判断是否在内接圆内，统计内接圆内点的个数，最后根据公式求出π。

准确度基本上随着样本个数增加而增加，下面是测试的一组统计：

测试个数 算出的π的值104.01003.210003.192100003.13881000003.1372810000003.142496100000003.14147268    

再测下去计算太慢，也没有意义了，毕竟我不是真的想知道π是多少。

根据统计结果可以看出来，测试样本越大，精确度也就越高，频率也就越接近概率值。

所以这套蒙特卡洛方法可以用来求出积分，也是使用积分区域的面积和外接长方形的面积比的方式估算。

例如：

求

 ${\int }_0^1{e}^{-x^2/2}/\sqrt{2\pi }dx$

可以知道这是一个单调递减函数，最大值就是x=0的时候也就是1/√2π

π，最小值是x=1的时候，也就是1/e^(0.5)/√2π

那么外接长方形面积就应该是1/√2π-1/e^(0.5)/√2π

所以积分面积=长方形面积*积分内点的个数/样本个数

import numpy as npimport math#calculate the functiondef f(x):    return math.exp(-1*x*x/2)/math.sqrt(2*math.pi)test_num=10000maxy=1/math.exp(0.5)/math.sqrt(2*math.pi)#calculate the top boundrectangle_y=np.random.random_sample(test_num)*(1/math.sqrt(2)-maxy)rectangle_x=np.random.random_sample(test_num)inarea_point_num=0#count how many points in the areafor point_count in range(len(rectangle_x)):    if(rectangle_y[point_count]<f(rectangle_x[point_count])):        inarea_point_num+=1print 'area: '+str(1.0*inarea_point_num/test_num*(1/math.sqrt(2)-maxy)

求出来的结果是0.340433079875，与准确值相差无几