--------------
#   Blatt5   #
--------------


# 1a)

n = 3
N = 10000
mu = 15
sig = 2
x = rnorm(n*N,mean=mu,sd=sig)
Phi = matrix(x,nrow = N, ncol = n)


# 1b)

MeanPhi = rowSums(Phi)/n


# 1c)

rm(x)        # just to be sure, oben gab es schon ein x
x = c(1,3)   # apparently mean = 2
x
sd(x)
var(x)


# 1d)

SdPhi = apply(Phi,1,sd)     # 1 bedeutet: apply to rows, 
                            # 2 bedeutet: apply to columns

# 1e)

Y = (MeanPhi - mu)/(SdPhi/sqrt(n))


# 1f)

hist(Y,prob=TRUE,breaks=50,xlim=c(-15,15),ylim=c(0,0.5))

# trotz breaks=50 hat das Histogramm sehr wenige Rechtecke. 
# Das liegt daran, dass mit breaks der gesamte Bereich der 
# erzeugten Zufallszahlen ueberdeckt werden soll. Wenn also 
# die meisten Zahlen etwa zwischen -10 und 10 liegen, aber 
# es gibt einige wenige Ausreisser etwa bei -238 und bei 
# +425, schauen wir eben:

summary(Y)

# dann wird das gesamte Intervall [-238,+425] in breaks=50 
# gleichbreite Rechecke zerlegt, der interessante Bereich 
# [-10,10] hat dann eben nur eine sehr grobe Aufloesung. 
# Das koennen wir etwa folgendermassen fixen:

Deltax = 0.25          # soll die Breite eines Rechtecks 
                       # werden

ymax = max(Y)
ymin = min(Y)

# wir brauchen Nbins Rechtecke mit Nbins*Deltax = ymax-ymin
# also:

Nbins = (ymax-ymin)/Deltax
Nbins = round(Nbins)   # sollte eine ganze Zahl sein

# und jetzt probieren wir einfach:

hist(Y, prob=TRUE, breaks=Nbins, xlim=c(-10,10), ylim=c(0,0.5) )

# das sieht gut aus. Jetzt noch die theoretische Dichte:

curve( dt(x,df=(n-1)), add=TRUE, col="red" )


# 1g) 

N = 10000
mu = 15
sig = 2

Deltax = 0.25          # Breite der Rechtecke im Histogramm

for(n in 2:10)
{
  x = rnorm(n*N,mean=mu,sd=sig)
  Phi = matrix(x,nrow = N, ncol = n)
  MeanPhi = rowSums(Phi)/n
  SdPhi = apply(Phi,1,sd)  
  Y = (MeanPhi - mu)/(SdPhi/sqrt(n))
  ymin = min(Y)
  ymax = max(Y)
  print(c(ymin,ymax))                         # for information only
  Nbins = round((ymax-ymin)/Deltax)

  # Ueberschrift fuer das Histogramm:
  info = paste("Histogramm von Y_n und Dichte t_{n-1},   n = ",n,sep="")
  hist(Y, prob=TRUE, breaks=Nbins, xlim=c(-10,10), ylim=c(0,0.5), main=info)
  curve(dt(x,df=(n-1)),add=TRUE,col="red")

  # naechstes Bild kommt erst, wenn wir Enter druecken:
  readline(prompt = "press Enter to continue..")
}