Simulation von Gleichverteilungen

Ein (Pseudo-) Zufallszahlengenerator ist ein Algorithmus, der eine deterministische Folge von ganzen Zahlen (Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),...  mit Werten zwischen 0 und einem Maximalwert m-1 erzeugt, welche durch eine vorgegebene Klasse statistischer Tests nicht von einer Folge von Stichproben unabhängiger, auf {0,1,2,...,m-1} gleichverteilter Zufallsgrößen unterscheidbar ist. Ein Zufallszahlengenerator erzeugt also nicht wirklich zufällige Zahlen. Die von "guten" Zufallszahlengeneratoren erzeugten Zahlen haben aber statistische Eigenschaften, die denen von echten Zufallszahlen in vielerlei (aber nicht in jeder) Hinsicht sehr ähnlich sind.

Konkret werden die Pseudozufallszahlen üblicherweise über eine deterministische Rekurrenzrelation vom Typ

(Subscript[x, n+1])=f((Subscript[x, n-k+1]),(Subscript[x, n-k+2]),...,(Subscript[x, n])),       n=k,k+1,k+2,...,

aus "Saatwerten"  (Subscript[x, 1]),(Subscript[x, 2]),...,(Subscript[x, k]) erzeugt. In vielen Fällen hängt die Funktion f nur von der letzten erzeugten Zufallszahl (Subscript[x, n]) ab. Wir betrachten einige Beispiele:

Lineare Kongruenzgeneratoren (LCG)

Bei linearen Kongruenzgeneratoren ist die Rekurrenzrelation vom Typ

(Subscript[x, n+1])=(a (Subscript[x, n])+c) mod m,       n=0,1,2,... .

Hierbei sind a, c und m gegeignet zu wählende positive ganze Zahlen, z.B.:

ZX81 Generator:               m=(2^16)+1, a=75, c=0

RANDU, IBM 360/370:     m=(2^31), a=65539, c=0

Marsaglia Generator:          m=(2^32), a=69069, c=1

Langlands Generator:         m=(2^48), a=142412240584757, c=11

Um einen ersten Eindruck zu erhalten, wie die Qualität der erzeugten Pseudozufallszahlen von a, c und m abhängt, implementieren wir die Generatoren mit Mathematica:

f[x_]:=Mod[a x +c,m]

Beispiel 1

Wir beginnen zur Demonstration mit dem Beispiel eines ganz schlechten LCG:

a=11;c=0;m=63;
pseudorandomdata=NestList[f,1,300];
ListPlot[pseudorandomdata]

Die Folge von Zufallszahlen ist in diesem Fall periodisch mit einer Periode, die viel kleiner ist als die maximal mögliche (63). Dies rechnet man auch leicht nach.

Periodizität mit Periode kleiner als m kann man leicht ausschliessen. Es gilt nämlich:

Satz (Knuth).  Die Periode eines LCG ist gleich m genau dann, wenn
                       (i) c und m teilerfremd sind,
                       (ii) jeder Primfaktor von m ein Teiler von a-1 ist, und
                       (iii) falls 4 ein Teiler von m ist, dann auch von a-1.

Beweis siehe Knuth: The art of computer programming, Vol. 2.

Beispiel 2:  ZX 81 Generator

Hier ergibt sich ein besseres Bild, solange wir nur die Verteilung der einzelnen Zufallszahlen betrachten:

a=75;c=0;m=2^16+1;
pseudorandomdata=NestList[f,1,30000];
ListPlot[pseudorandomdata]

Fassen wir jedoch Paare ((Subscript[x, i]),(Subscript[x, i+1])) von aufeinanderfolgenden Pseudozufallszahlen als Koordinaten eines zweidimensionalen Pseudozufallsvektors auf, und betrachten die empirische Verteilung dieser Vektoren, so ergibt sich keine besonders gute Approximation einer zweidimensionalen Gleichverteilung:

blocks=Partition[pseudorandomdata,2]; ListPlot[blocks]

Beispiel 3:  RANDU

Hier scheinen sowohl die einzelnen Pseudozufallszahlen (Subscript[x, i]) als auch die Vektoren ((Subscript[x, i]),(Subscript[x, i+1]))näherungsweise gleichverteilt zu sein:

a=65539;c=0;m=2^31;
pseudorandomdata=NestList[f,1,30000];
ListPlot[pseudorandomdata]

blocks=Partition[pseudorandomdata,2]; ListPlot[blocks]

