[PYTHON] Estimating π by Monte Carlo method

Estimate π by the Monte Carlo method.

Method

Random numbers with uniform x and y coordinates are plotted as points in the square. By counting the number of points in the circle, the probability of the points in the circle can be calculated.

The area can be obtained from this probability.

Specifically, we assume a square whose x and y coordinates are (-1, -1) and (1,1) and a circle that is in close contact with it.

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

fig,ax = plt.subplots(figsize=(5, 5))
rect = patches.Rectangle((-1,-1),2,2,linewidth=1, fill=False, color='blue')
ax.add_patch(rect)
ax.add_artist(plt.Circle((0, 0), 1, fill=False, color='r'))

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)

Count the number of points that fall inside this red circle.

Try to hit 1000 dots

df = pd.DataFrame(
    (1 + 1) * np.random.rand(1000, 2) - 1,
    columns=['x', 'y'])

fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(df.x, df.y, s=1, color='black')
rect = patches.Rectangle((-1,-1),2,2,linewidth=1, fill=False, color='blue')
ax.add_patch(rect)
circle = plt.Circle((0, 0), 1, fill=False, color='r')
ax.add_artist(circle)

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)

Calculate how many of these point clouds are in the circle.

df["is_inside_circle"] = df.x**2 + df.y**2 <= 1
len(df[df.is_inside_circle]) # 798

How to count

I'm counting the number of points that become.

As a result, 798 out of 1000 are in the circle.

Considering that 79.8% of the squares are made up of points, the area of the circle is

In the circle formula, if π = 3.14

The value is closer.

The way 1000 points are arranged comes from the way numpy random numbers are generated. Therefore, there is a bias that happens to have many dots in the circle.

Therefore, by increasing the number of trials (number of simulations), the probability distribution of the calculated π value can be searched.

Now, let's simulate this 1000 times.

Hit 1000 points in the rectangle and try to calculate π 1000 times

results = [] #Put the estimated value of π in this variable
for i in range(1000):
    rand_num = 1000
    xy = (1 + 1) * np.random.rand(rand_num, 2) - 1
    df = pd.DataFrame(xy, columns=['x', 'y'])
    df["is_inside_circle"] = df.x**2 + df.y**2 <= 1
    ans = 4*len(df[df.is_inside_circle == True])/rand_num
    results.append(ans)

#Formatting for drawing
r_df = pd.DataFrame(results, columns=['area'])
plt.hist(r_df.area, bins=50)

print(f"mean: {np.mean(results)}") #average: 3.1412359999999997
print(f"variance: {np.var(results)}") #Distributed: 0.0026137523040000036

By trying 1000 times

It was found that "in the Monte Carlo method using 1000 points, π can take values in the range of 3.1412 on average and 0.0026 in variance."

The standard deviation can be found by taking the square root of the variance.

np.sqrt(np.var(results)) # 0.051124869721105436

Assuming that the probability distribution of the estimated value of π follows a normal distribution, the following can be found.

68% of the values are within ± 1σ from the average, 95% of the values are within ± 2σ from the average, and 99.7% are within ± 3σ from the average. (Σ = standard deviation = 0.051 here)

Let's actually calculate.

Number of points in the range ± 1σ from the average

result_mean = np.mean(results)
sigma = np.sqrt(np.var(results))

#Number of points in the range ± 1σ from the average
len([result for result in results if result >result_mean + (-1) * sigma and result < result_mean + 1 * sigma])
# 686

#Number of points in the range ± 2σ from the average
len([result for result in results if result >result_mean + (-2) * sigma and result < result_mean + 2 * sigma])
# 958

#Number of points in the range ± 3σ from the average
len([result for result in results if result >result_mean + (-3) * sigma and result < result_mean + 3 * sigma])
# 998

From this, it can be seen that a number close to the normal distribution can be obtained.

Summary

In this article, I said, "I'll hit 1000 points in a rectangle and try to calculate π 1000 times." It has two variables: the number of points to hit in the rectangle and the number of simulations of π.

By changing this variable

--Increase the number of points to hit in the rectangle => ** The value of π per trial becomes close to the true value **. --Increase the number of simulations of π => ** The mean of the estimated values of π becomes closer to the true value, and the variance becomes smaller **.

That is realized.