Konwersja rozkładu jednorodnego na rozkład normalny

106

Jak mogę przekształcić rozkład równomierny (jak generuje większość generatorów liczb losowych, np. Między 0,0 a 1,0) na rozkład normalny? A jeśli chcę mieć wybraną średnią i odchylenie standardowe?

Terhorst
źródło
3
Czy masz specyfikację języka, czy jest to tylko ogólne pytanie dotyczące algorytmu?
Bill the Lizard,
3
Ogólne pytanie o algorytm. Nie obchodzi mnie, który język. Wolałbym jednak, aby odpowiedź nie opierała się na konkretnych funkcjach, które zapewnia tylko ten język.
Terhorst,

Odpowiedzi:

47

Algorytm Ziggurat jest dość skuteczny w tym, choć Transformacja Boxa-Mullera jest łatwiejszy do wdrożenia od zera (a nie szalone powolny).

Tyler
źródło
7
Zwykłe ostrzeżenia dotyczące generatorów przystających liniowych dotyczą obu tych metod, więc użyj przyzwoitego generatora podrzędnego. Twoje zdrowie.
dmckee --- kociak byłego moderatora
3
Na przykład Mersenee Twister, czy masz inne sugestie?
Gregg Lind
47

Istnieje wiele metod:

  • Czy nie używać Box Muller. Zwłaszcza jeśli narysujesz wiele liczb gaussowskich. Box Muller daje wynik, który jest zaciskany między -6 a 6 (zakładając podwójną precyzję. Sytuacja pogarsza się w przypadku pływaków). I jest naprawdę mniej skuteczny niż inne dostępne metody.
  • Ziggurat jest w porządku, ale wymaga wyszukiwania w tabeli (i pewnych poprawek specyficznych dla platformy ze względu na problemy z rozmiarem pamięci podręcznej)
  • Moim ulubionym jest Ratio-of-uniforms, tylko kilka dodawania / mnożenia i log 1/50 czasu (np. Spójrz tam ).
  • Odwrócenie CDF jest wydajne (i przeoczone, dlaczego?), Masz dostępne szybkie implementacje, jeśli przeszukujesz google. Jest to obowiązkowe w przypadku liczb quasi-losowych.
Alexandre C.
źródło
2
Czy jesteś pewien co do mocowania [-6,6]? To dość znacząca kwestia, jeśli jest prawdziwa (i warta uwagi na stronie wikipedii).
redcalx
1
@locster: tak powiedział mi mój nauczyciel (studiował takie generatory i ufam jego słowu). Może znajdę dla ciebie odniesienie.
Alexandre C.,
7
@locster: ta niepożądana właściwość jest również współdzielona przez odwrotną metodę CDF. Zobacz cimat.mx/~src/prope08/randomgauss.pdf . Można to złagodzić stosując jednolity RNG, który ma niezerowe prawdopodobieństwo uzyskania liczby zmiennoprzecinkowej bardzo bliskiej zeru. Większość RNG tego nie robi, ponieważ generują (zazwyczaj 64-bitową) liczbę całkowitą, która jest następnie odwzorowywana na [0,1]. To sprawia, że ​​metody te nie nadają się do próbkowania ogonów zmiennych gaussowskich (pomyśl o wycenie opcji o niskim / wysokim strajku w finansach obliczeniowych).
Alexandre C.,
6
@AlexandreC. Żeby wyjaśnić dwa punkty, używając liczb 64-bitowych, ogony wychodzą do 8,57 lub 9,41 (niższa wartość odpowiada konwersji na [0,1) przed zrobieniem log). Nawet przy ograniczeniu do [-6, 6] szanse znalezienia się poza tym zakresem wynoszą około 1,98e-9, co jest wystarczające dla większości ludzi, nawet w nauce. Dla liczb 8,57 i 9,41 jest to 1,04e-17 i 4,97e-21. Liczby te są tak małe, że różnica między próbkowaniem Box Mullera a prawdziwym próbkowaniem gaussowskim pod względem wspomnianego limitu jest prawie czysto akademicka. Jeśli potrzebujesz
czegoś
6
Myślę, że sugestia, aby nie używać transformacji Boxa Müllera, jest myląca dla dużego odsetka użytkowników. Dobrze jest wiedzieć o ograniczeniu, ale jak wskazuje CrazyCasta, w przypadku większości aplikacji, które nie są w dużym stopniu zależne od wartości odstających, prawdopodobnie nie musisz się tym martwić. Na przykład, jeśli kiedykolwiek polegałeś na próbkowaniu z normalnego przy użyciu numpy, polegałeś na transformacji Boxa Mullera (forma współrzędnych biegunowych) github.com/numpy/numpy/blob/… .
Andreas Grivas
30

Zmiana rozkładu dowolnej funkcji na inną wymaga użycia odwrotności żądanej funkcji.

