[PYTHON] Sprinkle rice grains to find the circumference (Monte Carlo method)

Postscript

This is an article that tried to implement the Monte Carlo method. I completely forgot.

1.First of all

While reading a children's statistics book [1], I saw an article about sprinkling rice grains and calculating pi. It seemed interesting, so I calculated the pi using the randomly generated x and y coordinate data.

[procedure] It's almost like this

  1. Prepare a square and a circle inscribed in that square.
  2. Plot there randomly (randomly). The point is that this ** crap **!
  3. Find the pi from the number of coordinates inside and outside the circle. image.png [Purpose] Let's see how many coordinates (sample, that is, rice grains) are needed to calculate a value close to the true pi, or a value close to the pi.

2. Implement

It was difficult to actually sow rice grains and experiment, so I tried it on a computer.

2-1. Import modules.

import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.patches as patches

Next, define some lists and constants. This time, we will use a square with one side and a circle with a radius of 0.5.

The purpose of this time is to investigate how many grains of rice should be used to calculate a value close to the true pi, so an example of the number of samples was tentatively defined in n_list.

#Number of coordinates generated
n_list = [10, 20, 50, 70, 100, 500, 1000, 5000, 10000, 100000,1000000]

#Radius of circle
radius = 0.5

#Initialization of the list of pi
pi_list = []

2-2. Generate random coordinates

Defines a function that generates random coordinate data (corresponding to a grain of rice wrapped in data).

def gen_random_coordinate(n):
    #Initialize the generated coordinate list
    x=[]
    y=[]

    #Generate coordinates with random numbers
    x = np.random.rand(n)
    y = np.random.rand(n)
        
    return x,y

2-3. Visualize

Visualize how you actually plot random coordinates.

#Graph drawing function
def make_graph(x,y):
    
    #Create a Figure object and one Axes object belonging to it at the same time
    fig, ax = plt.subplots()
    #Creating a circle
    circle = patches.Circle(xy=(0.5,0.5), radius=0.5, ec='r', fill=False)

    #Adjust to the correct aspect ratio
    ax.set_aspect('equal')
    ax.add_patch(circle)
    ax.scatter(x,y)
    
    plt.axis('scaled')
    plt.grid()
    plt.xlim(0, 1.0)
    plt.ylim(0, 1.0)

    plt.show()

2-4. Calculate pi from coordinate data

Calculate pi using the number of points inside and the number of points outside the circle.

def calc_circular_constant(x,y):
    #Initialization of the number of points inside the circle
    inner=0
    
    #Initialization of the number of points outside the circle
    outer=0
    
    #Count the number of points, whether each point is inside or outside the circle
    for xx , yy in zip(x,y):
        
        #Calculate the distance from the center of the circle
        dist=np.sqrt(np.square(xx-radius)+np.square(yy-radius))
        
        if dist > radius:
            outer +=1
        else:
            inner +=1
            
    print('inner:outer = ', str(inner), ':', str(outer))
    
    #Calculate pi
    pi = inner /(outer+inner) *4
    pi_list.append(pi)
     
    return pi

2-5. Execute

Execute each of the functions defined above. You checked how many samples you needed to get closer to the true pi. So I will try them in order using a list of sample numbers n_list.

for n in n_list:
    x,y = gen_random_coordinate(n)
    make_graph(x,y)
    print(calc_circular_constant(x,y))

3. Experimental results

The following is the output result. The number of samples starts from 10 and the number of points increases as you go down.

image.png Points inside the circle: Points outside the circle = 8: 2 Calculated pi: 3.2

image.png Points inside the circle: Points outside the circle = 16: 4 Calculated pi: 3.2

image.png Points inside the circle: Points outside the circle = 39: 11 Calculated pi: 3.12

image.png Points inside the circle: Points outside the circle = 53: 17 Calculated pi: 3.0285714285714285

image.png Points inside the circle: Points outside the circle = 79: 21 Calculated pi: 3.16

image.png Points inside the circle: Points outside the circle = 399: 101 Calculated pi: 3.192

image.png Points inside the circle: Points outside the circle = 793: 207 Calculated pi: 3.172

image.png Points inside the circle: Points outside the circle = 3894: 1106 Calculated pi: 3.1152

image.png Points inside the circle: Points outside the circle = 7818: 2182 Calculated pi: 3.1272

image.png Points inside the circle: Points outside the circle = 78611: 21389 Calculated pi: 3.14444

image.png Points inside the circle: Points outside the circle = 785814: 214186 Calculated pi: 3.143256

4. Summary

In this experiment, it was confirmed that as the number of samples increased, it approached the true pi. The graph turned deep blue when the number of data exceeded 5000. In the end, it was 3.143256, which was close to the true pi (3.141592653589 ..).

If the number of data is about 100, it is still far from 3.16. (That's right, but I wrote the reason in 5.)

The data to be generated must also be uniformly random data. With the above method, if you generate data that follows a normal distribution, you will not be able to calculate a decent pi. (I wonder if I can correct the outside using the number of data inside the circle, I will try it at another time.)

5. Finally

Even if I write it now, no matter how well it goes with 100, it will be calculated in 0.04 steps, so after 3.12 it will be 3.16. (I didn't really touch on this point because I wanted to do the above experiment.)

If you want to achieve results like the pi you learn in elementary school, you need 400 data. (Because the step of calculated pi is 0.01). If you do it over and over again, the number 3.14 is likely to come out.

Even if you count the rice grains only on the outside of the circle, sprinkling 400 grains of rice and trying many times seems to be a tough experiment. The size of the rice grains themselves is also large, so it seems that you have to prepare paper of the appropriate size.

If you have children, how about doing a free study of the suddenly long "spring break"?

References

[1] What is statistics useful for? (Author) Yoshiyuki Wakui (Edit) Children's Science Editorial Department / Seibundo Shinkosha https://www.seibundo-shinkosha.net/book/kids/20689/

Recommended Posts

Sprinkle rice grains to find the circumference (Monte Carlo method)
#Monte Carlo method to find pi using Python
Introduction to Monte Carlo Method
Find the ratio of the area of Lake Biwa by the Monte Carlo method
Monte Carlo method
Try implementing the Monte Carlo method in Python
The first Markov chain Monte Carlo method by PyStan
Calculation of the shortest path using the Monte Carlo method
Increase the speed of the Monte Carlo method of Cython cut-out implementation.
[Statistics] Visualize and understand the Hamiltonian Monte Carlo method with animation.
Simulate Monte Carlo method in Python
Estimating π by Monte Carlo method