Generieren Sie Zufallszahlen nach einer Normalverteilung in C / C ++

114

Wie kann ich nach einer Normalverteilung in C oder C ++ einfach Zufallszahlen generieren?

Ich möchte Boost nicht verwenden.

Ich weiß, dass Knuth ausführlich darüber spricht, aber ich habe seine Bücher momentan nicht zur Hand.

Damien
quelle
2
Duplikat des einen oder anderen von stackoverflow.com/questions/75677/… und stackoverflow.com/questions/1109446/…
dmckee --- Ex-Moderator Kätzchen

Antworten:

92

Es gibt viele Methoden, um aus einem regulären RNG Gauß-verteilte Zahlen zu generieren .

Die Box-Muller-Transformation wird üblicherweise verwendet. Es werden korrekt Werte mit einer Normalverteilung erzeugt. Die Mathematik ist einfach. Sie generieren zwei (einheitliche) Zufallszahlen und erhalten durch Anwenden einer Formel zwei normalverteilte Zufallszahlen. Geben Sie eine zurück und speichern Sie die andere für die nächste Anforderung einer Zufallszahl.

S.Lott
quelle
10
Wenn Sie Geschwindigkeit benötigen, ist die polare Methode jedoch schneller. Und der Ziggurat-Algorithmus noch mehr (wenn auch viel komplexer zu schreiben).
Joey
2
Ich habe hier eine Implementierung der Ziggurat gefunden. people.sc.fsu.edu/~jburkardt/c_src/ziggurat/ziggurat.html Es ist ziemlich vollständig.
Dwbrito
24
Beachten Sie, dass C ++ 11 std::normal_distributiongenau das hinzufügt , was Sie verlangen, ohne sich mit mathematischen Details zu befassen.
3
Es ist nicht garantiert, dass std :: normal_distribution auf allen Plattformen konsistent ist. Ich mache jetzt die Tests und MSVC bietet einen anderen Wertesatz als beispielsweise Clang. Die C ++ 11-Engines scheinen die gleichen Sequenzen zu generieren (bei gleichem Startwert), aber die C ++ 11-Distributionen scheinen mit unterschiedlichen Algorithmen auf verschiedenen Plattformen implementiert zu sein.
Arno Duvenhage
47

C ++ 11

C ++ 11 bietet std::normal_distribution, so würde ich heute vorgehen.

C oder älter C ++

Hier sind einige Lösungen in aufsteigender Reihenfolge der Komplexität:

  1. Addiere 12 einheitliche Zufallszahlen von 0 bis 1 und subtrahiere 6. Dies entspricht dem Mittelwert und der Standardabweichung einer normalen Variablen. Ein offensichtlicher Nachteil ist, dass der Bereich - im Gegensatz zu einer echten Normalverteilung - auf ± 6 begrenzt ist.

  2. Die Box-Muller-Transformation. Dies ist oben aufgeführt und relativ einfach zu implementieren. Wenn Sie jedoch sehr genaue Stichproben benötigen, beachten Sie, dass die Box-Muller-Transformation in Kombination mit einigen einheitlichen Generatoren unter einer Anomalie namens Neave-Effekt 1 leidet .

  3. Für beste Präzision empfehle ich, Uniformen zu zeichnen und die inverse kumulative Normalverteilung anzuwenden, um zu normalverteilten Variablen zu gelangen. Hier ist ein sehr guter Algorithmus für inverse kumulative Normalverteilungen.

1. HR Neave, "Zur Verwendung der Box-Muller-Transformation mit multiplikativen kongruenten Pseudozufallszahlengeneratoren", Applied Statistics, 22, 92-97, 1973

Peter G.
quelle
Hätten Sie vielleicht noch einen Link zum PDF über den Neave-Effekt? oder die Originalreferenz des Zeitschriftenartikels? danke
pyCthon
2
@stonybrooknick Die ursprüngliche Referenz wurde hinzugefügt. Coole Bemerkung: Während "Box Muller Neave" gegoogelt wurde, um die Referenz zu finden, tauchte diese Frage zum Stapelüberlauf auf der ersten Ergebnisseite auf!
Peter G.
Ja, es ist nicht jeder außerhalb bestimmter kleiner Gemeinschaften und Interessengruppen bekannt
pyCthon
@ Peter G. Warum sollte jemand Ihre Antwort ablehnen? - Möglicherweise hat dieselbe Person auch meinen Kommentar unten gemacht, mit dem ich einverstanden bin, aber ich fand Ihre Antwort sehr gut. Es wäre gut, wenn SO gemachte Downvotes einen echten Kommentar erzwingen würden. Ich vermute, dass die meisten Downvotes alter Themen nur frivol und trolly sind.
Pete855217
"Addiere 12 einheitliche Zahlen von 0-1 und subtrahiere 6." - Verteilung dieser Variablen hat Normalverteilung? Können Sie einen Zusammenhang mit der Ableitung herstellen, da während des zentralen Grenzwertsatzes der Ableitung n -> + inf eine sehr notwendige Annahme ist?
Bruziuz
31

