How to plot the Monte Carlo pi histogram?

Question

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')

Ouput:

What am I missing or doing wrong?

Answer 1

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

1. Plotting an histogram may not be a useful visualization. Indeed, the bins are quite large, and it's hard to distinguish between a correct approximation (like 3.14) and a bad one (like 3.2)
2. This algorithm takes a really long time to converge to pi.

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

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.

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)

Now for the explanation:

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.

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.