Simulation von Zügen aus einer Gleichverteilung mit Zügen aus einer Normalverteilung

Ich habe kürzlich eine Data Science-Interviewressource gekauft, in der eine der Wahrscheinlichkeitsfragen wie folgt lautete:

Wie kann man bei gegebenen Ziehungen aus einer Normalverteilung mit bekannten Parametern Ziehungen aus einer Gleichverteilung simulieren?

Mein ursprünglicher Denkprozess war, dass wir für eine diskrete Zufallsvariable die Normalverteilung in K eindeutige Unterabschnitte aufteilen könnten, wobei jeder Unterabschnitt eine gleiche Fläche unter der Normalkurve hat. Dann können wir bestimmen, welchen der K-Werte die Variable annimmt, indem wir erkennen, in welchen Bereich der Normalkurve die Variable letztendlich fällt.

Dies würde jedoch nur für diskrete Zufallsvariablen funktionieren. Ich habe nachgeforscht, wie wir dasselbe für kontinuierliche Zufallsvariablen tun könnten, aber leider konnte ich nur Techniken wie Inverse Transformation Sampling finden, die als Eingabe einer einheitlichen Zufallsvariablen dienen und Zufallsvariablen aus einer anderen Verteilung ausgeben . Ich dachte, dass wir diesen Prozess vielleicht in umgekehrter Reihenfolge durchführen könnten, um einheitliche Zufallsvariablen zu erhalten?

Ich habe auch darüber nachgedacht, die normalen Zufallsvariablen möglicherweise als Eingaben in einen linearen Kongruenzgenerator zu verwenden, bin mir aber nicht sicher, ob dies funktionieren würde.

Irgendwelche Gedanken darüber, wie ich diese Frage angehen könnte?

Im Sinne der Verwendung einfacher algebraischer Berechnungen, die nichts mit der Berechnung der Normalverteilung zu tun haben , möchte ich auf Folgendes eingehen. Sie sind so angeordnet, wie ich es mir vorgestellt habe (und mussten daher immer kreativer werden), aber ich habe das Beste - und das Überraschendste - für die Ewigkeit aufgehoben.