Eine schnelle und einfache Methode besteht darin, eine Anzahl gleichmäßig verteilter Zufallszahlen zu summieren und ihren Durchschnitt zu ermitteln. Im zentralen Grenzwertsatz finden Sie eine vollständige Erklärung, warum dies funktioniert.

Paul R.
quelle
+1 Sehr interessanter Ansatz. Ist es verifiziert, wirklich normalverteilte Unterensembles für kleinere Gruppen zu geben?
Morlock
4
@Morlock Je größer die Anzahl der gemittelten Samples ist, desto näher kommt man einer Gaußschen Verteilung. Wenn Ihre Anwendung strenge Anforderungen an die Genauigkeit der Verteilung stellt, ist es möglicherweise besser, etwas Strengeres wie Box-Muller zu verwenden. Bei vielen Anwendungen, z. B. der Erzeugung von weißem Rauschen für Audioanwendungen, können Sie jedoch mit einer relativ geringen Anzahl davonkommen von gemittelten Proben (zB 16).
Paul R
2
Wie können Sie dies parametrisieren, um eine bestimmte Varianz zu erhalten? Angenommen, Sie möchten einen Mittelwert von 10 mit einer Standardabweichung von 1?
Morlock
1
@ Ben: Könntest du mich auf ein effizientes Algo dafür hinweisen? Ich habe die Mittelungstechnik bisher nur verwendet, um ungefähr Gaußsches Rauschen für die Audio- und Bildverarbeitung mit Echtzeitbeschränkungen zu erzeugen. Wenn es eine Möglichkeit gibt, dies in weniger Taktzyklen zu erreichen, kann dies sehr nützlich sein.
Paul R
1
@Petter: Sie haben wahrscheinlich im allgemeinen Fall Recht, wenn es um Gleitkommawerte geht. Es gibt jedoch immer noch Anwendungsbereiche wie Audio, in denen Sie schnelles ganzzahliges (oder festes Punkt) Gaußsches Rauschen wünschen, und Genauigkeit ist nicht allzu wichtig, in denen die einfache Mittelungsmethode effizienter und nützlicher ist (insbesondere für eingebettete Anwendungen, in denen dies möglicherweise nicht einmal der Fall ist Hardware-Gleitkomma-Unterstützung sein).
Paul R
24

Ich habe ein C ++ - Open Source-Projekt für den Benchmark zur Generierung normalverteilter Zufallszahlen erstellt .

Es vergleicht mehrere Algorithmen, einschließlich

  • Methode des zentralen Grenzwertsatzes
  • Box-Muller-Transformation
  • Marsaglia polare Methode
  • Ziggurat-Algorithmus
  • Inverse Transformations-Abtastmethode.
  • cpp11randomverwendet C ++ 11 std::normal_distributionmit std::minstd_rand(es ist eigentlich eine Box-Muller-Transformation in Clang).

Die Ergebnisse der floatVersion mit einfacher Genauigkeit ( ) auf iMac [email protected], Clang 6.1, 64-Bit:

normaldistf

Auf Richtigkeit überprüft das Programm den Mittelwert, die Standardabweichung, die Schiefe und die Kurtosis der Proben. Es wurde festgestellt, dass die CLT-Methode durch Summieren von 4, 8 oder 16 einheitlichen Zahlen keine gute Kurtosis aufweist wie die anderen Methoden.

Der Ziggurat-Algorithmus bietet eine bessere Leistung als die anderen. Es ist jedoch nicht für die SIMD-Parallelität geeignet, da es Tabellensuche und Verzweigungen benötigt. Box-Muller mit SSE2 / AVX-Befehlssatz ist viel schneller (x1.79, x2.99) als die Nicht-SIMD-Version des Zikkurat-Algorithmus.

