Konvertieren einer Gleichverteilung in eine Normalverteilung

106

Wie kann ich eine gleichmäßige Verteilung (wie die meisten Zufallszahlengeneratoren beispielsweise zwischen 0,0 und 1,0 erzeugen) in eine Normalverteilung umwandeln? Was ist, wenn ich einen Mittelwert und eine Standardabweichung meiner Wahl möchte?

Terhorst
quelle
3
Haben Sie eine Sprachspezifikation oder ist dies nur eine allgemeine Algorithmusfrage?
Bill the Lizard
3
Allgemeine Algorithmusfrage. Es ist mir egal, welche Sprache. Ich würde es jedoch vorziehen, wenn die Antwort nicht auf bestimmten Funktionen beruht, die nur diese Sprache bietet.
Terhorst

Antworten:

47

Der Ziggurat-Algorithmus ist dafür ziemlich effizient, obwohl die Box-Muller-Transformation einfacher von Grund auf neu zu implementieren ist (und nicht verrückt langsam).

Tyler
quelle
7
Die üblichen Warnungen zu linearen kongruenten Generatoren gelten für beide Methoden. Verwenden Sie daher einen anständigen Underling-Generator. Prost.
dmckee --- Ex-Moderator Kätzchen
3
Wie Mersenee Twister, oder haben Sie andere Vorschläge?
Gregg Lind
47

Es gibt viele Methoden:

  • Sie nicht verwenden Box Muller. Besonders wenn Sie viele Gaußsche Zahlen zeichnen. Box Muller liefert ein Ergebnis, das zwischen -6 und 6 eingeklemmt ist (unter der Annahme einer doppelten Präzision. Mit Schwimmern verschlechtern sich die Dinge.). Und es ist wirklich weniger effizient als andere verfügbare Methoden.
  • Ziggurat ist in Ordnung, benötigt jedoch eine Tabellensuche (und einige plattformspezifische Optimierungen aufgrund von Problemen mit der Cache-Größe).
  • Das Verhältnis der Uniformen ist mein Favorit, nur ein paar Additionen / Multiplikationen und ein Log 1/50 der Zeit (z. B. dort schauen ).
  • Das Invertieren der CDF ist effizient (und wird übersehen, warum?). Wenn Sie Google durchsuchen, stehen Ihnen schnelle Implementierungen zur Verfügung. Dies ist für Quasi-Zufallszahlen obligatorisch.
Alexandre C.
quelle
2
Sind Sie sicher über die [-6,6] Klemmung? Dies ist ein ziemlich wichtiger Punkt, wenn er wahr ist (und verdient einen Hinweis auf der Wikipedia-Seite).
Redcalx
1
@locster: Das hat mir ein Lehrer gesagt (er hat solche Generatoren studiert und ich vertraue seinem Wort). Möglicherweise kann ich Ihnen eine Referenz finden.
Alexandre C.
7
@locster: Diese unerwünschte Eigenschaft wird auch von der inversen CDF-Methode gemeinsam genutzt. Siehe cimat.mx/~src/prope08/randomgauss.pdf . Dies kann durch Verwendung eines einheitlichen RNG gemildert werden, das eine Wahrscheinlichkeit ungleich Null aufweist, um eine Gleitkommazahl sehr nahe Null zu ergeben. Die meisten RNG tun dies nicht, da sie eine (normalerweise 64-Bit-) Ganzzahl erzeugen, die dann auf [0,1] abgebildet wird. Dies macht diese Methoden ungeeignet, um Schwänze von Gaußschen Variablen abzutasten (denken Sie an die Preisgestaltung für Optionen mit niedrigem / hohem Streik in der Computerfinanzierung).
Alexandre C.
6
@AlexandreC. Um in zwei Punkten klar zu sein: Bei Verwendung von 64-Bit-Zahlen gehen die Schwänze entweder auf 8,57 oder 9,41 (der niedrigere Wert entspricht der Konvertierung in [0,1), bevor das Protokoll erstellt wird). Selbst wenn auf [-6, 6] geklemmt, liegen die Chancen, außerhalb dieses Bereichs zu liegen, bei etwa 1,98e-9, was für die meisten Menschen selbst in der Wissenschaft gut genug ist. Für die Zahlen 8,57 und 9,41 werden dies 1,04e-17 und 4,97e-21. Diese Zahlen sind so gering, dass der Unterschied zwischen einer Box-Muller-Stichprobe und einer echten Gauß-Stichprobe in Bezug auf diese Grenze fast rein akademisch ist. Wenn Sie etwas Besseres brauchen, addieren Sie einfach vier davon und dividieren Sie durch 2.
CrazyCasta
6
Ich denke, der Vorschlag, die Box Muller-Transformation nicht zu verwenden, ist für einen großen Prozentsatz der Benutzer irreführend. Es ist großartig, über die Einschränkung Bescheid zu wissen, aber wie CrazyCasta betont, müssen Sie sich für die meisten Anwendungen, die nicht stark von Ausreißern abhängig sind, wahrscheinlich keine Sorgen machen. Wenn Sie beispielsweise jemals von der Abtastung einer Normalen mit Numpy abhängig waren, haben Sie sich auf die Box-Muller-Transformation (Polarkoordinatenform) github.com/numpy/numpy/blob/… verlassen .
Andreas Grivas
29

