# # Sitzung des Stochastik-Praktikums, WS 2013/2014, JGU Mainz # 25.11.2013 ## ################################################ ## ## Nochmal zur 2-dim Normalverteilung (und zu 2d/3d-Grafiken), ## vgl. auch Uebungsblatt 3, Aufg. 2 und Sitzung_4.11.2013.R ## z.B.: schrw <- 0.1 x <- seq(from=-3, to=3, by=schrw) y <- seq(from=-2, to=2, by=schrw) f <- function(x, y) { ## die 2-dimensionale Standardnormaldichte exp(-(x^2+y^2)/2)/(2*pi) } 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: werte2 <- outer(x, y, f) ## (s.a. ?outer ) max(abs(werte-werte2)) ## s.a. ?image, ?contour, etc. image(x,y,werte) contour(x,y,werte) image(x,y,werte, main="2-dim. Normaldichte"); contour(x,y,werte, add=TRUE, nlevels=20) persp(x,y,werte) persp(x,y,werte, phi=30, theta=-30) ## Um die 3d-Grafik einzufaerben: ## die "Hoehenwerte" der Facettenmittelpunkte: wm <- (werte2[-1,-1]+werte2[-1,-length(y)]+werte2[-length(x),-1]+werte2[-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,werte2,col=farbw, phi=30, theta=-30) persp(x,y,werte2,col=farbw, phi=30, theta=-40) for (i in seq(0, 360, by=5)) { persp(x,y,werte2,col=farbw, phi=30, theta=i) } ## ## vergleichen wir mit simulierten Daten: anz.beob <- 1000 ## eigentlich zu klein: ersetze durch anz.beob <- 1e6 z1 <- rnorm(anz.beob) z2 <- rnorm(anz.beob) ## zaehle, wie oft Beobachtungen in [x[i] +/- schrw/2] x [y[j] +/- schrw/2] liegen wbeob <- matrix(0, nrow=length(x), ncol=length(y)) for (i in 1:length(x)) { for (j in 1:length(y)) { wbeob[i,j] <- sum( (z1>=x[i]- schrw/2) & (z1=y[j] - schrw/2) & (z2 < y[j] + schrw/2)) ## beachte: "&", nicht "&&" } } sum(wbeob) wbeob <- wbeob/sum(wbeob) ## betrachten wir die empirische Verteilung graphisch: image(x,y,wbeob) contour(x,y,wbeob) ## schoener: image(x,y,wbeob, main=c("empirische Verteilung von",paste(anz.beob,"2-dim normalvert. ZVn"))) contour(x,y,wbeob, add=TRUE, nlevels=20) ## und als 3d-Grafik ## die "Hoehenwerte" der Facettenmittelpunkte (zum Einfaerben): wm <- (wbeob[-1,-1]+wbeob[-1,-length(y)]+wbeob[-length(x),-1]+wbeob[-length(x),-length(y)])/4 anz.farben <- 40 br <- c(-0.01, unique(quantile(wm, seq(0,1, length.out=anz.farben+1), include.lowest=TRUE))) ## um wiederholten Wert 0 auszuschliessen a <- cut(wm, br) #farbw <- terrain.colors(anz.farben)[a] farbw <- heat.colors(anz.farben)[a] persp(x,y,wbeob,col=farbw, phi=30, theta=-30, zlab="Empir. Dichte") persp(x,y,wbeob,col=farbw, phi=30, theta=-40, zlab="Empir. Dichte") for (i in seq(0, 360, by=5)) { persp(x,y,wbeob,col=farbw, phi=30, theta=i) } ## save(wbeob, file='wbeob1e6.RData') ## Daten speichern und ## load('wbeob1e6.RData') ## wieder laden ## ################################################## ## ## Anmerkung: ## aehnliche Befehle sind (siehe auch jeweils die Online-Hilfe) filled.contour(x,y,werte) # (zeichnet auch eine Legende) # weitere Alternativen (aus dem R-Paket lattice): library(lattice) contourplot(werte) levelplot(werte) contourplot(werte, row.values=x, column.values=y, xlab='x', ylab='y') levelplot(werte, row.values=x, column.values=y, xlab='x', ylab='y') wireframe(werte) wireframe(werte, row.values=x, column.values=y, xlab='x', ylab='y', main="Titel") wireframe(werte, row.values=x, column.values=y, xlab='x', ylab='y', drape=TRUE, main="Titel") cloud(werte) cloud(werte) ## ############################################################################## ## Poissonverteilung, Poissonapproximation der Binomialverteilung, Poissonprozess # # 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 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) # # Beispieldatensätze, die gut zur Poissonverteilung passen # # 1. Radioaktive Zerfälle, siehe # W. Feller, An introduction to probability theory and # its applications Vol I (3ed., Wiley, 1968), S. 160 rzerfall <- c(57, 203, 383, 525, 532, 408, 273, 139, 45, 27, 16) Zerf <- data.frame(Zerfaelle=0:10, AnzBeobInt=rzerfall) # Bem.: die 16 in Zeile 10 steht eigentlich für # Intervalle mit 10 oder mehr Zerfallsereignissen attach(Zerf) barplot(height=Zerf$AnzBeobInt, names.arg=Zerf$Zerfaelle, main="Anz. beobachteter Intervalle mit geg. Anz. Zerfälle", xlab="Zerfälle") # Passen diese Daten zu einer Poissonverteilung # (d.h. zur Hypothese, dass die Anz. beobachteter Zerfälle in # einem Beobachtungsintervall Poisson-verteilt ist), # zu welchem Parameter lambda? gesAnzBeobInt <- sum(AnzBeobInt) gesAnzBeobInt # =2608 gesAnzZerfaelle <- sum(AnzBeobInt*Zerfaelle) gesAnzZerfaelle # ein naheliegender Wert (sog. "Momentenschätzer") lambdahut <- gesAnzZerfaelle/gesAnzBeobInt lambdahut erw <- c(gesAnzBeobInt*dpois(0:9, lambda=lambdahut), gesAnzBeobInt*ppois(9, lambda=lambdahut, lower.tail=F)) Zerf2 <- cbind(Zerf, Erwartet=erw) Zerf2 # Füge erwartete Anzahlen dem Plot hinzu: points(1.2*(0:10)+0.75, erw, col="red") # Um die Güte der Anpassung einzuschätzen bietet sich die # sog. chi-Quadrat-Statistik an: chiqua <- sum(((AnzBeobInt-erw)^2)/erw) chiqua # pchisq(chiqua, df=10, lower.tail=FALSE) # Übrigens: Der Momentenschätzer ist hier zugleich der ML-Schätzer. detach(Zerf) # 2. von Bortkiewicz' # Angaben über in der preußischen Armee in den Jahren 1875-1894 # (in insgesamt 10 Armeekorps) durch "Schlag eines Pferdes getötete", # siehe Ladislaus von Bortkiewicz, # Das Gesetz der kleinen Zahlen, Teubner, 1898, S. 24f. # Bem.: Das Urheberrecht für dieses Buch ist ausgelaufen, es ist z.B. # frei auf Google books erhältlich. pfunf <- c(109, 65, 22, 3, 1, 0) sum(pfunf) # =200 (insgesamt 200 "Korpsjahre") # wir passen lambda wie oben an: lambdahut <- sum(pfunf*(0:5))/sum(pfunf) lambdahut Bortk <- data.frame(Tote=0:5, AnzBeob=pfunf, Erw=200*c(dpois(0:4, lambda=lambdahut), ppois(4, lambda=lambdahut, lower.tail=F))) Bortk # ############################################################## # # Zum Poissonprozess # # ############################################################## # Vorüberlegung: Wartezeit auf den ersten Erfolg # in einer (reskalierten) Münzwurffolge # (oder: eine reskalierte geometrische Verteilung approximiert die # Exponentialverteilung) # n, lambda Parameter (n >> 1) # X_1, X_2, ... u.a. Bernoulli(lambda/n), # S_n := min { i : X_i = 1 } - 1, # dann konvergiert S_n/n in Verteilung gegen Exp(lambda): # Wir vergleichen die Verteilungsfunktionen für verschiedene n: lambda <- 1.5 n.werte <- c(10, 20, 50, 100) curve(pexp(x, lambda), xlim=c(-0.01,2), xlab="x", ylab="Ws für Wert <= x", main=paste("Verteilungsfkt. Exp(",lambda,") und Approximationen")) farben <- c("blue","darkred","darkgreen","gray") for (i in 1:length(n.werte)) { n <- n.werte[i] curve(pgeom(n*x, prob=lambda/n), add=T ,col=farben[i], n=500) } legend("bottomright", legend=c("Exp",n.werte), col=c("black",farben), lty=1) # lambda <- 1.5 n.werte <- c(10, 20, 50, 100) curve(pexp(x, lambda), xlim=c(-0.01,2), xlab="x", ylab="Ws für Wert <= x", main=paste("Verteilungsfkt. Exp(",lambda,") und Approximationen")) farben <- c("blue","darkred","darkgreen","gray") for (i in 1:length(n.werte)) { n <- n.werte[i] curve(pgeom(n*x, prob=lambda/n), add=T ,col=farben[i], n=500) } legend("bottomright", legend=c("Exp",n.werte), col=c("black",farben), lty=1) # # Poissonprozess definiert aus unabhaengig exponentialverteilten Wartezeiten, # Simulation von Pfaden lambda <- 0.4 # beobachtungsdauer <- 100 # Einen Pfad eines Poissonprozesses für geg. Anz. Schritte simulieren simPPpfad <- function(lambda, schritte) { data.frame(x=c(0, cumsum(rexp(schritte, rate=lambda))), y=0:schritte) } # zeichne einen solchen Pfad (Parameter ... werden "durchgereicht"): malePPpfad <- function(pf, ...) { points(pf$x, pf$y, type="s", ...) } beobachtungsdauer <- 100 # Neues Zeichenfenster plot(0, xlim=c(0,beobachtungsdauer), ylim=c(0,1.1*beobachtungsdauer*lambda), type="n", xlab="Zeit", ylab="") # ein ad hoc-Ansatz für die Anzahl notwendiger Schritte: pfad <- simPPpfad(lambda, round(1.1*beobachtungsdauer*lambda)) malePPpfad(pfad) # füge einige Pfade hinzu: farben <- c("blue","darkred","darkgreen","gray") for (i in 1:4) { pfad <- simPPpfad(lambda, round(1.1*beobachtungsdauer*lambda)) malePPpfad(pfad, col=farben[i]) } # alternativ simuliere einen Pfad, bis eine gewisse Zeit verstrichen ist: simPPpfad.mitzeit <- function(lambda, maxzeit) { x <- 0 y <- 0 ort <- 0; zeit <- 0 while(zeit