Daher werde ich vorschlagen, Box-Muller für die Architektur mit SIMD-Befehlssätzen zu verwenden, und kann ansonsten Zikkurat sein.


PS Der Benchmark verwendet ein einfachstes LCG-PRNG zur Erzeugung gleichmäßig verteilter Zufallszahlen. Daher ist es für einige Anwendungen möglicherweise nicht ausreichend. Der Leistungsvergleich sollte jedoch fair sein, da alle Implementierungen dasselbe PRNG verwenden, sodass der Benchmark hauptsächlich die Leistung der Transformation testet.

Milo Yip
quelle
2
"Der Leistungsvergleich sollte jedoch fair sein, da alle Implementierungen dasselbe PRNG verwenden". Abgesehen davon, dass BM einen Eingangs-RN pro Ausgang verwendet, während CLT viel mehr usw. verwendet, ist die Zeit zum Generieren eines einheitlichen Zufalls # von Bedeutung.
Greggo
14

Hier ist ein C ++ - Beispiel, das auf einigen Referenzen basiert. Dies ist schnell und schmutzig. Sie sollten die Boost-Bibliothek nicht neu erfinden und verwenden.

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

Sie können ein QQ-Diagramm verwenden, um die Ergebnisse zu untersuchen und festzustellen, wie gut es einer realen Normalverteilung entspricht (ordnen Sie Ihre Stichproben 1..x, wandeln Sie die Ränge in Proportionen der Gesamtzahl von x um, dh wie viele Stichproben, erhalten Sie die z-Werte und zeichnen Sie sie auf. Eine gerade Linie nach oben ist das gewünschte Ergebnis.

Pete855217
quelle
1
Was ist sampleNormalManual ()?
Lösen von Rätseln
@solvingPuzzles - Entschuldigung, der Code wurde korrigiert. Es ist ein rekursiver Aufruf.
Pete855217
1
Dies kann bei einem seltenen Ereignis zum Absturz führen (wenn Sie Ihrem Chef die Anwendung zeigen, klingelt es?). Dies sollte mithilfe einer Schleife und nicht mithilfe einer Rekursion implementiert werden. Die Methode sieht ungewohnt aus. Wie ist die Quelle / wie heißt sie?
das Schwein
Box-Muller transkribiert aus einer Java-Implementierung. Wie gesagt, es ist schnell und schmutzig, zögern Sie nicht, es zu reparieren.
Pete855217
1
FWIW, viele Compiler werden in der Lage sein, diesen bestimmten rekursiven Aufruf in einen "Sprung an die Spitze der Funktion" umzuwandeln. Die Frage ist, ob Sie sich darauf verlassen möchten :-) Die Wahrscheinlichkeit, dass> 10 Iterationen erforderlich sind, beträgt 1 zu 4,8 Millionen. p (> 20) ist das Quadrat davon usw.
Greggo
12

Verwenden std::tr1::normal_distribution .

Der Namespace std :: tr1 ist kein Teil von boost. Es ist der Namespace, der die Bibliothekszusätze aus dem C ++ Technical Report 1 enthält und unabhängig von Boost in aktuellen Microsoft-Compilern und gcc verfügbar ist.

JoeG
quelle
25
Er hat nicht nach Standard gefragt, er hat nach "nicht Boost" gefragt.
JoeG
12

So generieren Sie die Beispiele auf einem modernen C ++ - Compiler.

#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
quelle
das generatorsollte wirklich ausgesät werden.
Walter
Es wird immer ausgesät. Es gibt einen Standard-Startwert.
Petter
4

Wenn Sie C ++ 11 verwenden, können Sie Folgendes verwenden 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);

Es gibt viele andere Verteilungen, mit denen Sie die Ausgabe der Zufallszahlen-Engine transformieren können.

Drew Noakes
quelle
Das wurde bereits von Ben erwähnt ( stackoverflow.com/a/11977979/635608 )
Mat
3

Ich habe die Definition des PDF unter http://www.mathworks.com/help/stats/normal-distribution.html befolgt und mir Folgendes ausgedacht :

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

Es ist vielleicht nicht der beste Ansatz, aber es ist ganz einfach.