Um die Verteilung einer Funktion auf eine andere zu ändern, müssen Sie die Umkehrung der gewünschten Funktion verwenden.

Mit anderen Worten, wenn Sie eine bestimmte Wahrscheinlichkeitsfunktion p (x) anstreben, erhalten Sie die Verteilung, indem Sie darüber integrieren -> d (x) = Integral (p (x)) und deren Umkehrung verwenden: Inv (d (x)) . Verwenden Sie nun die Zufallswahrscheinlichkeitsfunktion (die gleichmäßig verteilt ist) und wandeln Sie den Ergebniswert durch die Funktion Inv (d (x)). Sie sollten zufällige Werte erhalten, die entsprechend der von Ihnen gewählten Funktion verteilt werden.

Dies ist der generische mathematische Ansatz. Wenn Sie ihn verwenden, können Sie jetzt jede Wahrscheinlichkeits- oder Verteilungsfunktion auswählen, die Sie haben, solange sie eine inverse oder gute inverse Approximation hat.

Hoffe das hat geholfen und danke für die kleine Bemerkung zur Verwendung der Distribution und nicht der Wahrscheinlichkeit selbst.

Adi
quelle
4
+1 Dies ist eine übersehene Methode zum Generieren von Gaußschen Variablen, die sehr gut funktioniert. Inverse CDF kann in diesem Fall effizient mit der Newton-Methode berechnet werden (Ableitung ist e ^ {- t ^ 2}). Eine anfängliche Näherung ist als rationaler Bruch leicht zu erhalten, sodass Sie 3-4 Auswertungen von erf und exp benötigen. Es ist obligatorisch, wenn Sie Quasi-Zufallszahlen verwenden. In diesem Fall müssen Sie genau eine einheitliche Zahl verwenden, um eine Gaußsche Zahl zu erhalten.
Alexandre C.
9
Beachten Sie, dass Sie die kumulative Verteilungsfunktion und nicht die Wahrscheinlichkeitsverteilungsfunktion invertieren müssen. Alexandre impliziert dies, aber ich dachte, es könnte nicht schaden, es expliziter zu erwähnen - da die Antwort das PDF zu suggerieren scheint
ltjax
Sie können das PDF verwenden, wenn Sie bereit sind, eine Richtung relativ zum Mittelwert zufällig auszuwählen. verstehe ich das richtig
Mark McKenna
2
Dies wird als inverse Transformationsabtastung
Bindestrich bezeichnet
Hier ist eine verwandte Frage in SE mit einer allgemeineren Antwort mit einer schönen Erklärung.
Bindestrich
23

Hier ist eine Javascript-Implementierung unter Verwendung der polaren Form der Box-Muller-Transformation.

/*
 * Returns member of set with a given mean and standard deviation
 * mean: mean
 * standard deviation: std_dev 
 */
function createMemberInNormalDistribution(mean,std_dev){
    return mean + (gaussRandom()*std_dev);
}

/*
 * Returns random number in normal distribution centering on 0.
 * ~95% of numbers returned should fall between -2 and 2
 * ie within two standard deviations
 */
function gaussRandom() {
    var u = 2*Math.random()-1;
    var v = 2*Math.random()-1;
    var r = u*u + v*v;
    /*if outside interval [0,1] start over*/
    if(r == 0 || r >= 1) return gaussRandom();

    var c = Math.sqrt(-2*Math.log(r)/r);
    return u*c;

    /* todo: optimize this algorithm by caching (v*c) 
     * and returning next time gaussRandom() is called.
     * left out for simplicity */
}
user5084
quelle
5

Verwenden Sie den zentralen Grenzwertsatz Wikipedia-Eintrag Mathworld-Eintrag zu Ihrem Vorteil.

Generieren Sie n der gleichmäßig verteilten Zahlen, summieren Sie sie, subtrahieren Sie n * 0,5 und Sie haben die Ausgabe einer ungefähr normalen Verteilung mit einem Mittelwert von 0 und einer Varianz von (1/12) * (1/sqrt(N))(siehe Wikipedia zu gleichmäßigen Verteilungen für die letzte).