Innymi słowy, jeśli dążysz do określonej funkcji prawdopodobieństwa p (x), otrzymasz rozkład przez całkowanie po niej -> d (x) = całka (p (x)) i użycie jej odwrotności: Inv (d (x)) . Teraz użyj funkcji prawdopodobieństwa losowego (które mają rozkład równomierny) i rzuć wartość wyniku za pomocą funkcji Inv (d (x)). Powinieneś otrzymać losowe wartości rzutowane z rozkładem zgodnie z wybraną funkcją.

To jest ogólne podejście matematyczne - używając go możesz teraz wybrać dowolną funkcję prawdopodobieństwa lub rozkładu, o ile ma ona odwrotne lub dobre odwrotne przybliżenie.

Mam nadzieję, że to pomogło i dziękuję za małą uwagę na temat korzystania z rozkładu, a nie samego prawdopodobieństwa.

Adi
źródło
4
+1 Jest to przeoczona metoda generowania zmiennych gaussowskich, która działa bardzo dobrze. Odwrotność CDF można w tym przypadku efektywnie obliczyć metodą Newtona (pochodna to e ^ {- t ^ 2}), początkowe przybliżenie jest łatwe do uzyskania jako ułamek wymierny, więc potrzebne są 3-4 oceny erf i exp. Jest to obowiązkowe, jeśli używasz liczb quasi-losowych, przypadek, w którym musisz użyć dokładnie jednej liczby jednolitej, aby uzyskać liczbę gaussowską.
Alexandre C.
9
Zwróć uwagę, że musisz odwrócić dystrybuantę, a nie rozkład prawdopodobieństwa. Alexandre sugeruje to, ale pomyślałem, że wspomnienie o tym bardziej otwarcie może nie zaszkodzić - ponieważ odpowiedź wydaje się sugerować plik PDF
ltjax
Możesz użyć pliku PDF, jeśli jesteś przygotowany do losowego wyboru kierunku względem średniej; rozumiem to, prawda?
Mark McKenna
2
Nazywa się to próbkowaniem z odwrotną transformacją
dashesy
1
Tutaj jest powiązane pytanie w SE z bardziej uogólnioną odpowiedzią z ładnym wyjaśnieniem.
dashesy
23

Oto implementacja javascript wykorzystująca polarną postać transformacji Boxa-Mullera.