MJVC
quelle
-1 Funktioniert nicht für zB RANDN2 (0.0, d + 1.0). Makros sind dafür berüchtigt.
Petter
Das Makro schlägt fehl, wenn rand()von RANDUeine Null zurückgibt, da Ln (0) undefiniert ist.
InterDist
Haben Sie diesen Code tatsächlich ausprobiert? Es sieht so aus, als hätten Sie eine Funktion erstellt, die Rayleigh-verteilte Zahlen generiert . Vergleichen Sie mit der Box-Muller-Transformation , bei der sie sich multiplizieren cos(2*pi*rand/RAND_MAX), während Sie sich mit multiplizieren (rand()%2 ? -1.0 : 1.0).
HelloGoodbye
1

Die FAQ-Liste von comp.lang.c bietet drei verschiedene Möglichkeiten, um auf einfache Weise Zufallszahlen mit einer Gaußschen Verteilung zu generieren.

Sie können einen Blick darauf werfen: http://c-faq.com/lib/gaussian.html

Delgan
quelle
1

Box-Muller-Implementierung:

#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
quelle
1

Es gibt verschiedene Algorithmen für die inverse kumulative Normalverteilung. Die beliebtesten in der quantitativen Finanzierung werden auf http://chasethedevil.github.io/post/monte-carlo--inverse-cumulative-normal-distribution/ getestet.

Meiner Meinung nach gibt es keinen großen Anreiz, etwas anderes als den Algorithmus AS241 von Wichura zu verwenden : Er ist maschinenpräzise, ​​zuverlässig und schnell. Engpässe treten bei der Gaußschen Zufallszahlengenerierung selten auf.

Darüber hinaus zeigt es den Nachteil von Ziggurat-ähnlichen Ansätzen.

Die Top-Antwort hier befürwortet Box-Müller, Sie sollten sich bewusst sein, dass es bekannte Mängel gibt. Ich zitiere https://www.sciencedirect.com/science/article/pii/S0895717710005935 :

In der Literatur wird Box-Muller manchmal als etwas minderwertig angesehen, hauptsächlich aus zwei Gründen. Erstens, wenn man die Box-Muller-Methode auf Zahlen eines schlechten linearen Kongruenzgenerators anwendet, bieten die transformierten Zahlen eine extrem schlechte Abdeckung des Raums. Diagramme transformierter Zahlen mit spiralförmigen Schwänzen finden sich in vielen Büchern, insbesondere im klassischen Buch von Ripley, der wahrscheinlich der erste war, der diese Beobachtung machte. "

jherek
quelle
0

1) Die grafisch intuitive Möglichkeit, Gaußsche Zufallszahlen zu generieren, besteht in der Verwendung der Monte-Carlo-Methode. Sie würden mit Ihrem Pseudozufallszahlengenerator in C einen zufälligen Punkt in einem Feld um die Gaußsche Kurve erzeugen. Mit der Verteilungsgleichung können Sie berechnen, ob dieser Punkt innerhalb oder unterhalb der Gaußschen Verteilung liegt. Wenn dieser Punkt innerhalb der Gaußschen Verteilung liegt, haben Sie Ihre Gaußsche Zufallszahl als x-Wert des Punktes.

Diese Methode ist nicht perfekt, da die Gaußsche Kurve technisch gegen unendlich geht und Sie keine Box erstellen können, die sich in der x-Dimension der Unendlichkeit nähert. Aber die Guass'sche Kurve nähert sich in der y-Dimension ziemlich schnell 0, also würde ich mir darüber keine Sorgen machen. Die Einschränkung der Größe Ihrer Variablen in C kann Ihre Genauigkeit eher einschränken.

2) Eine andere Möglichkeit wäre die Verwendung des zentralen Grenzwertsatzes, der besagt, dass unabhängige Zufallsvariablen beim Hinzufügen eine Normalverteilung bilden. Unter Berücksichtigung dieses Theorems können Sie eine Gaußsche Zufallszahl approximieren, indem Sie eine große Anzahl unabhängiger Zufallsvariablen hinzufügen.

Diese Methoden sind nicht die praktischsten, aber das ist zu erwarten, wenn Sie keine bereits vorhandene Bibliothek verwenden möchten. Denken Sie daran, dass diese Antwort von jemandem stammt, der wenig oder keine Erfahrung mit Kalkül oder Statistik hat.

dan dan
quelle
0

