Generuj liczby losowe zgodnie z rozkładem normalnym w C / C ++

Odpowiedzi:

92

Istnieje wiele metod generowania liczb o rozkładzie Gaussa na podstawie zwykłego RNG .

Transformacja Boxa-Mullera jest powszechnie używany. Prawidłowo generuje wartości z rozkładem normalnym. Matematyka jest łatwa. Generujesz dwie (jednolite) liczby losowe, a stosując do nich wzór, otrzymujesz dwie liczby losowe o normalnym rozkładzie. Zwróć jeden, a drugi zachowaj na następne żądanie losowej liczby.

S.Lott
źródło
10
Jeśli potrzebujesz szybkości, metoda biegunowa jest jednak szybsza. A algorytm Ziggurat jeszcze bardziej (choć znacznie bardziej skomplikowany do napisania).
Joey,
2
znalazłem implementację Ziggurata tutaj people.sc.fsu.edu/~jburkardt/c_src/ziggurat/ziggurat.html Całkiem kompletne.
dwbrito
24
Uwaga, C ++ 11 dodaje, std::normal_distributionktóry robi dokładnie to, o co prosisz, bez zagłębiania się w szczegóły matematyczne.
3
Nie gwarantuje się spójności std :: normal_distribution na wszystkich platformach. Robię teraz testy, a MSVC zapewnia inny zestaw wartości niż na przykład Clang. Wydaje się, że silniki C ++ 11 generują te same sekwencje (biorąc pod uwagę to samo ziarno), ale dystrybucje C ++ 11 wydają się być implementowane przy użyciu różnych algorytmów na różnych platformach.
Arno Duvenhage
47

C ++ 11

C ++ 11 oferuje std::normal_distribution, tak bym dzisiaj poszedł.

C lub starszy C ++

Oto kilka rozwiązań w kolejności rosnącej złożoności:

  1. Dodaj 12 jednakowych liczb losowych od 0 do 1 i odejmij 6. To dopasuje średnią i odchylenie standardowe normalnej zmiennej. Oczywistą wadą jest to, że zakres jest ograniczony do ± 6 - w przeciwieństwie do prawdziwego rozkładu normalnego.

  2. Transformacja Boxa-Mullera. Jest to wymienione powyżej i jest stosunkowo proste do wdrożenia. Jeśli jednak potrzebujesz bardzo precyzyjnych próbek, pamiętaj, że transformata Box-Mullera w połączeniu z niektórymi jednorodnymi generatorami cierpi na anomalię zwaną Neave Effect 1 .

  3. Aby uzyskać najlepszą precyzję, sugeruję rysowanie mundurów i stosowanie odwrotnego skumulowanego rozkładu normalnego, aby uzyskać rozkład normalny. Oto bardzo dobry algorytm odwrotnych skumulowanych rozkładów normalnych.

1. HR Neave, „On using the Box-Muller Transformation with multiplicative congruential pseudolandom number generators”, Applied Statistics, 22, 92-97, 1973

Peter G.
źródło
Czy przypadkiem miałbyś inny link do pliku PDF na temat efektu Neave? lub odniesienie do oryginalnego artykułu w czasopiśmie? dziękuję
pyCthon
2
@stonybrooknick Oryginalne odniesienie zostało dodane. Fajna uwaga: podczas wyszukiwania w Google „box muller neave” w celu znalezienia odniesienia, to samo pytanie o przepełnienie stosu pojawiło się na pierwszej stronie wyników!
Peter G.
tak, to nie jest dobrze znane poza niektórymi małymi społecznościami i grupami interesu
pyCthon
@Peter G. Dlaczego ktoś miałby negatywnie oceniać Twoją odpowiedź? - prawdopodobnie ta sama osoba też zrobiła mój komentarz poniżej, z czym nie mam problemu, ale uznałem, że twoja odpowiedź była bardzo dobra. Byłoby dobrze, gdyby tak skierowane głosy przeciwne wymuszały prawdziwy komentarz. Podejrzewam, że większość głosów przeciw starych tematów jest po prostu frywolna i trollowa.
Pete855217
„Dodaj 12 jednolitych liczb od 0 do 1 i odejmij 6”. - rozkład tej zmiennej będzie miał rozkład normalny? Czy możesz podać link z wyprowadzeniem, ponieważ podczas wyprowadzania centralne twierdzenie graniczne n -> + inf jest bardzo potrzebne.
bruziuz
31

