# Stochastik-Praktikum, WS 2010/2011 # Sitzung vom 31.1.2011 # Themen: # Monte-Carlo-Integration, # zur Varianzreduktiom: antithetische Variable, # Importance sampling, Kontrollvariable # # MCMC, Beispiel Ising-Modell ######################################################## # # 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) sqrt(1-x^2) N <- 100000 # Anz. verwendete ZVn # 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) thetahut <- mean(gwerte) # Schätzwert für das Integral thetahut pi/2 # Wie genau ist der Schätzer? sd(gwerte) sd(gwerte)/sqrt(N) alpha <- 0.05 # auf asymptotischer Normalität beruhendes (1-alpha)-Konfidenzintervall # anhand dieser Simulationswerte ist: thetahut+c(-1,1)*qnorm(1-alpha/2)*(sd(gwerte)/sqrt(N)) # Monte-Carlo-Integration hat Konvergenzordnung 1/2: Nwerte <- round(seq(from=1000, to=N, length.out=10)) #thetahut.ausw <- sapply(Nwerte, function(k) sum(gwerte[1:k])/k) thetahut.ausw <- sapply(Nwerte, function(k) mean(sapply(runif(k, min=-1, max=1), g))) plot(Nwerte, thetahut.ausw, xlab="Anz. ZVn", ylab="Schätzwert") abline(h=pi/2) plot(Nwerte, sqrt(Nwerte)*(thetahut.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))) thetahut <- mean(gwerte); thetahut # Schätzwert für das Integral 4*pi/3 # der exakte Wert sd(gwerte) sd(gwerte)/sqrt(N) # 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 # 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) stdfehler.direkt <- sd(direkt)/sqrt(N) curve(sqrt(1-x^2),xlim=c(0,1),ylim=c(0,1.5)) curve(1.5*sqrt(1-x), add=TRUE, col="red") legend("topright", lty=1, col="red", legend="Beta(1,.1.5)-Dichte") curve(sqrt(1-x^2)/(1.5*sqrt(1-x)),xlim=c(0,1), 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) # "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 ##################################################### # # "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) 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 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 alpha <- 0.05 # auf asymptotischer Normalität beruhendes (1-alpha)-Konfidenzintervall # anhand dieser Simulationswerte ist: MC.schaetzer1 + c(-1,1)*qnorm(1-alpha/2)*MC.stdfehler1 4*pi/3 # der wahre Wert # prüfen wir anhand eines Meta-Experiments und QQ-Plots, ob die # Normalitätsannahme plausibel scheint: Nmeta <- 50 MC.werte1replikate <- replicate(Nmeta, mean(replicate(N, 8*f(runif(3,-1,1))))) qqnorm(MC.werte1replikate) 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 das 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 Integral cos(8*x)*sign(x)/(abs(x)^0.75) für x=-1 ... 1 # per Monte-Carlo 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="Monte-Carlo-Schätzung eines Integrals") 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) } # 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) } ############################################### ## ## MCMC ("Markov chain Monte Carlo") : ## Bsp. Gibbs-Sampler für das Curie-Weiss-Modell L <- 30 # ggf. variieren #zustCW <- sample(c(-1,1), L, replace=T) zustCW <- rep(1, times=L) #zustCW <- 2*rbinom(L, size=1, prob=0.65)-1 beta <- 0.7 h <- 0.0 Gibbs.sampler1.CW <- function(i) { # Koordinate i neu ziehen m <- sum(zustCW) mn.plus <- m+ifelse(zustCW[i]==1,0,2); mn.minus <- mn.plus-2 p <- 1/(1+exp(beta*((mn.minus^2-mn.plus^2)/(2*L)-h*(mn.minus-mn.plus)))) 2*rbinom(1,size=1,prob=p)-1 } Gibbs.sampler.CW <- function() { Gibbs.sampler1.CW(sample(1:L,1)) } # replicate(25, Gibbs.sampler.CW()) schr <- 500 #1e3 m.werte <- numeric(schr) for (t in 1:schr) { i<-sample(1:L,1) zustCW[i]<-Gibbs.sampler1.CW(i) m.werte[t]<-sum(zustCW) } plot(1:schr, m.werte[1:schr], type="l") # eine längere Beobachtungsreihe: zustCW <- rep(1, times=L) schr <- 1e4 for (t in 1:schr) { i<-sample(1:L,1) zustCW[i]<-Gibbs.sampler1.CW(i) m.werte[t]<-sum(zustCW) } plot(1:schr, m.werte[1:schr], type="l") burnin <- 1500 # schätzen wir die Korrelation von Werten im Abstand a akl <- 100 a.werte <- 1:akl mw <- mean(m.werte[burnin:schr]) mvar <- var(m.werte[burnin:schr]) emp.autokorr <- numeric(akl) for (i in 1:akl) { emp.autokorr[i] <- (mean(m.werte[burnin:(schr-akl)]*m.werte[(i-1)+(burnin:(schr-akl))])-mw^2)/mvar } x11() plot(0:(akl-1), emp.autokorr, xlab="Zeitdiff.", main=paste("Empir. Autokorr. (beta=",beta,"h=",h,")")) # sim.Gibbs.Kette <- function(schritte) { for (schr in 1:schritte) { i <- sample(1:L,1); zustCW[i]<<-Gibbs.sampler1.CW(i) # beachte "<--" für Zuweisung an globale Variable, s.a. ?"<<-" } } # Schätze mittlere Magnetisierung via MCMC: burnin <- 1500; lag <- 500 ss <- 500 beob <- numeric(ss) #zustCW <- sample(c(-1,1), L, replace=T) zustCW <- rep(1, times=L) #zustCW <- 2*rbinom(L, size=1, prob=0.65)-1 sim.Gibbs.Kette(burnin) beob[1] <- sum(zustCW) for (j in 2:ss) { sim.Gibbs.Kette(lag) beob[j] <- sum(zustCW) } mean(beob) hist(beob, breaks=(-L:(L+1))-0.5, prob=T) ## # für kleines L kann man in diesem Fall auch # numerisch "exakt" rechnen: L <- 30 beta <- 0.7; h <- 0.0 Gibbsgew.unnorm <- function(m) { exp(-beta*((1-m^2/L)/2-h*m)) } gg <- sapply(seq(from=-L, to=L, by=2), Gibbsgew.unnorm) Gibbsgew <- gg/sum(gg) Gibbsgew # mittlere Magnetisierung: sum(seq(from=-L, to=L, by=2)*Gibbsgew) # konstruiere ein Diagramm: h.werte <- c(0, 0.001, 0.01, 0.05, 0.1) beta.werte <- seq(from=0, to=2, length.out=25) mm <- matrix(0, ncol=length(beta.werte), nrow=length(h.werte)) for (i in 1:length(h.werte)) { h <- h.werte[i] for (j in 1:length(beta.werte)) { beta <- beta.werte[j] gg <- sapply(seq(from=-L, to=L, by=2), function(m) exp(-beta*((1-m^2/L)/2-h*m))) mm[i,j]<-sum(seq(from=-L, to=L, by=2)*gg)/sum(gg)/L } } farben <- rainbow(length(h.werte)) plot(0, xlim=c(min(beta.werte), max(beta.werte)), ylim=c(min(mm), max(mm)), xlab="beta", ylab="mittlere Magnetisierung", main=paste("Curie-Weiss-Modell, L=",L), type="n") for (i in 1:length(h.werte)) { points(beta.werte, mm[i,], col=farben[i], type="l") } ############################################### ## ## Zum (2d-)Ising-Modell # ## Größe des zugrundeliegenden Gitters: #Lx <- 32 #Ly <- 32 Lx <- 80 Ly <- 80 # Lx <- Ly <- 256 # Konstant +1 als Startzustand # zust <- matrix(1, nrow=Lx, ncol=Ly) # oder auch zust <- matrix(sample(c(-1,1), size=Lx*Ly, replace=T), nrow=Lx, ncol=Ly) J <- 1.0 beta <- 1.0 # 0.03 #0.1 # 1.0 h <- 0.0 ## ## Zur Rand-Periodisierung: ## (i %% Lx)+1 verwandelt i aus {1,...,Lx} in seinen "rechten" Nachbarn, ## ((i-2) %% Lx)+1 verwandelt i in seinen "linken" Nachbarn, ## wobei 1 und Lx als benachbart gelten. Gibbs.local.proposal.pR <- function(i,j) { ## Gibbs-Aktualisierung bei zust[i,j] vorschlagen, ## periodischer Rand summe.nachbarn <- zust[(i %% Lx)+1,j]+zust[((i-2) %% Lx)+1,j]+zust[i,(j %% Ly)+1]+zust[i,((j-2)%%Ly)+1] 2*rbinom(1,size=1,prob=1/(1+exp(-2*beta*(J*summe.nachbarn+h))))-1 } ## image(zust) male.zust <- function(...) { ## Aktuellen Zustand zeichnen: +1 = weiß, -1 = schwarz. ## Beachte: weitere Grafikparameter werden an die Funktion ## image "durchgereicht" image(1:Lx, 1:Ly, zust, col=c("black","white"), breaks=c(-2,0,2), xlab="x", ylab="y", ...) } male.zust() anz.schr <- 100000 zeichne.nach.jeweils <- Lx*Ly # 10*Lx*Ly male.zust(main="Konfiguration am Anfang") for (schr in 1:anz.schr) { i <- sample(1:Lx,1); j <- sample(1:Ly,1) zz<-Gibbs.local.proposal.pR(i,j) zust[i,j]<-zz if (schr %% zeichne.nach.jeweils == 0) { male.zust(main=paste("Konfiguration nach", schr, "Schritten", "(beta=",beta," h=",h,")")) } } ## die Kette für eine Anzahl Schritte simulieren: sim.Gibbs.Kette <- function(schritte) { for (schr in 1:schritte) { i <- sample(1:Lx,1); j <- sample(1:Ly,1) zust[i,j]<<-Gibbs.local.proposal.pR(i,j) # beachte "<--", s.a. ?"<<-" } } ## Bem: Mit ## system.time(sim.Gibbs.Kette(1e5)) ## kann man seine "Geduld kalibrieren". wdh <- 20 mittlere.magn <- numeric(wdh) for (w in 1:wdh) { print(w) zust <- matrix(sample(c(-1,1), size=Lx*Ly, replace=T), nrow=Lx, ncol=Ly) sim.Gibbs.Kette(1e5) mittlere.magn[w]<-mean(zust) } mittlere.magn plot(mittlere.magn) mean(mittlere.magn) ## Bem.: Alles im Speicher löschen: rm(list = ls())