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)

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

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

Selon la formule du cercle, si π = 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

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é.