# # Sitzung des Stochastik-Praktikums, WS 2013/2014, JGU Mainz # 2.12.2013 # Illustrationen des Gesetzes der großen Zahlen # "Gegenbeispiele": Cauchy-Verteilung, Polya-Urne # Konvergenz der empirischen Verteilungen (Satz von Glivenko-Cantelli) # ################################################################# # # Illustrationen des Gesetzes der großen Zahlen # # Zur Einstimmung: # Konvergenz der empirischen Häufigkeiten beim Münzwurf beobachten AnzWuerfe <- 1000 p <- 0.5 # Erfolgsw'keit der Münze # Folge.MA <- rbinom(AnzWuerfe, 1, p) Folge.MB <- rbinom(AnzWuerfe, 1, p) # Wir betrachten die relativen Häufigkeiten der # Erfolge nach n Würfen (als Folge in n): # relH.MA[n] = (Anz. Erfolge unter Würfen 1,...,n)/n relH.MA <- cumsum(Folge.MA)/(1:AnzWuerfe) relH.MB <- cumsum(Folge.MB)/(1:AnzWuerfe) par(mfrow=c(1,2)) # betrachten wir zunächst Anfangsstücke der Länge 50: plot(1:50, relH.MA[1:50], type="p", xlab="Anz. Würfe", ylab="rel. Häufigk. d. Erfolge", ylim=c(0,1), sub="Versuchsreihe A") abline(h=p, lty=2) plot(1:50, relH.MB[1:50], type="p", xlab="Anz. Würfe", ylab="rel. Häufigk. d. Erfolge", ylim=c(0,1), sub="Versuchsreihe B") abline(h=p, lty=2) par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", "Zwei Münzwurf-Versuchsreihen (jew. die ersten 50 Würfe)")) # Um allzu repetitiven Code zu vermeiden lohnt es sich, # eine Funktion für unsere Plots zu definieren: male.haeufigk <- function(Folge, Untertitel="") { plot(Folge, type="l", xlab="Anz. Würfe", ylab="rel. Häufigk. d. Erfolge", ylim=c(0,1), sub=Untertitel) } # nun jeweils die ersten 200 Würfe: x11() # öffne ein neues Grafikfenster par(mfrow=c(1,2)) male.haeufigk(relH.MA[1:200], "Versuchsreihe A") abline(h=p, lty=2) male.haeufigk(relH.MB[1:200], "Versuchsreihe B") abline(h=p, lty=2) par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", "Zwei Münzwurf-Versuchsreihen (jew. die ersten 200 Würfe)")) # Schliesslich alle Würfe: x11() # öffne ein neues Grafikfenster par(mfrow=c(1,2)) male.haeufigk(relH.MA, "Versuchsreihe A") abline(h=p, lty=2) male.haeufigk(relH.MB, "Versuchsreihe B") abline(h=p, lty=2) par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", paste("Zwei Münzwurf-Versuchsreihen (jew.", AnzWuerfe, "Würfe)"))) # Die beobachteten relativen Häufigkeiten nach AnzWuerfe Würfen waren jeweils: relH.MA[AnzWuerfe]; relH.MB[AnzWuerfe] # das alles gilt natürlich auch für unfaire Münzen: p <- 0.4 # Folge.MA <- rbinom(AnzWuerfe, 1, p) Folge.MB <- rbinom(AnzWuerfe, 1, p) # Die relativen Häufigkeiten der Erfolge nach n Würfen (als Folge in n): relH.MA <- cumsum(Folge.MA)/(1:AnzWuerfe) relH.MB <- cumsum(Folge.MB)/(1:AnzWuerfe) par(mfrow=c(1,2)) male.haeufigk(relH.MA, "Versuchsreihe A") abline(h=p, lty=2) male.haeufigk(relH.MB, "Versuchsreihe B") abline(h=p, lty=2) par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", paste("Zwei (unfaire) Münzwurf-Versuchsreihen (jew.", AnzWuerfe, "Würfe)"))) # ################################################################# # Intermezzo/Wiederholung: die Normalverteilung # mu reelle Zahl, sigma>0: mu <- 0.0; sigma <- 1.0 # Die Funktion normaldichte <- function(x) (1/sqrt(2*pi)*sigma)*exp(-(x-mu)^2/(2*sigma^2)) # ist die Dichte der Normalverteilung # (mit Mittelwert mu und Varianz sigma^2, "Gauß-Glocke") curve(normaldichte, xlim=c(-3.5,3.5), xlab="x", ylab="Dichte", col="red", lwd=2, main="Dichte der (Standard-)Normalverteilung") abline(v=mu); abline(v=mu-sigma, lty=2); abline(v=mu+sigma, lty=2) text(mu, 0.01, "mu", pos=4); text(mu+sigma, 0.01, "mu+sigma", pos=4); text(mu-sigma, 0.01, "mu-sigma", pos=2); # Übrigens: R kennt diese Funktion als dnorm # (das "d" steht in der R-Funktionsnamensphilosophie für "density") normaldichte(2.0)-dnorm(2.0) ?dnorm # Rs Online-Hilfe über die Normalverteilung # allgemeiner also mu <- 5.5; sigma <- 0.7 normaldichte(1.0) - dnorm(1.0, mean=mu, sd=sigma) # Schauen wir uns an, was die Parameter mu und sigma (für die Dichte) bewirken: plot(0, xlim=c(-3,7), ylim=c(0,1), xlab="x", ylab="Dichte", main="Verschiedene Normaldichten", type="n") # dieser "leere" Plot # erzeugt das Koordinatensystem curve(dnorm(x, mean=0, sd=1.0), col="blue", lwd=2, add=TRUE) curve(dnorm(x, mean=0, sd=0.4), col="red", lwd=2, add=TRUE, n=250) #(siehe Online-Hilfe für den "n-Parameter", er regelt, wie glatt das Bild wirkt curve(dnorm(x, mean=3, sd=1.0), col="darkgreen", lwd=2, add=TRUE) legend("topright", legend=c("mu=0, sigma=1", "mu=0, sigma=0.4", "mu=3, sigma=1"), lwd=2, col=c("blue", "red", "darkgreen")) # R generiert normalverteilte (Pseudo-)Zufallszahlen mittels rnorm # (das "r" steht in der R-Funktionsnamensphilosophie für "random"): rnorm(1) rnorm(1, mean=5, sd=2) # man kann gewünschten Erwartungswert u. Varianz angeben # 12 normalverteilte (Pseudo-)Zufallszahlen rnorm(12, mean=5, sd=2) # ################################################################# # Zurück zur Illustration des Gesetzes der großen Zahlen: # (Standard-)Normalverteilte Summanden Anz <- 10000 X.A <- rnorm(Anz) X.B <- rnorm(Anz) # Schauen wir die Werte an # (und sehen, dass wir keine Struktur erkennen ...) par(mfrow=c(1,2)) plot(X.A, xlab="Index", ylab="Wert", ylim=c(-4,4), type="p", pch=20, cex=0.5, sub="Folge A") abline(h=0, lty=2, col="red") plot(X.B, xlab="Index", ylab="Wert", ylim=c(-4,4), type="p", pch=20, cex=0.5, sub="Folge B") abline(h=0, lty=2, col="red") par(mfrow=c(1,1)) title(main=paste("Zweimal jeweils",Anz,"normalverteilte Zufallswerte")) # normierte (Partial-)Summen: # S.A[n] = (X.A[1]+...+X.A[n])/n, etc. S.A <- cumsum(X.A)/(1:Anz) S.B <- cumsum(X.B)/(1:Anz) # Wir definieren eine Funktion für unsere Plots: male.partialsummen <- function(Folge, Untertitel="", typ="l") { plot(Folge, type=typ, xlab="Anz. Summanden", ylab="Summe/Anzahl", sub=Untertitel, ylim=c(-2,2)) } # Summe/Anzahl: die ersten 50 Werte x11() # öffne ein neues Grafikfenster par(mfrow=c(1,2)) male.partialsummen(S.A[1:50], "Versuchsreihe A", "p") abline(h=0, lty=2) male.partialsummen(S.B[1:50], "Versuchsreihe B", "p") abline(h=0, lty=2) par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", "Summe/Anzahl für zwei u.i. normalverteilte Versuchsreihen", "(jew. die ersten 50 Werte)")) # dasselbe für die ersten 200 Werte: x11() par(mfrow=c(1,2)) male.partialsummen(S.A[1:200], "Versuchsreihe A") abline(h=0, lty=2) male.partialsummen(S.B[1:200], "Versuchsreihe B") abline(h=0, lty=2) par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", "Summe/Anzahl für zwei u.i. normalverteilte Versuchsreihen", "(jew. die ersten 200 Werte)")) # und für alle Werte: x11() par(mfrow=c(1,2)) male.partialsummen(S.A[1:1000], "Versuchsreihe A") abline(h=0, lty=2) male.partialsummen(S.B[1:1000], "Versuchsreihe B") abline(h=0, lty=2) par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", "Summe/Anzahl für zwei u.i. normalverteilte Versuchsreihen", "(jew. die ersten 1000 Werte)")) # ######################################################## # # Gesetz der großen Zahlen und die empirische Verteilung: a <- 0.5; b <- 2.0 # Vorbetrachtung: Relativer Anteil der Werte in [a,b] als Funktion von n. # Die Folge 1(a < X[i] <= b), i=1,2,3,... ist eine Münzwurffolge mit # Erfolgsw'keit P(a < X <= b) = pnorm(a) - pnorm(b). # In R ist pnorm die Verteilungsfunktion der (Standard-)Normalverteilung # (das "p" steht in der R-Funktionsnamensphilosophie # für "probability distribution function") pnorm(0) pnorm(1.96) pnorm(b)-pnorm(a) zwischenaundb <- ifelse((a < X.A) & (X.A <=b), 1, 0) # beachte: (a < X.A) & (X.A <=b) erzeugt einen _Vektor_ von Wahrheitswerten, # indem es eintragsweise auf X.A wirkt # (während (a < X.A) && (X.A <=b) nur einen Wahrheitswert erzeugte) haeufigk <- cumsum(zwischenaundb)/(1:Anz) plot(haeufigk[1:1000], type="l", ylim=c(0,1), xlab="n", ylab=paste("rel. Anteil zwischen",a,"und",b,"bis n"), main="Versuchsreihe A (u.i.v., standard-normal)") abline(h=pnorm(b)-pnorm(a), lty=2) # das funktioniert natürlich auch für jede andere Wahl von a und b, z.B. a <- -1; b <- 1 # (und damit ggf. obigen Code nochmals durchlaufen lassen) # Nebenbei bemerkt: Bei der Gaußglocke liegt ca. 2/3 des # Wahrscheinlichkeitsgewichts zwischen mu-sigma und mu+sigma: curve(dnorm, from=-3, to=3, main="(Standard-)Normaldichte"); abline(v=c(-1,1)) pnorm(1)-pnorm(-1) # # Konvergenz der empirischen Verteilungen: Histogramme # betrachte die ersten 250 Werte: x11() par(mfrow=c(1,2)) br=seq(from=-4,to=4, by=0.2) hist(X.A[1:250], prob=TRUE, xlim=c(-4,4), breaks=br, xlab="Versuchsreihe A", main="") curve(dnorm, col="red", add=TRUE) hist(X.B[1:250], prob=TRUE, xlim=c(-4,4), breaks=br, xlab="Versuchsreihe B", main="") curve(dnorm, col="red", add=TRUE) par(mfrow=c(1,1)) title(main="Histogramm der ersten 250 Werte und Normaldichte") # betrachte die ersten 1000 Werte: x11() par(mfrow=c(1,2)) br=seq(from=-4,to=4, by=0.2) hist(X.A[1:1000], prob=TRUE, xlim=c(-4,4), breaks=br, xlab="Versuchsreihe A", main="") curve(dnorm, col="red", add=TRUE) hist(X.B[1:1000], prob=TRUE, xlim=c(-4,4), breaks=br, xlab="Versuchsreihe B", main="") curve(dnorm, col="red", add=TRUE) par(mfrow=c(1,1)) title(main="Histogramm der ersten 1000 Werte und Normaldichte") # und schliesslich alle Werte: x11() par(mfrow=c(1,2)) br=seq(from=-4.6,to=4.6, by=0.2) hist(X.A, prob=TRUE, xlim=c(-4,4), breaks=br, xlab="Versuchsreihe A", main="") curve(dnorm, col="red", add=TRUE) hist(X.B, prob=TRUE, xlim=c(-4,4), breaks=br, xlab="Versuchsreihe B", main="") curve(dnorm, col="red", add=TRUE) par(mfrow=c(1,1)) title(main="Histogramm und Normaldichte") # # Das stimmt natürlich nicht nur für normalverteilte Zufallsgrößen. # Beispiel: Exponentialverteilte Daten # # (Die Konvergenz der empirischen Verteilung stimmt stets für # unabhängige, identisch verteilte Zufallsgrößen, dies ist eine # der theoretischen Grundlagen der Statistik: # wenn man ein Experiment beliebig oft unabhängig wiederholen kann, kann # man die Verteilung des Ergebnisses prinzipiell beliebig genau kennen lernen). lambda <- 2.0 daten <- rexp(10000, rate=lambda) x11() par(mfrow=c(1,3)) br <- seq(from=0, to=max(5, max(daten+0.1)), by=0.1) hist(daten[1:250], prob=T, main="", breaks=br, xlim=c(0,4), ylim=c(0,2.4)) # die Dichte der Exponentialverteilung: dichtefkt <- function(x) lambda*exp(-lambda*x) # (R kennt diese Funktion natürlich als dexp) curve(dichtefkt,col="red",add=TRUE) hist(daten[1:1000], prob=T, main="", breaks=br, xlim=c(0,4), ylim=c(0,2.4)) curve(dexp(x,rate=lambda),col="red",add=TRUE) hist(daten[1:10000], prob=T, main="", breaks=br, xlim=c(0,4), ylim=c(0,2.4)) curve(dexp(x,rate=lambda),col="red",add=TRUE) par(mfrow=c(1,1)) title(main=c(paste("U.i. Exp(",lambda,")-verteilte Daten: Histogramm"), "und theoretische Dichte")) # ################################################################### # Ein Beispiel einer Verteilung(sklasse) mit endlichem Erwartungswert, # aber moeglicherweise unendlicher Varianz: # Sei a > 1, # 1) f_a(x) := a x^(-1-a), x>=1 ist eine W'dichte, # 2) die zugehoerige Verteilungsfunktion ist # F_a(x) = 0 fuer x<=1, = 1-x^(-a) fuer x>0, # 3) die inverse Verteilungsfunktion ist # F_a^{-1}(u) = (1-u)^(-1/a). # 4) Fuer a>1 ist # integrate(x f_a(x), 1, infinity) = # integrate(a x^(-a), 1, infinity) = a/(a-1) < infinity, # d.h. der Erwartungswert ist a/(a-1) # (fuer a <= 1 gibt es keinen Erwartungswert) rfa <- function(n, a=1.5) { # simuliere n ZVn mit Dichte f_a, # zentriert so dass der EW = 0 (1-runif(n, min=0, max=1))^(-1/a)-a/(a-1) } replikate <- 5000 a <- 1.9 ## auch andere Werte probieren, z.B. a <- 1.25 oder a <- 1.6 x <- rfa(replikate, a) mean(x) hist(x) ## Wir sehen: es gibt einige wenige sehr extreme Werte max(x) hist(x, xlim=c(-a/(a-1)-0.1,5), breaks=c(-Inf,seq(-a/(a-1)-0.1,5,by=0.1),Inf)) abline(v=mean(x), lwd=2, col='red') mean(x[x>0]) # Schauen wir uns an, wie sich empirische Mittelwerte # (als Funktion der Anzahl summierter Kopien) hier verhalten: schr <- 1000; mul <- 1000 replikate <- schr*mul x <- rfa(replikate, a) ausw <- mul*(1:schr) y <- cumsum(x)[ausw]/ausw plot(ausw, y, type="l", xlab=paste("Anz. Summanden n"), ylab="Empir. Mittelwert (X_1+...+X_n)/n", main=paste("Kopien von X mit Dichte f_",a,sep="")) abline(h=0, col='red') ## ################################################################### ## ################################################################### # # Ohne Erwartungswert scheitert das Gesetz der großen Zahlen: # Beispiel Cauchy-Verteilung # Bem.: Die (Standard-)Cauchy-Verteilung hat Dichte (1/pi)*1/(1+x^2) # (beachte: (d/dx) arctan(x) = 1/(1+x^2), # d.h. die Verteilungsfunktion ist arctan(x)/pi+1/2) # R kennt diese als dcauchy bzw. pcauchy # Die Cauchy-Dichte ist wesentlich "schwerschwänziger" als die Normaldichte: curve(dnorm, from=-5, to=5, col="red", xlab="x", ylab="Dichte") curve(dcauchy, col="blue", add=TRUE) legend("topright", lwd=1, col=c("red","blue"), legend=c("Normal","Cauchy")) # ... was man auf der log-Skale extrem deutlich sieht: curve(log(dnorm(x)), from=-20, to=20, col="red", xlab="x", ylab="log(Dichte)") curve(log(dcauchy(x)), col="blue", add=TRUE) legend("topright", lwd=1, col=c("red","blue"), legend=c("Normal","Cauchy")) Anz <- 10000 daten.Cauchy <- rcauchy(Anz) min(daten.Cauchy); max(daten.Cauchy) hist(daten.Cauchy, probability=TRUE) boxplot(daten.Cauchy) hist(daten.Cauchy[daten.Cauchy > -10 & daten.Cauchy < 10], probability=TRUE, ylim=c(0,.4)) curve(dcauchy, col="red", add=TRUE) norm.summen.C <- cumsum(daten.Cauchy)/(1:Anz) plot(norm.summen.C[1:100]) abline(h=0, lty=2) plot(norm.summen.C[1:1000], type="l") abline(h=0, lty=2) plot(norm.summen.C[1:10000], type="l") abline(h=0, lty=2) # Die Unabhängigkeit der Summanden kann man nicht (völlig) # aufgeben: # Beispiel Polya-Urne # (hier zeigt sich sozusagen ein "zufälliges Gesetz der großen Zahlen") # vergl. auch Aufg. 1 (Übungsblatt 5) 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)