蒙特卡罗模拟（python实现）

心已赠人 2022-10-14 00:44 155阅读 0赞

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

## 蒙特卡罗模拟 ##

蒙特卡罗法是一种以抽样和随机数的产生为基础的随机性方法，因此也称为随机抽样法、计算机随机模拟法等。蒙特卡罗方法的基本原理是通过数字模拟试验，得到所要求解的出现某种事件的概率，作为问题的近似解。很久之前人们就已经开始使用蒙特卡罗方法来解决问题了，早在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()

![20210604103936889.png_pic_center][]

算例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()

![20210604105313979.png_pic_center][]

![在这里插入图片描述][20210604105552160.jpg_pic_center]  
作者：电气-余登武  
![20210604110050545.jpg_pic_center][]

[Link 1]: https://blog.csdn.net/KOBEYU652453/article/details/117459596
[20210604103936889.png_pic_center]: /images/20221014/8d5d834baaa047ccaebb67c2f5dfae6a.png
[20210604105313979.png_pic_center]: /images/20221014/414bb490a17e4c93bfb56c2a44518c94.png
[20210604105552160.jpg_pic_center]: /images/20221014/58f1906285fd4b8192356a37e421b9a9.png
[20210604110050545.jpg_pic_center]: /images/20221014/e30d990b4e6d46e1b60a450a324a56b6.png