如何绘制蒙特卡洛 pi 直方图？答案

【问题标题】：How to plot the Monte Carlo pi histogram?如何绘制蒙特卡洛 pi 直方图？
【发布时间】：2020-08-27 16:26:24
【问题描述】：

我正在尝试从 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')

输出：

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

【问题讨论】：

标签： python python-3.x for-loop matplotlib montecarlo

【解决方案1】：

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

绘制直方图可能不是一种有用的可视化方法。实际上，这些 bin 非常大，很难区分正确的近似值（如 3.14）和错误的近似值（如 3.2）
此算法需要很长时间才能收敛到 pi。

作为参考，我使用了相同的方法并得到了这些结果（顺便说一下，您的代码中有一个错字，您应该将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 之间进行采样，但这不会改变估计值。
我们使用 numpy.linalg.norm 和参数 axis=1 计算每个 2 坐标向量的范数。
如果采样点在圆内，我们考虑包含 True 的布尔数组，否则考虑 False
我们在这个数组上应用一个累积和。由于在 Python 中，True 在被视为整数时被视为1，因此累积总和在索引i 处包含仅考虑i 第一个样本时圆圈中的点数。
我们除以np.arange(1, N + 1)，它在索引i处包含相应的采样点数。
我们乘以 4 得到 pi 的近似值。

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

【讨论】：