蒙特卡洛算法

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

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

直接看程序:

import numpy as np
 
radius=1
test_num=100000
#random generate 10000 points in a sqare width is 2,and the center is in origin
sqare_x=2*np.random.random_sample(test_num)-1
sqare_y=2*np.random.random_sample(test_num)-1
incircle_point_num=0
#count how many points in the incircle
for 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+=1
print 'PI: '+str(4.0*incircle_point_num/test_num)

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

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

测试个数     算出的π的值
10    4.0
100    3.2
1000    3.192
10000    3.1388
100000    3.13728
1000000    3.142496
10000000    3.14147268
      
      
 

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

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

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

例如：

求

 ∫10e−x2/2/√2πdx∫01e−x2/2/2πdx

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

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

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

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

import numpy as np
import math
#calculate the function
def f(x):
    return math.exp(-1*x*x/2)/math.sqrt(2*math.pi)
 
test_num=10000
maxy=1/math.exp(0.5)/math.sqrt(2*math.pi)#calculate the top bound
rectangle_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 area
for point_count in range(len(rectangle_x)):
    if(rectangle_y[point_count]<f(rectangle_x[point_count])):
        inarea_point_num+=1
print 'area: '+str(1.0*inarea_point_num/test_num*(1/math.sqrt(2)-maxy) )

求出来的结果是0.340433079875，与准确值相差无几
————————————————
版权声明：本文为CSDN博主「蒙特卡洛家的树」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/u011021773/article/details/78301557