# # Stochastik-Praktikum, WS 2020/2021, JGU Mainz # 23.11.2020 # # ######################################################### # heute: # # (*) Reprise: mehr zu Grafiken in R # # (*) Wann ist die Poissonapproximation der # Binomialverteilung gut? (Und wie gut ist sie dann?) # # ######################################################### ## ####################################### ## Mehr zu Rs Grafikfaehigkeiten ## ####################################### ## Ein (von R mitgeliefertes) Beispiel eines ## 2-dim Datensatzes: ## die old faithful-Daten (s.a. ?faithful) ## ["Old Faithful" ist ein bekannter Geysir im Yellowstone-Nationalpark] data(faithful) faithful ## enthaelt 272 Datenpaare der Art (Dauer einer Eruption, Wartezeit) plot(faithful$waiting, faithful$eruptions, xlab="waiting time needed for eruption [mn]", ylab="duration of eruption [mn]") title(main = c("Old faithful geyser", "Daten aus Azzalini+Bowman (1990)")) range(faithful$waiting) hist(faithful$waiting, breaks=seq(40.5, 110.5, 2), col=8 ) range(faithful$eruptions) hist( faithful$eruptions, breaks=seq(0.5, 6.5, length.out=35), col=8 ) ## ############################################################ ## 2d- und 3d-Grafiken ## Farben, Farbbeschreibungen, ## Farbskalen: heat.colors, rainbow, ... ?rgb ?heat.colors ?topo.colors ?colorRampPalette ## z.B. rgb(0.7,0.4,0.2) heat.colors(5) # Rs benannte Farben: colors() # fuer grafische Darstellungen sind parametrisierbare # Farbverlaeufe oft praktisch, z.B. AnzFarben <- 25 Farbe <- colorRampPalette(c("blue", "green"))(AnzFarben) # interpoliert von blau nach gruen plot(-1, main="Farbspiele", xlim=c(0,AnzFarben), ylim=c(0,1), xlab="", ylab="",type="n") for (i in 1:AnzFarben) { rect(i-1,0,i,1,col=Farbe[i]) } ## z.B. auch ausprobieren: Farbe <- heat.colors(AnzFarben) # oder Farbe <- colorRampPalette(c("blue", "yellow", "green"))(AnzFarben) ## 2d-Grafiken mit image, contour # Wir moechten eine auf R^2 (oder einem Auschnitt) # definierte Funktion darstellen (etwa eine W'dichte), # z.B.: f <- function(x, y) { # (dies ist die 2-dim. Standardnormaldichte) exp(-x^2/2 - y^2/2)/(2*pi) } x <- seq(from=-3, to=3, by=0.1) y <- seq(from=-2, to=2, by=0.1) # Schreibe Funktionswerte werte <- matrix(0, nrow=length(x), ncol=length(y)) for (i in 1:length(x)) { for (j in 1:length(y)) { werte[i,j] <- f(x[i],y[j]) } } ## Bemerke: dasselbe auf einmal leistet werte2 <- outer(x, y, f) ## (s.a. ?outer ) max(abs(werte-werte2)) # (ist tatsaechlich = 0) ## s.a. ?image, ?contour, etc. image(x,y,werte) contour(x,y,werte) image(x,y,werte, main="Glocke") contour(x,y,werte, add=TRUE, nlevels=20) # wir koennen dem Plot nun weitere Informationen hinzufuegen, z.B. points(0,0,pch=6,col="green") text(0.1,0.1,"Maximalwert",col="white") ## (uebrigens: Rs Punkt-Symbole sind) plot(0:25, rep(1, times=26), pch=0:25, xlab="Index",ylab="",main="Punktsymbol-Nummern in R") ## fuer "schoenere" Plots mit Legende gibt es Zusatzpakete: # aehnliches kann library(lattice) filled.contour(x,y,werte) filled.contour(x,y,werte,color.palette=terrain.colors) # s.a. die Hilfe: ?filled.contour # aber: filled.contour "verstellt" das Koordinatensystem, # wir koennen nun nicht mehr (ohne Zusatzaufwand) # weitere Informationen hinzufuegen. Beispiel: points(0,0,pch=6) # uebrigens: # Zeilenindizes werden auf der x-Achse, # Spaltenindizes auf der y-Achse aufgetragen: m <- matrix(1:12,nrow=3) m image(m) # weitere Funktionen zum Experimentieren, # auch mit anderen Gittern x <- seq(from=-3, to=3, by=0.1) # auch mit verschiedenen Gittern experimentieren ... y <- seq(from=-2, to=2, by=0.1) f <- function(x, y) { z <- x+1i*y; Re(cos(z)/(abs(z)+0.1)) } # ein unsymmetrischeres Beispiel (besser mit anderem Gitter: y-Bereich vergroessern): f <- function(x, y) { 0.35*exp(-x^2/2 - y^2/2)/(2*pi) + 0.65*exp(-((x-2.5)^2-2*sqrt(0.7)*(x-2.5)*(y-3.7)+(y-3.7)^2)/(2*(1-0.7)))/(2*pi*(1-0.7)) } werte <-outer(x,y,f) # dann Bilder wie oben generieren, z.B. ... AnzFarben <- 30 # image(x,y,werte,nlevel=AnzFarben, col=rainbow(AnzFarben)) image(x,y,werte, col=rainbow(AnzFarben)) contour(x,y,werte, add=TRUE) # ########################################################## # # Erinnerung: Drehungen in R^2 t <- 0.1*pi M <- matrix(c(cos(t),sin(t),-sin(t),cos(t)),nrow=2) ## Drehmatrix, Drehwinkel t M x <- c(1,2) M %*% x # (strengenommen muesste man M %*% matrix(x, ncol=1) schreiben, da # x hier als Spaltenvektor = 2 x 1-Matrix aufgefasst wird) M %*% matrix(x, ncol=1) ## (das Resultat ist dasselbe) # P.S.: as.vector(M %*% x) # nimmt die "Zeilendimensions-Information" wieder weg # Bsp: (Bunte Rechtecke drehen) p1 <- c(-3,-0.5); p2 <- c(3,-0.5) p3 <- c(3,0.5); p4 <- c(-3,0.5) # (Vier Eckpunkte) farben <- rainbow(10) plot(0, xlim=c(-4,4),ylim=c(-4,4),type='n') for (i in 1:9) { polygon(x=c(p1[1],p2[1],p3[1],p4[1]),y=c(p1[2],p2[2],p3[2],p4[2]),col=farben[i]) p1 <- as.vector(M %*% p1); p2 <- as.vector(M %*% p2) p3 <- as.vector(M %*% p3); p4 <- as.vector(M %*% p4) } # Uebrigens: R kann auch "Transparenz/Deckkraft" in Farben verwenden: plot(0, xlim=c(-4,4),ylim=c(-4,4),type='n') for (i in 1:9) { polygon(x=c(p1[1],p2[1],p3[1],p4[1]),y=c(p1[2],p2[2],p3[2],p4[2]),col=rgb(0.97,0.3,0.2,alpha=0.3)) p1 <- as.vector(M %*% p1); p2 <- as.vector(M %*% p2) p3 <- as.vector(M %*% p3); p4 <- as.vector(M %*% p4) } ## ####################################################################### ## ## 3d-Grafiken mit persp persp(x,y,werte) persp(x,y,werte, phi=30, theta=-30) # Sichtwinkel variieren, siehe auch ?persp persp(x,y,werte, phi=10) # "Kobreite" phi persp(x,y,werte, phi=20) persp(x,y,werte, phi=30) persp(x,y,werte, theta=10) # "Azimuthwinkel theta persp(x,y,werte, theta=20) persp(x,y,werte, theta=30) # siehe auch ?persp # Achsen-Details, z.B. persp(x,y,werte,ticktype="detailed",zlab="f(x,y)") persp(x,y,werte,ticktype="detailed",nticks=3, zlab="f(x,y)") persp(x,y,werte,ticktype="detailed",nticks=5, cex=0.2, zlab="f(x,y)") # persp-Plots einfaerben: # die "Hoehenwerte" der Facettenmittelpunkte: wm <- (werte[-1,-1]+werte[-1,-length(y)]+werte[-length(x),-1]+werte[-length(x),-length(y)])/4 anz.farben <- 22 a <- cut(wm, quantile(wm, seq(0,1, length.out=anz.farben+1), include.lowest=TRUE)) farbw <- terrain.colors(anz.farben)[a] # farbw <- heat.colors(anz.farben)[a] persp(x,y,werte,col=farbw, phi=30, theta=-30) persp(x,y,werte,col=farbw, phi=30, theta=-40) for (i in seq(0, 360, by=5)) { persp(x,y,werte,col=farbw, phi=30, theta=i) } ## Beispiel: der von R mitgeliefert "volcano"-Datensatz, siehe ?volcano # (topograhische Hoehen-Information (in m) auf 10m x 10m-Gitter fuer Maunga Whau-Gipfel) data("volcano") v <- volcano dim(v) image(v) image(v, col=topo.colors(20)) contour(v, add=TRUE) persp(v) # einfaerben: ## die "Hoehenwerte" der Facettenmittelpunkte: vm <- (v[-1,-1]+v[-1,-dim(v)[2]]+v[-dim(v)[1],-1]+v[-dim(v)[1],-dim(v)[2]])/4 anz.farben <- 25 a <- cut(vm, quantile(vm, seq(0,1, length.out=anz.farben+1), include.lowest=TRUE)) farbw <- terrain.colors(anz.farben)[a] # farbw <- heat.colors(anz.farben)[a] persp(v,col=farbw, border=NA, phi=30, theta=-30) persp(v,col=farbw, border=NA, phi=60, theta=0, ltheta=30, lphi=60, shade=0.9) # oder "esoterischer" eingefaerbt: farbfunktion <- colorRampPalette(c("darkblue","yellow","red")) farbw <- colorRampPalette(c("darkblue","red"))(anz.farben)[a] farbw <- farbfunktion(anz.farben)[a] persp(v,col=farbw, border=NA, phi=30, theta=-30) ## Zusatz: Dreh- und zoombare 3d-Plots mit Zusatzpaket rgl # (siehe auch ?persp3d) library(rgl) persp3d(v) # einfaerben: anz.farben2 <- 11 a2 <- cut(v, quantile(v, seq(0,1, length.out=anz.farben2+1), include.lowest=TRUE)) farbw3d <- matrix(terrain.colors(anz.farben2)[a2], nrow=dim(v)[1]) persp3d(v,col=farbw3d,smooth=FALSE) # (beachte, dass die Farbbeschreibung der Facetten bei persp3d # etwas anders funktioniert als bei persp, siehe ?persp3d ) ## ####################################################################### ## Poisson-Approximation der Binomialverteilungen illustrieren ## ## ####################################################################### # Für n groß und p klein so dass lambda:=n*p nicht "zu groß" ist, # werden die Gewichte der Bin(n,p)-Verteilung gut durch # die Gewichte der Poisson(lambda)-Verteilung approximiert # (vgl. auch Proposition 1.26 der Vorlesung) n <- 20; p <- 0.1 # ggf variieren, z.B. # n <- 100; p <- 0.02 # n <- 100; p <- 0.1 bgew <- dbinom(0:n, size=n, prob=p) pgew <- dpois(0:n, lambda=n*p) titel <- c(paste("Bin(",n,",",p,")-Gewichte und Pois(",n*p,")-Gewichte",sep=""), paste("Summe der (absoluten) Differenzen:", sum(abs(bgew-pgew)))) plot(0:n, bgew, col="blue", pch=1, main=titel, ylab="Gewicht") points(0:n, pgew, col="red", pch=2) legend("topright", legend=c("Binomial","Poisson"), col=c("blue","red"), pch=c(1,2)) # # Für welche Wertekombinationen von n und p ist die # Poissonappoximation gut? # Wir verschaffen uns grafisch einen Eindruck. # # Folgende Funktion bestimmt # die Summe der Absolutbeträge der Differenzen der Binomial(n,p)- und der # Poisson(n*p)-Gewichte, den sog. "Totalvariationsabstand" abstandBinPoi <- function(n,p) { #n <- n[1]; p <- p[1] sum(abs(dbinom(0:n, size=n, prob=p)-dpois(0:n, lambda=n*p))+ppois(n,lambda=n*p,lower.tail=F)) } nwerte <- 1:100 pwerte <- seq(from=0.01, to=0.1, by=0.005) # Spass mit Rs Vektoroperationen (siehe ?outer, ?Vectorize) abstaende <- outer(nwerte, pwerte, FUN=Vectorize(abstandBinPoi,c("n","p"))) # abstaende[i,j] enthält nun das Ergebnis von # abstandBinPoi(nwerte[i], pwerte[j]) # Dasselbe (aber langsamer) liefert abstaende2 <- matrix(0, nrow=length(nwerte),ncol=length(pwerte)) for (i in 1:length(nwerte)) { for (j in 1:length(pwerte)) abstaende2[i,j] <- abstandBinPoi(nwerte[i], pwerte[j]) } # die Abstände als farbkodierte 2d-Grafik: image(nwerte, pwerte, abstaende, xlab="n", ylab="p", main="Güte der Poissonapproximation als Fkt. von n und p") # füge Konturlinien hinzu: contour(nwerte, pwerte, abstaende, add=T, labcex=0.8) # alternative Darstellung: filled.contour(nwerte, pwerte, abstaende, xlab="n", ylab="p", main="Güte der Poissonapproximation als Fkt. von n und p", col=heat.colors(12)) # Die Poissonapproximation ist sehr gut, wenn n*p^2 sehr klein ist: # wiederhole obiges für x11() nwerte <- 1:100; pwerte <- seq(from=0.01, to=0.02, by=0.0005) # # Robustheit der Poisson-Approximation illustrieren # # A[i] unabhängige Ereignisse, die jeweils mit W'keit p[i] # eintreten, i=1,...,n, (die p[i] müssen nicht identisch sein), # S:=1(A[1])+...+1(A[n]), # lambda:=p[1]+...+p[n] # Dann gilt grob gesprochen: # Wenn die individuellen p[i] klein sind, so ist S approximativ # Poisson(lambda)-verteilt # ( als mathematischer Satz formuliert: # \sum_{k=0}^\infty | P(S=k) - dpois(k, lambda) | # <= 2 \sum_{i=1}^n p[i]^2 , # siehe z.B. Georgii, Stochastik, 3. Aufl., de Gruyter, 2007, # Satz 5.23; man spricht hier von "Konvergenz bezüglich dem # Totalvariationsabstand") n <- 200 # Wir wählen "chaotische" p[i], beispielsweise zufällig: p <- abs(rnorm(n, mean=2/n, sd=1/n)) # ggf n und mean, sd variieren # sieht chaotisch aus: plot(p) # weniger chaotische Darstellung plot(sort(p, decreasing=TRUE)) lambda <- sum(p) lambda # zur Kontrolle: ist \sum_{i=1}^n p[i]^2 klein? sum(p^2) # Simuliere S=1(A[1])+...+1(A[n]) simS <- function() sum(runif(n) <= p) simS() # Wir studieren die Güte der Poissonapproximation per Simulation: wdh <- 10000 Swerte <- replicate(wdh, simS()) # und betrachten das Histogramm: maxSwert <- max(Swerte) hi <- hist(Swerte, prob=T, breaks=(-0.5):(maxSwert+0.5), main="Vert. der simulierten S-Werte") # und die zugehörigen Poissongewichte: points(0:maxSwert, dpois(0:maxSwert, lambda), col="red") # die empirischen Häufigkeiten: hi$density # zum Vergleich: die Poissongewichte dpois(0:maxSwert, lambda) # der Abstand der beiden Verteilungen (in der Totalvariationsnorm) ist: sum(abs(hi$density-dpois(0:maxSwert, lambda)))+ppois(maxSwert, lambda, lower.tail=F) # ######################################################## # # Ein Beispieldatensatz, der (ganz) gut zur Poissonverteilung passt: # # "bacteria and blood counts", siehe # W. Feller, An introduction to probability theory and # its applications Vol I (3ed., Wiley, 1968), Ch. VI.7 (e), S. 163 beob <- c(3, 7, 14, 21, 20, 19, 7, 9) Bact <- data.frame(AnzKolonien=0:7, AnzQuadrate=beob) Bact attach(Bact) barplot(height=AnzQuadrate, names.arg=AnzKolonien, main="Anz. beobachteter Quadrate mit geg. Anz. Kolonien", xlab="Bakterienkolonien") # Passen diese Daten zu einer Poissonverteilung # (d.h. zur Idee, dass die Anz. beobachteter # Bakterienkolonien in einem Quadrat der Petrischale # Poisson-verteilt ist), zu welchem Parameter lambda? gesAnzQuadrate <- sum(AnzQuadrate) gesAnzQuadrate # == 100 gesAnzKolonien <- sum(AnzQuadrate*AnzKolonien) gesAnzKolonien # ein naheliegender Wert (sog. "Momentenschätzer") # (vgl. auch Bsp. 1.26 der Vorlesung ("die durch Schlag eines # Pferdes ... getoeteten" nach L. von Bortkewitsch)) lambdahut <- gesAnzKolonien/gesAnzQuadrate lambdahut erw <- c(gesAnzQuadrate*dpois(0:7, lambda=lambdahut)) Bact2 <- cbind(Bact, Erwartet=erw) Bact2 # Füge erwartete Anzahlen dem Plot hinzu: points(1.2*(0:7)+0.75, erw, col="red") # Um die Güte der Anpassung einzuschätzen bietet sich die # sog. chi-Quadrat-Statistik an: chiqua <- sum(((AnzQuadrate-erw)^2)/erw) chiqua # pchisq(chiqua, df=7, lower.tail=FALSE) ## (fuer Experten: feste Nullhypothese) # pchisq(chiqua, df=6, lower.tail=FALSE) ## (fuer Experten: 1 Parameter angepasst) ## ##############################################