1. Kehre die Box-Mueller-Technik um : Aus jedem Normalenpaar können zwei unabhängige Uniformen als (im Intervall ) und (im Intervall ).atan2 ( Y , X ) [ - π , π ] exp ( - ( X 2 + Y 2 ) / 2 ) [ 0 , 1 $(X,Y)$$\text{atan2}(Y,X)$$[-\pi, \pi]$$\exp(-(X^2+Y^2)/2)$$[0,1]$

2. Nehmen Sie die Normalen in Zweiergruppen und Sie ihre Quadrate, um eine Folge von Variablen . Die Ausdrücke aus den Paaren erhalten Y 1 , Y 2 , … , Y i , $\chi^2_2$$Y_1, Y_2, \ldots, Y_i, \ldots$

   $$X_i = \frac{Y_{2i}}{Y_{2i-1}+Y_{2i}}$$

   wird eine -Verteilung haben, die einheitlich ist.$\text{Beta}(1,1)$

   Dass dies nur grundlegende, einfache Arithmetik erfordert, sollte klar sein.

3. Da die exakte Verteilung des Pearson - Korrelationskoeffizienten einer 4 - Paar - Stichprobe aus einer bivariaten Standardnormalverteilung gleichmäßig auf verteilt ist , können wir die Normalen einfach in Gruppen von 4 Paaren (dh 8 Werten in jede Menge) und geben den Korrelationskoeffizienten dieser Paare zurück. (Dies beinhaltet einfache Arithmetik plus zwei Quadratwurzeloperationen.)$[-1,1]$

4. Seit der Antike ist bekannt, dass eine zylindrische Projektion der Kugel (eine Fläche im Dreiraum) flächengleich ist . Dies impliziert, dass bei der Projektion einer gleichmäßigen Verteilung auf die Kugel sowohl die horizontale Koordinate (entsprechend dem Längengrad) als auch die vertikale Koordinate (entsprechend dem Breitengrad) gleichmäßige Verteilungen aufweisen. Da die trivariate Normalverteilung sphärisch symmetrisch ist, ist ihre Projektion auf die Kugel gleichmäßig. Das Erhalten des Längengrads ist im Wesentlichen dieselbe Berechnung wie der Winkel in der Box-Müller-Methode ( siehe ), aber der projizierte Breitengrad ist neu. Die Projektion auf die Kugel normalisiert lediglich ein Dreifach der Koordinaten und an diesem Punktz X 3 i - 2 , X 3 i - 1 , X 3 $(x,y,z)$$z$ist der projizierte Breitengrad. Nehmen Sie also die Normalvariablen in Dreiergruppen, , und berechnen Sie$X_{3i-2}, X_{3i-1}, X_{3i}$

   $$\frac{X_{3i}}{\sqrt{X_{3i-2}^2 + X_{3i-1}^2 + X_{3i}^2}}$$

   für .$i=1, 2, 3, \ldots$

5. Da die meisten Computersysteme Zahlen in repräsentieren binären , in der Regel einheitliche Zahlenerzeugung beginnt durch gleichmäßig verteilte Erzeugung ganze Zahlen zwischen und (oder einer hohen Leistung von bis Computer - Wortlänge bezogen) und Neuskalierung sie nach Bedarf. Solche Ganzzahlen werden intern als Zeichenfolgen mit Binärziffern dargestellt. Wir können unabhängige Zufallsbits erhalten, indem wir eine Normalvariable mit ihrem Median vergleichen. Es reicht also aus, die Normal-Variablen in Gruppen zu unterteilen, die der gewünschten Anzahl von Bits entsprechen, diese mit ihrem Mittelwert zu vergleichen und die resultierenden Folgen von Wahr / Falsch-Ergebnissen zu einer Binärzahl zusammenzusetzen. Schreiben2 32 - 1 2 32 k H H ( x ) = 1 x > 0 H ( x ) = 0 [ 0 , 1 $0$$2^{32}-1$$2$$32$$k$für die Anzahl der Bits und für das Vorzeichen (d. h. wenn und andernfalls) können wir den resultierenden normierten einheitlichen Wert in mit ausdrücken Formel$H$$H(x)=1$$x\gt 0$$H(x)=0$$[0, 1)$

   $$\sum_{j=0}^{k-1} H(X_{ki - j})2^{-j-1}.$$

   Die Variablen können aus jeder stetigen Verteilung gezogen werden, deren Median (z. B. eine Standardnormalverteilung). Sie werden in Gruppen von wobei jede Gruppe einen solchen pseudoeinheitlichen Wert erzeugt. 0 $X_n$$0$$k$

6. Die Rückweisungsabtastung ist eine standardmäßige, flexible und leistungsstarke Methode, um Zufallsvariablen aus beliebigen Verteilungen zu ziehen. Angenommen, die Zielverteilung hat PDF . Ein Wert wird gemäß einer anderen Verteilung mit PDF gezeichnet . Im Zurückweisungsschritt wird unabhängig von ein einheitlicher Wert , der zwischen und liegt, und mit verglichen : Wenn er kleiner ist, wird beibehalten, andernfalls wird der Vorgang wiederholt. Dieser Ansatz scheint jedoch kreisförmig zu sein: Wie erzeugen wir eine einheitliche Variable mit einem Prozess, der zunächst eine einheitliche Variable benötigt?Y g U 0 g ( Y ) Y f ( Y ) $f$$Y$$g$$U$$0$$g(Y)$$Y$$f(Y)$$Y$

   Die Antwort ist, dass wir eigentlich keine einheitliche Variable benötigen, um den Ablehnungsschritt auszuführen. Stattdessen können wir (unter der Annahme von ) eine faire Münze werfen, um zufällig eine oder zu erhalten . Dies wird als erstes Bit in der Binärdarstellung einer einheitlichen Variablen im Intervall ] interpretiert . Wenn das Ergebnis ist, bedeutet dies ; ansonsten . Die Hälfte der Zeit reicht aus, um den Ablehnungsschritt zu entscheiden: Wenn , die Münze jedoch , sollte akzeptiert werden; wenn0 1 U [ 0 , 1 ) 0 0 ≤ U < 1 / 2 1 / 2 ≤ U < 1 f ( Y ) / g ( Y ) ≥ 1 / 2 0 Y f ( Y ) / g ( Y ) < 1 / 2 1 Y U $g(Y)\ne 0$$0$$1$$U$$[0,1)$$0$$0 \le U \lt 1/2$$1/2\le U \lt 1$$f(Y)/g(Y) \ge 1/2$$0$$Y$$f(Y)/g(Y) \lt 1/2$ aber die Münze ist , sollte abgelehnt werden; Andernfalls müssen wir die Münze erneut werfen, um das nächste Bit von . Weil - egal welchen Wert hat - es eine Chance gibt, nach jedem anzuhalten, beträgt die erwartete Anzahl von nur .$1$$Y$$U$1 / 2 1 / 2 ( 1 ) + 1 / 4 ( 2 ) + 1 / 8 ( 3 ) + ⋯ + 2 - n ( n ) + ⋯ = $f(Y)/g(Y)$$1/2$$1/2(1)+1/4(2)+1/8(3)+\cdots+2^{-n}(n)+\cdots=2$

   Die Auswahl der Ablehnungen kann sich lohnen (und effizient sein), vorausgesetzt, die erwartete Anzahl von Ablehnungen ist gering. Dies erreichen wir, indem wir das größtmögliche Rechteck (das eine gleichmäßige Verteilung darstellt) unter eine normale PDF-Datei einfügen.

   

   Wenn Sie Calculus verwenden, um die Fläche des Rechtecks ​​zu optimieren, sollten die Endpunkte bei , wobei die Höhe , wodurch die Fläche etwas kleiner wird größer als . Indem Sie diese Standardnormaldichte als und alle Werte außerhalb des Intervalls automatisch verwerfen und ansonsten das Verwerfungsverfahren anwenden, erhalten Sie gleichmäßige Variationen in effizient:exp ( - 1 / 2 ) / $\pm 1$0,48g[-1,1][-1,1$\exp(-1/2)/\sqrt{2\pi}\approx 0.241971$$0.48$$g$$[-1,1]$$[-1,1]$

   • In einem Bruchteil von der Zeit liegt die jenseits von und wird sofort verworfen. ( ist die normale Standard-CDF.)[ - 1 , 1 ] $2\Phi(-1) \approx 0.317$$[-1,1]$$\Phi$

   • In dem verbleibenden Teil der Zeit muss die binäre Zurückweisungsprozedur befolgt werden, wobei durchschnittlich zwei weitere Normalvariablen erforderlich sind.

   • Die Gesamtprozedur erfordert durchschnittlich Schritte.$1/(2\exp(-1/2)/\sqrt{2\pi}) \approx 2.07$

   Die erwartete Anzahl von Normalvariablen, die benötigt werden, um jedes einheitliche Ergebnis zu erzielen, wird berechnet

   $$\sqrt{2 e \pi}\left(1-2\Phi(-1)\right) \approx 2.82137.$$

   Obwohl dies ziemlich effizient ist, ist zu beachten, dass (1) für die Berechnung des normalen PDF- Dokuments die Berechnung eines Exponentials erforderlich ist und (2) der Wert ein für alle Mal vorberechnet werden muss. Es ist immer noch etwas weniger berechnend als die Box-Müller-Methode ( vgl. ).$\Phi(-1)$

7. Die Ordnungsstatistik einer Gleichverteilung weist exponentielle Lücken auf. Da die Summe der Quadrate von zwei Normalen (mit dem Mittelwert Null) exponentiell ist, können wir eine Realisierung von unabhängigen Uniformen erzeugen, indem wir die Quadrate von Paaren solcher Normalen summieren, die kumulative Summe davon berechnen und die Ergebnisse neu skalieren, damit sie in das Intervall fallen und den letzten fallen lassen (der immer gleich ). Dies ist ein erfreulicher Ansatz, da nur Quadrieren, Summieren und (am Ende) eine einzige Division erforderlich ist.[ 0 , 1 ] $n$$[0,1]$$1$

   Die Werte werden automatisch in aufsteigender Reihenfolge angezeigt. Wenn eine solche Sortierung gewünscht wird, ist diese Methode allen anderen rechnerisch überlegen , da sie die -Kosten einer Sortierung vermeidet . Wenn jedoch eine Folge unabhängiger Uniformen benötigt wird, reicht es aus, diese Werte zufällig zu sortieren . Da (wie in der Box-Müller-Methode zu sehen ist, siehe ) die Verhältnisse jedes Normalenpaares unabhängig von der Summe der Quadrate jedes Paares sind, haben wir bereits die Möglichkeit, diese zufällige Permutation zu erhalten: Ordnen Sie die kumulativen Summen nach den entsprechenden Verhältnissen . (Wenn sehr groß ist, könnte dieser Prozess in kleineren Gruppen vonO ( n log ( n ) ) n n k 2 ( k + 1 ) k k O ( n log ( k ) ) , O ( n ) 2 N ( 1 + 1 / k ) $n$$O(n\log(n))$$n$$n$$k$mit geringem Wirkungsgradverlust, da jede Gruppe nur Normalen benötigt, um einheitliche Werte zu erzeugen . Für festes ist der asymptotische Rechenaufwand daher = , wobei Normalvariablen erforderlich sind, um einheitliche Werte zu erzeugen .)$2(k+1)$$k$$k$$O(n\log(k))$$O(n)$$2n(1+1/k)$$n$

8. In hervorragender Näherung sieht jede normale Variation mit einer großen Standardabweichung über Bereiche mit viel kleineren Werten gleichmäßig aus. Wenn wir diese Verteilung in den Bereich rollen (indem wir nur die Bruchteile der Werte nehmen), erhalten wir eine Verteilung, die für alle praktischen Zwecke einheitlich ist. Dies ist äußerst effizient und erfordert eine der einfachsten Rechenoperationen von allen: Runden Sie einfach jede Normalvariable auf die nächste ganze Zahl ab und behalten Sie den Überschuss bei. Die Einfachheit dieses Ansatzes wird überzeugend, wenn wir eine praktische Implementierung untersuchen:$[0,1]$R

   rnorm(n, sd=10) %% 1
   

   erzeugt zuverlässig neinheitliche Werte im Bereich auf Kosten von nur normalen Variationen und fast keiner Berechnung.$[0,1]$n

   (Selbst wenn die Standardabweichung , weicht das PDF dieser Annäherung von einem einheitlichen PDF, wie in der folgenden Abbildung gezeigt, um weniger als einen Teil von . Um es zuverlässig zu erkennen, wäre eine Stichprobe von erforderlich values ​​- das übersteigt bereits die Möglichkeiten eines Standard-Zufallstests. Bei einer größeren Standardabweichung ist die Ungleichmäßigkeit so gering, dass sie nicht einmal berechnet werden kann. Zum Beispiel bei einer SD von wie im Code gezeigt, das Maximum Abweichung von einem einheitlichen PDF beträgt nur .)10 8 10 16 10 10 - $1$$10^8$$10^{16}$$10$$10^{-857}$

   

In jedem Fall können normale Variablen "mit bekannten Parametern" leicht neu zentriert und in die oben angenommenen Standardnormalen skaliert werden. Anschließend können die resultierenden gleichmäßig verteilten Werte neu zentriert und skaliert werden, um jedes gewünschte Intervall abzudecken. Diese erfordern nur Grundrechenarten.

Die Leichtigkeit dieser Konstruktionen wird durch den folgenden RCode verdeutlicht, der für die meisten nur ein oder zwei Zeilen verwendet. Ihre Richtigkeit wird durch die resultierenden nahezu einheitlichen Histogramme bestätigt, die jeweils auf unabhängigen Werten basieren (dies erfordert für alle sieben Simulationen etwa 12 Sekunden). Falls Sie sich Sorgen über das Ausmaß der Abweichungen in einem dieser Diagramme machen, finden Sie am Ende ein Histogramm mit einheitlichen Werten, die mit dem Generator für einheitliche Zufallszahlen simuliert wurden .$100,000$R

Alle diese Simulationen wurden unter Verwendung eines auf Behältern basierenden Tests auf Gleichförmigkeit getestet ; Keiner konnte als signifikant ungleichmäßig angesehen werden (der niedrigste p-Wert betrug für die Ergebnisse, die mit dem tatsächlichen Generator für gleichmäßige Zahlen generiert wurden !). 1000 3 $\chi^2$$1000$$3\%$R

set.seed(17)
n <- 1e5
y <- matrix(rnorm(floor(n/2)*2), nrow=2)
x <- c(atan2(y[2,], y[1,])/(2*pi) + 1/2, exp(-(y[1,]^2+y[2,]^2)/2))
hist(x, main="Box-Mueller")

y <- apply(array(rnorm(4*n), c(2,2,n)), c(3,2), function(z) sum(z^2))
x <- y[,2] / (y[,1]+y[,2])
hist(x, main="Beta")

x <- apply(array(rnorm(8*n), c(4,2,n)), 3, function(y) cor(y[,1], y[,2]))
hist(x, main="Correlation")

n.bits <- 32; x <-  (2^-(1:n.bits)) %*% matrix(rnorm(n*n.bits) > 0, n.bits)
hist(x, main="Binary")

y <- matrix(rnorm(n*3), 3)
x <- y[1, ] / sqrt(apply(y, 2, function(x) sum(x^2)))
hist(x, main="Equal area")

accept <- function(p) { # Using random normals, return TRUE with chance `p`
  p.bit <- x <- 0
  while(p.bit == x) {
    p.bit <- p >= 1/2
    x <- rnorm(1) >= 0
    p <- (2*p) %% 1
  }
  return(x == 0)
}
y <- rnorm(ceiling(n * sqrt(exp(1)*pi/2))) # This aims to produce `n` uniforms
y <- y[abs(y) < 1]
x <- y[sapply(y, function(x) accept(exp((x^2-1)/2)))]
hist(x, main="Rejection")

y <- matrix(rnorm(2*(n+1))^2, 2)
x <- cumsum(y)[seq(2, 2*(n+1), 2)]
x <- x[-(n+1)] / x[n+1]
x <- x[order(y[2,-(n+1)]/y[1,-(n+1)])] 
hist(x, main="Ordered")

x <- rnorm(n) %% 1 # (Use SD of 5 or greater in practice)
hist(x, main="Modular")

x <- runif(n)      # Reference distribution
hist(x, main="Uniform")

Sie können einen Trick verwenden, der dem sehr ähnlich ist, was Sie erwähnen. Angenommen, ist eine normale Zufallsvariable mit bekannten Parametern. Dann wissen wir, dass seine Verteilungsfunktion und auf gleichmäßig verteilt sind . Beachten Sie, dass wir dies für sehen, um dies zu beweisen≤ μ , σ 2 ≤ μ , σ 2 ( X ) ( 0 , 1 ) d ≤ ( 0 , 1 $X \sim N(\mu, \sigma^2)$$\Phi_{\mu,\sigma^2}$$\Phi_{\mu,\sigma^2}(X)$$(0,1)$$d \in (0,1)$

$P(\Phi_{\mu,\sigma^2}(X) \leq d) = P(X \leq \Phi_{\mu,\sigma^2}^{-1}(d)) = d$ .

Die obige Wahrscheinlichkeit ist eindeutig Null für nicht positives und für . Dies reicht aus, um zu zeigen, dass eine gleichmäßige Verteilung auf da wir gezeigt haben, dass die entsprechenden Maße für einen Generator der Borel Algebra gleich sind on . Auf diese Weise können Sie einfach die normalverteilten Daten mit der Verteilungsfunktion umwandeln und Sie erhalten gleichmäßig verteilte Daten.1 d ≥ 1 ≤ μ , σ 2 ( X ) ( 0 , 1 ) σ $d$$1$$d \geq 1$$\Phi_{\mu,\sigma^2}(X)$$(0,1)$$\sigma$$\mathbb{R}$