#------------------------------# # # # week6a: t-, F- und # # chi^2-Verteilung # # # #----------------------------- # # Wir wollen zunaechst die analytischen Ausdruecke # fuer die Dichten aus der Vorlesung checken, dazu # definieren wir eigene Funktionen: Dichte_t = function( y , n ) { res = 1/(1+y^2/n)^((n+1)/2) cn = gamma((n+1)/2)/(gamma(n/2)*sqrt(pi*n)) return(cn*res) } Dichte_F = function( y , k , ell ) { res = y^(k/2-1)/(1+k/ell*y)^((k+ell)/2) ckell = (k/ell)^(k/2) * gamma((k+ell)/2)/( gamma(k/2)*gamma(ell/2) ) return( ckell*res ) } Dichte_chi = function( y , n ) { res = y^(n/2-1) * exp(-y/2) cn = 1/( 2^(n/2) * gamma(n/2) ) return(cn*res) } # Wir schauen uns die Sachen an und vergleichen mit den # in R vorimplementierten Dichten: # t-Verteilung: x = seq(from=-10,to=10,by=0.01) for(n in 1:10) { info = paste("n =",n) plot(x,Dichte_t(x,n),main=info) points(x,dt(x,n),cex = 0.5,col="red") # cex = character expansion, cex = 1 normale Groesse, # cex = 0.5 halbe Groesse, cex = 2 doppelte Groesse usw. readline("press enter to continue..") } # F-Verteilung: x = seq(from=0,to=4,by=0.01) for(k in 1:10) { for(ell in 1:10) { info = paste("(k,l) = (",k,",",ell,")", sep="" ) plot(x,Dichte_F(x,k,ell),main=info) points(x,df(x,k,ell),cex = 0.5,col="red") readline("press enter to continue..") } } # chi^2 - Verteilung: x = seq(from=0,to=15,by=0.01) for(n in 1:10) { info = paste("n =",n) plot(x,Dichte_chi(x,n),main=info) points(x,dchisq(x,n),cex = 0.5,col="red") readline("press enter to continue..") } #------------------------------------# # Simulation der Verteilungen aus # # normalverteilten Zufallszahlen # #------------------------------------# # t-Verteilung: N = 10000 # Anzahl der simulierten y's # fuer das Histogramm for(n in 1:10) # das n fuer die t_n-Verteilung { phi0 = rnorm(N) phi = rnorm(n*N) Phi = matrix(phi,N,n) PhiQuad = Phi^2 sn = rowSums(PhiQuad) y = phi0/sqrt(sn/n) info = paste("n =",n) # fuer n=1,2,3 sieht das sehr komisch aus, # aber scheint schon zu passen -> fixen wir # weiter unten hist(y,prob=TRUE,breaks=50,main=info) curve(dt(x,n),add=TRUE,col="red") readline("press enter to continue..") } # F-Verteilung: N = 10000 # Anzahl der simulierten y's # fuer das Histogramm for(k in 1:10) { for(ell in 2:10) # ell=1 kann sehr, sehr grosse y's generieren, { # error Nbreaks too large weiter unten.. # Z=Zaehler, N=Nenner phiZ = rnorm(k*N) PhiZ = matrix(phiZ,N,k) PhiZQuad = PhiZ^2 sZ = rowSums(PhiZQuad) phiN = rnorm(ell*N) PhiN = matrix(phiN,N,ell) PhiNQuad = PhiN^2 sN = rowSums(PhiNQuad) y = sZ/sN * ell/k info = paste("(k,l) = (",k,",",ell,")", sep="" ) Nbreaks = (max(y)-min(y))/0.1 hist(y, ,xlim=c(0,5), prob=TRUE, breaks=Nbreaks, main=info ) curve(df(x,k,ell),add=TRUE,col="red") readline("press enter to continue..") } } # chi^2 - Verteilung: N = 10000 # Anzahl der simulierten y's # fuer das Histogramm for(n in 1:10) # das n fuer die t_n-Verteilung { phi = rnorm(n*N) Phi = matrix(phi,N,n) PhiQuad = Phi^2 y = rowSums(PhiQuad) info = paste("n =",n) hist(y,prob=TRUE,breaks=50,main=info) curve(dchisq(x,n),add=TRUE,col="red") readline("press enter to continue..") } # t-Verteilung, jetzt auch hier mit schoenen Histogrammen # fuer kleine Werte von n: N = 10000 for(n in 1:10) { phi0 = rnorm(N) phi = rnorm(n*N) Phi = matrix(phi,N,n) PhiQuad = Phi^2 sn = rowSums(PhiQuad) y = phi0/sqrt(sn/n) info = paste("n =",n) Nbreaks = (max(y)-min(y))/0.2 hist(y, xlim=c(-10,10), prob=TRUE, breaks=Nbreaks, main=info) curve(dt(x,n),add=TRUE,col="red") readline("press enter to continue..") }