蒙特卡罗(Monte Carlo)方法:简单来说,蒙特卡洛的基本原理简单描述是先大量模拟,然后计算一个事件发生的次数,再通过这个发生次数除以总模拟次数,得到想要的结果,精髓就是:用统计结果去计算频率,从而得到真实值的近似值。蒙特卡洛方法可以应用在很多场合,但求的是近似解,在模拟样本数越大的情况下,越接近与真实值,但样本数增加会带来计算量的大幅上升。
1.求圆周率pi的近似值:
(1)正方形内部有一个相切的圆,它们的面积之比是π/4。现在,在这个正方形内部,随机产生10000个点,计算它们与中心点的距离,从而判断是否落在圆的内部,若这些点均匀分布,则圆周率 。
其中res:表示落到圆内投点数 n:表示总的投点数
(2)代码
import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.patches import Circle# 投点次数n = 10000# 圆的半径、圆心r = 1.0a,b = (0.,0.)# 正方形区域x_min, x_max = a-r, a+ry_min, y_max = b-r, b+r# 在正方形区域内随机投点x = np.random.uniform(x_min, x_max, n) #均匀分布y = np.random.uniform(y_min, y_max, n)#计算点到圆心的距离d = np.sqrt((x-a)**2 + (y-b)**2)#统计落在圆内点的数目res = sum(np.where(d < r, 1, 0))#计算pi的近似值(Monte Carlo:用统计值去近似真实值)pi = 4 * res / nprint('pi: ',pi)#画个图fig = plt.figure()axes = fig.add_subplot(111)axes.plot(x, y,'ro',markersize = 1)plt.axis('equal') #防止图像变形circle = Circle(xy=(a,b), radius=r ,alpha=0.5)axes.add_patch(circle)plt.show()
(3)运行结果:
n = 30000#矩形边界区域x_min, x_max = 0.0, 1.0y_min, y_max = 0.0, 1.0#在矩形区域内随机投点xx = np.random.uniform(x_min, x_max, n)#均匀分布y = np.random.uniform(y_min, y_max, n)#统计落在函数y=x^2下方的点res = sum(np.where(y < f(x), 1 ,0))#计算定积分的近似值integral = res / nprint('integeral: ', integral)# 画个图fig = plt.figure()axes = fig.add_subplot(111)axes.plot(x, y,'ro',markersize = 1)plt.axis('equal') # 防止图像变形axes.plot(np.linspace(x_min, x_max, 10), f(np.linspace(x_min, x_max, 10)), 'b-') # 函数图像#plt.xlim(x_min, x_max)plt.show()
(3)运行结果
联系客服