如何绘制 Monte Carlo 圆周率直方图？

Question

我正在尝试通过 Monte Carlo 方法绘制 pi 的直方图分布，但每次我运行模拟而不是近似对称的直方图，峰值在 3.14 左右。输出直方图也有一些差距，我认为我正确地近似了 pi。我的代码如下：

   [...(importing relevant modules)]
    N = 1000 #total number of random points

    circlex = []
    circley = []
    squarex = []
    squarey = []

    pis = []

    for i in range(1, M + 1):
        x1 = random.uniform(-1,1)
        y1 = random.uniform(-1,1)
        if x1**2 + y1**2 <=1:
            circlex.append(x1)
            circley.append(y1)
        else:
            squarex.append(x1)
            squarey.append(y1)

        pi = (4.0*len(circlex))/i
        pis.append(pi)
    print(pi)
    print(pis)
    plt.hist(pis, color = 'g')

输出：

我错过了什么或做错了什么？

Answer 1

您的代码实际上是正确的。但是有两件事你忘了考虑：

绘制直方图可能不是有用的可视化。实际上，bins 很大，很难区分正确的近似值（如 3.14）和错误的近似值（如 3.2）
这个算法需要很长时间才能收敛到圆周率。

作为参考，我用同样的方法得到了那些结果（顺便说一句，你的代码中有一个错字，你应该把range(1, M + 1)转换成range(1, N + 1)）：

approx_pi(N=100)         # 3.2
approx_pi(N=1_000)       # 3.188
approx_pi(N=10_000)      # 3.1372
approx_pi(N=100_000)     # 3.145
approx_pi(N=1_000_000)   # 3.14378
approx_pi(N=10_000_000)  # 3.141584

因此，不要害怕为 N 取更大的值以获得更准确的结果。另外，考虑绘制 pi 近似值的演变而不是直方图来可视化你的近似值是如何收敛的。

最后，根据您的目标，使用 numpy 可以获得更快的代码：

import numpy as np

pis = 4 * np.cumsum(np.linalg.norm(np.random.random(size=(N, 2)), axis=1)<= 1) / np.arange(1, N + 1)

现在进行说明：

我们定义了一个形状为 (N,2) 的随机数组。这对应于 x 和 y 的 N 个样本。请注意，它们是在 0 和 1 之间采样的，但这不会改变估计值。
我们使用带有参数 axis=1 的 numpy.linalg.norm 计算每个 2 坐标向量的范数。
如果采样点在圆内，我们认为布尔数组包含True，否则False
我们对该数组应用累加和。由于在 Python 中，True 在被视为整数时被视为 1，因此累积和包含在索引 i 处仅考虑时圆中的点数i 第一个样本。
我们除以 np.arange(1, N + 1)，索引 i 处包含相应的采样点数。
我们乘以 4 得到圆周率的近似值。

这段代码真的很难看，但速度也更快（比迭代版本快大约 10 倍）。我认为根据您的需要，您可能会对此感兴趣。

如何绘制 Monte Carlo 圆周率直方图？

How to plot the Monte Carlo pi histogram?

python

python-3.x

for-loop

matplotlib

montecarlo