n = 10 gibt dir schnell etwas halbwegs Anständiges. Wenn Sie etwas mehr als halbwegs Anständiges wollen, entscheiden Sie sich für eine Tyler-Lösung (wie im Wikipedia-Eintrag zu Normalverteilungen angegeben ).

jilles de wit
quelle
1
Dies ergibt keine besonders enge Normalität (die "Schwänze" oder Endpunkte liegen nicht nahe an der tatsächlichen Normalverteilung). Box-Muller ist besser, wie andere vorgeschlagen haben.
Peter K.
1
Box Muller hat auch falsche Schwänze (es gibt eine Zahl zwischen -6 und 6 in doppelter Genauigkeit zurück)
Alexandre C.
n = 12 (summiere 12 Zufallszahlen im Bereich von 0 bis 1 und subtrahiere 6) ergibt stddev = 1 und mean = 0. Dies kann dann verwendet werden, um eine beliebige Normalverteilung zu erzeugen. Multiplizieren Sie einfach das Ergebnis mit dem gewünschten Standard und addieren Sie den Mittelwert.
JerryM
2

Ich würde Box-Muller benutzen. Zwei Dinge dazu:

  1. Am Ende erhalten Sie zwei Werte pro Iteration. In der
    Regel wird ein Wert zwischengespeichert und der andere zurückgegeben. Beim nächsten Aufruf eines Beispiels geben Sie den zwischengespeicherten Wert zurück.
  2. Box-Muller gibt einen Z-Score an
    Sie müssen dann den Z-Score um die Standardabweichung skalieren und den Mittelwert addieren, um den vollen Wert in der Normalverteilung zu erhalten.
hughdbrown
quelle
Wie skaliert man den Z-Score?
Terhorst
2
skaliert = Mittelwert + stdDev * zScore // gibt Ihnen normal (Mittelwert, stdDev ^ 2)
yoyoyoyosef
2

Wobei R1, R2 zufällige einheitliche Zahlen sind:

NORMALE VERTEILUNG mit SD von 1: sqrt (-2 * log (R1)) * cos (2 * pi * R2)

Das ist genau ... Sie müssen nicht all diese langsamen Schleifen machen!

Erik Aronesty
quelle
Bevor mich jemand korrigierte ... hier ist die Annäherung, die ich mir ausgedacht habe: (1,5- (R1 + R2 + R3)) * 1,88. Ich mag es auch.
Erik Aronesty
2

Es scheint unglaublich, dass ich nach acht Jahren noch etwas hinzufügen könnte, aber für den Fall von Java möchte ich die Leser auf die Random.nextGaussian () -Methode verweisen , die für Sie eine Gaußsche Verteilung mit einem Mittelwert von 0,0 und einer Standardabweichung von 1,0 generiert.

Eine einfache Addition und / oder Multiplikation ändert den Mittelwert und die Standardabweichung an Ihre Bedürfnisse.

Pepijn Schmitz
quelle
1

Dies ist meine JavaScript-Implementierung von Algorithmus P ( Polar-Methode für normale Abweichungen ) aus Abschnitt 3.4.1 von Donald Knuths Buch The Art of Computer Programming :

function normal_random(mean,stddev)
{
    var V1
    var V2
    var S
    do{
        var U1 = Math.random() // return uniform distributed in [0,1[
        var U2 = Math.random()
        V1 = 2*U1-1
        V2 = 2*U2-1
        S = V1*V1+V2*V2
    }while(S >= 1)
    if(S===0) return 0
    return mean+stddev*(V1*Math.sqrt(-2*Math.log(S)/S))
}
Alessandro Jacopson
quelle
0

Das Standard-Python-Bibliotheksmodul random hat das, was Sie wollen:

Normalvariable (mu, Sigma)
Normalverteilung. mu ist der Mittelwert und Sigma ist die Standardabweichung.

Schauen Sie sich für den Algorithmus selbst die Funktion in random.py in der Python-Bibliothek an.

Der manuelle Eintrag ist hier

Brent.Longborough
quelle
1
Leider verwendet die Python-Bibliothek Kinderman, AJ und Monahan, JF, "Computergenerierung von Zufallsvariablen unter Verwendung des Verhältnisses einheitlicher Abweichungen", ACM Trans Math Software, 3, (1977), S. 257-260. Dabei werden zwei einheitliche Zufallsvariablen verwendet, um den Normalwert zu generieren, und nicht nur eine einzige. Daher ist es nicht offensichtlich, wie dieser als das vom OP gewünschte Mapping verwendet werden soll.
Ian
0

