Zum #PiTag: Ein Weg zu pi //2074

Heute geht es wieder rund.
Anlässlich des diesjährigen e-Tages hatte ich auf die Schnelle über eine numersche Berechnung von e mit einem Monte-Carlo-Verfahren breakgeplaint.
Zur Feier des heutigen pi-Tages möchte ich ein bisschen ausführlicher darüber bloggen, wie man auf ähnliche Weise pi mit dem Computer berechnen kann. Das ist übrigens eine ganz einfache Übung, die man mit so ziemlich jeder Programmiersprache machen kann, und die ich Anfängern einfach mal empfehlen möchte.

Für den Einheitskreis gilt in einem kartesischen Koordinatensystem: x^2 + y^2 = 1. Seine Fläche ist gleich pi.
Betrachten wir nur den 1. Quadranten, so liegen die Werte von x und y jeweils zwischen 0 und 1. Die Fläche ist ein Viertelkreis, also pi/4.
Wenn wir Paare von Zufallszahlen (x,y) mit x,y in [0,1[ erzeugen, so liegen diese Punkte alle innerhalb eines Quadrates mit der Fläche 1.
Ein Teil davon liegt innerhalb der Kreiskurve, ein anderer Teil außerhalb. Wenn wir davon ausgehen, dass all diese Punkte völlig zufällig und gleichverteilt sind, so wird die Anzahl der Punkte innerhalb des Kreises N_i proportional zu dessen Fläche sein, während die Gesamtanzahl N_ges proportional zur Fläche des Quadrates ist. Für große Zahlen wird also N_i/N_ges gegen pi/4 gehen.
Um zu prüfen, ob sich ein Punkt innerhalb des Kreises befindet, berechnen wir einfach die Summe der Komponentenquadrate x^2 + y^2. Ist dies kleiner (bzw. kleiner oder gleich – das ist Geschmackssache, und dürfte eigentlich auf lange Sicht nichts ändern, d.h. der Gleichheitsfall dürfte praktisch nicht vorkommen, hängt auch davon ab, ob man mit Single- oder Double-Gleitkommazahlen rechnet) als 1, so inkrementieren wir den Zähler N_i. Der Zähler N_ges wird sowieso bei jedem Durchlauf erhöht.

Weiter unten beim Bonusmaterial seht ihr ein paar Beispiel-Durchläufe.
Die erste Spalte der Tabelle gibt den Index, bzw. die bisherige Gesamtanzahl an. Dann folgen die Zufallszahlen x und y. Wenn x*x + y*y kleiner als 1 ist, so wird der Zähler N_i in der vierten Spalte inkrementiert. Der vierfache Quotient aus N_i und N_ges gibt dann in der letzten Spalte eine Näherung für pi an.
Für mehrere Durchläufe habe ich das für die ersten paar Hundert Samples dann aufgeplottet. Man sieht, dass der Näherungswert nur sehr langsam gegen pi geht. Die numerische Integration am e-Tag konvergierte dagegen wesentlich besser.
Wenn man es sich recht überlegt, ist das nicht überraschend. Bei der Integration kam auf jede Zufallszahl ein Funktionswert. Diese Funktionswerte wurden im Prinzip gemittelt und trugen so wesentlich zum Ergebnis bei. Das geht natürlich am besten, wenn die Funktion einigermaßen konstant und monoton ist, also nicht volatil und auch nicht stärker gekrümmt. Eventuelle Singularitäten im Integrationsbereich könnten auch durch das Monte-Carlo-Verfahren nicht abgefangen werden.
Beim hier benutzen Algorithmus dagegen braucht man einen Punkt aus zwei Zufallszahlen, bei dem man erst abprüfen muss, in welchem Bereich er liegt. Dies ist aber nur das Kriterium um einen zugeordneten Zähler zu inkrementieren oder nicht. Da verwundert es nicht, dass man nur recht langsam an einen guten Wert für pi herankommt (auch wenn man zufällig nach 14 Würfen durchaus 11 Treffer haben kann). Das ist ähnlich ineffizient, als wolle man den Wert für ein Halb durch häufige Münzwürfe ermitteln, indem man die Anzahl von „Zahl“ zur Gesamtanzahl der Würfe ins Verhältnis setzt.

Und wer sich wundert, dass die Ergebnisse manchmal längere Zeit weit über oder weit unter dem exakten Zahlenwert für pi liegen: tja, so ist das halt mit Statistiken. Erst für sehr große Zahlen haben sie eine Aussagekraft. Bei kleineren Ensembles überwiegt der Einfluss des individuellen Falls.

---

N_ges	x	y	N_i	4*N_i/N_ges
1	0.998967118794098	0.595742204692215	0	0
2	0.975805654888973	0.467658712528646	0	0
3	0.187581549631432	0.855376057792455	1	1.33333333333333
4	0.609713102458045	0.0805352982133627	2	2
5	0.291903404751793	0.692890904378146	3	2.4
6	0.957869877805933	0.529505324549973	3	2
7	0.604051471920684	0.221956855151802	4	2.28571428571429
8	0.604007393820211	0.560130152851343	5	2.5
9	0.736354042077437	0.676822786685079	5	2.22222222222222
10	0.332407142734155	0.909002791158855	6	2.4
11	0.197703889803961	0.481591243762523	7	2.54545454545455
12	0.411775207845494	0.410620452836156	8	2.66666666666667
13	0.365421066293493	0.29703240422532	9	2.76923076923077
14	0.766812087735161	0.542733677662909	10	2.85714285714286
15	0.649498504819348	0.0293121063150465	11	2.93333333333333
16	0.959352832986042	0.0195470154285431	12	3
17	0.896105438703671	0.0350089385174215	13	3.05882352941176
18	0.150952495401725	0.292146272026002	14	3.11111111111111
19	0.327223631786183	0.00679539190605283	15	3.15789473684211
20	0.468791889725253	0.0655273180454969	16	3.2