Szybką i łatwą metodą jest po prostu zsumowanie liczby równomiernie rozłożonych liczb losowych i obliczenie ich średniej. Zobacz centralne twierdzenie graniczne, aby uzyskać pełne wyjaśnienie, dlaczego to działa.

Paul R.
źródło
+1 Bardzo ciekawe podejście. Czy zweryfikowano, że rzeczywiście daje się podzestawy o normalnym rozkładzie dla mniejszych grup?
Morlock,
4
@Morlock Im większa liczba uśrednionych próbek, tym bardziej zbliżasz się do rozkładu Gaussa. Jeśli twoja aplikacja ma ścisłe wymagania co do dokładności dystrybucji, może być lepiej, jeśli użyjesz czegoś bardziej rygorystycznego, takiego jak Box-Muller, ale w przypadku wielu aplikacji, np. Generowanie białego szumu dla aplikacji audio, możesz uciec z dość małą liczbą uśrednionych próbek (np. 16).
Paul R
2
Poza tym, jak sparametryzować to, aby uzyskać pewną wariancję, powiedzmy, że chcesz uzyskać średnią 10 z odchyleniem standardowym 1?
Morlock
1
@Ben: czy możesz mi wskazać wydajne algorytmy? Używałem techniki uśredniania tylko do generowania szumu około Gaussa do przetwarzania dźwięku i obrazu z ograniczeniami w czasie rzeczywistym - jeśli istnieje sposób na osiągnięcie tego w mniejszej liczbie cykli zegara, może to być bardzo przydatne.
Paul R
1
@Petter: prawdopodobnie masz rację w ogólnym przypadku, jeśli chodzi o wartości zmiennoprzecinkowe. Nadal istnieją obszary zastosowań, takie jak audio, w których potrzebujesz szybkiego szumu gaussowskiego (lub punktu stałego), a dokładność nie jest zbyt ważna, gdzie prosta metoda uśredniania jest bardziej wydajna i użyteczna (szczególnie w przypadku aplikacji osadzonych, gdzie może nawet nie być być sprzętową obsługą zmiennoprzecinkowych).
Paul R
24

Stworzyłem projekt open source w C ++ dla standardowego testu porównawczego generowania liczb losowych .

Porównuje kilka algorytmów, w tym

  • Metoda centralnego twierdzenia granicznego
  • Transformacja Boxa-Mullera
  • Metoda polarna Marsaglia
  • Algorytm ziggurata
  • Metoda próbkowania z odwrotną transformacją.
  • cpp11randomużywa C ++ 11 std::normal_distributionz std::minstd_rand(w rzeczywistości jest to transformacja Boxa-Mullera w clang).

Wyniki wersji z pojedynczą precyzją ( float) na iMac Corei5-3330S@2,70GHz, clang 6.1, 64-bit:

normaldistf

Dla poprawności program weryfikuje średnią, odchylenie standardowe, skośność i kurtoozę próbek. Stwierdzono, że metoda CLT polegająca na sumowaniu 4, 8 lub 16 liczb jednolitych nie ma dobrej kurtozy, tak jak inne metody.

Algorytm Ziggurat ma lepszą wydajność niż inne. Jednak nie nadaje się do równoległości SIMD, ponieważ wymaga wyszukiwania w tabeli i rozgałęzień. Box-Muller z zestawem instrukcji SSE2 / AVX jest znacznie szybszy (x1,79, x2,99) niż wersja algorytmu ziggurat bez SIMD.

Dlatego zasugeruję użycie Box-Mullera dla architektury z zestawami instrukcji SIMD, a w przeciwnym razie może być zigguratem.


PS benchmark wykorzystuje najprostszy LCG PRNG do generowania równomiernie rozłożonych liczb losowych. W przypadku niektórych zastosowań może to nie wystarczyć. Ale porównanie wydajności powinno być uczciwe, ponieważ wszystkie implementacje używają tego samego PRNG, więc test porównawczy testuje głównie wydajność transformacji.

Milo Yip
źródło
2
„Ale porównanie wydajności powinno być uczciwe, ponieważ wszystkie implementacje używają tego samego PRNG” .. Z wyjątkiem tego, że BM używa jednego wejściowego RN na wyjście, podczas gdy CLT używa ich o wiele więcej itd., Więc czas na wygenerowanie jednolitego, losowego # ma znaczenie.
greggo
14

