# # Sitzung des Stochastik-Praktikums, WS 2013/2014, JGU Mainz # 6.1.2014 ## Themen: ## Monte-Carlo-Integration, ## Methoden zur Varianzreduktion: ## Importance sampling, antithetische Variable, Kontrollvariable ## Simulation vorgegebener Verteilungen mittels ## Verwerfungsmethode ("Rejection sampling") ######################################################## # # Monte-Carlo-Integration # Beispiel: Integriere sqrt(1-x^2) von -1 bis 1 # (wir wissen: dies ergibt pi/2, die halbe Fläche des Einheitskreises) # (ausserdem: (d/dx 0.5*x*sqrt(1-x^2)+0.5*arcsin(x)=sqrt(1-x^2) ) f <- function(x) ifelse(abs(x)<1, sqrt(1-x^2), 0) curve(f, xlim=c(-1.2,1.2)) ## Plot von f N <- 100000 # Anz. verwendete ZVn, ggfs. erhöhen # schreibe f=2*f*0.5 und interpretiere als Integral bzgl. unif([-1,1]) g <- function(x) 2*f(x) gwerte <- sapply(runif(N, min=-1, max=1), g) I.hut <- mean(gwerte) # Schätzwert für das Integral I.hut pi/2 # Wie genau ist der Schätzer? sd(gwerte) sd(gwerte)/sqrt(N) # der Standardfehler, d.h. die geschaetzte Varianz # von I.hut # Monte-Carlo-Integration hat Konvergenzordnung 1/2, # schauen wir es uns an: Nwerte <- round(seq(from=1000, to=N, length.out=10)) Nwerte I.hut.ausw <- sapply(Nwerte, function(k) mean(sapply(runif(k, min=-1, max=1), g))) I.hut.ausw plot(Nwerte, I.hut.ausw, xlab="Anz. Replikate", ylab="Schätzwert") abline(h=pi/2, lty=2) plot(Nwerte, sqrt(Nwerte)*(I.hut.ausw-pi/2), xlab="Anz. ZVn", ylab="(Schätzwert-wahrer Wert)*sqrt(Auswertungen)") abline(h=0) # # Beispiel: Volumen der Einheitskugel im R^3 # N <- 100000 f <- function(x,y,z) ifelse(x*x+y*y+z*z<=1,1,0) # schreibe f=8*f*(0.5)^3 und interpretiere als Integral bzgl. unif([-1,1]^3) gwerte <- replicate(N, 8*f(runif(1,min=-1,max=1),runif(1,min=-1,max=1),runif(1,min=-1,max=1))) I.hut <- mean(gwerte); I.hut # Schätzwert für das Integral 4*pi/3 # der exakte Wert sd(gwerte) sd(gwerte)/sqrt(N) # geschätzte Varianz von I.hut # Zum Vergleich: Numerische Integration auf einem festen Rechtecksgitter: N1 <- round(N^(1/3)) N1 gitter <- seq(from=-1, to=1, length.out=N1) masche <- gitter[2]-gitter[1] gittervol <- masche^3 summe <- 0.0 for (x in gitter) { for (y in gitter) { for (z in gitter) { summe <- summe+f(x,y,z)*gittervol } } } summe; 4*pi/3 ################################################ ## ## "Tricks" zur Varianzreduktion bei Monte-Carlo-Integration ## 1) Importance sampling # Bsp.: Integal sqrt(1-x^2), x=0..1 N <- 10000 # ggf. variieren # zunächst einfache Monte-Carlo-Integration direkt <- sapply(runif(N), function(x) sqrt(1-x^2)) mean(direkt) sd(direkt) sd(direkt)/sqrt(N) # Standardfehler des "direkten Schätzers" stdfehler.direkt <- sd(direkt)/sqrt(N) # Vergleichen wir die "Zielfunktion" f curve(sqrt(1-x^2),xlim=c(0,1),ylim=c(0,1.5)) # ... und eine Wahrscheinlichkeitsdichte, gemäß der wir leicht # simulieren können (die Dichte von Beta(1.0,1.5)) curve(1.5*sqrt(1-x), add=TRUE, col="red") legend("topright", lty=1, col="red", legend="Beta(1.0,1.5)-Dichte") curve(sqrt(1-x^2)/(1.5*sqrt(1-x)),xlim=c(0,1), lty=2, add=TRUE) # Benutze Beta(1,1.5) als Vorschlagsverteilung h <- function(x) sqrt(1-x^2)/(1.5*sqrt(1-x)) gewichtet <- sapply(rbeta(N,shape1=1,shape2=1.5), h) mean(gewichtet); pi/4 sd(gewichtet) sd(gewichtet)/sqrt(N) stdfehler.gewichtet <- sd(gewichtet)/sqrt(N) stdfehler.gewichtet sd(direkt)/sd(gewichtet) ## "Tricks" zur Varianzreduktion bei Monte-Carlo-Integration ## 2) "Antithetische" Variable: benutze U und 1-U zugleich # zu Vergleichszwecken: gleiche Anzahl Funktionsauswertungen N2 <- round(N/2) f <- function(x) sqrt(1-x^2) antithet <- sapply(runif(N2), function(u) (f(u)+f(1-u))/2) mean(antithet) sd(antithet) sd(antithet)/sqrt(N2) stdfehler.antithet <- sd(antithet)/sqrt(N2) # Schliesslich: Vergleich der Genauigkeit der 3 Methoden: stdfehler.direkt; stdfehler.gewichtet; stdfehler.antithet ##################################################### ## "Tricks" zur Varianzreduktion bei Monte-Carlo-Integration ## 3) "Kontrollvariable" # nochmal: Volumen der Einheitskugel im R^3 per # Monte-Carlo-Integration # (als geeignetes Integral bezgl. der uniformen Vert. auf [-1,1]^3 auffassen) f <- function(v) ifelse(v[1]*v[1]+v[2]*v[2]+v[3]*v[3]<=1,1,0) MC.werte1 <- replicate(N, 8*f(runif(3,-1,1))) MC.schaetzer1 <- mean(MC.werte1) MC.stdfehler1 <- sd(MC.werte1)/sqrt(N) MC.schaetzer1; MC.stdfehler1 4*pi/3 # der wahre Wert ## Exkurs: ## prüfen wir anhand eines Meta-Experiments und QQ-Plots, ob eine ## Normalitätsannahme plausibel scheint: Nmeta <- 50 MC.werte1replikate <- replicate(Nmeta, mean(replicate(N, 8*f(runif(3,-1,1))))) qqnorm(MC.werte1replikate) h <- function(v) 1-(v[1]*v[1]+v[2]*v[2]+v[3]*v[3]) # h hat Erwartungswert 0 unter der uniformen Vert. auf [-1,1]^3 mean(replicate(1000, h(runif(3,-1,1)))) # ... was wir empirisch "verifizieren" N <- 1000 # Anz. Simulationen, ggf. variieren, z.B. 10000 oder 1e5 konst <- 6.0 # variieren, z.B. 1.0, 3.0, 6.0, 10 fmod <- function(v) 8*f(v)-konst*h(v) # Schätzen wir das Volumen der Einheitskugel im R^3 mit Hilfe von fmod: MC.werte2 <- replicate(N, fmod(runif(3,-1,1))) MC.schaetzer2 <- mean(MC.werte2); MC.stdfehler2 <- sd(MC.werte2)/sqrt(N) MC.schaetzer2; MC.stdfehler2 # Vergleichen wir: cat(paste("Schätzer1:", MC.schaetzer1, " Standardfehler", MC.stdfehler1,"\n", "mit Wahl konst=", konst, " : Schätzer2:", MC.schaetzer2, " Standardfehler", MC.stdfehler2,"\n", "Std.fehler1/Std.fehler2=", MC.stdfehler1/MC.stdfehler2,"\n")) # Schätzen wir die optimale konst: hwerte <- numeric(N); f8werte <- numeric(N) for (i in 1:N) { v<-runif(3,-1,1) hwerte[i]<-h(v); f8werte[i]<-8*f(v) } var(hwerte) cov(hwerte, f8werte) # cov(hwerte, f8werte)/var(hwerte) ######################################################### # # Bem./Hinweis zur Endlichkeit der Varianz eines # Monte-Carlo-Schätzers # Ang., wir sollten das Integral über # cos(8*x)*sign(x)/(abs(x)^0.75) für x von -1 bis 1 # per Monte-Carlo-Methode berechnen: # (Wir wissen natürlich bereits aus Symmetrie, dass der wahre Wert = 0 ist.) f <- function(x) cos(8*x)*sign(x)/(abs(x)^0.75) # Graph anschauen: xwerte <- c(seq(-1,0.01,length.out=100),seq(0.01,1,length.out=100)) plot(xwerte, sapply(xwerte, f), ylim=c(-6,6), pch=20, xlab="x", ylab="f") N <- 1e4 werte <- replicate(N, 2*f(runif(1,-1,1))) Nwerte <- seq(from=100, to=N, length.out=12) schaetzer <- numeric(length(Nwerte)) schaetzer.str <- numeric(length(Nwerte)) for (i in 1:length(Nwerte)) { werte <- replicate(Nwerte[i], 2*f(runif(1,-1,1))) schaetzer[i]<-mean(werte) schaetzer.str[i]<-sd(werte) } plot(Nwerte, schaetzer, xlab="n", ylab="Schätzer, basierend auf n ZVn", main=c("Monte-Carlo-Schätzung eines Integrals","Schätzwert +/+ Standardfehler")) abline(h=0, lty=2) for (i in 1:length(Nwerte)) { lines(c(Nwerte[i], Nwerte[i]), c(schaetzer[i]-schaetzer.str[i]/sqrt(Nwerte[i]), schaetzer[i]+schaetzer.str[i]/sqrt(Nwerte[i])), col="blue", lwd=2) } # Ggf. wiederholen. # Bem.: Tatsächlich hat unser Schätzer hier unendliche Varianz. # # Wir können die Singularität durch eine geeignete # Vorschlagsverteilung "heilen", z.B. # Beta(0.25, 1) mit "zufälligem Vorzeichen" # (mit Dichte 0.5*dbeta(abs(x),shape1=0.25, shape2=1) auf [-1,1]) N <- 1e4 w <- function(x) f(x)/(0.5*dbeta(abs(x),shape1=0.25, shape2=1)) schaetzer <- numeric(length(Nwerte)) schaetzer.str <- numeric(length(Nwerte)) for (i in 1:length(Nwerte)) { werte <- replicate(Nwerte[i], w(sample(c(-1,1),1)*rbeta(1,shape1=0.25,shape2=1))) schaetzer[i]<-mean(werte) schaetzer.str[i]<-sd(werte) } x11() plot(Nwerte, schaetzer, xlab="n", ylab="Schätzer, basierend auf n ZVn", main="Monte-Carlo-Schätzung eines Integrals (mit Importance sampling)") abline(h=0, lty=2) for (i in 1:length(Nwerte)) { lines(c(Nwerte[i], Nwerte[i]), c(schaetzer[i]-1.96*schaetzer.str[i]/sqrt(Nwerte[i]), schaetzer[i]+1.96*schaetzer.str[i]/sqrt(Nwerte[i])), col="blue", lwd=2) } ##################################################### ## ## Simulation vorgegebener Verteilungen mittels ## 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") ## Bem./Vorgriff: Ein statistisches Standard-Verfahren, um zu ## prüfen, ob die simulierten Werte tatsächlich aus der ## behaupteten hypergeometrischen Verteilung stammen, ist ## der 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")