Zufallszahlen manuell generieren

Wie kann ich aus einer gegebenen Verteilung manuell eine Zufallszahl generieren, beispielsweise 10 Realisierungen aus der Standardnormalverteilung?

Wenn "manuell" "mechanisch" einschließt, stehen Ihnen viele Optionen zur Verfügung. Um eine Bernoulli-Variable mit der Hälfte der Wahrscheinlichkeit zu simulieren, können wir eine Münze werfen: für Schwänze, 1 für Köpfe. Um eine geometrische Verteilung zu simulieren, können wir zählen, wie viele Münzwürfe benötigt werden, bevor wir Köpfe erhalten. Um eine Binomialverteilung zu simulieren, können wir unsere Münze n- mal werfen (oder einfach n Münzen werfen ) und die Köpfe zählen. Die "quincunx" oder "bean machine" oder "Galton box" ist eine kinetischere Alternative - warum nicht eine in Aktion setzen und selbst sehen ? Es scheint, als gäbe es keine "gewichtete Münze"$0$$1$$n$$n$Wenn wir jedoch den Wahrscheinlichkeitsparameter unserer Bernoulli- oder Binomialvariablen auf andere Werte als ändern möchten, können wir dies mit der Nadel von Georges-Louis Leclerc, Comte de Buffon , tun. Um die diskrete Gleichverteilung auf { 1 , 2 , 3 , 4 , 5 , 6 } zu simulieren , werfen wir einen sechsseitigen Würfel. Fans von Rollenspielen werden auf exotischere Würfel gestoßen sein , zum Beispiel tetraedrische Würfel, die gleichmäßig aus { 1 , 2 , 3 , 4 } entnommen werden können.$p = 0.5$$\{1, 2, 3, 4, 5, 6\}$$\{1,2,3,4\}$Während mit einem Spinner oder Roulette-Rad kann man noch weiter gehen. ( Bildnachweis )

Müssten wir heute verrückt sein, um auf diese Weise Zufallszahlen zu generieren, wenn es nur ein Befehl auf einer Computerkonsole ist - oder wenn wir eine geeignete Tabelle mit Zufallszahlen zur Verfügung haben, einen Streifzug zu den staubigeren Ecken des Bücherregals? Nun vielleicht, obwohl ein physikalisches Experiment etwas angenehm Tastbares hat. Für Menschen, die vor dem Computerzeitalter arbeiteten, war es jedoch praktischer, Zufallsvariablen manuell zu simulieren, und zwar vor weit verbreiteten Zufallszahltabellen (von denen später mehr zur Verfügung stehen). Als Buffon das St. Petersburg-Paradoxon untersuchte- das berühmte Münzwurfspiel, bei dem sich der Gewinn jedes Mal verdoppelt, wenn ein Kopf geworfen wird, der Spieler bei den ersten Schwänzen verliert und dessen erwartete Auszahlung kontraintuitiv unendlich ist - er musste die geometrische Verteilung mit simulieren = 0,5 . Zu diesem Zweck hat er anscheinend ein Kind angeheuert, das eine Münze wirft, um 2048 Spiele des St. Petersburg-Spiels zu simulieren. Diese simulierte geometrische Verteilung ist in Stigler (1991) wiedergegeben :$p=0.5$

Tosses Frequency
1      1061
2      494
3      232
4      137
5      56
6      29
7      25
8      8
9      6

In demselben Aufsatz, in dem er diese empirische Untersuchung des St. Petersburg-Paradoxons veröffentlichte, stellte Buffon auch die berühmte " Buffon's Needle " vor. Wenn eine Ebene in einem Abstand durch parallele Linien in Streifen unterteilt ist und eine Nadel der Länge l ≤ d darauf fällt, beträgt die Wahrscheinlichkeit, dass die Nadel eine der Linien kreuzt, 2 $d$$l \leq d$ .$\frac{2l}{\pi d}$