Monte-Carlo-Methode Der intuitivste Weg, dies zu tun, wäre die Verwendung einer Monte-Carlo- Methode. Nehmen Sie einen geeigneten Bereich -X, + X. Größere Werte von X führen zu einer genaueren Normalverteilung, die Konvergenz dauert jedoch länger. ein. Wählen Sie eine Zufallszahl z zwischen -X bis X. b. Halten Sie mit einer Wahrscheinlichkeit fest, N(z, mean, variance)wo N die Gaußsche Verteilung ist. Andernfalls fallen lassen und zu Schritt (a) zurückkehren.

Jagat
quelle
-1

Schau dir an, was ich gefunden habe.

Diese Bibliothek verwendet den Ziggurat-Algorithmus.

Dwbrito
quelle
-3

Computer ist ein deterministisches Gerät. Es gibt keine Zufälligkeit bei der Berechnung. Darüber hinaus kann das arithmetische Gerät in der CPU die Summierung über einen endlichen Satz von Ganzzahlen (Durchführen einer Auswertung im endlichen Feld) und einen endlichen Satz von reellen rationalen Zahlen auswerten. Und auch bitweise Operationen durchgeführt. Mathe macht einen Deal mit großartigeren Mengen wie [0.0, 1.0] mit unendlich vielen Punkten.

Sie können mit einem Controller einen Draht im Computer hören, aber würde er gleichmäßige Verteilungen haben? Ich weiß es nicht. Wenn jedoch angenommen wird, dass das Signal das Ergebnis einer großen Menge unabhängiger Zufallsvariablen ist, erhalten Sie eine ungefähr normalverteilte Zufallsvariable (dies wurde in der Wahrscheinlichkeitstheorie bewiesen).

Es gibt Algorithmen, die als Pseudozufallsgenerator bezeichnet werden. Wie ich dachte, besteht der Zweck des Pseudozufallsgenerators darin, die Zufälligkeit zu emulieren. Und die Kriterien für Goodnes sind: - Die empirische Verteilung wird (in gewissem Sinne - punktuell, einheitlich, L2) zu theoretischen Werten konvergiert. - Werte, die Sie vom Zufallsgenerator erhalten, scheinen ideenabhängig zu sein. Natürlich ist es aus "realer Sicht" nicht wahr, aber wir gehen davon aus, dass es wahr ist.

Eine der beliebtesten Methoden - Sie können 12 irv mit gleichmäßigen Verteilungen summieren ... Aber um ehrlich zu sein, während der Ableitung des zentralen Grenzwertsatzes mit Hilfe der Fourier-Transformation, Taylor-Reihe, müssen einige Male n -> + inf-Annahmen getroffen werden. Zum Beispiel theoretisch - Ich persönlich verstehe nicht, wie Leute eine Summierung von 12 irv mit gleichmäßiger Verteilung durchführen.

Ich hatte Wahrscheinlichkeitstheorie in der Universität. Und besonders für mich ist es nur eine mathematische Frage. In der Universität habe ich folgendes Modell gesehen:


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

So wie es zu tun war, war es nur ein Beispiel. Ich denke, es gibt andere Möglichkeiten, es zu implementieren.

Der Beweis, dass es richtig ist, findet sich in diesem Buch "Moskau, BMSTU, 2004: XVI Wahrscheinlichkeitstheorie, Beispiel 6.12, S.246-247" von Krishchenko Alexander Petrovich ISBN 5-7038-2485-0

Leider weiß ich nicht, ob es eine Übersetzung dieses Buches ins Englische gibt.

Bruziuz
quelle
Ich habe mehrere Abstimmungen. Lassen Sie mich wissen, was hier schlecht ist?
Bruziuz
Die Frage ist, wie man Pseudozufallszahlen im Computer erzeugt (ich weiß, die Sprache ist hier locker), es ist keine Frage der mathematischen Existenz.
user2820579
Ja, du hast Recht. Und die Antwort ist, wie man eine Pseudozufallszahl mit Normalverteilung basierend auf einem Generator mit gleichmäßiger Verteilung erzeugt. Wenn der Quellcode bereitgestellt wurde, können Sie ihn in einer beliebigen Sprache umschreiben.
Bruziuz
Klar, ich denke der Typ sucht zB "Numerische Rezepte in C / C ++". Übrigens, nur um unsere Diskussion zu ergänzen, geben die Autoren dieses letzten Buches interessante Referenzen für ein paar Pseudozufallsgeneratoren, die die Standards für "anständige" Generatoren erfüllen.
user2820579
1
Ich habe hier ein Backup erstellt: sites.google.com/site/burlachenkok/download
bruziuz