Oto przykład C ++, oparty na niektórych odniesieniach. Jest to szybkie i brudne, lepiej nie wymyślać ponownie i nie używać biblioteki boost.

#include "math.h" // for RAND, and rand
double sampleNormal() {
    double u = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double v = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double r = u * u + v * v;
    if (r == 0 || r > 1) return sampleNormal();
    double c = sqrt(-2 * log(r) / r);
    return u * c;
}

Możesz użyć wykresu QQ, aby zbadać wyniki i zobaczyć, jak dobrze przybliża on rzeczywisty rozkład normalny (uszereguj próbki 1..x, zamień rangi na proporcje całkowitej liczby x tj. Ile próbek, uzyskaj wartości z i wykreśl je. Prosta w górę jest pożądanym wynikiem).

Pete855217
źródło
1
Co to jest sampleNormalManual ()?
rozwiązywanie łamigłówek
@solvingPuzzles - przepraszam, poprawiłem kod. To połączenie rekurencyjne.
Pete855217
1
To musi się zawiesić podczas jakiegoś rzadkiego wydarzenia (pokazanie aplikacji swojemu szefowi dzwoni dzwonkiem?). Powinno to być realizowane przy użyciu pętli, a nie rekursji. Ta metoda wygląda na nieznaną. Jakie jest źródło / jak się nazywa?
świnia
Box-Muller przepisał z implementacji Java. Jak powiedziałem, jest szybki i brudny, możesz to naprawić.
Pete855217
1
FWIW, wiele kompilatorów będzie w stanie zamienić to konkretne wywołanie rekurencyjne w „skok na górę funkcji”. Pytanie brzmi, czy chcesz na to liczyć :-) Ponadto prawdopodobieństwo, że zajmie to> 10 iteracji, wynosi 1 na 4,8 miliona. p (> 20) jest kwadratem tego itd.
greggo
12

Użyj std::tr1::normal_distribution.

Przestrzeń nazw std :: tr1 nie jest częścią boost. Jest to przestrzeń nazw, która zawiera dodatki do bibliotek z C ++ Technical Report 1 i jest dostępna w aktualnych kompilatorach Microsoft i gcc, niezależnie od boost.

JoeG
źródło
25
Nie prosił o standard, prosił o „nie doładowanie”.
JoeG
12

W ten sposób generujesz próbki na nowoczesnym kompilatorze C ++.

#include <random>
...
std::mt19937 generator;
double mean = 0.0;
double stddev  = 1.0;
std::normal_distribution<double> normal(mean, stddev);
cerr << "Normal: " << normal(generator) << endl;
Petter
źródło
generatorpowinien być naprawdę zaszczepiono.
Walter
To jest zawsze zaszczepione. Istnieje domyślne ziarno.
Petter
4

Jeśli używasz C ++ 11, możesz użyć std::normal_distribution:

#include <random>

std::default_random_engine generator;
std::normal_distribution<double> distribution(/*mean=*/0.0, /*stddev=*/1.0);

double randomNumber = distribution(generator);

Istnieje wiele innych dystrybucji, których można użyć do przekształcenia danych wyjściowych silnika liczb losowych.

Drew Noakes
źródło
Wspomniał o tym już Ben ( stackoverflow.com/a/11977979/635608 )
Mat
3

Postępowałem zgodnie z definicją pliku PDF podaną w http://www.mathworks.com/help/stats/normal-distribution.html i wymyśliłem to:

const double DBL_EPS_COMP = 1 - DBL_EPSILON; // DBL_EPSILON is defined in <limits.h>.
inline double RandU() {
    return DBL_EPSILON + ((double) rand()/RAND_MAX);
}
inline double RandN2(double mu, double sigma) {
    return mu + (rand()%2 ? -1.0 : 1.0)*sigma*pow(-log(DBL_EPS_COMP*RandU()), 0.5);
}
inline double RandN() {
    return RandN2(0, 1.0);
}

To może nie jest najlepsze podejście, ale jest dość proste.

MJVC
źródło
-1 Nie działa np. Dla RANDN2 (0.0, d + 1.0). Makra są z tego znane.
Petter
Makro nie powiedzie się, jeśli rand()of RANDUzwróci zero, ponieważ Ln (0) jest niezdefiniowane.
interDist
Czy faktycznie wypróbowałeś ten kod? Wygląda na to, że utworzyłeś funkcję, która generuje liczby o rozkładzie Rayleigha . Porównaj z transformacją Boxa-Mullera , gdzie mnożą się przez cos(2*pi*rand/RAND_MAX), a ty mnożysz przez (rand()%2 ? -1.0 : 1.0).
HelloGoodbye,
1

