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)
スクリーンショット 2020-04-26 16.38.04.png

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

x^{2} + y^{2} \leq 1

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

2 \times 2 \times 0.798 = 3.192

In the circle formula, if π = 3.14

1 \times 1 \times 3.14 = 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
スクリーンショット 2020-04-26 16.48.24.png

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.

Recommended Posts

Estimating π by Monte Carlo method
Monte Carlo method
Dominion compression play analyzed by Monte Carlo method
Introduction to Monte Carlo Method
The first Markov chain Monte Carlo method by PyStan
Simulate Monte Carlo method in Python
Find the ratio of the area of Lake Biwa by the Monte Carlo method
Gomoku AI by Monte Carlo tree search
#Monte Carlo method to find pi using Python
Try implementing the Monte Carlo method in Python
[Statistics] Let's explain sampling by Markov chain Monte Carlo method (MCMC) with animation.
Calculation of the shortest path using the Monte Carlo method
Sprinkle rice grains to find the circumference (Monte Carlo method)
[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy
Finding pi with a three-line function [Python / Monte Carlo method]
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method
Classify data by k-means method
[Statistics] Visualize and understand the Hamiltonian Monte Carlo method with animation.
I couldn't install pypy3.6-7.3.1 on macOS + pyenv, but I could install pypy3.6-7.3.0. I felt the wind of pypy by the Monte Carlo method.