基于python实现Monte Carlo

一、背景介绍

蒙特卡洛模型（Monte Carlo Model），又称统计模拟法或随机抽样技术，是一种基于概率和统计理论的随机模拟方法。这种方法通过将所求解的问题与一定的概率模型相联系，使用电子计算机实现统计模拟或抽样，以获得问题的近似解。蒙特卡洛模型的基本思想主要包括以下几点：

• 概率模型构建：为了求解数学、物理、工程技术以及管理等方面的问题，首先建立一个概率模型或随机过程。这个模型或过程的参数，如概率分布或数学期望，与问题的解相关联。
• 随机抽样：通过对模型或过程的观察或抽样试验，收集大量的随机样本数据。这些样本数据用于计算所求参数的统计特征。
• 近似解求解：利用算术平均值或其他统计量作为所求解的近似值。通过增加样本数量，可以提高解的精度。

蒙特卡洛模型的应用领域十分广泛，包括但不限于：

• 物理领域：在计算物理、物理化学和相关应用领域，蒙特卡洛方法被用于处理从复杂的量子色动力学计算到设计隔热罩和空气动力学形式等各种问题。
• 工程领域：在工程设计中的敏感性分析和定量概率分析方面，蒙特卡罗方法发挥着重要作用。例如，在微电子工程中，它被用于分析模拟和数字集成电路中的相关和非相关变化。
• 计算生物学：在计算生物学的各个领域，如系统演化的贝叶斯推断和基因组、蛋白质或膜等生物系统的研究中，蒙特卡洛方法也发挥着重要作用。
• 虚拟现实技术：在虚拟现实技术中，蒙特卡洛方法被用于光线追踪、碰撞检测和随机走样等方面，以提高渲染效果和计算效率。

蒙特卡洛模型的发展历史悠久，早在1777年，法国Buffon就提出了用投针实验的方法求圆周率π，这被认为是蒙特卡罗方法的起源。到了20世纪40年代，S.M.乌拉姆和J.冯·诺伊曼为研制核武器而首先提出了蒙特卡罗模拟的概念。随着电子计算机的出现和计算能力的提升，蒙特卡罗方法在各个领域的应用得到了广泛的拓展和深入。

二、利用Monte Carlo计算$\pi$

在一个边长为2的正方形区域内随机投点，并检查每个点是否落在单位圆内（即，到原点的距离是否小于或等于1）。最后，我们根据落在圆内的点数与总点数的比例来估计π的值。
请注意，由于蒙特卡洛模拟是基于随机抽样的，因此每次运行程序得到的π的近似值可能会有所不同。但是，随着投点数量的增加，估计值的精度也会逐渐提高。

具体实施步骤：
（1）投点数量$counts$、计数器$count$；
（2）在$[-2,2]$区间内，随机生成一个坐标($x$, $y$)；
（3）如果满足条件$count$ = $count$ + 1；
（4）圆的面积：正方形面积 = $\pi :4=count:counts$，模拟结果$\pi$ = $4count/counts$。

三、Monte Carlo优缺点

1.优点：

• 预测精确；
• 连续和非连续都能处理。

2.缺点：

• 计确时间通常较长；
• 误差是概率误差。

四、以计算半径为1的圆面积为例

完整代码如下

import random
import math
import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
counts = 10000#循环次数

#作一个美丽的圆、正方形
def plot_ax(counts):
    plt.xlim(xmax=1.1,xmin=-1.1)
    plt.ylim(ymax=1.1,ymin=-1.1)
    x1=np.linspace(-1,1,counts)#-1到1中画counts个点
    y1=[-1 for i in range(counts)]
    plt.plot(x1,y1,color='green')
    x2=np.linspace(-1,1,counts)#-1到1中画counts个点
    y2=[1 for i in range(counts)]
    plt.plot(x2,y2,color='green')
    x3=[1 for i in range(counts)]
    y3=np.linspace(-1,1,counts)#-1到1中画counts个点
    plt.plot(x3,y3,color='green')
    x4=[-1 for i in range(counts)]
    y4=np.linspace(-1,1,counts)#-1到1中画counts个点
    plt.plot(x4,y4,color='green')
    x1=np.linspace(-1,1,counts)#-1到1中画counts个点
    y1 = (1-x1**2)**0.5
    plt.plot(x1,y1,color='red')
    x1=np.linspace(-1,1,counts)#-1到1中画counts个点
    y1 = -(1-x1**2)**0.5
    plt.plot(x1,y1,color='red')
    x1=np.linspace(-1,1,counts)#-1到1中画counts个点
    y1 = [0 for i in range(counts)]
    plt.plot(x1,y1,color='black')
    x1=[0 for i in range(counts)]
    y1 = np.linspace(-1,1,counts)#-1到1中画counts个点
    plt.plot(x1,y1,color='black')
    plt.tick_params(axis='x',colors='black')
    plt.tick_params(axis='y',colors='black')
    ax.set_aspect('equal', adjustable='box')
def π(counts):
    count = 0  # 计数器
    for i in range(counts):
        #生成随机坐标(x,y)
        x=random.uniform(-1,1)#生成0-1之间的随机数
        #print(x)
        y=random.uniform(-1,1)#生成0-1之间的随机数
        #print(y)
        distance=x**2+y**2
        #print(y_predict)
        if distance<=1:#坐标点distance小于等于1，在圆内，即才符合条件，计数器加1
            count = count+1
            plt.plot(x,y,'y.')
        else:
            plt.plot(x, y, 'b.')
    ax.set_aspect('equal', adjustable='box')
    plt.show()
    predict_value=4*count/counts
    return predict_value
plot_ax(counts)
print(π(counts))