# # Stochastik-Praktikums, WS 2020/2021, JGU Mainz # 8.2.2021 ## ############################################################# ## ## Thema: ## Markovketten, speziell: Konvergenz und Gleichgewichte ## ##eine Anwendung (exemplarisch): Markov chain Monte Carlo ## am Beispiel des Gibbs-Samplers fuer das ## Curie-Weiss- und das Ising-Modell ## ## ############################################################# # (Erinnerung: Simulation eines Uebergangsschritts gemaess # Uebergangsmatrix A) sim.MKschritt <- function(A, x) sample(1:ncol(A), 1, prob=A[x,]) # (und ein "naives" Stueck Code fuer Potenzen von A:) matrixpotenz <- function(A, n) { B <- diag(x=1,nrow=dim(A)[1]) while(n>0) { B <- B %*% A; n <- n-1 } B } # #################################### # Ein simples "Wetter-Beispiel" # (1="Sonne", 2="Bedeckt", 3="Regen") A <- matrix(c(0.6,0.3,0.1, 0.3,0.2,0.5, 0,0.1,0.9), ncol=3, byrow=TRUE) A ############################################### # # numerische "Beobachtung" der Konvergenz in Gleichgewicht # fuer die "Wetter"-Kette mit Uebergangsmatrix A # Verteilung der Kette nach 30, 100 bzw. 150 Schritten # bei Start im Zustand 1 c(1,0,0) %*% matrixpotenz(A, 30) c(1,0,0) %*% matrixpotenz(A, 100) c(1,0,0) %*% matrixpotenz(A, 150) # und bei Start im Zustand 2 c(0,1,0) %*% matrixpotenz(A, 30) c(0,1,0) %*% matrixpotenz(A, 100) c(0,1,0) %*% matrixpotenz(A, 150) # (Wir sehen: auf lange Sicht haengt hier # die Verteilung nicht vom Startzustand ab.) # # Das ist anders, wenn Periodeneffekte auftreten, # betrachte dazu das Beispiel der Irrfahrt auf dem Kreis: ############################## # Uebergangsmatrix der symmetrischen Irrfahrt # auf dem Kreis aus L Knoten uebgsmat.Irrfahrt.auf.Kreis <- function(L) { A <- matrix(0, nrow=L, ncol=L) for (i in 2:(L-1)) { A[i,i-1] <- 1/2; A[i,i+1]<-1/2 } A[1,2] <- 1/2; A[1,L] <- 1/2 A[L,L-1] <- 1/2; A[L,1] <- 1/2 A } B <- uebgsmat.Irrfahrt.auf.Kreis(6) B # Man sieht: Nur Eintraege mit Zeilennummer+Spaltennummer # gerade sind > 0 matrixpotenz(B, 50) matrixpotenz(B, 51) # Dies ist allgemein so für die Irrfahrt auf einem Kreis mit # gerader Knotenanzahl, nicht aber bei ungerader Knotenanzahl, # wie folgendes Beispiel zeigt: C <- uebgsmat.Irrfahrt.auf.Kreis(7) C matrixpotenz(C, 50) matrixpotenz(C, 100) # Die Eintraege sind alle sehr nahe an 1/7: 1/7 ## ##################################################### # Uebrigens: Gleichgewichte zu A sind Links-Eigenvektoren # der Uebergangsmatrix A zum Eigenwert 1, man koennte dazu # auch (numerische) lineare Algebra verwenden: eigen(t(A)) # beachte t(A), d.h. die Transponierte verwenden, # denn R berechnet Rechts-Eigenvektoren # (zur Kontrolle:) ez <- eigen(t(A)) ez$vectors[,1] (ez$vectors[,1]) %*% A ######################################################## ######################################################## # # Das "Pagerank"-Beispiel (vgl. auch die Folien), # hier der dort verwendete R-Code # Adjazenzmatrix L <- matrix(c(0,1,1,0,1,0,0, 1,0,0,1,0,0,0, 0,1,0,0,1,0,0, 0,0,1,0,1,0,0, 0,1,1,0,0,0,0, 0,0,0,0,0,0,1, 0,0,0,0,0,1,0), ncol=7) N <- 7 alpha <- 0.9 # Bau der Uebergangsmatrix "von Hand": A <- matrix(0, nrow=N, ncol=N) for (i in 1:N) { c <- sum(L[i,]) for (j in 1:N) { A[i,j] <- alpha*ifelse(c>0, L[i,j]/c, 1/N)+(1-alpha)/N } } # (Alternativ als Zweizeiler:) C <- matrix(rep(rowSums(L),each=N), nrow=N, byrow=TRUE) Aalt <- alpha*ifelse(C>0, L/C, 1/N)+(1-alpha)*matrix(1/N,nrow=N,ncol=N) # Einen Pfad der Markovkette mit Uebergangsmatrix A fuer eine gewisse Anzahl Schrite simulieren schritte <- 100000 pfad <- integer(schritte) # Startpunkt pfad[1] <- sample(1:N) for (i in 2:schritte) pfad[i] <- sim.MKschritt(A, pfad[i-1]) # relative Zeitanteile in verschiedenen Zustaenden verfolgen: sapply(1:N, function(i) sum(pfad==i))/schritte anteil2 <- cumsum(pfad==2)/(1:schritte) anteil4 <- cumsum(pfad==4)/(1:schritte) # und zeichnen: zeitgitter <- seq(from=1, to=schritte, by=100) plot(zeitgitter, anteil2[zeitgitter], type="l", xlab="Schritte", col="black", ylab="relativer Anteil") points(zeitgitter, anteil4[zeitgitter], type="l", col="darkred") legend("bottomright",legend=c("in Zustand 2", "in Zustand 4"), lty=1, col=c("black","darkred")) # # Wie schnell konvergiert die Kette ins Gleichgewicht? Dazu: # 2^n-te Potenz einer Matrix (schnell) berechnen: matrix2erpotenz <- function(M, n) { B <- M for (i in 1:n) B <- B %*% B B } # Verteilung nach 2^15 Schritten ist numerisch vom # Gleichgewicht ununterscheidbar: pi <- rep(1/N, times=N) %*% matrix2erpotenz(A,15) pi pi %*% A pi - pi %*% A # Abstand vom Gleichgewicht bestimmen und zeichnen: schritte <- 30 stvert <- rep(1/N, times=N) abst <- numeric(schritte) B <- A for (i in 1:schritte) { abst[i] <- sum(abs(stvert %*% B - pi)) B <- B %*% A } plot(1:schritte, log10(abst), xlab="Schritte") ######################################################## ######################################################## # # MCMC ("Markov chain Monte Carlo") : # Bsp. Gibbs-Sampler fuer das Curie-Weiss-Modell # Pierre Curie (1859-1906), frz. Physiker # Pierre-Ernest Weiss (1865-1940), frz. Physiker # Josiah Willard Gibbs (1839-1903), US-amerik. Physiker, Mathematiker, Chemiker # Anz. "Elementarmagnete" im Modell: N <- 30 # ggf. variieren #zustCW <- sample(c(-1,1), N, replace=T) zustCW <- rep(1, times=N) #zustCW <- 2*rbinom(N, size=1, prob=0.65)-1 beta <- 0.7 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*N)))) 2*rbinom(1,size=1,prob=p)-1 } Gibbs.sampler.CW <- function() { Gibbs.sampler1.CW(sample(1:N,1)) } # replicate(25, Gibbs.sampler.CW()) schr <- 500 #1e3 m.werte <- numeric(schr) for (t in 1:schr) { i<-sample(1:N,1) zustCW[i]<-Gibbs.sampler1.CW(i) m.werte[t]<-sum(zustCW) } plot(1:schr, m.werte[1:schr], type="l") # eine laengere Beobachtungsreihe: zustCW <- rep(1, times=N) schr <- 1e4 for (t in 1:schr) { i<-sample(1:N,1) zustCW[i]<-Gibbs.sampler1.CW(i) m.werte[t]<-sum(zustCW) } plot(1:schr, m.werte[1:schr], type="l") burnin <- 1500 # schaetzen 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:N,1); zustCW[i]<<-Gibbs.sampler1.CW(i) # beachte "<<-" für Zuweisung an globale Variable, s.a. ?"<<-" } } # Schaetze mittlere Magnetisierung via MCMC: burnin <- 1500; lag <- 500 ss <- 500 beob <- numeric(ss) #zustCW <- sample(c(-1,1), N, replace=T) zustCW <- rep(1, times=N) #zustCW <- 2*rbinom(N, 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=(-N:(N+1))-0.5, prob=T) mean(abs(beob)) ## # für kleines N kann man in diesem Fall auch # numerisch "exakt" rechnen: N <- 70 # (sollte hier eine gerade Zahl sein, damit es woertlich klappt ...) beta <- 0.7 ## ggfs. variieren, z.B. : beta <- 1.3 Gibbsgew.unnorm <- function(m) { exp(-beta*((1-m^2/N)/2))*choose(N,(N-m)/2) } gg <- sapply(seq(from=-N, to=N, by=2), Gibbsgew.unnorm) Gibbsgew <- gg/sum(gg) Gibbsgew # mittlere (absolute) Magnetisierung: sum(abs(seq(from=-N, to=N, by=2))*Gibbsgew)/N # konstruiere ein Diagramm: beta.werte <- seq(from=0, to=2, length.out=25) mam <- numeric(length(beta.werte)) for (j in 1:length(beta.werte)) { beta <- beta.werte[j] gg <- sapply(seq(from=-N, to=N, by=2), function(m) exp(-beta*((1-m^2/N)/2))) mam[j] <- sum(abs(seq(from=-N, to=N, by=2))*gg)/sum(gg)/N } plot(beta.werte, mam, xlim=c(min(beta.werte), max(beta.werte)), xlab="beta", ylab="mittlere absolute Magnetisierung", main=paste("Curie-Weiss-Modell, N=",N), type="l") ############################################### ## ## Zum (2d-)Ising-Modell ## ## (nach Ernst Ising, 1900-1998 benannt) ## Groesse des zugrundeliegenden Gitters: #Lx <- 32 #Ly <- 32 Lx <- 80 Ly <- 80 # Lx <- Ly <- 256 Lx <- Ly <- 128 # 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 # "Kopplungskonstante" 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="", ylab="", ...) } male.zust() anz.schr <- 100000 zeichne.nach.jeweils <- round(Lx*Ly/10) # 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,")",sep="")) } } ## die Kette fuer 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())