#-------------------------------------# # week4b: R-Simulation des # # zentralen Grenzwertsatzes # #-------------------------------------# Wir wollen den zentralen Grenzwertsatz jetzt mit Hilfe einer R-Simulation verifizieren. Erinnern wir uns eben an die Aussage des Satzes: Zentraler Grenzwertsatz: Es sei X_1, X_2, X_3, ... eine Folge von unabhaengigen, identisch verteilten Zufalls- variablen mit E[ X_k ] = mu V[ X_k ] = sigma^2 fuer alle k. Wir betrachten die Summe S_n = X_1 + X_2 + . . . + X_n die einen Erwartungswert E[ S_n ] = n*mu und eine Varianz V[ S_n ] = n*sigma^2 hat. Schliesslich definieren wir die normiertn Groessen Z_n := ( S_n - E[S_n] ) / sqrt( V[S_n] ) = ( S_n - n*mu ) / ( sqrt(n)*sigma ) Dann gilt: Die Verteilung von Z_n konvergiert gegen die Standard-Normalverteilung. Das heisst lim_{n to infty} Prob[ Z_n in [x,x+dx) ] = 1/sqrt(2 pi) * exp[ - x^2/2 ] * dx Bemerkung: Diese Formulierung des zentralen Grenzwert- satzes bekommen wir, wenn wir fuer das allgemeine f im Theorem 2.2.4 die Funktion f(z):= chi( z in [x,x+dx) ) waehlen. Schauen wir uns jetzt die Sachen in R an: Ueberpruefen wir das zunaechst vielleicht mal fuer gleichverteilte Zufallszahlen, sagen wir, gleich- verteilt auf dem Intervall [10,20]. Wir haben hier im wesentlichen 2 Parameter, einmal das kleine n, das ist das n aus dem zentralen Grenz- wertsatz, was dann also nach unendlich gehen soll. Und dann wollen wir ja folgendes machen: Wir simu- lieren die Z_n's, fuer, etwa, n=50 oder n=100, und wollen dann verifizieren, dass die Verteilung der Z_n's naeherungsweise aussieht wie eine Gauss'sche Glockenkurve. Wenn wir uns die Verteilung einer zufaelligen Groesse anschauen wollen, bedeuted das, dass wir moeglichst viele von den Z_n's simulieren wollen, und uns dann ein Histogramm der Z_n's an- schauen wollen. Also brauchen wir noch einen zweiten Parameter fuer die Anzahl der Z_n's, die wir simulieren wollen, nehmen wir dafuer etwa den Buchstaben gross N. Fuer N = 1000 sehen die Histogramme mitunter immer noch ein bischen ruckelig aus, fuer N = 10000 sind die in der Regel immer schoen smoothig. Also haben wir die beiden Parameter n und N, mit denen man dann natuerlich ein bischen herumexperi- mentieren kann. ### Start R-Session ### n = 100 N = 10000 # Wir brauchen insgesamt n*N Zufallszahlen: x = runif(n*N, min=10, max=20) x # der Erwartungswert ist offensichtlich mu = 15, aber # fuer die Varianz sigma^2 ist das schon nicht mehr # so unmittelbar klar, was das ist, deshalb nehmen wir # einfach das Stichprobenmittel und die Stichproben- # Varianz oder Standardabweichung. Wir brauchen ja # Werte fuer mu und sigma, um die Z_n's berechnen zu # koennen: mu = mean(x) mu # sehr nahe bei 15 sigma = sd(x) sigma # etwa 2.88 (20-10)/sqrt(12) # etwa dasselbe # Um die Z_n's zu berechnen, schreiben wir die n*N # Zufallszahlen in eine N x n - Matrix: X = matrix(x, nrow=N, ncol=n) head(X) tail(X) # sieht ok aus # Wir berechnen die S_n's, das sind dann die Zeilen- # summen der Matrix X: Sn = rowSums(X) Sn # ist ein Vektor der Lange N=10000 mit # Werten um die 1500 = n*mu, das passt # Wir berechnen die Z_n's: Zn = (Sn - n*mu) / ( sqrt(n)*sigma ) Zn hist(Zn) hist(Zn,breaks=50) # Um das mit einer W'keitsdichte vergleichen zu koennen, # brauchen wir eine relative, keine absolute, Skalierung, # so dass der Flaecheninhalt von allen Rechtecken = 1 ist; # von der Form her bleibt das Bild dasselbe, aber die # Zahlen an der y-Achse, an der vertikalen Achse, aendern # sich: hist(Zn, breaks=50, prob=TRUE) # in dieses Bild soll jetzt ebenfalls die Gauss'sche # Glockenkurve, etwa in rot, mit rein: curve(dnorm(x), col="red", add=TRUE) # ok, das passt.. # Damit haben wir also den zentralen Grenzwertsatz fuer # gleichverteilte Zufallszahlen durch Simulation veri- # fizieren koennen. # JETZT: Wir wollen folgendes machen: In R gibt es ja eine # Reihe von eingebauten W'keitsverteilungen Vert, nach denen # wir dann auch mit dem rVert()-Befehl Zufallszahlen erzeugen # koennen. Wir moechten jetzt die Verteilung oder den Ver- # teilungsnamen (und die Parameter n und N) beliebig vorgeben # koennen, und dann soll automatisch das entsprechende # Histogramm dazu erzeugt werden. Die in R eingebauten Ver- # teilungen koennen Sie sich hier http://hsrm-mathematik.de/SS2023/semester4/Stochastik2/W'keitsverteilungen-in-R.pdf # noch einmal anschauen. Wir wollen jetzt also den Verteilungs- # name in eine Textvariable schreiben, etwa Vert = "unif" # eine Text-Variable Vert # und wollen damit dann Zufallszahlen erzeugen: rVert = paste("r",Vert) rVert # ok, das Leerzeichen muss weg.. rVert = paste("r",Vert,sep="") rVert # passt jetzt.. x = runif(100) x y = rVert(100) # geht so nicht.. # Ok, das kann so nicht klappen, weil rVert einfach eine Text- # Variable ist, aber nicht als Funktionsaufruf erkannt wird. # Fuer sowas gibt es typischerweise einen "Call"-Befehl oder # aehnliches. Schauen wir mal in die Hilfe, ?call # ok, da gibt es was.. # versuchen wir mal y = call(rVert,100) y # okay.. ?call # es gibt ein "do.call", nehmen wir das mal: y = do.call(rVert,100) # geht immer noch nicht.. # das 2.Argument muss eine Liste sein, die aus den # Variablen besteht, die an rVert = runif uebergeben # werden sollen: im einfachsten Fall wollen wir nur # eine 100 uebergeben, fuer 100 Zufallszahlen: plist = list(100) # "plist" etwa fuer "parameter list" #just for info: mode(plist) # ist vom Datentyp "Liste", ok str(plist) y = do.call(rVert,plist) y # ok, das klappt # Das funktioniert jetzt, aber wir haben gleichverteilte # Zufallszahlen auf dem Intervall [0,1], das sind die # Default-Werte wenn man min und max nicht angibt; wir # wollen ja Zahlen zwischen 10 und 20 haben. # # Allgemein: Die Parameter fuer die gewaehlte W'keits- # verteilung (also hier min und max) muss man mit in die # plist packen: rm(plist) # loeschen plist = list(100, 10, 20) plist # Listen-Elemente werden immer mit doppelten # eckigen Klammern [[k]] angegeben und aufgerufen y = do.call(rVert,plist) y # okay, das tut's jetzt # Zurueck zum zentralen Grenzwertsatz: Wir wollen jetzt also # eine Funktion schreiben, die als Argumente die Parameter # n, N, Vert, plist bekommt, und dann das entsprechende # Histgramm generiert. Das geht jetzt folgendermassen: # # Waehlen wir einen Funktionsnamen, etwa SimuCLT # ("CLT" fuer "central limit theorem"), dann: # # wir wollen N mal die Zufallsvariable Z_n simulieren, wir # waehlen etwa n = 50 und N = 10000 als Default-Werte: SimuCLT = function( N = 10000 , n = 50 , Vert , plist ) { # Zufallszahlen: rVert = paste("r", Vert, sep="") plist = c(N*n, plist) # wir brauchen n*N Zufallszahlen x = do.call(rVert, plist) mu = mean(x) # wir schaetzen E[X_k] und V[X_k] anstatt sigma = sd(x) # die theoretischen Werte zu nehmen # Z_n's: X = matrix(x, ncol=n, nrow=N) Sn = rowSums(X) Zn = (Sn - n*mu)/(sqrt(n)*sigma) # graphics: par(mfrow=c(1,2)) # "plot array": 2 plots in einem Fenster, # zeilenweise angeordnet info1 = "Verteilung der X_1,...,X_n" info2 = "Verteilung von Z_n" # der erste plot: hist(x, breaks=50, xlim=c(mu-5*sigma,mu+5*sigma), prob=TRUE, main=info1 ) # der zweite plot: hist(Zn, breaks=50, xlim=c(-5,5), prob=TRUE, main=info2 ) curve(dnorm, add=TRUE, col="red") } # Probieren wir das aus: myVert = "unif" mylist = list(10,20) # gleichverteilt auf [10,20] SimuCLT( 10000 , 50 , myVert , mylist ) # ok SimuCLT( Vert=myVert , plist=mylist ) # auch ok SimuCLT( Vert="unif" , plist=list(20,30) ) # auch ok SimuCLT( "unif" , plist=list(20,30) ) # falscher Aufruf SimuCLT( Vert="exp" , plist=list(1) ) SimuCLT( Vert="binom" , plist=list(12,1/6) ) SimuCLT( n=500 , Vert="binom" , plist=list(12,1/6) ) SimuCLT( Vert="binom" , plist=list(1,0.5) ) # ok, man braucht hier offensichtlich SimuCLT( n=500 , Vert="binom" , plist=list(1,0.5) ) # groessere n's als nur 50.. SimuCLT( n=5000, Vert="binom" , plist=list(1,0.5) ) # dauert ein bischen.. # Die SimuCLT-Funktion generiert nur ein Histogramm, aber gibt # keine Werte zurueck. Wenn eine Funktion eine Zahl, einen Vektor, # eine Matrix oder alles gleichzeitig berechnen und zurueckgeben # soll, sieht die Syntax so aus: test1 = function( x ) { result = x^2 return(result) } test1(7) test2 = function( n ) { res1 = rnorm(n) # ein Vektor res2 = matrix(rnorm(n*n),nrow=n,ncol=n) # eine Matrix reslist = list(res1,res2) # Datentyp Liste # just for info: names(reslist) = c("Zufallsvektor","Zufallsmatrix") # verschiedene Datentypen kann man immer # in einer Liste zusammenfassen # und dann alles gleichzeitig zurueckgeben: return(reslist) } test2(4) res = test2(5) res[[1]] res[[2]] # oder: res$Zufallsvektor res$Zufallsmatrix