Jak pobierać próbki z rozkładu normalnego ze znaną średnią i wariancją przy użyciu konwencjonalnego języka programowania?

36

Nigdy nie miałem kursu statystyki, więc mam nadzieję, że pytam w odpowiednim miejscu.

Załóżmy, że mam tylko dwa dane opisujące rozkład normalny: średnią i wariancję . Chcę użyć komputera do losowego pobierania próbek z tej dystrybucji, tak aby uszanować te dwie statystyki. $\mu$ $\sigma^2$

To całkiem oczywiste, że mogę poradzić sobie ze średnią, po prostu normalizując około 0: po prostu dodaj do każdej próbki przed wysłaniem próbki. Ale nie widzę, jak programowo generować próbki, aby uszanować . $\mu$ $\sigma^2$

Mój program będzie w konwencjonalnym języku programowania; Nie mam dostępu do żadnych pakietów statystycznych.

normal-distribution sampling computing

— Fixee
źródło

Czy twój język ma generator liczb losowych? Czy ten generator ma tylko rozkład jednolity, czy może też generować z rozkładu normalnego?

— ttnphns

@ttnphns: Niemal każdy język komputerowy ma generator liczb losowych. Są to w przeważającej mierze jednolite generatory w niektórych skończonych domenach.

— Fixee

33

Jeśli możesz próbkować z danego rozkładu ze średnią 0 i wariancją 1, możesz łatwo próbować z transformacji położenia w skali tego rozkładu, który ma średnią i wariancję . Jeśli jest próbką ze średniego rozkładu 0 i wariancji 1, to jest próbką o średniej i wariancji . Wystarczy więc przeskalować zmienną o odchylenie standardowe (pierwiastek kwadratowy wariancji) przed dodaniem średniej . $\mu$ $\sigma^2$ $x$

σ x + μ

$\sigma x + \mu$

μ

$\mu$

σ^{2}

$\sigma^2$

σ

$\sigma$

μ

$\mu$

Sposób uzyskania symulacji z rozkładu normalnego ze średnią 0 i wariancją 1 to inna historia. To zabawne i interesujące wiedzieć, jak zaimplementować takie rzeczy, ale niezależnie od tego, czy korzystasz z pakietu statystycznego, czy języka programowania, zalecam uzyskanie i użycie odpowiedniej funkcji lub biblioteki do generowania liczb losowych. Jeśli potrzebujesz porady na temat używanej biblioteki, możesz dodać szczegółowe informacje na temat używanych języków programowania.

Edycja: W świetle komentarzy, kilku innych odpowiedzi oraz faktu, że Fixee zaakceptował tę odpowiedź, podam więcej szczegółów na temat tego, w jaki sposób można wykorzystać transformacje zmiennych jednolitych do uzyskania normalnych zmiennych.

Jedną z metod, wspomnianą już w komentarzu VitalStatistix , jest metoda Boxa-Mullera, która przyjmuje dwie niezależne jednolite zmienne losowe i wytwarza dwie niezależne normalne zmienne losowe. Podobna metoda, która pozwala uniknąć obliczeń dwóch funkcji transcendentalnych sin i cos kosztem kilku kolejnych symulacji, została opublikowana jako odpowiedź przez francogrex .
Całkowicie ogólną metodą jest transformacja jednolitej zmiennej losowej za pomocą funkcji odwrotnego rozkładu. Jeśli jest równomiernie rozłożone na wówczas ma standardowy rozkład normalny. Chociaż nie ma jednoznacznego wzoru analitycznego dla , można go obliczyć na podstawie dokładnych przybliżeń liczbowych. Obecna implementacja w R (ostatnio sprawdziłem) korzysta z tego pomysłu. Metoda jest koncepcyjnie bardzo prosta, ale wymaga dokładnej implementacji , co prawdopodobnie nie jest tak rozpowszechnione, jak (inne) funkcje transcendentalne $U$ $[0,1]$ $Φ^{- 1} (U)$ $\Phi^{-1}(U)$ $\Phi^{-1}$ $\Phi^{-1}$ log , grzech i cos .
Kilka odpowiedzi wspomina o możliwości zastosowania centralnego twierdzenia granicznego do przybliżenia rozkładu normalnego jako średniej jednolitych zmiennych losowych. Zasadniczo nie jest to zalecane. Przedstawione argumenty, takie jak dopasowanie średniej 0 i wariancji 1, oraz względy poparcia rozkładu nie są przekonujące. W ćwiczeniu 2.3 w „Wprowadzaniu metod Monte Carlo z R.” autorstwa Christiana P. Roberta i George'a Caselli generator ten nazywa się przestarzały, a przybliżenie nazywa się bardzo słabym .
Istnieje wiele innych pomysłów. Rozdział 3, a w szczególności sekcja 3.4, w „The Art of Computer Programming” Vol. 2 autorstwa Donalda E. Knutha to klasyczne odniesienie do generowania liczb losowych. Brian Ripley napisał Computer Generation of Random Variables: A Tutorial , który może być przydatny. Zalecana jest także książka wspomniana przez Roberta i Casellę, a może rozdział 2 w innej książce „Metody statystyczne Monte Carlo”.

