# Stochastik-Praktikum, WS 2010/2011 # Sitzung vom 24.1.2011 # Themen # (Pseudo-)Zufallsgeneratoren # Tests für (Pseudo-)Zufallszahlen # Verwerfungsmethode # Pseudozufallsgeneratoren (lin. Kongruenzen) # ######################################################### # Simulation von (Pseudo-)Zufallszahlen # Für sehr viele Simulationsprobleme benötigt man # sog. "Pseudozufallszahlen": # Beobachtungen oder Werte x[1], x[2], x[3], ... die keine # erkennbare Regelmäßigkeit haben und für die Zwecke der # Berechnung als unabhängig und zufällig generiert angenommen werden # (dürfen), d.h. wir tun beispielsweise so, als ob die x[1], x[2], ... # durch wiederholtes, unabhängiges Drehen eines Glücksrad gewonnen # worden wären -- obwohl sie aus Praktikabilitätsgründen durch einen # (deterministischen) Algorithmus (einen sog. (Pseudo-)Zufallsgenerator) # im Computer erzeugt worden sind. # # Bemerkung: # Eine "Standardreferenz" zu Pseudozufallszahlen ist Chapter 3 in # Donald E. Knuth, The art of computer programming, # Vol. 2 / Seminumerical algorithms, 3rd Ed., Addison-Wesley, 1998. # Das dort behandelte Material geht weit über den hier vorgestellten # Stoff hinaus. # Ein einfaches Beispiel sind die # linearen Kongruenzengeneratoren: # Beginne mit einem (ganzzahligen) Startwert ("Zufallssame", "random seed"). # Wenn der aktuelle Wert x ist, so ist der nächste Wert # (a*x+c) mod M # z.B. M <- 2048; a <- 65; c <- 1 # "Selbstbau"-linearer Kongruenzengenerator: linKonGen <- function() { neu <- (a*x + c) %% M # folgendes setzt den Wert der (globalen) Variable x auf den von neu assign("x", neu, .GlobalEnv) # gebe "uniform" auf (0,1] verteilten Wert aus: (neu+1)/M } # Startwert: x <- 5 N <- 5000 werte <- numeric(N) for (i in 1:N) werte[i]<-linKonGen() # Übrigens: Dasselbe leistet (wesentlich schneller) # werte <- replicate(N, linKonGen()) # Anschauen: plot(werte) lines(werte,add=TRUE) hist(werte, prob=TRUE) werte <- replicate(500, linKonGen()) # Zeichne Wert gegen Nachfolger-Wert: plot(werte, werte[c(2:500,1)]) # zum Vergleich: x11() werteR <- runif(500) plot(werteR, werteR[c(2:500,1)]) # Perioden können ein Problem sein (speziell bei kleinem M): M <- 16 #a <- 5; c<-1 a <- 4; c <-0 x <- 5 replicate(30, linKonGen()) # (in diesem Beispiel: Wert 4 = Wert 20, etc.) # Zwei Beispiele: M <- 256 a <- 17; c <- 1 # gibt relativ gleichmäßige Verteilung N <- M werte <- numeric(N) for (i in 1:N) werte[i]<-linKonGen() plot(M*werte, M*werte[c(2:N,1)]) x11() M <- 256 a <- 129; c <- 1 # gibt wenig gleichmäßige Verteilung N <- M werte <- numeric(N) for (i in 1:N) werte[i]<-linKonGen() plot(M*werte, M*werte[c(2:N,1)]) # Rs Zufallsgenerator-Hilfe ?set.seed ############################################ # # Eine Illustration zu Knuths Warnung # "... random numbers should not be generated with a method chosen # at random. Some theory should be used." # (loc. cit., S. 6, siehe auch Exercise 11, S.8) m <- 100000 f <- sample(1:m, m, replace=T) # Definiere "Zufallsfolge" durch y[n+1]=f[y[n]] # mit zufälligem Startwert y[0]. # Diese wird offenbar schließlich zyklisch, wie lang ist der # resultierende Zyklus? gesehen <- logical(m) # initialisiert mit m-mal FALSE y <- sample(1:m,1) y # der Startwert anfangslaenge <- 0 while(gesehen[y]==FALSE) { anfangslaenge <- anfangslaenge+1 z<-y gesehen[z]<-TRUE y<-f[y] } anfangslaenge # Wie lang war das Stück vor dem Zyklus? # Wir sind auf einen Zyklus gestoßen, die Werte von z und y sind # nun so, dass f[z]=y, d.h. wir haben einen Zyklus # y_0=y, y_1=f[y_0], y_2=f[y_1], ..., y_{n-1}=f[y_{n-2}]=z, y_n=f[z]=y_0. # Bestimme die Länge des Zyklus: zyklaenge <- 1 yy <- f[y] while (yy != y) { zyklaenge <- zyklaenge+1 yy <- f[yy] # cat(yy, " ") # ggf. ent-kommentieren, um Zykel anzuschauen } zyklaenge # Ggf. einige Male wiederholen. # Beobachtung: auch bei großem m läuft die Folge typischerweise # auf recht kurze Zyklen. ############################################ # # (empirische) Tests mit Zufallsgeneratoren # hier ggf. eigenen Zufallsgenerator einsetzen zufgen <- function() runif(1) # zufgen <- linKonGen # Beispielwerte: # M <- 2048; a <- 65; c <- 1 # (das Mini-Beispiel vom Anfang, versagt kläglich) # M <- 10^10; a<-3141592621; c <- 1 # M <- 2^31; a <- 65539; c <- 0 # Dies ist "RANDU", s.a. Knuth, S. 107 # M <- 2^35; a <- 2^18+1; c <- 1 # "Generator F" aus Knuth, S. 47 # ggf. auch Startwert variieren: # x <- 5; x <- 314159265 # 1. Sind die Werte gleichmäßig verteilt? # Wir benutzen den chi^2-Test: # Zerlege (0,1] in kl Klassen (0,1/kl], (1/kl,2/kl],...,((kl-1)/k1,1] # wenn wir n Werte generieren lassen, erwarten wir in jeder Klasse n/kl # Werte (mit "Zufallsfluktuationen") kl <- 5 # variieren n <- 100000 beob <- rep(0, times=kl) for (i in 1:n) { w <- ceiling(kl*zufgen()) # verwandelt Wert aus (0,1] in Wert aus {1,...,kl} beob[w] <- beob[w]+1 } beob beob/n beob/n-rep(1/kl, times=kl) # Einschätzung der Größe der Abweichung vom "theoretischen Mittelwert" # via chi^2-Statistik chiquadratwert <- sum((beob-n/kl)^2/(n/kl)) chiquadratwert # Wie wahrscheinlich wäre eine solche Abweichung für # "echte" Zufallszahlen (p-Wert des Tests)? pchisq(chiquadratwert, df=kl-1, lower.tail=F) # Dasselbe mit einem R-Befehl: chisq.test(beob) # # Analog für d-Tupel von Werten: # (sollten bei Vergröberung uniform verteilt sein auf kl^d Klassen) d <- 3 # Andere Werte einsetzen, z.B. 2, 3, 5 kl <- 6 # Andere Werte einsetzen, z.B. 2, 3, 4, 5, 6 beob <- rep(0, times=kl^d) n <- 10000 for (i in 1:n) { # gewinne eine ganze Zahl aus {0,1,...,kl^d-1} # durch Darstellung mit d Ziffern im kl-System stelle <- 1; w <- 0 for (j in 1:d) { w <- w+(ceiling(kl*zufgen())-1)*stelle stelle <- stelle*kl } w <- w+1 # (Rs array-Indizes beginnen bei 1, nicht bei 0) beob[w]<-beob[w]+1 } beob beob-n/(kl^d) beob/n-1/(kl^d) # Wie gut passt die empirische Verteilung von d-Tupeln? chisq.test(beob) # # "Lückentest": Wie lange muss man warten # (genauer: wieviele Fehlversuche), bis der # nächste Wert in [s,t] kommt? s <- 0.2; t <- 0.5 sim.wartezeit <- function() { i <- 0 repeat { u <- zufgen() if (u >=s && u <=t) break i<-i+1 } i } L <- 50000 wz <- replicate(L, sim.wartezeit()) # vergleiche solche Wartezeiten mit der geometrischen Verteilung # mit Erfolgsparameter t-s: # grafisch: if (max(wz)>=10) { klassen <- c(seq(from=-0.5, to=10.5,by=1), max(wz)+0.5) geom.gewichte <- c(dgeom(0:10, prob=t-s), pgeom(10, prob=t-s, lower.tail=F)) } else { klassen <- seq(from=-0.5, to=10.5,by=1) geom.gewichte <- dgeom(0:10, prob=t-s) } hi <- hist(wz, prob=T, breaks=klassen) points(hi$mids, geom.gewichte/(hi$breaks[-1]-hi$breaks[-length(hi$breaks)]), col="red") # und mit dem chi^2-Test: chisq.test(hi$counts, p=geom.gewichte) # # Maximalwert von t Beobachtungen: # Für U_1,..., U_t u.a. unif([0,1]) ist # P(max(U1,...,U_t)<=x)=x^t = int_0^t t*x^{t-1} dt, # d.h. max(U1,...,U_t) sollte Beta(3,1)-verteilt sein: t <- 3 # ggf. variieren max(replicate(t, zufgen())) L <- 10000 mwerte <- replicate(L, max(replicate(t, zufgen()))) # Vergleiche empirische und theoretische Verteilungsfkt. plot(ecdf(mwerte)) curve(pbeta(x, shape1=3, shape2=1), add=T, col="red") # # Befunde: # Rs eingebauter Zufallsgenerator "besteht" diese Tests, # die Beispiel-LKG haben z.T. Schwierigkeiten. ##################################################### # # Verwerfungsmethode (rejection sampling) # Zugrundeliegende Theorie # (wir betrachten der Einfachheit halber hier den diskreten Fall): # E endliche Menge, f(x), g(x) für x aus E Wahrscheinlichkeitsgewichte, # es gelte f(x)/g(x) <= C für eine Konstante C und # X sei gemäß den Gewichten g verteilt (d.h. P(X=x)=g(x)) und wir # nehmen an, wir besitzen einen Simulationsalgorithmus für X. # Sei U unabhängig von X, uniform auf [0,1] verteilt. # Generiere einen Zufallswert Y folgendermaßen: # Simuliere ein Paar (X, U) (unabhängig), # falls U <= f(X)/(C*g(X)), akzeptiere diesen Wert und setze Y:=X, # ansonsten simuliere ein neues, unabhängiges Paar (X,U), etc. # Warum funktioniert dieses Verfahren? # Es ist P(X=x, U <= f(X)/(C*g(X))) = g(x) P(U <= f(x)/(C*g(x))) # = (g(x) f(x))/(C g(x)) = f(x)/C, also # P(U <= f(X)/(C*g(X))) = \sum_x P(X=x, U <= f(X)/(C*g(X))) = 1/C # und folglich # P(Y=x) = P(X=x, U <= f(X)/(C*g(X)) | U <= f(X)/(C*g(X))) = f(x). # Wir sehen zugleich: Die Wahrscheinlichkeit, dass ein Vorschlag # erfolgreich ist, ist P(U <= f(X)/(C*g(X))=1/C, d.h. wir # müssen geometrisch mit Erfolgsparameter 1/C oft vorschlagen, # und die Methode ist umso effizienter, je näher C an 1 liegt, d.h. # je ähnlicher f und g sich sind. # Bem.: # 1. Im kontinuierlichen Fall gilt ein analoges Argument, in dem # die Wahrscheinlichkeitsgewichte durch Wahrscheinlichkeitsdichten # ersetzt werden. # 2. Wir sehen: Es kommt in obigem Argument nicht darauf an, dass # f(x) Wahrscheinlichkeitsgewichte sind, sondern nur, dass # f(x) >= 0 und \sum_x f(x) endlich: # Sei \sum_x f(x) = K != 1, so finden wir # P(U <= f(X)/(C*g(X))) = K/C und somit hat Y die # (Wahrscheinlichkeits-)Gewichte f(x)/K. # Beispiel: Ziehen ohne Zurücklegen via Ziehen mit Zurücklegen simulieren, # d.h. aus einem Zufallsgenerator für die Binomialverteilung einen # für die hypergeometrische Verteilung bauen. # M weiße, N schwarze Kugeln in Urne, ziehe n ohne Zurücklegen, # beobachte, wieviele weiße Kugeln in der Stichprobe M <- 30; N <- 20; n <- 15 p <- M/(M+N) # Die Gewichte wären: #dhyper(0:n, M, N, n) #dbinom(0:n, n, p) C <- max(dhyper(0:n, M, N, n)/dbinom(0:n, n, p)) C sim.hypgeom.VM <- function() { repeat { x <- rbinom(1, n, p) u <- runif(1) if (u <= dhyper(x, M, N, n)/(C*dbinom(x, n, p))) break } x } # Überzeugen wir uns, dass das funktioniert: wdh <- 10000 beob <- replicate(wdh, sim.hypgeom.VM()) # zunächst per Auge: hi <- hist(beob, breaks=0:(n+1)-0.5, prob=T, main="Empir. Häufigkeiten der Ergebn. d. Verwerfungsroutine") points(0:n, dhyper(0:n, M, N, n), col="red") # und per chi^2-Test: chisq.test(hi$counts, p=dhyper(0:n, M, N, n)) # und falls wir der chi^2-Approximation hier nicht trauen: chisq.test(hi$counts, p=dhyper(0:n, M, N, n), simulate.p.value=T, B=500) # Das gibt jedenfalls keinen Grund, an der Verwerfungsmethode zu zweifeln. # # Ein kontinuierliches Beispiel: # Beta-Verteilung (mit Parametern a, b >= 1) aus einem # unif([0,1])-verteilten Vorschlag generieren # Die Beta(a,b)-Verteilung (a, b > 0) hat Dichte # B(a,b) x^(a-1)*(1-x)^(b-1) auf [0, 1], wo # B(a,b) = gamma(a+b)/(gamma(a)*gamma(b)) ("Beta-Funktion", s.a. ?beta) a <- 1.5; b <- 2 curve(x^(a-1)*(1-x)^(b-1), xlim=c(0,1)) # nimmt ihr Maximum bei x=(a-1)/(a+b-2), # wie man leicht durch Ableiten prüft abline(v=(a-1)/(a+b-2), col="blue") dbeta.unnorm <- function(x) x^(a-1)*(1-x)^(b-1) C <- dbeta.unnorm((a-1)/(a+b-2))/1 sim.beta.VM0 <- function() { repeat { x <- runif(1); u <- runif(1) if (u <= dbeta.unnorm(x)/C) # bemerke: Es genügt, die "Zieldichte" # bis auf Normierung zu kennen break } x } # Überzeugen wir uns wiederum, dass das funktioniert: wdh <- 10000 beob <- replicate(wdh, sim.beta.VM0()) # zunächst per Auge: klassen <- seq(from=0, to=1, by=0.05) hi <- hist(beob, breaks=klassen, prob=T, main="Empir. Häufigkeiten der Ergebn. d. Verwerfungsroutine") curve(dbeta(x,a,b), add=T, col="red") klassengewichte <- (pbeta(klassen[-1],a,b)-pbeta(klassen[-length(klassen)],a,b)) points(hi$mids, klassengewichte/(klassen[-1]-klassen[-length(klassen)]), col="blue") # und per chi^2-Test: chisq.test(hi$counts, p=klassengewichte) # und falls wir der chi^2-Approximation hier nicht trauten: chisq.test(hi$counts, p=klassengewichte, simulate.p.value=T, B=500) # # Beta-Verteilung (mit Parametern 0 1) aus einem # Beta(a,1)-Vorschlag gewinnen # (Beta(a,1) hat Dichte a x^(a-1) auf [0,1], also Verteilungsfunktion # F(x)=x^a und inverse Vert.fkt. F^{-1}(y)=y^(1/a) (für y in [0,1])) a <- 0.4; b <- 2 curve(x^(a-1)*(1-x)^(b-1), xlim=c(0.001,1), ylim=c(0,20)) dbeta.unnorm <- function(x) x^(a-1)*(1-x)^(b-1) quotient <- function(x) (1-x)^(b-1) # dies ist (bis auf Normierug) # die Beta(a,b)-Dichte geteilt durch # die Beta(a,1)-Dichte C <- 1 sim.beta.VM1 <- function() { repeat { x <- (runif(1))^(1/a) # dies ist Beta(a,1)-verteilt: # Methode der Inversion der Verteilungsfunktion u <- runif(1) if (u <= quotient(x)/C) # bemerke: Es genügt, die "Zieldichte" # bis auf Normierung zu kennen break } x } # Überzeugen wir uns wiederum, dass das funktioniert: wdh <- 10000 beob <- replicate(wdh, sim.beta.VM1()) # zunächst per Auge: klassen <- seq(from=0, to=1, by=0.05) hi <- hist(beob, breaks=klassen, prob=T, main="Empir. Häufigkeiten der Ergebn. d. Verwerfungsroutine") curve(dbeta(x,a,b), add=T, col="red") klassengewichte <- (pbeta(klassen[-1],a,b)-pbeta(klassen[-length(klassen)],a,b)) points(hi$mids, klassengewichte/(klassen[-1]-klassen[-length(klassen)]), col="blue") # und per chi^2-Test: chisq.test(hi$counts, p=klassengewichte) # und falls wir der chi^2-Approximation hier nicht trauten: chisq.test(hi$counts, p=klassengewichte, simulate.p.value=T, B=500) # # Verwerfungsmethode und Zerlegung des Wertebereichs: # # Sei f gewünschte "Ziel-Dichte" (bzw. -Gewichte). # Annahme: Wir können den Wertebereich in disjunkte Teile B1, B2 zerlegen # und können X1 mit Werten in B1 mit Dichte g1 bzw. # X2 mit Werten in B2 mit Dichte g2 simulieren, und es gilt # f(x)/g1(x) <= C1 für x in B1, f(x)/g2(x) <= C2 für x in B2. # Verfahren: # Wähle I=1 mit W'keit C1/(C1+C2), I=2 mit W'keit C2/(C1+C2) # Wenn I=1: Simuliere (X1,U), akzeptiere Y:=X1, wenn U <= f(X)/(C1*g1(X)), # sonst beginne das Verfahren komplett neu, # analog wenn I=2: Simuliere (X2,U), akzeptiere # Y:=X2, wenn U <= f(X)/(C2*g2(X)), etc. # # Man kann analog zu obigem Argument im "unzerteilten Fall" zeigen, # dass Y die Verteilung(sgewichte|dichte) f hat. # Ein entsprechend angepasstes Verfahren erlaubt, den Raum in mehr # als 2 Stücke zu zerlegen. # Beispiel: Gamma-Verteilung mit Formparameter ga < 1. # (Dichte ist x^(ga-1)*exp(-x)/gamma(ga).) ga <- 0.45 curve(x^(ga-1)*exp(-x), xlim=c(0,3), ylim=c(0,5), n=1000) abline(v=1, lty=2) curve(x^(ga-1), xlim=c(0.001,1), add=T, col="red") curve(exp(-x), xlim=c(1,5), add=T, col="blue") # Wir benutzen links der 1 als Vorschlagsverteilung Beta(ga, 1), # rechts der 1 eine um 1 verschobene Exp(1): ga.dichte.unnormiert <- function(x) x^(ga-1)*exp(-x) vorschl.dichte1 <- function(x) ga*x^(ga-1) C1 <- 1/ga # ga.dichte.unnormiert/vorschl.dichte1 ist <= 1/ga auf (0,1] vorschl.dichte2 <- function(x) exp(x-1) C2 <- exp(-1) # ga.dichte.unnormiert/vorschl.dichte2 ist <= 1/e auf (1,\infty) sim.gamma.VM <- function() { repeat { if (runif(1) <=C1/(C1+C2)) { x <- (runif(1))^(1/ga) # dies ist Beta(ga,1)-verteilt: u <- runif(1) if (u <= ga.dichte.unnormiert(x)/(C1*vorschl.dichte1(x))) break } else { x <- -log(runif(1))+1 # dies ist eine um 1 nach rechts verschobene Exp(1) u <- runif(1) if (u <= ga.dichte.unnormiert(x)/(C2*vorschl.dichte2(x))) break } } x } # Überzeugen wir uns wiederum, dass das funktioniert: wdh <- 10000 beob <- replicate(wdh, sim.gamma.VM()) # zunächst per Auge: klassen <- seq(from=0, to=max(beob)+0.1, by=0.05) hi <- hist(beob, breaks=klassen, prob=T, main="Empir. Häufigkeiten der Ergebn. d. Verwerfungsroutine") curve(dgamma(x,shape=ga), add=T, col="red") klassengewichte <- (pgamma(klassen[-1],shape=ga)-pgamma(klassen[-length(klassen)],shape=ga)) points(hi$mids, klassengewichte/(klassen[-1]-klassen[-length(klassen)]), col="blue") # und per chi^2-Test: chisq.test(c(hi$counts, 0), p=c(klassengewichte,1-sum(klassengewichte))) # oder per Simulation des p-Werts: chisq.test(c(hi$counts, 0), p=c(klassengewichte,1-sum(klassengewichte)), simulate.p.value=T, B=500)