如何绘制蒙特卡罗 pi 直方图? [英] How to plot the Monte Carlo pi histogram?

问题描述

我试图从蒙特卡罗方法绘制 pi 的直方图分布，但每次运行模拟时我得到的直方图要么向左倾斜，要么向右倾斜，而不是近似对称的直方图，峰值约为3.14.输出直方图也有一些差距，我认为我正确地逼近了 pi.我的代码如下:

I am trying to plot a histogram distribution of pi from the Monte Carlo method but I am getting histograms that are either skewed to the left or right each time I run the simulation instead of a histogram that is approximately symmetric and peaks at around 3.14. The output histograms also have some gaps and I think I am approximating pi correctly. My code is below:

   [...(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')

输出:

我遗漏了什么或做错了什么?

What am I missing or doing wrong?

推荐答案

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

Your code is actually correct. But there are two things you forgot to take into account:

1. 绘制直方图可能不是一种有用的可视化方法.实际上，bins 非常大，很难区分正确的近似值(如 3.14)和错误的近似值(如 3.2)
2. 这个算法需要很长时间才能收敛到 pi.

作为参考，我使用了相同的方法并得到了那些结果(顺便说一下，您的代码中有一个错字，您应该将 range(1, M + 1) 转换为 范围(1, N + 1)):

As a reference, I used the same method and get those results (by the way, there's a typo in your code, you should transform range(1, M + 1) into 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 近似值的演变而不是直方图，以可视化您的近似值是如何收敛的.

Hence, don't be afraid to take larger values for N to get more accurate results. Plus, consider plotting the evolution of the approximation of pi rather than an histogram to visualize how your approximation is converging.

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

Finally, depending on what your goal is, it is possible to get a faster code using 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)

现在解释:

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

1. We define a random array of shape (N,2). This corresponds to the N samples of x and y. Note that they are sampled here between 0 and 1, but this does not change the estimation.
2. We compute the norm of each of these 2-coordinates vectors using numpy.linalg.norm with argument axis=1.
3. We consider the boolean array containing True if the sampled point lies within the circle, and False otherwise
4. We apply a cumulative sum on this array. Since in Python, True is considered as 1 when considered as an integer, the cumulative sum contains at index i the number of points that are in the circle when considering only the i first samples.
5. We divide by np.arange(1, N + 1) which contains at index i the corresponding number of sampled points.
6. We multiply by 4 to get an approximation of pi.

这段代码真的很丑，但也快得多(比迭代版本快 10 倍左右).我认为根据您的需要，您可能会对此感兴趣.

This code is really ugly, but also way faster (around 10 times faster than the iterative version). I thought that depending on your needs, you may find this of interest.

这篇关于如何绘制蒙特卡罗 pi 直方图?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！