Fast jede Programmiersprache bietet einen eingebauten
Zufallszahlengenerator, der als Funktionsaufruf verfügbar ist.
Was ist das eigentlich? Wozu benötigt man ihn, und wie
funktioniert er? Diese Fragen sollen wir in diesem Artikel
beantwortet werden.
Das wohl jedem bekannte Vorbild ist der Würfel, der bei jeder
Benutzung ein nicht vorhersagbares Ergebnis liefert, natürlich in
gewissen Grenzen. Das bedeutet, seine Werte sind ganze Zahlen
zwischen 1 und 6. In Verallgemeinerung können wir also sagen: Ein
Zufallszahlengenerator ist eine Funktion, die bei jedem Aufruf ein
nicht vorhersagbares Ergebnis liefert. Lediglich den Bereich, aus
dem die Werte stammen, können wir festlegen.
Ein Anwendungsschwerpunkt sind sicher die ungezählten Computer-
spiele. Der Zufallszahlengenerator sorgt dafür, daß beim Action-
Spiel der Gegner immer an anderer Stelle und zu einem anderen
Zeitpunkt auftaucht. Und er ist dafür verantwortlich, daß ein
Schachprogramm bei jedem Spiel eine andere Eröffnung wählt.
Auch in Wissenschaft und Technik werden Zufallszahlengeneratoren
angewandt. Denn Simulationsprogramme erweisen sich als genauso
zuverlässig wie Realtests, sind aber wesentlich billiger. Eine
Autofirma kostet z.B. jeder einzelne Crashtest ein Auto. Ein
Computerprogramm, das die Versuchsbedingungen genau simuliert,
bedeutet hingegen eine einmalige Ausgabe. Zum selben Preis kann
das Experiment dann beliebig oft wiederholt werden. Zufallszahlen-
generatoren bewirken kleine Unregelmäßigkeiten, die in der Natur
auch auftreten. Ohne sie würde jede Simulation exakt gleich
ablaufen, und das Programm wäre wertlos.
Unter dem Stichwort "Monte-Carlo-Methoden" bedient sich nun auch
die moderne Mathematik der Zufallszahlengeneratoren. In der
Primzahltheorie besteht eine Aufgabe darin, festzustellen, ob eine
gegebene Zahl eine Primzahl ist. Hier wurden in letzter Zeit
Methoden entwickelt, die auf der Benutzung von Zufallszahlen
basieren.
Was ist ein Zufallszahlengenerator?
Fragen wir uns nun, was wir unter "Zufallszahlen" überhaupt
verstehen. Einer einzelnen Zahl können wir sicherlich nicht die
Eigenschaft "zufällig" zuschreiben. Oder wer wollte behaupten, die
Zahl 2 wäre zufällig, die Zahl 3 hingegen nicht? Sinnvoll wird der
Begriff des Zufalls wohl erst, wenn wir mehrere Zahlen betrachten.
Die Zahlenfolge 1, 2, 3, 4, 5... würden wir in dieser Reihenfolge nicht
als zufällig annehmen. Denn offensichtlich käme als nächstes die 6
an die Reihe. Und wir hatten ja bereits festgestellt, daß alle
Zufallszahlen "unvorhersagbar" sein sollen. Die Zahlenfolge 4, 7, 1, 2,
6... hingegen werden wir als zufällig annehmen. Denn eine Regel,
wie diese zustandekommen, ist nicht abzulesen. So könnte als
nächstes die 3 kommen, aber genauso die 5. Möglich wäre auch eine
der bereits dagewesenen Zahlen.
Betrachten wir also Folgen von Zahlen und unter diesen die
"Zufallsfolgen". Welche Forderungen werden wir an eine beliebige
Folge stellen, um sie als Zufallsfolge zu akzeptieren? Wichtig
ist, daß die jeweils nächste Zahl unvorhersagbar ist. Genauso
wichtig ist aber, daß alle möglichen Werte mit derselben
Wahrscheinlichkeit auftreten. Auf das Beispiel des Würfels
übertragen: Die 6 muß genauso wahrscheinlich sein wie die 1.
Würfeln wir sehr oft, müssen alle Zahlen etwa gleich oft fallen.
Dies ist eine Eigenschaft, die wir sehr einfach testen können. Der
Ausdruck "etwa gleich oft" ist aber nicht genau bestimmt. Würfeln
wir 1200-mal, müßte jede Zahl "etwa" 200-mal auftreten. Das wird
aber wohl nie exakt der Fall sein. Doch auch wenn eine Zahl nur
180-mal vorkommt, eine andere 220-mal, können wir von einem
"fairen" Würfel sprechen. Beim nächsten Mal kann das ja genau
umgekehrt aussehen. Bei einem Verhältnis von 380 zu 20 hingegen
wären Zweifel angebracht.
Woher kommen Zufallszahlen?
Unsere nächste Frage gilt der Art und Weise, wie wir zu
vernünftigen Zufallsfolgen gelangen können. Neben dem Werfen eines
Würfels gibt es ähnliche gleichwertige Methoden, etwa das Werfen
einer Münze oder die Benutzung eines Roulettes. Andere
Möglichkeiten wären z.B. das Entnehmen einer Folge von Zahlen aus
dem Telefonbuch. Hierzu müßten wir dieses jeweils mit
geschlossenen Augen aufschlagen und mit dem Finger an eine
beliebige Stelle tippen. Die dort stehende Zahl wäre in die Folge
zu übernehmen. Oder wir könnten uns die Nummern der
vorbeifahrenden Autos notieren.
Während die ersten drei Methoden zu einem brauchbaren Resultat
führen, sind die Telefon- und Auto-Zufallszahlen mit Vorsicht zu
genießen. So gibt es in Großstädten oft nur Telefonnummern mit 6
oder mehr Stellen. Kleinere Zahlen würden in dieser Folge somit
nicht auftreten. Bei den Autonummern wäre es genau umgekehrt: Hier
werden große Zahlen seltener ausgegeben.
Aufgrund derartiger Schwierigkeiten veröffentlichte L.H.C.Tippett
bereits 1927 eine Tafel mit 40000 Zufallszahlen, die er auf
irgendeine Weise aus den Daten von Umfragen ermittelt hatte.
Später wurden Maschinen konstruiert, deren Aufgabe lediglich im
Berechnen von Zufallsfolgen bestand. Resultate waren ähnliche
Tafeln von teilweise erheblich größerem Umfang.
Heute berechnen wir unsere Zufallsfolgen direkt mit dem Rechner.
Wie gesagt, stellen die meisten Programmiersprachen bereits Generatoren
zur Verfügung. Das ist jedoch nicht immer der Fall. Und dann sind
wir gezwungen, Zufallsfolgen selbst nachzubilden. Wie das geht, wollen
wir uns nun ansehen.
Erzeugung von Zufallsfolgen
Meist werden Zufallsfolgen rekursiv definiert. Das bedeutet, wir
suchen uns einen beliebigen Anfangswert. Das zweite Element der
Folge wird durch einen bestimmten Ausdruck aus diesem Anfangswert
berechnet, das dritte Element aus dem zweiten, das vierte aus dem
dritten, usw. Wir müssen uns somit immer den aktuellen Wert
merken.
Jedoch ergibt sich hierbei ein Problem: Die Folge wird irgendwann
zyklisch. D.h. irgendwann erreichen wir einen Wert, der vorher
schon einmal da war. Und dann wiederholt sich dieser gesamte
Abschnitt der Folge beliebig oft. In der Praxis versucht man
daher, die Periodenlänge zu maximieren. Genau untersucht hat man
die sog. "lineare Kongruenz". Hier wird der Folgenwert s(n+1) wie
folgt berechnet: s(n+1)=(a*s(n)+c) mod m mit verschiedenen
Konstanten a, c und m. Durch die Berechnung "mod m" können wir im
Höchstfall eine Periodenlänge von m erreichen. Dies ist auch
möglich, wie der folgende Satz aussagt:
Die lineare Kongruenz s(n+1)=(a*s(n)+c) mod m hat genau dann die
Periodenlänge m, wenn gilt: 1.) ggT(c,m)=1.
2.) b=a-1 ist Vielfaches jedes Primfaktors von m.
3.) Ist m Vielfaches von 4, dann auch b.
|
Als Beispiel wählen wir m=45. Um eine Periodenlänge von 45 zu
ereichen, muß der größte gemeinsame Teiler von m und c gleich 1
sein. Wir wählen für c also einfach eine Primzahl, etwa c=7. Nun
müssen wir die Primfaktoren von m bestimmen, das sind 3 und 5. m
ist nicht Vielfaches von 4, daher ist Punkt 3 des Satzes nicht zu
beachten. Die Zahl b muß somit Vielfaches von 3*5=15 sein, z.B.
b=30. Es ergibt sich a=31. Die Folge ist somit s(n+1)=(31*s(n)+7)
mod 45. Anfangswert darf jede beliebige Zahl von 0 bis 44 sein.
Wir wollen unsere Möglichkeiten natürlich voll ausschöpfen, eine
Periodenlänge von 45 ist uns daher viel zu wenig. Wir wollen mit
Long Integers rechnen. Diese werden meist mit 32 Bits dargestellt,
es gibt somit 2^32 verschiedene Werte. Wir können für unseren
Generator also eine Periodenlänge von 2^32 anstreben. Wenden wir
den Satz an, müssen wir m=2^32 setzen. c darf dann kein Vielfaches
von 2 sein, wir wählen 2001. Der einzige Primfaktor von m ist 2.
Punkt 3 muß nun berücksichtigt werden, da 4 Teiler von m ist.
Somit können wir für b jedes Vielfache von 4 nehmen, z.B. 9012.
Damit haben wir a=9013, und die gesuchte Rekursionsformel lautet:
s(n+1)=(9013*s(n)+2001) mod 2^32.
Die Realisierung
Um dies zu realisieren, müssen wir auf die Programmiersprache
eingehen, die wir verwenden wollen. Wir müssen uns im Programm
immer den letzten Wert der Folge merken. Dies ist in vielen
Sprachen nur durch eine globale Variable möglich, was eine
unschöne Lösung ist. Denn hierdurch kann es zu Konflikten mit den
Variablen im Hauptprogramm kommen. Die Sprache C hingegen bietet
mit den sog. "statischen" Variablen die Möglichkeit, die Werte
eines Unterprogramms zwischen zwei Aufrufen aufzubewahren.
Gleichzeitig mit der Vereinbarung als statische Variable muß noch
eine Initialisierung stattfinden. Wir können hierfür durchaus den
Wert 1 nehmen. Die Anweisung muß dann lauten: "static unsigned
long s=1". Beim ersten Aufruf der Funktion wird diese
Initialisierung ausgeführt, später nicht mehr.
Und so sieht das aus:
float linear()
{
static unsigned long s=1;
float r;
s=9013*s+2001;
r=s/4+0x20000000;
r=r/0x40000000;
return(r);
}
|
Nach der Vereinbarung von s wird eine float-Variable r erklärt, in der das
Resultat zwischengespeichert wird, welches wir zurückgeben wollen.
Nun kommt die Rekursionsformel zur Anwendung. Hier könnte man
überrascht sein, daß explizit keine Rechnung "mod 2^32" zu sehen
ist. Das liegt einfach daran, daß diese automatisch erfolgt. Denn
bei Überschreiten des zulässigen Rechenbereichs läßt C einfach die
höchstwertigen Bits weg, was der von uns gewünschten Modulo-
Rechnung gleichkommt. Aus diesem Grund ist es sinnvoll, gerade
m=2^32 zu setzen. Hieraus ergibt sich auch die Konsequenz, daß die
von uns gewählten Konstanten a, c und m nur dann sinnvoll sind,
wenn Long Integers tatsächlich in 32 Bits gespeichert werden. In
allen anderen Fällen müssen andere Konstanten zum Einsatz kommen.
Die nun folgenden Anweisungen werden notwendig, weil wir reelle
Zufallszahlen zwischen 0 und 1 zurückgeben wollen. Dies ist
günstiger, da ohnehin meist eine Umrechnung stattfinden muß. Dann
kommt man nicht um das Rechnen mit float-Zahlen herum. Es findet
also eine Transformation auf diesen Bereich statt. Schließlich
wird der Wert von r zurückgegeben.
Noch ein Wort zum Anfangswert s=1. Natürlich kann hier jeder
beliebige Wert stehen. Oft kann es auch erforderlich sein, den
Wert von außen zu setzen, um zu unterschiedlichen Folgen zu
gelangen. Dann müßte dieser Wert als Parameter übergeben werden.
Wir können die Qualität dieses Zufallszahlengenerators leicht
testen, indem wir eine größere Anzahl von Werten der Folge
ausrechnen lassen. Den Bereich zwischen 0 und 1 teilen wir in
gleich große Intervalle ein und zählen, wie oft ein Wert dieser
Folge in jedem Intervall liegt. Es werden nicht zu große
Abweichungen auftreten.
Hier ein entsprechendes Testprogramm:
#define ANZ 10000 /* Anzahl Tests */
#define GEN 100 /* Genauigkeit */
main()
{
int a[GEN], i, k, max, min;
long j;
float r, linear();
for ( i=0; imax) max=a[i];
}
printf("Maximum: %d, Minimum: %d\n",max,min);
printf("Abweichung: %f\n",(max-min)/2.0/ANZ*GEN);
}
|
Hier können zwei Parameter gesetzt werden: GEN gibt an, in wieviele
Intervalle der Wertebereich eingeteilt wird, ANZ ist die Anzahl der
durchzuführenden Tests. Am Ende werden Maximal- und Minimalwert
ausgegeben sowie die relative Abweichung.
Nun könnte man versucht sein, zu sagen: Die lineare Kongruenz ist
viel zu einfach aufgebaut, um einen guten Zufallszahlengenerator
abzugeben. Die verwendete Formel müßte viel komplizierter sein. Um
hier ein Vergleichsbeispiel zu haben, ist im folgenden ein anderer
Generator zu finden. Er arbeitet auch rekursiv, ermittelt aber den
jeweils nächsten Wert durch eine Mischung aus Textverschiebung und
Berechnungen.
float seltsam()
{
static unsigned long s=1;
float r;
char c[11], h;
sprintf(c,"%010ld",s);
h=c[1]; c[1]=c[8]; c[8]=h;
h=c[2]; c[2]=c[4]; c[4]=h;
h=c[5]; c[5]=c[7]; c[7]=h;
c[1]=48+c[2]-c[3]; if (c[1]<48) c[1]+=10;
c[4]=48+c[5]-c[6]; if (c[4]<48) c[4]+=10;
c[7]=48+c[8]-c[9]; if (c[7]<48) c[7]+=10;
c[3]++; if (c[3]>57) c[3]-=10;
c[6]++; if (c[6]>57) c[6]-=10;
c[9]++; if (c[9]>57) c[9]-=10;
c[0]=(c[3]+c[6]+c[9]-144)%10+48;
sscanf(c,"%10ld",&s);
r=s/4+0x20000000;
r=r/0x40000000;
return(r);
}
|
Wir arbeiten ebenfalls mit dem Datentyp unsigned long. Im
folgenden bezeichnen wir die höchstwertige Stelle als Stelle 0,
die niedrigste mit 9. In jedem Schritt vertauschen wir also die
zweite mit der vierten Stelle, die fünfte mit der siebten und die
erste mit der achten. Sodann addieren wir die achte und neunte
Stelle und schreiben das Ergebnis in die siebte Stelle. Gleiches
führen wir für die Stellen 6, 5 und 4 sowie 3, 2 und 1 durch. Die
Stelle 0 ist in jedem Schritt die Summe der Stellen 3, 6 und 9.
Weil es sich jeweils um Ziffern handelt, wird selbstverständlich
mod 10 gerechnet.
Realisiert ist dieser Algorithmus mit Hilfe der Funktionen sprintf
und sscanf, welche Formatumwandlungen im Speicher vornehmen
können. Es zeigt sich nun, daß trotz dieser "Kompliziertheit" die
lineare Kongruenz bessere Resultate liefert, was die gleichmäßige
Verteilung anbetrifft. Dies können wir mit dem Programm aus
Listing 2 ja leicht testen.
Veränderung des Werte-Intervalls
Wir verfügen also über einen zufriedenstellenden Generator, der
uns Zufallsfolgen mit gleichmäßiger Verteilung liefert. Diese sind
immer reelle Zahlen zwischen 0 und 1. Wir sollten nun noch
erwähnen, wie man den Bereich verändert, in dem die Werte der
Zufallsfolge liegen. Wir suchen also Werte zwischen a und b, wobei
a<b gelten muß. Dazu müssen wir das Intervall [0,1] auf die Länge
des Intervall [a,b] strecken. Das geschieht durch Multiplikation
mit b-a. Nach dieser Operation befinden sich die Werte im
Intervall [0,b-a]. Nun wird noch die Addition von a erforderlich.
Liefert ein Generator random() also Werte zwischen 0 und 1, so
ergibt der Ausdruck a+(b-a)*random() Werte zwischen a und b.
Suchen wir Integer-Werte von a bis b, müssen wir eine "cast"-
Umwandlung benutzen. Der notwendige Ausdruck hat dann die Form
(int)(a+(b-a)*random()). Die Simulation eines Würfels oder eines
Roulettes ist nun sehr leicht möglich. Für den Würfel würde der
Ausdruck (int)(1+6*random()) in Frage kommen, für das Roulette
(int)(36*random()).
Wir sind bis jetzt also in der Lage, Ereignisse, die in
gleichmäßiger Verteilung auftreten, zu simulieren. Hierzu haben
wir einen leistungsfähigen Zufallszahlengenerator. Leider tritt
eine gleichmäßige Verteilung in der Natur aber höchst selten auf.
Nehmen wir als Beispiel hierzu einen 100-Meter-Läufer, der im
Durchschnitt aller Läufe seine Strecke in 11.0 Sekunden bewältigt.
Seine Bestzeit liegt bei 10.7, während man als schlechtesten Wert
bei seiner momentanen Form 11.3 annehmen muß.
Dieser Mann könnte nun unseren Zufallszahlengenerator befragen,
wie wohl das Resultat seines nächsten Laufes lauten wird. Die
hierzu nötige Formel wäre 10.7+0.6*random(). Da wir aber eine
gleichmäßige Verteilung programmiert haben, wäre die Zeit von 10.7
genauso wahrscheinlich wie 11.0. Dies widerspricht jedoch unserer
Erfahrung. Die Zeit 11.0 wird er viel häufiger erzielen als seine
Bestzeit von 10.7. Statistiker haben mehrere Modelle entwickelt, um
derartige Verteilungen beschreiben zu können. Sie sprechen von
"binomialer" oder "normaler" Verteilung.
Natürliche Verteilung
Wir wollen uns im folgenden von diesen exakt definierten
mathematischen Begriffen freimachen. Wir werden von einer
"natürlichen" Verteilung sprechen und damit eine Verteilung
meinen, bei der Werte einer Zufallsfolge in einem festgelegten
Intervall auftreten. Im Unterschied zur gleichmäßigen Verteilung
gibt es irgendwo in diesem Intervall einen Punkt, in dem sich die
auftretenden Werte häufen, während sie an anderer Stelle seltener
erscheinen. Was das genau heißen soll, wollen wir gar nicht
festlegen.
Um hierfür einen Zufallszahlengenerator zu programmieren, wollen
wir unsere lineare Kongruenz benutzen, die erhaltenen Werte aber
so transformieren, daß die gewünschten Eigenschaften zu Tage
treten. Es wird nun etwas mathematisch. Wer das nicht mag, für den
ist der Artikel hier fast zu Ende. Im vorletzten und letzten
Absatz sind die Resultate angegeben.
Wir benötigen eine Funktion, die reelle Werte zwischen 0 und 1 auf
ebensolche abbildet, wobei an gewisser Stelle eine Häufung
auftritt. Eine Funktion wie in Bild 2 ist hierfür geeignet. Diese
sorgt für eine Häufung um den Punkt 0.5 herum.
Es handelt sich um ein Polynom dritten Grades. Es hat die Form
f(x)=a*x^3+b*x^2+c*x+d. Hierbei sind die Koeffizienten a, b, c und
d noch zu bestimmen. Um dies zu tun, müssen wir die notwendigen
Bedingungen angeben. Zunächst soll gelten: f(0)=0 und f(1)=1.
Damit sind aber erst zwei Koeffizienten festgelegt.
Zwei weitere Bedingungen haben wir somit zur freien Verfügung.
Diese wollen wir benutzen, um unsere "natürliche" Verteilung
variieren zu können. Wir bringen die Parameter p und q ins Spiel.
Sie sollen angeben, wie groß die Steigung der Kurve in den Punkten
0 und 1 ist. Durch Variation beider Werte lassen sich sowohl der
Ort wie das Ausmaß der Häufung beeinflussen. f'(0)=p und f'(1)=q
lauten die Bedingungen.
Aus den ersten zwei Bedingungen ergeben sich nun die Gleichungen
d=0, und a+b+c+d=1. Für die anderen Bedingungen müssen wir die
Ableitung berechnen. Diese ergibt sich als: f'(x)=3*a*x^2+2*b*x+c.
Somit erhalten wir hieraus die Gleichungen c=p und 3*a+2*b+c=q.
Insgesamt haben wir also ein Gleichungssystem mit vier
Unbekannten. Neben d=0 und c=p, welches sofort abzulesen ist,
ergibt sich b=3-q-2*p und a=q+p-2. Die gesuchte Funktion lautet
somit: f(x)=(q+p-2)*x^3+(3-q-2*p)*x^2+p*x.
Zur schnelleren Berechnung stellen wir dies schließlich um:
f(x)=(((q+p-2)*x+(3-q-2*p))*x+p)*x. Somit ersparen wir uns einige
Multiplikationen. Nun noch zur praktischen Auswirkung der beiden
Parameter p und q: Beide sollten im Bereich von 0 bis 3 benutzt
werden, da sie sonst Werte außerhalb des Intervalls liefern. Gilt
p=q, ergibt sich eine Häufung bei 0.5. Ist p<>q, ist die Häufung
zum kleineren Wert hin verschoben. Wählt man p=q=3, erhält man
eine maximale Häufung um den Punkt 0.5, bei p=3, q=0 eine maximale
Häufung um den Punkt 1.
Hier ist die Realisierung dieses Zufallszahlengenerators mit natürlicher
Verteilung:
float normal(p,q)
float p,q;
{
float x, linear();
x=linear();
return (((q+p-2)*x+(3-q-2*p))*x+p)*x;
}
|
Wie besprochen, benutzen wir hierbei unsere oben programmierte lineare Kongruenz.
Beim Aufruf sind p und q als Argumente anzugeben. Eine Verschiebung der Werte
in ein anderes Intervall wird hier genau wie oben durchgeführt.
Literatur
- Donald E.Knuth: The Art of Computer Programming (Volume 2), Chapter 3: Random Numbers