[PYTHON] Schätzung von π nach der Monte-Carlo-Methode

Schätzen Sie π nach der Monte-Carlo-Methode.

Methode

Zeichnen Sie eine Zufallszahl mit einheitlichen x- und y-Koordinaten als Punkte im Quadrat. Durch Zählen der Anzahl der Punkte im Kreis kann die Wahrscheinlichkeit von Punkten im Kreis berechnet werden.

Die Fläche kann aus dieser Wahrscheinlichkeit erhalten werden.

Nehmen Sie insbesondere ein Quadrat an, dessen x- und y-Koordinaten (-1, -1) und (1, 1) sind, und einen Kreis, der in engem Kontakt damit steht.

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)

Zählen Sie die Anzahl der Punkte innerhalb dieses roten Kreises.

Versuche 1000 Punkte zu treffen

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

Berechnen Sie, wie viele dieser Punkte sich im Kreis befinden.

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

Wie man zählt

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

Ich zähle die Anzahl der Punkte, die werden.

Infolgedessen befinden sich 798 von 1000 im Kreis.

Wenn man bedenkt, dass 79,8% der Quadrate aus Punkten bestehen, ist die Fläche des Kreises

2 \times 2 \times 0.798 = 3.192

Nach der Kreisformel ist π = 3,14

1 \times 1 \times 3.14 = 3.14

Der Wert ist näher.

Das Verfahren zum Anordnen von 1000 Punkten leitet sich vom Verfahren zum Erzeugen von numpy Zufallszahlen ab. Daher tritt eine Vorspannung auf, bei der zufällig viele Punkte im Kreis platziert werden.

Daher kann durch Erhöhen der Anzahl von Versuchen (Anzahl von Simulationen) die Wahrscheinlichkeitsverteilung des berechneten π-Werts gesucht werden.

Lassen Sie uns dies 1000 Mal simulieren.

Treffen Sie 1000 Punkte auf dem Quadrat und versuchen Sie, π 1000 Mal zu berechnen

results = [] #Geben Sie den geschätzten Wert von π in diese Variable ein
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)

#Formatierung zum Zeichnen
r_df = pd.DataFrame(results, columns=['area'])
plt.hist(r_df.area, bins=50)

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

Indem Sie es 1000 Mal versuchen

"Bei der Monte-Carlo-Methode mit 1000 Punkten kann π Werte im Bereich von durchschnittlich 3,1412 und einer Dispersion von 0,0026 annehmen."

Die Standardabweichung kann ermittelt werden, indem die Quadratwurzel der Varianz gezogen wird.

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

Unter der Annahme, dass die Wahrscheinlichkeitsverteilung des geschätzten Wertes von π einer Normalverteilung folgt, kann Folgendes gefunden werden.

68% der Werte liegen innerhalb von ± 1σ vom Durchschnitt, 95% der Werte liegen innerhalb von ± 2σ vom Durchschnitt und 99,7% liegen innerhalb von ± 3σ vom Durchschnitt. (Σ = Standardabweichung = 0,051 hier)

Lassen Sie uns tatsächlich berechnen.

Anzahl der Punkte im Bereich ± 1σ vom Durchschnitt

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

#Anzahl der Punkte im Bereich ± 1σ vom Durchschnitt
len([result for result in results if result >result_mean + (-1) * sigma and result < result_mean + 1 * sigma])
# 686

#Anzahl der Punkte im Bereich ± 2σ vom Durchschnitt
len([result for result in results if result >result_mean + (-2) * sigma and result < result_mean + 2 * sigma])
# 958

#Anzahl der Punkte im Bereich ± 3σ vom Durchschnitt
len([result for result in results if result >result_mean + (-3) * sigma and result < result_mean + 3 * sigma])
# 998

Daraus ist ersichtlich, dass eine Zahl nahe der Normalverteilung erhalten werden kann.

Zusammenfassung

In diesem Artikel sagte ich: "Ich werde 1000 Punkte in einem Quadrat treffen und versuchen, π 1000 Mal zu berechnen." Es gibt zwei Variablen: die Anzahl der Punkte im Viereck und die Anzahl der Simulationen von π.

Durch Ändern dieser Variablen

Das ist realisiert.

Recommended Posts

Schätzung von π nach der Monte-Carlo-Methode
Monte-Carlo-Methode
Dominion-Kompressionsspiel nach Monte-Carlo-Methode analysiert
Einführung in die Monte-Carlo-Methode
Die erste Markov-Ketten-Monte-Carlo-Methode von PyStan
Simulieren Sie die Monte-Carlo-Methode in Python
Finden Sie das Verhältnis der Fläche des Biwa-Sees nach der Monte-Carlo-Methode
Fünfäugige KI bei der Suche nach Monte-Carlo-Bäumen
#Monte Carlo-Methode zum Ermitteln des Umfangsverhältnisses mit Python
Versuchen Sie, die Monte-Carlo-Methode in Python zu implementieren
[Statistik] Ich werde die Abtastung nach der Markov-Ketten-Monte-Carlo-Methode (MCMC) mit Animation erläutern.
Berechnung der kürzesten Route nach der Monte-Carlo-Methode
Mit Reiskörnern bestreuen, um das Umfangsverhältnis zu ermitteln (Monte-Carlo-Methode)
[Wissenschaftlich-technische Berechnung mit Python] Monte-Carlo-Integration, numerische Berechnung, Numpy
Finden Sie das Umfangsverhältnis mit einer dreizeiligen Funktion [Python / Monte-Carlo-Methode]
[Wissenschaftlich-technische Berechnung nach Python] Monte-Carlo-Simulation nach der Metropolenmethode der Thermodynamik des 2D-Anstiegsspinsystems
Klassifizieren Sie Daten nach der k-means-Methode
[Statistik] Visualisieren und verstehen Sie die Hamiltonsche Monte-Carlo-Methode mit Animation.
Ich konnte pypy3.6-7.3.1 nicht mit macOS + pyenv installieren, aber ich konnte pypy3.6-7.3.0 installieren. Ich fühlte den Wind von Pypy nach der Monte-Carlo-Methode.