# # Stochastik-Praktikum, WS 2016/2017, JGU Mainz # 12.12.2016 ## ############################################################################## ## heute: ## Illustrationen des Gesetzes der großen Zahlen ## und ein "Gegenbeispiel" ## ## Monte-Carlo-Integration: "naiv" und ## mit Methoden zur Varianzreduktion: ## Importance sampling, Kontrollvariable # ################################################################# # # 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. #Dies machen wir insgesamt 4 mal und plotten alles auf eine Seite starkes_gesetz_binomial<-function(n,p){ x11() par(mfrow=c(2,2)) for(i in seq(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 starkes_gesetz_binomial(1000,0.5) #Nocheinmal mit 1000 Versuchen und p=0.75 starkes_gesetz_binomial(1000,0.75) #Wir sollten bei 1000 Versuchen die Konvergenz gegen p gut sehen können. #Aber bei 100 Versuchen können wir noch starke Abweichungen von p beobachten. starkes_gesetz_binomial(100,0.5) #Jetzt dasselbe mit der Normalverteilung #Wir simulieren n unabhängige normalverteilte Zufallsvariablen mit #Erwartungswert mu und Varianz sigma^2 #Dann plotten wir für m=1,2,...,n die relativen Häufigkeiten von Erfolgen #in den ersten m Versuchen. #Dies machen wir insgesamt 4 mal und plotten alles auf eine Seite starkes_gesetz_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 starkes_gesetz_normal(1000,mu=0,sigma = 1,FALSE) #Einmal mit 1000 Versuchen und mu=5, sigma=1 #Ohne Zentrieren starkes_gesetz_normal(1000,mu=5,sigma = 3,FALSE) #Vergleich von Konvergenzgeschwindkeit #100 Versuchen und mu=0, sigma=1 #Zum besseren Vergleich mit Zentrieren starkes_gesetz_normal(100,mu=0,sigma = 1,TRUE) #Einmal mit 100 Versuchen und mu=0, sigma=sqrt(10) starkes_gesetz_normal(100,mu=0,sigma = 3.162,TRUE) #Mit ein bisschen Glück sollten wir sehen, #dass 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 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 für 1