蒙特卡洛法模拟计算圆周率π

一、蒙特卡洛法介绍
蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是一种以概率统计理论为基础的数值计算方法，常用于特定条件下的概率计算问题。蒙特卡罗是摩纳哥的著名赌城，该法为表明其随机抽样的本质而命名。
算法思路简单也好理解：比如抛一枚硬币，假设我们开始不知道正面朝上的概率是多少，却有大量的时间来将硬币抛一万次，那么在一万次试验后，会发现正面朝上的次数接近一半，当然，抛的次数越多，概率越接近50%，蒙特卡洛方法是大数定律在实际应用问题上的体现。其优点十分明显，基本可以绕开问题本身的“黑盒”，不必考虑问题内部的结构而只关注问题的输入与输出，利用输出的结果来分析问题，适用于对离散系统进行计算仿真试验。
例如上例中，我们不用考虑硬币在空中停留多长时间，不用考虑抛出力度、硬币大小、空气阻力、风速等乱七八糟的问题，在大量的试验后只关注最后硬币哪面朝上，就能正确估算出硬币正面朝上的概率。同样的，例如我们不知道走到某个路口需要等红绿灯的概率，不知道某个产品线的合格率，蒙特卡洛法告诉你：模拟一万次试验后你就知道了O(∩_∩)O~
二、利用蒙特卡洛方法计算圆周率π
采用蒙特卡洛思想，首先在一个正方形区域内随机生成若干个均匀分布的点，随后判断哪些点在正方形的内切圆范围内。如果点的数量足够多，那么圆内点的数量与点的总数量的比值，就是圆的面积与正方形面积之比。利用点数量的比值与正方形面积就可以推出圆的面积，进而得出圆周率π。

上代码：

import numpy as np
import matplotlib.pyplot as plt
import math
from tqdm import tqdm
from numpy import mean

center_x = 0.0  #圆心横坐标
center_y = 0.0  #圆心纵坐标
r = 1.0         #半径

class Caculate_Pi:
    def __init__(self,dot_num):
        self.dot_num = dot_num   #初始化点的数量
        self.count = 0   #圆内点的计数器
        self.x = np.random.uniform(center_x - r,center_x + r,self.dot_num)   #点均匀分布的横轴范围
        self.y = np.random.uniform(center_y - r,center_y + r,self.dot_num)   #点均匀分布的纵轴范围
    
    def Judge(self,x,y,r):      #判断点是否在圆内
        return np.sqrt(x ** 2 + y ** 2) <= r
    
    def Count(self):            #计算圆内点的个数
        result = [self.Judge(self.x[i], self.y[i], r) for i in range(self.dot_num)]
        for i in result:
            if i:
                self.count += 1

dot_num = []
Pi = []
for i in tqdm(range(1,1000)):   #考虑从100到1000000个点模拟的π值如何变化
    dot_num.append(100 * i)
    demo = Caculate_Pi(100 *i)
    demo.Count()
    pi = demo.count * 4 /(100 * i)
    Pi.append(pi)
#可视化随着点数量增多模拟π的变化
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.xlabel('模拟次数')
plt.ylabel('模拟的Pi值')
plt.axhline(y = math.pi,c = "red")
plt.annotate('y=π',xy=(0,3.1415926),xytext=(-8000,3.14))
plt.plot(dot_num,Pi)
print('模拟得到的π均值:',mean(Pi[100:]))

运行结果：

可以看出模拟次数越多(从100到100000)，模拟的π值越来越趋近于真实的π（3.141592654）。本次计算最终模拟的平均值为

模拟得到的π均值:3.1414679762608873

已经比较接近了，以上就是利用蒙特卡洛法计算圆周率的小练习。
当然，计算的方法还包括经典的精确算法——割圆法、无穷级数法，其同样可以用代码实现，下次无聊的时候再写吧。