------------- # Aufgabe 1 # ------------- # 1a) mycol = rainbow(10) curve( dchisq(x,df=1), xlim=c(0,12), ylim=c(0,0.6), col=mycol[1], xlab="xi", ylab="densities", main="Dichte der chi^2-Verteilung mit n=1,2,...,10 Freiheitsgraden" ) for(i in 2:10) { curve(dchisq(x,df=i),add=TRUE,col=mycol[i]) } # 1b-d) n = 5 N = 10000 x = rnorm(n*N) Phi = matrix(x,nrow=N,ncol=n,byrow=TRUE) PhiSquared = Phi*Phi xi = rowSums(PhiSquared) hist(xi,prob=TRUE,breaks=20) curve(dchisq(x,df=n),add=TRUE,col="red") # 1e) # let's check for all n=1,2,...,10: N=10000 for(n in 1:10) { x = rnorm(n*N) Phi = matrix(x,nrow=N,ncol=n,byrow=TRUE) PhiSquared = Phi*Phi xi = rowSums(PhiSquared) info = paste("n =",n) hist(xi,prob=TRUE,breaks=50,xlim=c(0,30),ylim=c(0,0.5),main=info) curve(dchisq(x,df=n),add=TRUE,col="red") # wait a bit until next plot is shown: readline(prompt = "press Enter to continue..") } ------------- # Aufgabe 2 # ------------- # 2a) N = 10000 n = 5 phi0 = rnorm(N) phi = matrix(rnorm(n*N),ncol=n,nrow=N) y = phi0 / sqrt(rowSums(phi*phi)/n) hist(y,prob=TRUE,breaks=50) curve(dt(x,df=n),col="red",add=TRUE) # 2b) N = 10000 k = 40 l = 80 phi1 = matrix(rnorm(k*N),ncol=k,nrow=N) phi2 = matrix(rnorm(l*N),ncol=l,nrow=N) y = (rowSums(phi1*phi1)/k) / (rowSums(phi2*phi2)/l) hist(y,prob=TRUE,breaks=50) curve(df(x,df1=k,df2=l),col=2,add=TRUE) # Ergaenzung zu 2a: # wenn wir uns das wieder fuer die ersten 10 n's # anschauen wollen: N = 10000 for(n in 1:10) { phi0 = rnorm(N) phi = matrix(rnorm(n*N),ncol=n,nrow=N) y = phi0 / sqrt(rowSums(phi*phi)/n) info = paste("n =",n) hist(y,prob=TRUE,breaks=50,xlim=c(-15,15),ylim=c(0,0.5),main=info) curve(dt(x,df=n),add=TRUE,col="red") # wait a bit until next plot is shown: readline(prompt = "press Enter to continue..") } # die ersten Histogramme sehen etwas komisch aus: # das liegt daran, dass sehr grosse y's auftreten # koennen und das gesamte Intervall [min(y),max(y)] # wird in breaks=50 Rechtecke eingeteilt; man sieht # dann nur ein oder zwei Rechtecke, weil in den an- # deren kaum Zufallszahlen liegen: N = 10000 n = 1 phi0 = rnorm(N) phi = matrix(rnorm(n*N),ncol=n,nrow=N) y = phi0 / sqrt(rowSums(phi*phi)/n) summary(y) # etwa -2000 oder -10000 bis +1000 oder +12000 hist(y) ymin=min(y) ymax=max(y) Delta = 0.25 # ymax - ymin = Nbreaks * Delta wollen wir haben, # also: Nbreaks = (ymax-ymin)/Delta hist(y,breaks=Nbreaks) # jetzt lassen wir die grossen y's weg: hist(y,breaks=Nbreaks,xlim=c(-20,20)) # ok, das passt jetzt.. # Damit: Delta = 0.25 for(n in 1:10) { phi0 = rnorm(N) phi = matrix(rnorm(n*N),ncol=n,nrow=N) y = phi0 / sqrt(rowSums(phi*phi)/n) Nbreaks = (max(y)-min(y))/Delta info = paste("n =",n) hist(y,prob=TRUE,breaks=Nbreaks,xlim=c(-15,15),ylim=c(0,0.5),main=info) curve(dt(x,df=n),add=TRUE,col="red") # wait a bit until next plot is shown: readline(prompt = "press Enter to continue..") }