# Stochastik-Praktikum, WS 2019/20, JGU Mainz # 2.12.2019 ## ############################################################################## ## heute: ## Markovketten, speziell ## Beispiel Pólya-Urne ## Beispiel Galton-Watson-Prozess ######################################################## # Simulation von Markovketten (mit endlichem Zustandsraum) # ################################################ # (Erinnerung:) Matrizen in R: # M <- matrix(1:12, nrow=3, ncol=4) M # wenn man lieber "zeilenweise" auffuellen moechte: M2 <- matrix(1:12, nrow=3, ncol=4, byrow=T) M2 # t(...) transponiert eine Matrix: t(M2) # %*% ist das Matrixprodukt: A <- matrix(1:4, nrow=2) B <- matrix(c(1,0,1,1),nrow=2) A; B A %*% B # %*% berechnet auch Matrix mal Vektor v <- c(1,-1) A %*% v v %*% A # Wenn A eine (m x n)-Matrix (d.h. m Zeilen, n Spalten) # und v ein n-Vektor ist, so ist A %*% v eine (n x 1)-Matrix # (d.h. ein (n-)Spaltenvektor), analog wenn # u ein m-Vektor ist, so ist u %*% A eine (1 x m)-Matrix # (d.h. ein (m-)Zeilenvektor). M v <- c(1,3,-2,0) M %*% v # ist also "Kurzform" von M %*% matrix(v, ncol=1) u <- c(0,2,-1) u %*% M # ist also "Kurzform" von matrix(u, nrow=1) %*% M # as.vector verwandelt diese Dinge bei Bedarf wieder in # einen Vektor zurueck: as.vector(u %*% M) ################################################## # # Simulieren eines Uebergangsschritts # gemaess Uebergangsmatrix A von Startpunkt x aus: # (es bleibt allerdings dem Aufrufer ueberlassen, sicherzustellen, # dass A tatsaechlich eine stochastische Matrix und x eine pos. # ganze Zahl <= Anz. Zeilen von A ist) sim.MKschritt.v0 <- function(A, x) { u <- runif(1) y <- 0 s<-0 while(s0) { B <- B %*% A n <- n-1 } B } ### # Nebenbei: Eine ziemlich ineffiziente (aber sehr simple) rekursive # Implementierung von positiv-ganzzahligen Matrixpotenzen: "%^%" <- function(A, n) if(n == 1) A else A %*% (A %^% (n-1)) ############################## # # numerische "Beobachtung" der Konvergenz in Gleichgewicht # fuer die "Wetter"-Kette mit Übergangsmatrix 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 # Auftreffwahrscheinlichkeiten betrachten: # Wie wahrscheinlich ist es, bei Start in 5 die 14 vor # der Null zu treffen? A <- uebgsmat.Irrfahrt.auf.Kreis(15) # mache "Randzustaende" absorbierend: A[1,] <- 0; A[1,1] <- 1.0 A[15,] <- 0; A[15,15] <- 1.0 # Wie sieht es nach vielen Schritten aus? B <- matrixpotenz(A, 500) ## ggfs groessere Werte einsetzen, z.B.: B <- matrixpotenz(A, 1000) B[6,] # Verteilung bei Start in 5 plot(B[6,]) ## ############################################## ## Zur Pólya-Urne ## (F. Eggenberger, G. Pólya (Polya), ## Über die Statistik verketteter Vorgänge. ## Zs. f. angew. Math. u. Mech. 3, 279-290 (1923). ## Urne mit Anfangs w weißen und s schwarzen Kugeln, ## ziehe eine Kugel ("rein zufällig") heraus, ## lege zusammen mit Delta vielen Kugeln derselben Farbe zurück # Beispiel Polya-Urne (mit Delta=1) # polya.simulation <- function(w,s,zuege) { # simuliere Polya-Urne, startend mit w weißen und s schwarzen Kugeln # für zuege viele Züge resultat <- rep(0,times=zuege) for (i in 1:zuege) { if (runif(1)>1, # Startanteil schwarz: s/(s+w) = h/n, Delta/(s+w) = d/n d <- 14.2; h <- 5.5 ## ############################################### ## ## Ein kritischer (Bienayme-)Galton-Watson-Prozess ## (Francis Galton, 1822-1911; Henry William Watson, 1827-1903; ## Irénée-Jules Bienaymé, 1796-1878) ## Aus Jane Austen (1775-1817), "Persuasion": ## "Sir Walter Elliot, of Kellynch Hall, in Somersetshire, was a man who, ## for his own amusement, never took up any book but the Baronetage; ... ## there his faculties were roused into admiration and ## respect, by contemplating the limited remnant of the earliest patents ..." # Skalenparameter: N <- 20 # variieren, z.B. N <- 100; N <- 1000 startwert <- N #1 #N pfadsimGW <- function(n) { pf <- numeric(n+1) pf[1]<-startwert for (i in 1:n) { neue.gen <- 0 if (pf[i]>0) { for (j in 1:pf[i]) { # jedes Ind. hat 0 oder 2 Kinder, je mit W'keit 0.5 neue.gen <- neue.gen+2*rbinom(1,size=1,prob=0.5) } } pf[i+1]<-neue.gen } pf } pfad<-pfadsimGW(N) plot(0:N, pfad, xlab="Generation", ylab="Anzahl", type="l", main=paste("Galton-Watson-Prozess: Starte mit",N,"Ind."), ylim=c(0,2*N)) # ggf. mehrere Pfade einzeichnen: pfad2 <- pfadsimGW(N) points(0:N, pfad2, type="l") # Betrachte nun reskalierte Zeit- und Anzahl-Achsen: # eine Generation entspricht Zeit 1/N, # ein Individuum entspricht "Masse" 1/N x11() plot((0:N)/N, pfad/N, xlab="Zeit [ N Generationen]", ylab="Anzahl", type="l", main=paste("Galton-Watson-Prozess: Starte mit",N,"Ind."), ylim=c(0,2)) points((0:N)/N, pfad2/N, type="l") ## ggfs. mehrere Pfade simulieren und dann einzeichnen: pfad2 <- pfadsimGW(N); points((0:N)/N, pfad2/N, type="l")