# # Stochastik-Praktikum, WS 2019/20, JGU Mainz # 25.11.2019 ## ############################################################################## ## heute: ## Illustrationen des Gesetzes der großen Zahlen ## (und ein "Gegenbeispiel") ## ## Monte-Carlo-Integration: "naiv" und ## mit Kontrollvariable zur Varianzreduktion: ## (ggfs. auch: Importance sampling als Methode zur ## Varianzreduktion) # ################################################################# # # Illustrationen des Gesetzes der großen Zahlen # # Zur Einstimmung: # Konvergenz der empirischen Häufigkeiten beim Münzwurf beobachten # Wir simulieren n unabhängige binomialverteilte Zufallsvariablen/Münzwürfe # mit Erfolgsw'keit p. # Dann plotten wir für m=1,2,...,n die relativen Häufigkeiten von Erfolgen # in den ersten m Versuchen. n <- 100 p <- 0.5 x <- rbinom(n,size=1,prob=p) cumsum(x)/(1:n) plot(1:n,cumsum(x)/(1:n), xlab="Anz.Summanden m",ylab="Empir.Mittelwert (X_1+...+X_m)/m") abline(h=p, col='red') # Wir zeigen 4 Replikate dieses Experiments in einem Bild: GesdgrZ_binomial<-function(n,p){ x11() par(mfrow=c(2,2)) for(i in 1:4) { #Wir simulieren n Münzwürfe/Binomialverteilte ZV mit Erfolgsw'keit p Folge<-rbinom(n,size=1,prob=p) #Wir berechnen die empirische Häufigkeiten Frequenz<-cumsum(Folge)/seq(1,n) #nun erzeugen wir das passende Bild ylab<-"Empir.Mittelwert (X_1+...+X_m)/m" xlab<-"Anz.Summanden m" plot(Frequenz, ylab=ylab,xlab=xlab,type = "l",ylim=c(0,1)) abline(h=p,col="red") } par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", paste0("p=",p,",",n," Versuche"))) } # Einmal mit 1000 Versuchen und p=0.5 GesdgrZ_binomial(1000,0.5) # Noch einmal mit 1000 Versuchen und p=0.75 GesdgrZ_binomial(1000,0.75) # Wir sollten bei 1000 Versuchen die Konvergenz gegen p gut sehen koennen. # Bei 100 Versuchen koennen wir noch deutliche(re) Abweichungen von p beobachten. GesdgrZ_binomial(100,0.5) # Jetzt dasselbe mit der Normalverteilung: # Wir simulieren n unabhaengige normalverteilte Zufallsvariablen mit # Erwartungswert mu und Varianz sigma^2 # Dann plotten wir für m=1,2,...,n die relativen Haeufigkeiten von Erfolgen # in den ersten m Versuchen. # Wir zeigen (wieder) 4 Replikate dieses Experiments in einem Bild: GesdgrZ_normal<-function(n,mu,sigma,zentrieren){ x11() par(mfrow=c(2,2)) ylab<-"Empir.Mittelwert (X_1+...+X_m)/m" xlab<-"Anz.Summanden m" for(i in seq(1,4)){ Folge<-rnorm(n,mean=mu,sd=sigma) Frequenz<-cumsum(Folge)/seq(1,n) #Zentriere Bild, falls gewünscht if(zentrieren){ plot(Frequenz, ylab=ylab,xlab=xlab,type = "l",ylim=c(mu-1,mu+1)) } else{plot(Frequenz, ylab=ylab,xlab=xlab,type = "l")} abline(h=mu,col="red") } par(mfrow=c(1,1)) title(main=c("Illustration des Gesetzes der großen Zahlen:", paste0("mu=",mu,",sigma=",sigma,",",n," Versuche"))) } # Einmal mit 1000 Versuchen und mu=0, sigma=1 # Ohne Zentrieren GesdgrZ_normal(1000,mu=0,sigma = 1,FALSE) # Einmal mit 1000 Versuchen und mu=5, sigma=1 # Ohne Zentrieren GesdgrZ_normal(1000,mu=5,sigma = 3,FALSE) # Vergleich von Konvergenzgeschwindkeit # 100 Versuchen und mu=0, sigma=1 # Zum besseren Vergleich mit Zentrieren GesdgrZ_normal(100,mu=0,sigma = 1,TRUE) # Einmal mit 100 Versuchen und mu=0, sigma=sqrt(10) GesdgrZ_normal(100,mu=0,sigma = 3.162,TRUE) # Mit ein bisschen Glück sollten wir sehen, #d ass die Schwankungen für sigma=sqrt(10) deutlich stärker sind. # ################################################################### # Ohne Erwartungswert scheitert das Gesetz der großen Zahlen: # Ein Beispiel einer Verteilung(sklasse) mit unendlichem Erwartungswert # Sei 0 < 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>1, # 3) die inverse Verteilungsfunktion ist # F_a^{-1}(u) = (1-u)^(-1/a). # 4) Es ist # integrate(x f_a(x), 1, infinity) = # integrate(a x^(-a), 1, infinity) = infinity, # d.h. es gibt keinen Erwartungswert # Wir simuliere n ZVn mit Dichte f_a, indem wir F_a und die Inversionmethode verwenden. # Danach versehen wir diese mit einem zufaelligen Vorzeichen. # So erhalten wir eine symmetrische ZV ohne Erwartungswert rfa <- function(n, a=0.5) { (1-runif(n, min=0, max=1))^(-1/a)*sample(c(-1,1),size=n,replace=TRUE) } replikate <- 5000 a <- 0.5 ## auch andere Werte probieren, z.B. a <- 0.25 oder a <- 0.9 x <- rfa(replikate, a) mean(x) hist(x) ## Wir sehen: es gibt einige wenige sehr extreme Werte min(x); max(x) par(mfrow=c(1,2)) plot(x,xlab="Simulationen",ylab="Werte der Simulationen") plot(log(abs(x),base=10),xlab="Simulationen",ylab="Logarithmus der Werte der Simulationen zur Basis 10") par(mfrow=c(1,1)) 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: x <- rfa(replikate, a) y<-cumsum(x)/(1:replikate) par(mfrow=c(1,2)) plot(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') plot(x,xlab="Simulationen",ylab="Werte der Simulationen") abline(h=0, col='red') par(mfrow=c(1,1)) # Vergleich mit der Normalverteilung x<-rnorm(replikate) y<-cumsum(x)/(1:replikate) par(mfrow=c(1,2)) plot(y, type="l", xlab=paste("Anz. Summanden n"), ylab="Empir. Mittelwert (X_1+...+X_n)/n", title="Standardnormalverteilung") abline(h=0, col='red') plot(x,xlab="Simulationen",ylab="Werte der Simulationen") abline(h=0, col='red') par(mfrow=c(1,1)) # Was passiert fuer 1