/*
 * 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
źródło
5

Użyj centralnego twierdzenia granicznego wpisu wikipedii o mathworld na swoją korzyść.

Wygeneruj n równomiernie rozłożonych liczb, zsumuj je, odejmij n * 0,5 i otrzymasz wynik w przybliżeniu normalnego rozkładu ze średnią równą 0 i wariancją równą (1/12) * (1/sqrt(N))(patrz wikipedia o rozkładach jednorodnych dla tego ostatniego)

n = 10 daje coś w połowie przyzwoitego szybko. Jeśli chcesz czegoś więcej niż w połowie przyzwoitego, wybierz rozwiązanie Tylers (jak wspomniano we wpisie Wikipedii o normalnych dystrybucjach )

jilles de wit
źródło
1
Nie da to szczególnie bliskiej normy („ogony” lub punkty końcowe nie będą zbliżone do rzeczywistego rozkładu normalnego). Box-Muller jest lepszy, jak sugerowali inni.
Peter K.
1
Box Muller też ma złe ogony (zwraca liczbę od -6 do 6 z podwójną precyzją)
Alexandre C.
n = 12 (suma 12 liczb losowych z zakresu od 0 do 1 i odjęcie 6) daje odchylenie standardowe = 1 i średnią = 0. Można to następnie wykorzystać do wygenerowania dowolnego rozkładu normalnego. Po prostu pomnóż wynik przez żądane odchylenie standardowe i dodaj średnią.
JerryM,
3

Użyłbym Box-Mullera. Dwie rzeczy na ten temat:

  1. W rezultacie otrzymujesz dwie wartości na iterację.
    Zwykle buforujesz jedną wartość, a zwracasz drugą. Przy następnym wywołaniu próbki zwracasz zbuforowaną wartość.
  2. Box-Muller podaje wynik Z.
    Następnie należy wyskalować wynik Z za pomocą odchylenia standardowego i dodać średnią, aby uzyskać pełną wartość w rozkładzie normalnym.
hughdbrown
źródło
Jak skalujesz wynik Z?
Terhorst,
3
scaled = mean + stdDev * zScore // daje normalne (średnia, stdDev ^ 2)
yoyoyoyosef
2

Gdzie R1, R2 to losowe liczby jednolite:

ROZKŁAD NORMALNY, ze SD równym 1: sqrt (-2 * log (R1)) * cos (2 * pi * R2)

To jest dokładne ... nie musisz robić tych wszystkich wolnych pętli!

Erik Aronesty
źródło
Zanim ktoś mnie poprawił ... oto przybliżenie, które wymyśliłem: (1,5- (R1 + R2 + R3)) * 1,88. Też to lubię.
Erik Aronesty,
2

Wydaje się niewiarygodne, że mogłem coś do tego dodać po ośmiu latach, ale w przypadku Javy chciałbym zwrócić czytelnikom uwagę na metodę Random.nextGaussian () , która generuje rozkład Gaussa ze średnią 0,0 i odchyleniem standardowym 1,0.

Proste dodawanie i / lub mnożenie zmieni średnią i odchylenie standardowe zgodnie z Twoimi potrzebami.

Pepijn Schmitz
źródło
1

Standardowy moduł losowy biblioteki Pythona ma to, czego chcesz:

normalvariate (mu, sigma)
Rozkład normalny. mu to średnia, a sigma to odchylenie standardowe.

Jeśli chodzi o sam algorytm, spójrz na funkcję w random.py w bibliotece Pythona.

Ręczne wprowadzanie jest tutaj

Brent, Longborough
źródło
2
Niestety, biblioteka Pythona wykorzystuje Kinderman, AJ i Monahan, JF, "Computer generation of random variable using the ratio of uniform deviates", ACM Trans Math Software, 3, (1977), str. 257-260. Wykorzystuje to dwie jednolite zmienne losowe do wygenerowania wartości normalnej, a nie jedną, więc nie jest oczywiste, jak użyć jej jako odwzorowania, którego chciał PO.
Ian
1

Oto moja implementacja algorytmu P ( metoda biegunowa dla odchyleń normalnych ) z sekcji 3.4.1 książki Donalda Knutha 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
źródło
0

Myślę, że powinieneś spróbować tego w EXCEL: =norminv(rand();0;1) . Spowoduje to iloczyn liczb losowych, które powinny mieć rozkład normalny ze średnią zerową i jednoczącą wariancję. „0” można podać dowolną wartość, dzięki czemu liczby będą miały pożądaną średnią, a zmieniając „1”, uzyskasz wariancję równą kwadratowi wprowadzonego przez Ciebie tekstu.

Na przykład: =norminv(rand();50;3)ustąpi liczbom o rozkładzie normalnym z ŚREDNIA = 50 ODMIANA = 9.

Hipopotam
źródło
0

P Jak mogę przekształcić rozkład równomierny (jak generuje większość generatorów liczb losowych, np. Między 0,0 a 1,0) na rozkład normalny?

  1. Do implementacji oprogramowania znam kilka losowych nazw generatorów, które dają pseudojednorodną losową sekwencję w [0,1] (Mersenne Twister, Linear Congruate Generator). Nazwijmy to U (x)

  2. Istnieje obszar matematyczny, który nazywa się teorią prawdopodobieństwa. Pierwsza rzecz: jeśli chcesz modelować rv z rozkładem całkowym F, możesz spróbować po prostu obliczyć F ^ -1 (U (x)). W teorii pr udowodniono, że taki rv będzie miał rozkład całkowy F.

  3. Krok 2 można zastosować do wygenerowania rv ~ F bez użycia jakichkolwiek metod zliczania, gdy F ^ -1 można wyprowadzić analitycznie bez problemów. (np. dystrybucja eksp.)

  4. Aby zamodelować rozkład normalny, można obliczyć y1 * cos (y2), gdzie y1 ~ jest jednorodne w [0,2pi]. a y2 to dystrybucja releasei.

P: A jeśli chcę mieć wybrane średnie i odchylenie standardowe?

Możesz obliczyć sigma * N (0,1) + m.

Można wykazać, że takie przesunięcie i skalowanie prowadzą do N (m, sigma)

bruziuz
źródło
0

To jest implementacja Matlaba wykorzystująca polarną postać transformacji Boxa-Mullera :

Funkcja 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

A wywołanie histfit(randn_box_muller(10000000),100);tego jest wynikiem: Box-Muller Matlab Histfit

Oczywiście jest to naprawdę nieefektywne w porównaniu z randn wbudowanym w Matlab .

madx
źródło
0

Mam następujący kod, który może pomóc:

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)]
wielkie umysły myślą podobnie
źródło
0

Użycie zaimplementowanej funkcji rnorm () jest również łatwiejsze, ponieważ jest szybsze niż pisanie generatora liczb losowych dla rozkładu normalnego. Zobacz poniższy kod jako dowód

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

źródło
Nie ma jednak gwarancji powrotu, prawda? ;-)
Peter K.
5
Liczby losowe są zbyt ważne, aby pozostawić je przypadkowi.
Drew Noakes
Nie odpowiada na pytanie - normalna dystrybucja ma nieskończoną domenę.
Matt,