-------------------------------------
#                                   #
#      Die Normalverteilung und     #
#     der zentrale Grenzwertsatz    #
#                                   #
------------------------------------- 


Definition: Eine Zufallsvariable (oder Zufallszahl) phi 
heisst normalverteilt mit Mittelwert mu und Standard-
abweichung sigma, wenn 

  Prob[ phi in [x,x+dx) ] =                             (1)

    1/sqrt(2 pi sigma^2) * exp[ - 1/2 {(x-mu)/sigma}^2 ] dx

Eine normalverteilte Zufallsvariable heisst standard- 
normalverteilt, wenn mu=0 und sigma=1 ist. 


Die grosse Bedeutung der Normalverteilung wird im wesent-
lichen durch folgenden Sachverhalt verursacht:


Satz (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 normierte Summe 

  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


Schauen wir uns jetzt die Sachen in R an:



# Start R-Session:
#
# Wir wollen zu einer gegebenen W'keitsverteilung, die in R
# implementiert ist (diese Verteilungen sind in den drei Seiten 
# W'keitsverteilungen-in-R.pdf aufgelistet) den zentralen Grenz-
# wertsatz ueberpruefen. Wir informieren uns zunaechst ueber den 
# call()-Befehl:

?call

# probieren wir mal was aus: 

Vert = "unif"             # eine Text-Variable
Vert
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..

y = call(rVert,100)
y                         # we are getting closer..

?call    # -> do.call, let's try that:

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"
y = do.call(rVert,plist)
y                         # ok, das klappt

# fuer gleichverteilte Zufallszahlen zwischen a=-2 und b=2: 
# diese Information koennen wir in die plist packen:

rm(plist)                 # loeschen
plist = list(100,-2,2)
plist

y = do.call(rVert,plist)
y                         # klappt auch..


# Zurueck zum zentralen Grenzwertsatz (CLT = central 
# limit theorem): 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 )
{
  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 
  X = matrix(x,ncol=n,nrow=N)
  Sn = rowSums(X)
  Zn = (Sn - n*mu)/(sqrt(n)*sigma)
  par(mfrow=c(1,2))               # "plot array": 2 plots in einem Fenster
  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")
}


myVert = "unif"
mylist = list(0,4)        # gleichverteilt auf [0,4]

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) )  
        
SimuCLT( Vert="binom" , plist=list(1,0.1) ) 
SimuCLT( n=1000 , Vert="binom" , plist=list(1,0.1) ) 
SimuCLT( n=10000 , Vert="binom" , plist=list(1,0.1) )