Lista często zadawanych pytań dotyczących comp.lang.c zawiera trzy różne sposoby łatwego generowania liczb losowych z rozkładem Gaussa.

Możesz rzucić okiem: http://c-faq.com/lib/gaussian.html

Delgan
źródło
1

Wdrożenie Box-Mullera:

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
 // return a uniformly distributed random number
double RandomGenerator()
{
  return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );
}
 // return a normally distributed random number
double normalRandom()
{
  double y1=RandomGenerator();
  double y2=RandomGenerator();
  return cos(2*3.14*y2)*sqrt(-2.*log(y1));
}

int main(){
double sigma = 82.;
double Mi = 40.;
  for(int i=0;i<100;i++){
double x = normalRandom()*sigma+Mi;
    cout << " x = " << x << endl;
  }
  return 0;
}
Sysadmin
źródło
1

Istnieją różne algorytmy odwrotnego skumulowanego rozkładu normalnego. Najpopularniejsze w finansach ilościowych są testowane na http://chasethedevil.github.io/post/monte-carlo--inverse-cumulative-normal-distribution/

Moim zdaniem nie ma zbytniej zachęty do używania czegoś innego niż algorytm AS241 firmy Wichura : to precyzja maszyny, niezawodność i szybkość. Wąskie gardła rzadko występują w generowaniu liczb losowych Gaussa.

Ponadto pokazuje wady podejść podobnych do Zigguratu.

Najlepsza odpowiedź to zwolennicy Box-Müllera, należy mieć świadomość, że ma on znane wady. Cytuję https://www.sciencedirect.com/science/article/pii/S0895717710005935 :

w literaturze Box-Muller bywa uważany za nieco gorszego, głównie z dwóch powodów. Po pierwsze, jeśli zastosuje się metodę Boxa-Mullera do liczb ze złego generatora liniowego kongruencji, to przekształcone liczby zapewniają wyjątkowo słabe pokrycie przestrzeni. Wykresy przekształconych liczb ze spiralnymi ogonami można znaleźć w wielu książkach, zwłaszcza w klasycznej książce Ripleya, który prawdopodobnie był pierwszym, który dokonał tej obserwacji ”

jherek
źródło
0

1) Graficznie intuicyjny sposób generowania liczb losowych Gaussa polega na użyciu czegoś podobnego do metody Monte Carlo. Możesz wygenerować losowy punkt w ramce wokół krzywej Gaussa, używając swojego generatora liczb pseudolosowych w C. Możesz obliczyć, czy ten punkt znajduje się wewnątrz, czy pod rozkładem Gaussa, używając równania rozkładu. Jeśli ten punkt znajduje się w rozkładzie Gaussa, to masz swoją losową liczbę Gaussa jako wartość x punktu.

Ta metoda nie jest doskonała, ponieważ z technicznego punktu widzenia krzywa Gaussa ciągnie się w kierunku nieskończoności, a nie można było stworzyć prostokąta zbliżającego się do nieskończoności w wymiarze x. Ale krzywa Guassiana zbliża się do 0 w wymiarze y dość szybko, więc nie martwiłbym się tym. Ograniczenie rozmiaru twoich zmiennych w C może być czynnikiem ograniczającym dokładność.

2) Innym sposobem byłoby użycie Centralnego Twierdzenia Granicznego, które stwierdza, że ​​po dodaniu niezależnych zmiennych losowych tworzą one rozkład normalny. Pamiętając o tym twierdzeniu, można przybliżyć liczbę losową Gaussa, dodając dużą liczbę niezależnych zmiennych losowych.

Te metody nie są najbardziej praktyczne, ale należy się tego spodziewać, gdy nie chcesz korzystać z istniejącej biblioteki. Pamiętaj, że ta odpowiedź pochodzi od kogoś, kto ma niewielkie lub żadne doświadczenie w rachunku różniczkowym lub statystycznym.

dan dan
źródło
0

