蒙特卡洛方法计算圆周率的程序优化增强

发表于 更新于

关于蒙特卡洛方法，我已经连续发布了三篇博文：

在第一篇文章中我们探讨了蒙特卡洛方法使用面积法计算圆周率的经典案例，对MC方法有了一个初步的认识，在第二篇文章中我们在MC方法中使用不同的公式计算圆周率，得到了影响MC方法精度的条件，在第三篇文章中根据第二篇文章得到的条件实现了高精度计算圆周率的一个公式，从而验证了提高MC精度的条件。本文进一步，基于第三篇文章，在程序实现上将写死的程序代码以循环程序自动化。同时，引入mpmath库，实现迭代求根式的精度拓展，但是仍然使用numpy进行随机采样。

Python3 代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#! /usr/bin/env python3
# vim:fenc=utf-8
import numpy as np
from mpmath import mp, mpf
# 修改为 40 位十进制精度
mp.dps = 40

# 定义函数f(x)
def f(x):
    return 1/(1-x**2)**0.5

# Monte Carlo方法求积分
def monte_carlo_integral(f, m, b, num_samples=1000):
    # 根据半角公式用m设定嵌套次数
    count = 0
    while count < m:
        b = mp.sqrt(1/2-mp.sqrt(1-b**2)/2)
        count +=1
    # 先用半角公式根据b计算a
    a = mp.sqrt(1/2-mp.sqrt(1-b**2)/2)
    # 在区间[a, b]内随机生成num_samples个样本点
    samples = np.random.uniform(a, b, num_samples)
    
    # 计算这些样本点对应的函数值
    y_values = f(samples)
    
    # 根据平均高度估计积分值，乘以区间的宽度得到近似面积
    integral_estimate = 2**(m + 3)*(b - a) * np.mean(y_values)
    
    return integral_estimate

# 使用Monte Carlo方法估计积分
estimated_integral = monte_carlo_integral(f,42, mp.sqrt(2)/2)

print("Pi", estimated_integral)

计算结果

  • 运行程序得到的结果： 3.14159265358979313874736836854190094502
  • 资源消耗情况： 0.08s user 0.02s system 98% cpu 0.099 total
  • 按现有资料\(\pi\)的前20位精确结果为3.1415926535 8979323846, 对比可知我们的程序可以精确到小数点后15位, 用时仅0.08s, 所以这是一个不错的结果。

优化的细节

  • 动态精度控制：通过循环，使用参数m灵活调整精度(如m=42时，区间极小，方差极低）。
  • 高效性：样本需求低(1000), 计算速度快。
  • 可拓展性：增加m可进一步提高精度，代码结构清晰，易于维护。

在这四篇博客中，我们从初步了解蒙特卡罗（MC）方法开始，逐步进行优化改进，展示了三种能够计算圆周率的程序。其中，本文介绍的第三个程序以其高效性和灵活性脱颖而出，特别适合需要高精度结果并希望简化手动计算的场景。然而，值得注意的是，m的取值并非越大越好，因为随着m增大，随机点的取值区间(a,b)会减少，加之本程序使用了通常采用双精度浮点数的np.random.uniform 进行随机取样，为了确保精度，m需要保持在一个适当的范围内。此外，np.mean的使用也限制了精度的进一步提升。如果想要进一步提高利用蒙特卡罗方法计算圆周率的精度，可以从调整m的取值和改进均值计算方式这两个方面入手。尽管如此，蒙特卡罗方法的核心价值在于在一定精度下快速提供那些常规方法难以计算或计算量极大的问题的结果，因此本文在这些考量基础上未对精度做进一步优化，但指明了可能的改进路径，以期与读者共同探讨技术细节。