Ich denke, Sie sollten dies in EXCEL versuchen : =norminv(rand();0;1). Dies ergibt die Zufallszahlen, die normalerweise mit dem Mittelwert Null verteilt werden sollen, und vereint die Varianz. "0" kann mit einem beliebigen Wert angegeben werden, so dass die Zahlen den gewünschten Mittelwert haben. Durch Ändern von "1" erhalten Sie die Varianz, die dem Quadrat Ihrer Eingabe entspricht.

Zum Beispiel: =norminv(rand();50;3)ergibt die normalverteilten Zahlen mit MEAN = 50 VARIANCE = 9.

Nilpferd
quelle
0

F Wie kann ich eine Gleichverteilung (wie die meisten Zufallszahlengeneratoren, z. B. zwischen 0,0 und 1,0) in eine Normalverteilung umwandeln?

  1. Für die Software-Implementierung kenne ich einige zufällige Generatornamen, die Ihnen eine pseudo-einheitliche zufällige Sequenz in [0,1] geben (Mersenne Twister, Linear Congruate Generator). Nennen wir es U (x)

  2. Es gibt einen mathematischen Bereich, der Wahrscheinlichkeitstheorie genannt wird. Das erste: Wenn Sie rv mit der Integralverteilung F modellieren möchten, können Sie versuchen, nur F ^ -1 (U (x)) auszuwerten. In der Theorie wurde bewiesen, dass ein solches RV eine integrale Verteilung F haben wird.

  3. Schritt 2 kann angewendet werden, um rv ~ F ohne Verwendung von Zählmethoden zu erzeugen, wenn F ^ -1 ohne Probleme analytisch abgeleitet werden kann. (zB exp.distribution)

  4. Um die Normalverteilung zu modellieren, können Sie y1 * cos (y2) berechnen, wobei y1 ~ in [0,2pi] einheitlich ist. und y2 ist die relei Verteilung.

F: Was ist, wenn ich einen Mittelwert und eine Standardabweichung meiner Wahl möchte?

Sie können Sigma * N (0,1) + m berechnen.

Es kann gezeigt werden, dass eine solche Verschiebung und Skalierung zu N (m, Sigma) führt.

Bruziuz
quelle
0

Dies ist eine Matlab-Implementierung, die die polare Form der Box-Muller- Transformation verwendet:

Funktion randn_box_muller.m:

function [values] = randn_box_muller(n, mean, std_dev)
    if nargin == 1
       mean = 0;
       std_dev = 1;
    end

    r = gaussRandomN(n);
    values = r.*std_dev - mean;
end

function [values] = gaussRandomN(n)
    [u, v, r] = gaussRandomNValid(n);

    c = sqrt(-2*log(r)./r);
    values = u.*c;
end

function [u, v, r] = gaussRandomNValid(n)
    r = zeros(n, 1);
    u = zeros(n, 1);
    v = zeros(n, 1);

    filter = r==0 | r>=1;

    % if outside interval [0,1] start over
    while n ~= 0
        u(filter) = 2*rand(n, 1)-1;
        v(filter) = 2*rand(n, 1)-1;
        r(filter) = u(filter).*u(filter) + v(filter).*v(filter);

        filter = r==0 | r>=1;
        n = size(r(filter),1);
    end
end

Und dies aufzurufen histfit(randn_box_muller(10000000),100);ist das Ergebnis: Box-Muller Matlab Histfit

Offensichtlich ist es im Vergleich zum eingebauten Matlab- Randn wirklich ineffizient .

Madx
quelle
0

Ich habe den folgenden Code, der vielleicht helfen könnte:

set.seed(123)
n <- 1000
u <- runif(n) #creates U
x <- -log(u)
y <- runif(n, max=u*sqrt((2*exp(1))/pi)) #create Y
z <- ifelse (y < dnorm(x)/2, -x, NA)
z <- ifelse ((y > dnorm(x)/2) & (y < dnorm(x)), x, z)
z <- z[!is.na(z)]
Große Köpfe denken ähnlich
quelle
0

Es ist auch einfacher, die implementierte Funktion rnorm () zu verwenden, da sie schneller ist als das Schreiben eines Zufallszahlengenerators für die Normalverteilung. Siehe den folgenden Code als Beweis

n <- length(z)
t0 <- Sys.time()
z <- rnorm(n)
t1 <- Sys.time()
t1-t0
peterweethetbeter
quelle
-2
function distRandom(){
  do{
    x=random(DISTRIBUTION_DOMAIN);
  }while(random(DISTRIBUTION_RANGE)>=distributionFunction(x));
  return x;
}

quelle
Es ist jedoch nicht garantiert, dass wir zurückkehren, oder? ;-)
Peter K.
5
Zufallszahlen sind zu wichtig, um sie dem Zufall zu überlassen.
Drew Noakes
Beantwortet die Frage nicht - die Normalverteilung hat eine unendliche Domäne.
Matt