Mit Buffons Nadel kann daher eine Zufallsvariable X ∼ Bernoulli ( 2 l) simuliert werdenoderX∼Binomial(n,2$X \sim \text{Bernoulli}(\frac{2l}{\pi d})$, und wir können die Erfolgswahrscheinlichkeit anpassen, indem wir die Länge unserer Nadeln oder (vielleicht bequemer) den Abstand ändern, in dem wir die Linien anordnen. Eine alternative Verwendung von Buffon-Nadeln ist ein äußerst ineffizienter Weg, um eine probabilistische Approximation fürπ zu finden. Das Bild (Kredit) zeigt 17 Streichhölzer, von denen 11 eine Linie kreuzen. Wenn der Abstand zwischen den Linien gleich der Länge des Streichholzes eingestellt ist, wie hier, beträgt der erwartete Anteil an sich kreuzenden Streichhölzern$X \sim \text{Binomial}(n,\frac{2l}{\pi d})$$\pi$ und somit können wir abschätzen , π als doppelt reziproken Wert der beobachteten Fraktion: hier erhält man π =2⋅$\frac{2}{\pi}$$\hat \pi$. Verwendung2,5 cm Nadeln mit Linien 3 cm voneinander entfernt 1901 beanspruchte Mario Lazzarinium das Experiment durchgeführt wird, und nach 3408 Würfen erhielten π =$\hat \pi = 2 \cdot \frac{17}{11} \approx 3.1$ . Dies ist ein bekanntes Rational fürπ, das auf sechs Dezimalstellen genau ist. Badger (1994) liefert überzeugende Beweise dafür, dass dies betrügerisch war, nicht zuletzt, dass, um mit Lazzarinis Apparat eine Genauigkeit von sechs Dezimalstellen zu 95% zu erreichen, ein Geduldsspiel von 134 Billionen Nadeln geworfen werden muss! Mit Sicherheit ist Buffons Nadel als Zufallszahlengenerator nützlicher als als Methode zur Schätzung vonπ.$\hat \pi = \frac{355}{113}$$\pi$$\pi$

Unsere Generatoren waren bisher enttäuschend diskret. Was ist, wenn wir eine Normalverteilung simulieren wollen? Eine Möglichkeit besteht darin, zufällige Ziffern zu erhalten und daraus gute diskrete Näherungen für eine gleichmäßige Verteilung auf bilden. Anschließend werden einige Berechnungen durchgeführt, um diese in zufällige normale Abweichungen umzuwandeln. Ein Dreh- oder Rouletterad kann Dezimalstellen von null bis neun angeben. ein Würfel kann binäre Ziffern erzeugen; Wenn unsere arithmetischen Fähigkeiten mit einer funkigeren Basis fertig werden, würde sogar ein Standardsatz Würfel ausreichen. Andere Antworten haben diese Art von transformationsbasiertem Ansatz ausführlicher behandelt. Ich verschiebe jede weitere Diskussion bis zum Ende.$[0,1]$

Im späten neunzehnten Jahrhundert war die Nützlichkeit der Normalverteilung bekannt, und so gab es Statistiker, die zufällige normale Abweichungen simulieren wollten. Es erübrigt sich zu erwähnen, dass langwierige Handberechnungen nicht geeignet gewesen wären, außer den Simulationsprozess an erster Stelle einzurichten. Sobald dies festgestellt wurde, musste die Erzeugung der Zufallszahlen relativ schnell und einfach sein. Stigler (1991) listet die Methoden auf, die drei Statistiker dieser Ära anwenden. Alle untersuchten Glättungstechniken: Zufällige normale Abweichungen waren von offensichtlichem Interesse, z. B. um Messfehler zu simulieren, die geglättet werden mussten.

Der bemerkenswerte amerikanische Statistiker Erastus Lyman De Forest war an der Glättung von Lebenstabellen interessiert und stieß auf ein Problem, das die Simulation der absoluten Werte normaler Abweichungen erforderte. In einem Thema, das sich als laufend erweisen wird, hat De Forest tatsächlich eine Stichprobe aus einer Halbnormalverteilung gezogen . Außerdem, anstatt eine Standardabweichung von einem mit (dem wir Aufruf „Standard“ verwendet), wollte De Forest einen „wahrscheinlichen Fehler“ (mittlere Abweichung) von einem. Dies war die Form, die in der Tabelle "Wahrscheinlichkeit von Fehlern" in den Anhängen von "Ein Handbuch der sphärischen und praktischen Astronomie, Band II" von angegeben ist$Z \sim N(0, 1^2)$William Chauvenet . Aus dieser Tabelle interpolierte De Forest die Quantile einer Halbnormalverteilung von bis p = 0,995 , die er als "Fehler gleicher Häufigkeit" ansah .$p=0.005$$p=0.995$

Wenn Sie die Normalverteilung nach De Forest simulieren möchten, können Sie diese Tabelle ausdrucken und ausschneiden. De Forest (1876) schrieb, dass die Fehler "auf 100 gleichgroße Pappbits geschrieben wurden, die in einer Schachtel geschüttelt und alle einzeln herausgezogen wurden".

Der Astronom und Meteorologe Sir George Howard Darwin (Sohn des Naturforschers Charles) veränderte die Dinge, indem er ein sogenanntes "Roulette" zur Erzeugung zufälliger normaler Abweichungen entwickelte. Darwin (1877) beschreibt, wie:

$x$$\frac{720}{\sqrt \pi} \int_0^x e^{-x^2} dx$$+$$-$$+$$-$

"Index" sollte hier als "Zeiger" oder "Indikator" gelesen werden (vgl. "Zeigefinger"). Stigler weist darauf hin, dass Darwin wie De Forest eine halbnormale kumulative Verteilung um die Festplatte verwendete. Die anschließende Verwendung einer Münze zum Anbringen eines zufälligen Zeichens führt zu einer vollständigen Normalverteilung. Stigler merkt an, dass es unklar ist, wie fein die Skala abgestuft wurde, geht aber davon aus, dass die Anweisung zum manuellen Anhalten der Scheibe während der Drehung "die potentielle Vorspannung in Richtung eines Abschnitts der Scheibe zu verringern und den Vorgang zu beschleunigen" lautete.

Sir Francis Galton , übrigens ein Halbcousin von Charles Darwin, wurde bereits im Zusammenhang mit seinem Quincunx erwähnt. Während dies mechanisch eine Binomialverteilung simuliert, die nach dem De Moivre-Laplace-Theorem eine auffallende Ähnlichkeit mit der Normalverteilung aufweist (und gelegentlich als Lehrmittel für dieses Thema verwendet wird), hat Galton tatsächlich ein weitaus aufwändigeres Schema erstellt, als er es wünschte Probe aus einer Normalverteilung. Noch außergewöhnlicher als die unkonventionellen Beispiele oben in dieser Antwort, entwickelte Galton normalverteilte Würfel- oder genauer gesagt, ein Würfelsatz, der eine ausgezeichnete diskrete Annäherung an eine Normalverteilung mit einer mittleren Abweichung von eins ergibt. Diese Würfel aus dem Jahr 1890 werden in der Galton Collection am University College London aufbewahrt.

In einem Artikel von 1890 in Nature schrieb Galton:

Als Instrument zur zufälligen Auswahl habe ich nichts Besseres als Würfel gefunden. Es ist am mühsamsten, die Karten zwischen den aufeinanderfolgenden Zügen gründlich zu mischen, und die Methode, markierte Bälle in einem Beutel zu mischen und aufzurühren, ist noch mühsamer. Ein Teetotum oder irgendeine Form von Roulette ist diesen vorzuziehen, aber Würfel sind besser als alle. Wenn sie geschüttelt und in einen Korb geworfen werden, rasen sie so unterschiedlich gegeneinander und gegen die Rippen der Korbarbeit, dass sie wild umherstolpern, und ihre Positionen am Anfang geben keinen erkennbaren Hinweis darauf, was sie nach einer Weile sein werden einmal gut schütteln und werfen. Die Chancen, die ein Würfel bietet, sind vielfältiger als allgemein angenommen. Es gibt 24 gleiche Möglichkeiten und nicht nur 6, da jede Fläche vier Kanten hat, die verwendet werden können, wie ich zeigen werde.

Für Galton war es wichtig, in der Lage zu sein, schnell eine Folge normaler Abweichungen zu erzeugen . Nach jedem Wurf richtete Galton die Würfel nur durch Berühren aus und notierte die Punkte an den Vorderkanten. Er würfelte anfangs mehrere Würfel vom Typ I, deren Kanten halbnormale Abweichungen aufwiesen, ähnlich wie bei De Forest, aber mit 24 nicht 100 Quantilen. Für die größten Abweichungen (die auf den Würfeln vom Typ I als Leerzeichen markiert waren) würfelte er so viele der empfindlicheren Würfel vom Typ II (die nur bei einer feineren Abstufung große Abweichungen aufwiesen), wie er brauchte, um die Felder in seiner Sequenz auszufüllen . Um von halbnormalen zu normalen Abweichungen zu konvertieren, würfelt er mit Würfel III, der zuweist.$+$$-$$1 \frac 1 4$

Das Labor für mathematisch-statistische Experimente von Raazesh Sainudiin umfasst ein Studentenprojekt der Universität von Canterbury, Neuseeland, das Galtons Würfel reproduziert . Das Projekt umfasst empirische Untersuchungen, bei denen die Würfel viele Male gewürfelt wurden (einschließlich einer empirischen CDF, die beruhigend "normal" aussieht), und eine Anpassung der Würfelbewertungen, damit sie der Standardnormalverteilung folgen. Unter Verwendung von Galtons Originalwerten gibt es auch ein Diagramm der diskretisierten Normalverteilung, der die Würfelwerte tatsächlich folgen.

Wenn Sie im großen Stil bereit sind, das "Mechanische" zum Elektrischen zu dehnen, beachten Sie, dass RANDs Epos " A Million Random Digits" mit 100.000 normalen Abweichungen auf einer Art elektronischer Simulation eines Rouletterads basierte. Aus dem technischen Bericht (von George W. Brown, ursprünglich Juni 1949) entnehmen wir:

So motiviert entwickelten die RAND-Leute mit Unterstützung des Ingenieurs der Douglas Aircraft Company ein Elektro-Roulette-Rad, das auf einer Variation eines Vorschlags von Cecil Hastings basierte. Für die Zwecke dieses Vortrags ist eine kurze Beschreibung ausreichend. Eine Zufallsfrequenz-Impulsquelle wurde etwa einmal pro Sekunde durch einen Impuls mit konstanter Frequenz gesteuert, was durchschnittlich etwa 100.000 Impulse in einer Sekunde ergab. Impulsstandardisierungsschaltungen leiteten die Impulse an einen fünfstelligen Binärzähler weiter, so dass die Maschine im Prinzip wie ein Rouletterad mit 32 Positionen aussieht und durchschnittlich etwa 3000 Umdrehungen pro Umdrehung ausführt. Eine Binär-zu-Dezimal-Umwandlung wurde verwendet, wobei 12 der 32 Positionen weggeworfen wurden, und die resultierende Zufallsziffer wurde in einen IBM-Durchschlag eingespeist, wodurch Lochkartentabellen mit Zufallsziffern erhalten wurden.

$\chi^2$Tests der Häufigkeiten von geraden und ungeraden Ziffern ergaben, dass einige Chargen eine leichte Unwucht aufwiesen. Dies war in einigen Chargen schlimmer als in anderen, was darauf hindeutet, dass "die Maschine in dem Monat seit ihrer Einstellung ausgefallen ist ... Die Angaben zu dieser Maschine erfordern eine übermäßige Wartung, um sie in einwandfreiem Zustand zu halten". Es wurde jedoch ein statistischer Weg zur Lösung dieser Probleme gefunden:

Zu diesem Zeitpunkt hatten wir unsere ursprünglichen Millionen Ziffern, 20.000 IBM-Karten mit 50 Ziffern auf einer Karte, wobei die statistische Analyse eine geringe, aber wahrnehmbare ungerade Tendenz aufzeigte. Es wurde nun beschlossen, den Tisch neu zu ordnen oder zumindest zu verändern, indem mit ihm ein kleines Roulette gespielt wurde, um die ungerade-gerade-Tendenz zu beseitigen. Wir haben (Mod 10) die Ziffern jeder Karte ziffernweise zu den entsprechenden Ziffern der vorherigen Karte hinzugefügt. Die abgeleitete Tabelle mit einer Million Stellen wurde dann den verschiedenen Standardtests, Frequenztests, Serientests, Pokertests usw. unterzogen. Diese Millionen Stellen weisen eine saubere Gesundheitsbewertung auf und wurden als moderne Tabelle mit zufälligen Stellen von RAND übernommen.

Es gab natürlich guten Grund zu der Annahme, dass das Hinzufügen etwas Gutes bewirken würde. Im Allgemeinen ist der zugrunde liegende Mechanismus der begrenzende Ansatz von Summen von Zufallsvariablen, die das Einheitsintervall in der Rechteckverteilung auf die gleiche Weise modulieren, wie sich uneingeschränkte Summen von Zufallsvariablen der Normalität nähern. Diese Methode wurde von Horton und Smith von der Interstate Commerce Commission verwendet, um einige gute Chargen von scheinbar zufälligen Zahlen aus größeren Chargen von schlecht nicht zufälligen Zahlen zu erhalten.

$[0,1]$$u$$[0,1]$$F$$F^{-1}(u)$

Verweise

• Badger, L. (1994). " Lazzarinis Lucky Approximation von π ". Mathematik-Magazin . Mathematische Vereinigung von Amerika. 67 (2): 83–91.

• $(*)$

• Darwin, GH (1877). " Über fehlbare Maße variabler Mengen und über die Behandlung meteorologischer Beobachtungen. " Philosophical Magazine , 4 (22), 1–14

• De Forest, EL (1876). Interpolation und Anpassung von Serien . Tuttle, Morehouse und Taylor, New Haven, Conn.

• Galton, F. (1890). "Würfel für statistische Experimente". Nature , 42 , 13 & ndash ; 14

• Stigler, SM (1991). "Stochastische Simulation im neunzehnten Jahrhundert". Statistical Science , 6 (1), 89 & ndash; 97.

$(*)$"Wer arithmetische Methoden zur Erzeugung von Zufallszahlen in Betracht zieht, befindet sich natürlich in einem Zustand der Sünde. Denn es gibt, wie bereits mehrfach erwähnt, keine Zufallszahl - es gibt nur Methoden zur Erzeugung von Zufallszahlen und ein strenges Rechenverfahren ist natürlich keine solche Methode. "

Zum Beispiel kann ich auf meinem Linux-Betriebssystem überprüfen

$ date +%s.%N
1479733744.077762986
$ date +%s.%N
1479733980.615056616

> sqrt(-2*log(.077762986))*cos(2*pi*.615056616)
[1] -1.694815

Nachtrag: Da die Berechnung von Logarithmen und Cosinus als nicht manuell genug angesehen werden kann, gibt es für Box-Müller eine Variante, die die Verwendung dieser transzendentalen Funktionen vermeidet (siehe Aufgabe 2.9 in unserem Buch Monte Carlo Statistical Methods ):

Nun kann man wegen der Exponentialvariablen gegen diese Version argumentieren. Es gibt jedoch auch eine sehr clevere Möglichkeit, diese Variablen ohne Aufruf transzendenter Funktionen zu simulieren , wie in diesem Algorithmus aus Luc Devroyes uneinheitlichen Zufallsvarianten zusammengefasst :

Zugegeben, es ist die Berechnung von 1 / e erforderlich, jedoch nur einmal.

$(0,1)^2$$(0,1)$$(0,2\pi)$

Die Verwendung der CLT zur Annäherung an die Normalität ist sicherlich keine Methode, die ich jemals empfehlen würde, da (1) Sie immer noch andere Variablen benötigen, um den Durchschnitt zu ermitteln. Daher können Sie auch Uniformen im Box-Müller-Algorithmus verwenden, und (2) die Genauigkeit wächst erheblich langsam mit der Anzahl der Simulationen. Vor allem, wenn Sie eine diskrete Zufallsvariable wie das Ergebnis eines Würfels verwenden, auch bei mehr als sechs Gesichtern . Zitat von Thomas et al. (2007), eine Umfrage zum Für und Wider von Gaußschen Zufallsgeneratoren:

Der zentrale Grenzwertsatz ist natürlich ein Beispiel für eine "Näherungsmethode" - selbst wenn eine perfekte Arithmetik verwendet wird, ist die Ausgabe für endliches K nicht Gaußsch.

Hier ist ein kurzer Versuch, um das Problem zu veranschaulichen: Ich habe den 100-fachen Durchschnitt von 30 Würfelergebnissen generiert:

dies=apply(matrix(sample(1:6,30*100,rep=TRUE),ncol=30),1,mean)

dann normalisierte man diese Mittelwerte zu einem Mittelwert von Null - Varianz, die man variiert

stdies=(dies-3.5)/sqrt(35/12/30)

und betrachtete den normalen Sitz (oder das Fehlen davon) dieser Probe:

dies76/30122/30$D_i$

dies=matrix(apply(matrix(sample(0:5,15*200,rep=TRUE),nrow=15)/6^(1:15),2,sum),ncol=2) 
norma=sqrt(-2*log(dies[,1]))*c(cos(2*pi*dies[,2]),sin(2*pi*dies[,2]))

Die Passform ist so gut, wie es für eine normale Stichprobe der Größe 200 zu erwarten ist (zeichnen Sie eine weitere Stichprobe für eine echte normale Stichprobe auf norma=rnorm(100)):

wie ein Kolmogorov-Smirnov-Test weiter zeigt:

> ks.test(norma,pnorm)

        One-sample Kolmogorov-Smirnov test

data:  norma
D = 0.06439, p-value = 0.3783
alternative hypothesis: two-sided

Dies ist nicht genau zufällig, aber es sollte nahe genug sein, wie Sie scheinen ein grobes Experiment zu wollen.

Verwenden Sie Ihr Telefon, um einen Chronometer einzurichten. Nach gut 10 Sekunden hören Sie auf (je länger Sie warten, desto eher nähern Sie sich einem wirklich "zufälligen" Ergebnis, aber 10 Sekunden sind in Ordnung). Nehmen Sie die letzten Ziffern (z. B. 10,67 Sekunden ergeben 67). Wenden Sie die Perzentiltabelle für die Normalverteilung an. In diesem Beispiel müssen Sie nur nach 0,67 suchen, und Sie finden die Nummer. In diesem Fall beträgt Ihr Wert ungefähr 0,45. Dies ist nicht ganz genau, aber es gibt Ihnen eine solide Schätzung.

Wenn Sie unter 50 rollen, tun Sie einfach 100- [Ihr Ergebnis] und verwenden Sie die Tabelle. Ihr Ergebnis ist aufgrund der Symmetrie von N (0,1) dasselbe mit einem Minuszeichen.

$n$$+1$$-1$$n$$\sqrt n$$n$$N (0,1)$

Warum? Lassen

$\mathbb P (X_k = \pm 1) = \frac 12$

$Y := X_1 + X_2 + \cdots + X_n$

Normalisieren,

Wir erhalten eine Zufallsvariable mit Einheitsvarianz

Es ist erwähnenswert, dass Sie, sobald Sie eine Uniform (0,1) generieren können, eine beliebige Zufallsvariable generieren können, für die die inverse cdf berechnet werden kann, indem Sie einfach die einheitliche Zufallsvariable in die inverse CDF einfügen.

Wie kann man also eine Uniform (0,1) manuell berechnen? Nun, wie von @Silverfish erwähnt, gibt es eine Vielzahl von Würfeln, die von traditionellen RPG-Spielern verwendet werden. Eines davon ist ein zehnseitiger Würfel. Unter der Annahme, dass dies ein fairer Würfel ist, können wir jetzt eine diskrete Uniform erzeugen (0, 9).

$0.01, 0.02, ..., 0.99, 1.00$$0.001, 0.002, ..., 0.999, 1.000$

Wir können also einer durchgehenden Uniform (0,1) sehr nahe kommen, indem wir sie mit einer fein gitterten diskreten Uniformverteilung mit einigen 10-seitigen Würfeln approximieren. Dieses kann dann in eine inverse CDF eingesteckt werden, um die interessierende Zufallsvariable zu erzeugen.