Konwersja rozkładu jednorodnego na rozkład normalny

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?

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

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.

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.

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 */
}

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 )

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.

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!

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.

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

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))
}

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.

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)

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:

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

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)]

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

function distRandom(){
  do{
    x=random(DISTRIBUTION_DOMAIN);
  }while(random(DISTRIBUTION_RANGE)>=distributionFunction(x));
  return x;
}