Estimation de π par la méthode de Monte Carlo

Estimer π par la méthode de Monte Carlo.

Méthode

Tracez un nombre aléatoire avec des coordonnées x et y uniformes comme points dans le carré. En comptant le nombre de points dans le cercle, la probabilité des points dans le cercle peut être calculée.

La superficie peut être obtenue à partir de cette probabilité.

Plus précisément, nous supposons un carré dont les coordonnées x et y sont (-1, -1) et (1, 1) et un cercle qui est en contact étroit avec lui.

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)

Comptez le nombre de points à l'intérieur de ce cercle rouge.

Essayez de toucher 1000 points

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

Calculez le nombre de ces points dans le cercle.

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

Comment compter

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

Je compte le nombre de points qui deviennent.

En conséquence, 798 sur 1000 sont dans le cercle.

Considérant que 79,8% des carrés sont constitués de points, l'aire du cercle est

2 \times 2 \times 0.798 = 3.192

Selon la formule du cercle, si π = 3,14

1 \times 1 \times 3.14 = 3.14

La valeur est plus proche.

La méthode d'agencement de 1000 points est dérivée de la méthode de génération de nombres aléatoires numpy. Par conséquent, il se produit un biais dans lequel de nombreux points sont placés dans le cercle.

Par conséquent, en augmentant le nombre d'essais (nombre de simulations), la distribution de probabilité de la valeur π calculée peut être recherchée.

Simulons cela 1000 fois.

Frappez 1000 points dans le carré et essayez de calculer π 1000 fois

results = [] #Mettez la valeur estimée de π dans cette 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)

#Formatage pour le dessin
r_df = pd.DataFrame(results, columns=['area'])
plt.hist(r_df.area, bins=50)

print(f"mean: {np.mean(results)}") #moyenne: 3.1412359999999997
print(f"variance: {np.var(results)}") #Distribué: 0.0026137523040000036
スクリーンショット 2020-04-26 16.48.24.png

En essayant 1000 fois

"Dans la méthode de Monte Carlo utilisant 1000 points, π peut prendre des valeurs de l'ordre de 3,1412 en moyenne et 0,0026 en dispersion."

L'écart type peut être trouvé en prenant la racine carrée de la variance.

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

En supposant que la distribution de probabilité de la valeur estimée de π suit une distribution normale, on peut trouver ce qui suit.

68% des valeurs sont à ± 1σ de la moyenne, 95% des valeurs sont à ± 2σ de la moyenne et 99,7% sont à ± 3σ de la moyenne. (Σ = écart type = 0,051 ici)

Calculons réellement.

Nombre de points dans la plage ± 1σ de la moyenne

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

#Nombre de points dans la plage ± 1σ de la moyenne
len([result for result in results if result >result_mean + (-1) * sigma and result < result_mean + 1 * sigma])
# 686

#Nombre de points dans la plage ± 2σ de la moyenne
len([result for result in results if result >result_mean + (-2) * sigma and result < result_mean + 2 * sigma])
# 958

#Nombre de points dans la plage ± 3σ de la moyenne
len([result for result in results if result >result_mean + (-3) * sigma and result < result_mean + 3 * sigma])
# 998

À partir de là, on peut voir qu'un nombre proche de la distribution normale peut être obtenu.

Résumé

Dans cet article, j'ai dit: "Je vais toucher 1000 points dans un carré et essayer de calculer π 1000 fois". Il a deux variables: le nombre de points à atteindre dans le quadrilatère et le nombre de simulations de π.

En modifiant cette variable

--Augmentez le nombre de points à atteindre dans le carré => ** La valeur de π par essai devient proche de la vraie valeur **. --Augmenter le nombre de simulations de π => ** La moyenne des valeurs estimées de π se rapproche de la valeur réelle et la variance devient plus petite **.

Cela est réalisé.

Recommended Posts

Estimation de π par la méthode de Monte Carlo
Méthode de Monte Carlo
Jeu de compression Dominion analysé par la méthode de Monte Carlo
Introduction à la méthode Monte Carlo
La première méthode de Monte Carlo en chaîne de Markov par PyStan
Simuler la méthode Monte Carlo en Python
Trouvez le ratio de la superficie du lac Biwa par la méthode de Monte Carlo
IA à cinq yeux en recherchant les arbres de Monte Carlo
Méthode #Monte Carlo pour trouver le rapport de circonférence en utilisant Python
Essayez d'implémenter la méthode Monte Carlo en Python
[Statistiques] Je vais expliquer l'échantillonnage par la méthode de Monte Carlo en chaîne de Markov (MCMC) avec animation.
Calcul de l'itinéraire le plus court selon la méthode de Monte Carlo
Saupoudrer de grains de riz pour trouver le rapport de circonférence (méthode de Monte Carlo)
[Calcul scientifique / technique par Python] Intégration Monte Carlo, calcul numérique, numpy
Recherche du rapport de circonférence avec une fonction à 3 lignes [méthode Python / Monte Carlo]
[Calcul scientifique / technique par Python] Simulation de Monte Carlo par la méthode metropolis de la thermodynamique du système de spin ascendant 2D
Classer les données par la méthode k-means
[Statistiques] Visualisez et comprenez la méthode Hamiltonian Monte Carlo avec animation.
Je n'ai pas pu installer pypy3.6-7.3.1 avec macOS + pyenv, mais je pourrais installer pypy3.6-7.3.0. J'ai senti le vent du pypy par la méthode Monte Carlo.