Pseudozufallszahlengeneratoren nichtkonstanter Wahrscheinlichkeitsverteilung
mathematische Zufallszahlengeneratoren » Steuern pseudozufälliger Ereignisse

Voraussetzungen

Programme wie Spiele, Bildschirmschoner, technische oder naturwissenschaftliche Simulatoren, künstlerische und andere Applikationen werden oftmals von Zufallszahlen gesteuert. Dabei können folgende Anforderungen auftreten:

Folgende Rahmenbedingungen können die Realisation zusätzlich erschweren:

Grundlagen

In C[++,#] gibt es die Standardlibrary-Funktion rand(), die einen Wert zwischen 0 und RAND_MAX (0x7FFF=32767) – also 32768 verschiedene Werte liefert. (Compiler anderer Programmiersprachen sind zumeist in C++ geschrieben, daher kann dort analog vorgegangen werden.) Daraus kann folgendermaßen ein pseudozufälliges Doppelwort gebildet werden:

typedef unsigned long DWORD; // mapidbg.h
class Rnd { //……
enum { RMAX=RAND_MAX+1 };
static DWORD Rnd::DWord()
{ return (RMAX*RMAX*(DWORD)rand()+RMAX*rand()+rand()); }

Da rand()seine Werte immer in der gleichen Reihenfolge liefert, werden nur bestimmte Belegungen des Doppelwortes erzeugt. Soll bestmögliche Verteilung aller Bits von Rnd::DWord()erreicht werden, sollte man unterschiedliche rand()'s verwenden:

return (RMAX*RMAX*(DWORD)myrand1()+RMAX*myrand2()+myrand3());

Will man also ein Ereignis mit 17,3%iger Wahrscheinlichkeit auftreten lassen, genügt folgendes Konstrukt mit zumeist ausreichender Genauigkeit:

if((Rnd::DWord()%/*modulo*/1000)<173) { /*……*/ }

Sollte tatsächlich allerhöchste Präzision erforderlich sein, darf man von Rnd::DWord() nur einen Bereich auswerten, der durch den Modulokoeffizienten teilbar ist:

DWORD x; do { x=Rnd::DWord(); } while(x>=2147483000);
         // nächstkleineres unter 2^31 = (0x80000000/1000)*1000
if(x%100)<17) { /*……*/ }

Der eigene Generator

Wer – aus welchen Gründen auch immer (z.B. wegen Reproduzierbarkeit und Steuerbarkeit) – sich einen eigenen konstanten Zufallszahlenerzeuger programmieren will, kann das folgendermaßen durchführen:

static DWORD NextDW(DWORD&prev)
{ prev*=1103515245; prev+=12345; return prev; }

Der Multiplikator und der Summenbildner können dabei nach eigenem Gutdünken ausgewählt werden – der resultierende Generator muss aber unbedingt getestet werden!
Auch sollte man die niederwertigen Bits des erzeugten Doppelwortes nicht auswerten, da sie die geforderte Wahrscheinlichkeit, dass sie zu jeweils 50% (nicht) gesetzt sind, nicht in ausreichender Regelmäßigkeit erreichen:

static int my_rand(DWORD&prev) // prev ... Speicher mit altem Wert
{ return(int)((NextDW(prev)>>/*shift right*/17)&RAND_MAX); }

(siehe: Bjarne Stroustrup, The C++ Programming Language, Third Edition, page 685f)

Die Modifikationsfunktion f(x)

Will man nun nichtkonstante Wahrscheinlichkeit erreichen, erzeugt man folgendermaßen aus dem Wertebereich 0<=r<=RAND_MAX einen Wertebereich 0.0<=x<1.0:

double x=((double)r)/((double)(RAND_MAX+1));

Daraus bildet man den Eingangswertebereich zu einer Funktion f(x)

und vollzieht die Transformation in den gewünschten Wertebereich:

// MAX_CALC_VALS ... oberster zulässiger Wert des Ausgangsdoppelwortes
DWORD y = (DWORD)((double)MAX_CALC_VALS*f(x));

Je nach Ausgestaltung von y werden jetzt z.B. kleine Werte von y häufig, große Werte von y aber selten vorkommen.  Häufige Werte treten auf, wo f(x) flach ist (dy/dx«), seltene Werte an steilen Bereichen. (dy/dx»)
Eine Spiegelung der Eingangswerte um ihren Mittelwert läßt die sich ergebende Häufigkeit offensichtlich unverändert. Will man die Häufigkeit an ihrem Mittelwert spiegeln (z.B. für zunehmende statt abnehmende Häufigkeit), ist dabei der Unterschied zwischen Maximalwert und Anzahl der Werte (=Maximalwert+1, z.B: (RAND_MAX+1)) zu beachten!

Auswahl der Modifikationsfunktion f(x) mit rižArtë

Zur Unterstützung bei der Auswahl einer geeigneten Modifikationsfunktion wurde in dem Programm rižArtë ein eigenes Tool geschaffen. rižArtë ist Freeware und kann von der Eingangsseite dieser MJL-Homepage (links) downgeloaded werden.
Starten sie das Programm, quittieren sie die Eingangsmaske mit Esc und gehen sie in das rižArtë‑Menü ganz rechts.

Im nun erschienenen Werkzeugselektor wählen sie den Eintrag Randomizer Algorithms Analysis.

Nun erscheint ein Window, in dem Sie verschiedene Generatorenalgorithmen simulieren und testen können.

Im oberen Teil des Windows können Sie definieren, ob Sie den C-Standard- oder einen selbstdefinierten Generator verwenden wollen. Im Feld # %-Grid stellen Sie ein, auf wie viele unterschiedliche Werte die Ergebnisse projiziert werden. Mit dem Button Test wird eine Simulation gestartet. # of Generations bestimmt die Anzahl von Simulationsdurchläufen.
In dem Feld Modifier Function wird die Modifikationsfunktion definiert.
Betätigen sie den Button Select f(x) um zwischen den zur Verfügung gestellten Modifikationsfunktionen wählen zu können.

Im Feld Formula f(x) wird die nun eingestellte Modifikationsfunktion angezeigt.
Die Felder « und » ermöglichen schnelles Manövrieren ohne den Selektor zu betätigen.

Die Funktionen enthalten alle einen konstanten Parameter c, der den Kurvenverlauf modifiziert. Normalerweise kann sich c im Bereich zwischen 0.00001 und 10000.0 bewegen. Bei Funktionen mit eingeschränktem Wertebereich für c befindet sich dessen Einschränkungsdefinition in eckiger Klammer am Ende der Formel.

Das farbige Ergebnisfeld zeigt in Blau die Minima und Maxima aller Werte, die in der Darstellung auf diese Pixelbreite projizieren. Die roten Punkte sind die jeweiligen pixelbreiten Durchschnittswerte. In Gelb werden die Prozentzahlen entsprechend dem eingestellten Grid angezeigt. Im grünen Rand werden oben die Gesetzt-Häufigkeiten für einzelne Bits (vor der Modifikationsfunktion – für das Testen von my_rand()) und links das jeweilige Minimum, Maximum und der Durchschnitt über alle Werte angezeigt.

Nach dem Exekutieren eines einzelnen Testdurchlaufs erscheint die Parametrierung und das Ergebnis folgendermaßen:

Während diese Anzeige aktiv ist, befinden sich die Ergebnisse als Text am Clipboard und können in einen beliebigen Texteditor eingefügt werden.

Besonderheiten einzelner Funktionen

Die erste Formel liefert eine der Gauß'schen Normalverteilung ähnelnde Wahrscheinlichkeitskurve mit einer unsymmetrischen Glocke größter Häufigkeit. (Für die Erzeugung einer tatsächlichen Normalverteilung siehe http://de.wikipedia.org/wiki/Normalverteilung)
Wählen Sie den zweiten Eintrag x^c und stellen Sie den Wert von c auf 1.0, um einen selbstparametrierten my_rand()-Generator zu testen.
Funktion drei ersetzt die eingehenden Werte 0.0<=x<1.0 durch zwei Geradenstücke zwischen den Punkten P0 (0.0/0.0) und P1 (1.0/1.0) und erzeugt damit zwei konstante unterschiedliche Häufigkeiten.
Die Funktionen sind so normiert, dass sie den Wertebereich 0.0<=x<1.0 abbilden auf 0.0<=y<1.0. Für die hyperbolischen Funktionen mussten besondere Normierungen vorgenommen werden. Das Zeichen ¶ steht für die Zahl Pi π=3.14159…. Das Zeichen ^ steht für Potenzieren.

Sourcecode

Rückprojektion auf den DWORD-Wertebereich.

// MAX_CALC_VALS  ... oberster zulässiger Wert des Ausgangsdoppelwortes [0x7FFF]
DWORD Calcul(double y) { return MAX_CALC_VALS&(int)((double)MAX_CALC_VALS*y); }
DWORD CalcLi(DB x,DB c)
{ DB d=c/2.0; return Calcul( (x<d)? ((1-d)/d*x) : ((1-d) + d/(1-d)*(x-d)) ); }

Modifikationsfunktionen (in der gleichen Reihenfolge wie im Selektor):

// par_m_cst      ... Funktionsparametrierungskonstante
// fact           ... Eingangsvariable
#define PI    3.141592653589732384626433832795028841971e0 // PI
#define EULER 2.7182818284590452353602874713527           // Euler'sche Zahl
switch(methMode)
{ case METH_NORM: calc=Calcul(0.5*sqrt(-log((fact+exp(-4*par_m_cst)) /
                                             (1.0+exp(-4*par_m_cst)))/par_m_cst)); break;
  case METH_POTE: calc=Calcul(pow(fact,                  par_m_cst));              break;
  case METH_LINE: calc=CalcLi(fact,                      par_m_cst);               break;
  case METH_C_PY: calc=Calcul(((pow(par_m_cst,fact)-1.0)/(par_m_cst-1.0)));        break;
  case METH_LN_X: calc=Calcul( (log(par_m_cst+fact)-log(par_m_cst))/
                                (log(par_m_cst+1.0)-log(par_m_cst)));              break;
  case METH_SQR1: calc=Calcul((sqrt(1.0+pow(fact,       par_m_cst))-1.0) /
                              (sqrt(2.0)                           -1.0));         break;
  case METH_SQR2: calc=Calcul((sqrt(1.0+par_m_cst*fact*fact)-1.0) /
                              (sqrt(1.0+par_m_cst)          -1.0));                break;
  case METH_SINE: calc=Calcul(pow(sin(PI/2.0* fact),    par_m_cst));               break;
  case METH_COSI: calc=Calcul(pow(0.5+0.5*cos(PI*fact), par_m_cst));               break;
  case METH_TANG: calc=Calcul(pow(tan(PI/4.0* fact),    par_m_cst));               break;
  case METH_ASIN: calc=Calcul(pow(2.0/PI*asin(fact),    par_m_cst));               break;
  case METH_ACOS: calc=Calcul(pow(2.0/PI*acos(fact),    par_m_cst));               break;
  case METH_ATAN: calc=Calcul(pow(4.0/PI*atan(fact),    par_m_cst));               break;
  case METH_SINP: calc=Calcul(sin(PI/2.0* pow(fact,     par_m_cst)));              break;
  case METH_COSP: calc=Calcul(cos(PI/2.0* pow(fact,     par_m_cst)));              break;
  case METH_TANP: calc=Calcul(tan(PI/4.0* pow(fact,     par_m_cst)));              break;
  case METH_ASIP: calc=Calcul(2.0/PI*asin(pow(fact,     par_m_cst)));              break;
  case METH_ACOP: calc=Calcul(acos (  2.0*pow(fact,     par_m_cst)-1.0)/PI);       break;
  case METH_ATAP: calc=Calcul(4.0/PI*atan(pow(fact,     par_m_cst)));              break;
  case METH_SINH: calc=Calcul(sinh(par_m_cst*fact)/sinh(par_m_cst));               break;
  case METH_COSH: calc=Calcul((exp(par_m_cst*fact)+exp(-par_m_cst*fact)-2.0)/
                              (exp(par_m_cst     )+exp(-par_m_cst     )-2.0));     break;
  case METH_TANH: calc=Calcul(tanh(par_m_cst*fact)/tanh(par_m_cst));               break;
}

Bei Auswahl einer der obigen Funktionen können aus Gründen der Laufzeitersparnis Ausdrücke wie
(sqrt(1.0+par_m_cst)-1.0) oder
(exp(par_m_cst)+exp(-par_m_cst)-2.0) durch Konstante ersetzt werden.
log() bezeichnet nicht den dekadischen, sondern den natürlichen Logarithmus.

Falls Sie eigene Modifikationsfunktionen in das Testtool eingebracht haben wollen, wenden Sie sich bitte an rizarte@yahoo.de.