Metoda Monte Carlo Najbardziej intuicyjnym sposobem byłoby zastosowanie metody Monte Carlo. Weź odpowiedni zakres -X, + X. Większe wartości X spowodują dokładniejszy rozkład normalny, ale zbieżność zajmie więcej czasu. za. Wybierz losową liczbę z od -X do X. b. Zachowaj z prawdopodobieństwem, N(z, mean, variance)gdzie N jest rozkładem Gaussa. Upuść w przeciwnym razie i wróć do kroku (a).

Jagat
źródło
-3

Komputer jest urządzeniem deterministycznym. W obliczeniach nie ma przypadkowości. Ponadto urządzenie arytmetyczne w CPU może oceniać sumę po pewnym skończonym zbiorze liczb całkowitych (wykonując obliczenia w polu skończonym) i skończonym zbiorze rzeczywistych liczb wymiernych. A także wykonywał operacje bitowe. Matematyka radzi sobie z większymi zestawami, takimi jak [0.0, 1.0], z nieskończoną liczbą punktów.

Możesz posłuchać przewodu wewnątrz komputera z jakimś kontrolerem, ale czy miałby on jednolite dystrybucje? Nie wiem Ale jeśli przyjmiemy, że jego sygnał jest wynikiem akumulacji dużej ilości niezależnych zmiennych losowych, to otrzymamy zmienną losową o rozkładzie normalnym (zostało to udowodnione w teorii prawdopodobieństwa)

Istnieją algorytmy zwane - generatorem pseudolosowym. Uważam, że celem generatora pseudolosowego jest naśladowanie losowości. Kryteria dobrobytu są następujące: - rozkład empiryczny jest zbieżny (w pewnym sensie - punktowy, jednolity, L2) do teoretycznego - wartości, które otrzymujesz z generatora losowego, wydają się być niezależne. Oczywiście nie jest to prawdą z „prawdziwego punktu widzenia”, ale zakładamy, że to prawda.

Jedna z popularnych metod - można zsumować 12 irv z rozkładami jednorodnymi ... Ale szczerze mówiąc podczas wyprowadzania Centralne twierdzenie graniczne z pomocą transformaty Fouriera, szereg Taylora, trzeba mieć założenia n -> + inf. Na przykład teoretycznie - Osobiście nie rozumiem, jak ludzie wykonują zsumowanie 12 irv z równomiernym rozkładem.

Miałem teorię prawdopodobieństwa na uniwersytecie. A szczególnie dla mnie jest to tylko pytanie matematyczne. Na uniwersytecie widziałem następujący model:


double generateUniform(double a, double b)
{
  return uniformGen.generateReal(a, b);
}

double generateRelei(double sigma)
{
  return sigma * sqrt(-2 * log(1.0 - uniformGen.generateReal(0.0, 1.0 -kEps)));
}
double generateNorm(double m, double sigma)
{
  double y2 = generateUniform(0.0, 2 * kPi);
  double y1 = generateRelei(1.0);
  double x1 = y1 * cos(y2);
  return sigma*x1 + m;
}

Tak więc jak do zrobienia to był tylko przykład, myślę, że istnieją inne sposoby na jego realizację.

Dowód, że jest to poprawne, można znaleźć w tej książce „Moskwa, BMSTU, 2004: XVI Teoria prawdopodobieństwa, przykład 6.12, str. 246-247” autorstwa Krishchenko Aleksandra Pietrowicza ISBN 5-7038-2485-0

Niestety nie wiem o istnieniu tłumaczenia tej książki na język angielski.

bruziuz
źródło
Mam kilka głosów przeciw. Daj mi znać, co tu jest źle?
bruziuz
Pytanie brzmi, jak wygenerować liczby pseudolosowe w komputerze (wiem, język jest tu luźny), nie jest to kwestia matematycznego istnienia.
user2820579
Tak, masz rację. A odpowiedź brzmi: jak wygenerować liczbę pseudolosową z rozkładem normalnym na podstawie generatora o rozkładzie równomiernym. Kod źródłowy został dostarczony, możesz go przepisać w dowolnym języku.
bruziuz
Jasne, myślę, że facet szuka np. „Przepisów numerycznych w C / C ++”. Nawiasem mówiąc, aby uzupełnić naszą dyskusję, autorzy tej ostatniej książki podają interesujące referencje dotyczące kilku pseudolosowych generatorów, które spełniają standardy bycia „przyzwoitymi” generatorami.
user2820579
1
Zrobiłem kopię zapasową tutaj: sites.google.com/site/burlachenkok/download
bruziuz