Na koniec dnia poprawnie zaimplementowana metoda nie jest lepsza niż zastosowany jednolity generator liczb pseudolosowych. Osobiście wolę polegać na bibliotekach specjalnego przeznaczenia, które moim zdaniem są godne zaufania. Niemal zawsze polegam na metodach zaimplementowanych w języku R albo bezpośrednio w języku R, albo poprzez interfejs API w języku C / C ++. Oczywiście nie jest to rozwiązanie dla wszystkich, ale nie znam wystarczająco wielu innych bibliotek, aby polecać alternatywy.

— NRH
źródło

(+1) Dobra odpowiedź i porady dla PO.

— kardynał

18

Nie jestem pewien, czy robię tutaj niepotrzebny komentarz, ale jeśli masz dostęp tylko do jednolitego generatora liczb losowych, możesz użyć transformacji Boxa-Mullera do wygenerowania niezależnych liczb losowych N (0,1). W skrócie, jeśli U_1 i U_2 są niezależnymi losowaniami z rozkładu Uniform (0,1), to

i

\sqrt{- 2) \log (U_{1})} sałata (2) π U_{2)})

$\sqrt{-2 \log(U_1) } \cos(2\pi U_2)$

są rozdzielone jako niezależne zmienne losowe N (0,1). Podstawowy pomysł

\sqrt{- 2) \log (U_{1})} grzech (2) π U_{2)})

$\sqrt{-2 \log(U_1) }\sin(2\pi U_2)$

— VitalStatistix

2

@Vital: Nie jest to niepotrzebny komentarz; dobry. Transformacja Boxa-Mullera jest prawdopodobnie najłatwiejszym do zaprogramowania z minimalną szansą na nieumyślne zrobienie czegoś złego. Nie jest najszybszy , ale jest wystarczająco konkurencyjny. To powiedziawszy, korzystanie z ustalonej biblioteki kodów jest prawdopodobnie jeszcze bezpieczniejsze, zwłaszcza że miejscem, w którym można popełnić błąd, jest sposób generowania jednolitych zmiennych losowych !

— kardynał

@Vital: Dzięki, tego właśnie szukałem. Jeśli chcesz przekształcić swój komentarz w odpowiedź, chętnie go poprę.

— Fixee

1

@VitalStatistix, to świetny komentarz i wygląda na to, że tego właśnie szukał PO. Dlaczego nie przekształcić go w odpowiedź i być może rozwinąć ją nieco w ogólnej idei stosowania transformacji jednolitych zmiennych losowych. Wahałem się przed zrobieniem tego z tego powodu, o którym kardynał wspomina głównie dlatego, że nie wiem, czy domyślny generator jednolity z dowolnego języka jest dobrym generatorem.

— NRH

10

To jest naprawdę komentarz do odpowiedzi Michaela Lwa i komentarza Fixee, ale został opublikowany jako odpowiedź, ponieważ nie mam reputacji na tej stronie do komentowania.

$[0, 1]$ $6$ $1$

E [\sum_{i = 1}^{12} X_{i}] = \sum_{i = 1}^{12} E [X_{i}] = 12 \times \frac{1}{2} = 6

$E\left [\sum_{i=1}^{12} X_i\right ] = \sum_{i=1}^{12} E[X_i] = 12\times \frac{1}{2} = 6$

var [\sum_{i = 1}^{12} X_{i}] = \sum_{i = 1}^{12} var [X_{i}] = 12 \times \frac{1}{12} = 1.

$\text{var} \left [\sum_{i=1}^{12} X_i\right ] = \sum_{i=1}^{12} \text{var}[X_i] = 12\times \frac{1}{12} = 1.$

\sum_{i = 1}^{12} X_{i} - 6

$\sum_{i=1}^{12} X_i - 6$

\sqrt{10 / 12}

$\sqrt{10/12}$

\sum_{i = 1}^{12} X_{i} - 6

$\sum_{i=1}^{12} X_i - 6$

[- 6, 6]

$[-6, 6]$

6

$6$

— Dilip Sarwate
źródło

5

Oprócz odpowiedzi NRH, jeśli nadal nie masz możliwości wygenerowania losowych próbek ze „standardowego rozkładu normalnego” N (0,1), poniżej jest dobry i prosty sposób (ponieważ wspominasz, że nie masz danych statystycznych pakiet, poniższe funkcje powinny być dostępne w większości standardowych języków programowania).

1. Wygeneruj u i v jako dwie równomiernie rozmieszczone liczby losowe w zakresie od -1 do 1 przez
u = 2 r1 - 1iv = 2 r2 - 1

