Jak mogę łatwo wygenerować liczby losowe zgodnie z rozkładem normalnym w C lub C ++?

Nie chcę używać Boost.

Wiem, że Knuth długo o tym mówi, ale nie mam teraz pod ręką jego książek.

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.

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 ¹ .

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}

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.

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:

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.

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

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.

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;

Możesz użyć GSL . Podano kilka pełnych przykładów, aby zademonstrować, jak z niego korzystać.

Zajrzyj na: http://www.cplusplus.com/reference/random/normal_distribution/ . To najprostszy sposób tworzenia rozkładów normalnych.

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.

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.

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

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

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 ”

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.

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

Zobacz, co znalazłem.

Ta biblioteka używa algorytmu Ziggurat.

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.