先来个测试题:
用仿真的方法计算圆周率π,请写出计算步骤和简单的公式:
(思考十秒钟,想不出来可以继续往下看,后面有答案)
001 什么是蒙特卡洛模拟
蒙特卡洛模拟,又称为统计实验方法
以概率论和统计理论方法为基础的一种计算方法
通过随机数来解决很多计算问题
主要步骤是:
将实际问题转化为概率模型
通过计算机实现统计模拟,以获得问题的近似解
002 基本原理
蒙特卡洛模拟抓住事物运动的几何数量和几何特征,利用数学方法来模拟,是一种数字模拟实验。它是一个以概率模型为基础,按照这个模型所描绘的过程,通过模拟实验的结果,作为问题的近似解。
可以把蒙特卡洛解法分为三步
1.构造或描绘概率过程
2.实现从已知概率分布抽样
3.建立各种估计量
借助计算机技术,蒙特卡洛模拟方法有两大优点
1.简单,省去了繁复的数学推导和验算过程,使普通人能够理解
2.快速,确定了概率模型,后续运算完全用计算机实现
蒙特卡洛模拟的特点:
随机采样得到的近似解,随着随机采样数值增多,得到正确结果的概率越大
003 怎么用
回到最初的测试题,用蒙特卡洛模拟计算圆周率π
圆周率π的计算推导十分复杂,普通人难以明白,但是我们可以用数值模拟的方式计算出圆周率
在这个正方形内部有个内切圆,他们的面积之比是π/4
现在,在这个正方形内部随机生产10000个点(坐标是(x,y)),计算它们与中心点的距离,从而判断是否落在圆内。
园内的点数/正方形内的点数 = π/4
以此来计算得出π的值
直接上代码(python3.6)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
%matplotlib inline
n = 10000
r = 1.0
a,b = (0.0,0.0)
xmin, xmax = a-r, a+r
ymin, ymax = b-r, b+r
x = np.random.uniform(xmin,xmax,n)
y = np.random.uniform(ymin,ymax,n)
fig = plt.figure(figsize=(6,6))
axes = fig.add_subplot(1,1,1)#添加子图
#画子图
plt.plot(x,y,'ko',markersize = 1) #plot绘图 markersize表示点的大小;‘ro’r表示red,o表示圆圈
plt.axis('equal')
d = np.sqrt((x-a)**2 + (y-a)**2)
#res 得到圆中的点数
res = sum(np.where(d<r,1,0)) #numpy.where(conditon,x,y) 满足条件输出x,不满足输出y
pi = res/n*4
print('pi:',pi)
#计算pi的近似值,蒙特卡洛模拟方法,用统计值去近似真实值
#绘制圆形子图
circle = Circle(xy = (a,b), radius = r,alpha = 0.5, color = 'gray')
axes.add_patch(circle)#添加圆形子图
plt.grid(True,linestyle = '--',linewidth = '0.8')
plt.show()
#总结
#蒙特卡洛模拟是用统计值逼近真实值,展示了统计思想
JupyterNotebook运行结果如下
随着我们设置的点的数量n越多,计算的到的pi值就会越接近真实值3.1415926
好了,蒙特卡洛模拟就讲到这里。
这只是一个简单的例子,关于真实场景下运用蒙特卡洛模拟的例子,等我下一篇文章咯~
下一篇预告:用蒙特卡洛模拟人类财富分配