2. oblicz, w = u^2 + v^2jeśli w> 1 wróć do 1

3. zwróć u * zy = v * z z z= sqrt(-2ln(w)/w) Przykładowy kod wyglądałby tak:

u = 2 * random() - 1;
v = 2 * random() - 1;
w = pow(u, 2) + pow(v, 2);
if (w < 1) {
    z = sqrt((-2 * log(w)) / w);
    x = u * z;
    y = v * z;
    }

następnie użyj tego, co sugeruje MHR powyżej, aby uzyskać losowe odchylenia N(mu, sigma^2).

— francogrex
źródło

Kiedy opublikowałem swoją odpowiedź powyżej, nie zauważyłem, że @vitalStatistix dał ci algorytm transformacji Box-Mullera. Ten, który podałem powyżej, jest równie dobry, jak sądzę.

— francogrex,

2

Czy mógłbyś wyjaśnić powód generowania zmiennych normalnych z rozkładu jednorodnego (innego niż z perspektywy algorytmicznej), a nie tylko poprzez bezpośrednie użycie pdf rozkładu Gaussa / Normalnego? Czy jest to całkowicie błędne?

— Arun

4

@Arun Jeden powód: Metoda polarna Marsaglii jest przydatna, gdy masz tylko RNG, który generuje odchylenia mundurowe.

— chl

1

@Arun to najprostszy sposób. Możesz również wygenerować bezpośrednio z pliku PDF, używając na przykład metody „odrzucenia akceptacji”. Zamieściłem dla ciebie prosty przykład na mojej stronie (ponieważ w polu komentarza jest za mało miejsca).

— francogrex

4

Rozkład normalny pojawia się, gdy dodaje się wiele losowych wartości o podobnym rozkładzie (to znaczy podobnych do siebie). Jeśli dodasz do siebie dziesięć lub więcej równomiernie rozłożonych losowych wartości, suma jest prawie prawie rozkładem normalnym. (Dodaj więcej niż dziesięć, jeśli chcesz, aby było jeszcze bardziej normalne, ale dziesięć wystarczy na prawie wszystkie cele).

Powiedz, że twoje jednolite losowe wartości są równomiernie rozłożone między 0 a 1. Suma będzie wówczas wynosić od 0 do 10. Odejmij 5 od sumy, a średnia wynikowego rozkładu wyniesie 0. Teraz dzielisz wynik przez odchylenie standardowe (bliski) rozkład normalny i pomnożyć wynik przez pożądane odchylenie standardowe. Niestety nie jestem pewien, jakie jest standardowe odchylenie sumy dziesięciu jednolitych losowych odchyleń, ale jeśli będziemy mieli szczęście, ktoś powie nam w komentarzu!

Wolę rozmawiać ze studentami o rozkładzie normalnym w tych kategoriach, ponieważ użyteczność założenia rozkładu normalnego w wielu systemach wynika całkowicie z właściwości, że sumy wielu losowych wpływów prowadzą do rozkładu normalnego.

— Michael Lew
źródło

Używasz tutaj Central Limit Thm (to, że kilka losowych zmiennych sumuje się do normalnej zmiennej losowej). Nie zastanawiałem się nad tym, ponieważ myślałem, że będzie zbyt wolno, ale mówisz, że 10 wystarczy ?! Jest to lepsze niż obliczanie dziennika i sin / cos i sqrt!

— Fixee

Również średnia jednolitego rv dla [0,1] wynosi 0,5 przy wariancji 1/12. Jeśli zsumujesz 10 z nich, otrzymasz średnią 5 i wariancję 10/12 = 5/6.

— Fixee

1

Z pedagogicznego punktu widzenia metoda ta zapewnia miłą, przydatną dyskusję i demonstrację. Jednak zdecydowanie odradzam każdemu korzystanie z tego podejścia w praktyce.

— kardynał

1

@Fixee: Musisz się upewnić i zrównoważyć obliczenia

\log

$\log$ ,

\sin

$\sin$ ,

\cos

$\cos$ i pierwiastek kwadratowy z generowaniem dodatkowych jednorodnych zmiennych losowych. Na przykład procesory Intel mają wszystkie cztery te funkcje jako operacje wbudowane wykonywane sprzętowo. Pierwiastek kwadratowy jest podstawową „arytmetyczną” operacją zgodną ze standardami IEEE 754.

— kardynał

1

@Michael: Deklarowanie daje „prawo” dystrybucja jest nieco naciągane, zwłaszcza że dystrybucja zbliżenia posiada zwartą i wsparcie w wielu zastosowaniach, jeden robi opieki o tym, jak skutecznie się zmiennymi mogą być generowane. :) Chodzi o to, że dostępnych jest kilka znacznie lepszych opcji. Ale nadal uważam, że zapewnia coś przydatnego pedagogicznie.

— kardynał