蒙特卡罗模拟（python实现）

      博主专业是电气，电气科学实验中，有时用到蒙特卡罗模拟，今天讲解下。
      先附上个人电气博文目录传送门。见下。点击即可
学好电气全靠它，个人电气博文目录（持续更新中…）

蒙特卡罗模拟

      蒙特卡罗法是一种以抽样和随机数的产生为基础的随机性方法，因此也称为随机抽样法、计算机随机模拟法等。蒙特卡罗方法的基本原理是通过数字模拟试验，得到所要求解的出现某种事件的概率，作为问题的近似解。很久之前人们就已经开始使用蒙特卡罗方法来解决问题了，早在17世纪，发生事件的概率是用发生事件的频率来决定的。在20世纪计算机的使用，使得这样的实验可以大量快速的进行模拟。近年来，随着蒙特卡罗方法不仅被用于数学、物理、化学等领域，还应用于核科学与技术、机械工程、天文学和仪器科学与技术等。该法产生了不同的分支，有直接蒙特卡罗法，动力学蒙特卡罗法等 。

      本质就是通过模拟事件发生，然后通过计算某一特定事件发生的频率来计算得到某值（=概率）

      本文通过两个简单算例感受下蒙特卡罗模拟。

算例1：计算π

import numpy as np
import matplotlib.pyplot as plt
from pylab import *
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
from matplotlib.patches import Circle#导入圆
n = 100000
# 投点次数
r = 1.0           # 半径
a, b = (0., 0.)   # 圆心
# 圆的信息
x_min, x_max = a-r, a+r
y_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)) #np.where(d < r, 1, 0) 如果d小于r,返回1.否则返回0# 统计落在圆内的点的数目


pi = 4 * res / n #res / n 得到是 =圆面积/正方形面积  。4是本例中正方形面积。 
#公式 π*r*r=正方形面积（4）*（落在圆内的个数/落在正方形中的个数）
print('蒙特卡罗模拟得到的π值: ', pi) #得到的π值


#画图
fig = plt.figure(figsize = (6,6))
axes = fig.add_subplot(1,1,1)
plt.scatter(x,y,marker='*',linewidth = 0.1)
plt.axis('equal')
circle = Circle(xy = (a,b),radius = r, alpha = 0.5 ,color = 'red')
axes.add_patch(circle)
plt.title('蒙特卡罗模拟得到的π值:%f ' %pi)
plt.show()

算例2：计算积分 y = sin(x)

import numpy as np
import matplotlib.pyplot as plt
from pylab import *
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False

n = 10000#投点次数

#矩形边界
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0

x = np.random.uniform(x_min, x_max, n) # 均匀分布
y = np.random.uniform(y_min, y_max, n)

def f(x):
    return np.sin(x)
# 创建函数 y = sin(x)

res = sum(np.where(y < f(x), 1, 0)) #统计点在sin(x)下面的个数

integral = res / n
print('蒙特卡罗模拟得到的sin(x)积分值:%f: ', integral)

plt.scatter(x,y)#散点图

xi = np.linspace(0,1,100)
yi = np.sin(xi)
plt.plot(xi,yi,'--k')
plt.fill_between(xi, yi, 0, color ='red',alpha=0.5,label='area')
plt.grid()
plt.title('蒙特卡罗模拟得到的sin(x)积分值:%f: '%integral)
plt.show()

作者：电气-余登武