Fassen wir aber jeweils drei aufeinanderfolgende Pseudozufallszahlen als Koordinaten eines Vektors ((Subscript[x, i]),(Subscript[x, i+1]),(Subscript[x, i+2])) im (ℤ^3)auf, dann ist die empirische Verteilung dieser Pseudozufallsvektoren keine Gleichverteilung mehr, sondern konzentriert sich auf nur 15 zweidimensionalen Hyperebenen:

blocks3=Partition[pseudorandomdata,3]; ListPointPlot3D[blocks3]

Beispiel 4:  Marsaglia Generator

Der von Marsaglia 1972 vorgeschlagene LCG besteht dagegen alle obigen Tests (und einige weitere):

a=60069;c=1;m=2^32;
pseudorandomdata=NestList[f,1,30000];
ListPlot[pseudorandomdata]

blocks=Partition[pseudorandomdata,2]; ListPlot[blocks]

blocks3=Partition[pseudorandomdata,3]; ListPointPlot3D[blocks3]

Dies bedeutet natürlich nicht, daß die vom Marsaglia Generator erzeugte Folge eine für alle Zwecke akzeptable Approximation einer Folge von unabhängigen Stichproben von der Gleichverteilung ist. Da die Folge in Wirklichkeit deterministisch ist, kann man natürlich einen Test konstruieren, der sie von einer echten Zufallsfolge unterscheidet.

Shift-Register-Generatoren

Bei Shift-Register-Generatoren interpretiert man eine Zahl (Subscript[x, n])∈{0,1,...,(2^k)-1} zunächst als Binärzahl bzw. als Vektor aus (({0,1})^k), und wendet dann eine gegebene Matrix T darauf an, um (Subscript[x, n+1]) zu erhalten:

(Subscript[x, n+1])=T (Subscript[x, n]) ,       n=0,1,2,... .

Kombination von Zufallszahlengeneratoren

Zufallszahlengeneratoren lassen sich kombinieren, z.B. indem man die von mehreren Zufallszahlengeneratoren erzeugten Folgen von Pseudozufallszahlen aus {0,1,...,m-1} modulo m addiert. Auf diese Weise erhält man sehr leistungsfähige Zufallszahlengeneratoren, z.B. den Kiss Generator von Marsaglia, der einen LCG und zwei Shift-Register-Generatoren kombiniert, Periode (2^95) hat, und umfangreiche statistische Tests besteht.

Zufallszahlen aus [0,1)

Ein Zufallszahlengenerator kann natürlich nicht wirklich reelle Pseudozufallszahlen erzeugen, die die Gleichverteilung auf dem Intervall [0,1) simulieren, denn dazu würden unendlich viele "zufällige" Nachkommastellen benötigt. Stattdessen werden üblicherweise (pseudo-) zufällige Dezimalzahlen vom Typ

(Subscript[u, n])=((Subscript[x, n])/m) ,          (Subscript[x, n])ϵ{0,1,...,m-1},

erzeugt, wobei m vorgegeben ist (z.B. Darstellungsgenauigkeit des Computers), und (Subscript[x, n]) eine Folge ganzzahliger Pseudozufallszahlen aus {0, 1, ..., m - 1} ist. In Mathematica kann man mit

RandomReal[spec,WorkingPrecision->k]

pseudozufällige Dezimalzahlen mit einer beliebigen vorgegebenen Anzahl k von Nachkommastellen erzeugen.

Zufallspermutationen

Der folgende Algorithmus erzeugt eine (pseudo-)zufällige Permutation aus (Subscript[S, n]):

Algorithmus RPERM

rperm[n_]:=
Module[{x=Range[n],k,a},
   (* Beginn mit Liste {1,2,...,n} *)
Do[
k=RandomInteger[{i,n}]; 
a=x[[i]]; x[[i]]=x[[k]]; x[[k]]=a;  (* Vertausche x[[i]] und x[[k]] *)
,{i,n-1}];
(* Schleife, i läuft von 1 bis n-1 *) 
x    (* Ausgabe von x *)]

rperm[17]

Übung:

Sei

({{({{(Subscript[Ω, n])={1,2,... ,n}×{2,3,... ,n}×⋯×{n-1,n}.}})}})

a) Zeigen Sie, daß die Abbildung

({{({{X(ω)=(Subscript[τ, n-1,ω(n-1)])∘⋯∘(Subscript[τ, 2,ω(2)])∘(Subscript[τ, 1,ω(1)])}})}})

eine Bijektion von (Subscript[Ω, n]) nach (Subscript[S, n]) ist ((Subscript[τ, i,j]) bezeichnet die Transposition von i und j).

b) Folgern Sie, daß der Algorithmus oben tatsächlich eine Stichprobe einer gleichverteilten Zufallspermutation aus (Subscript[S, n]) simuliert.