remove(list=ls()) setwd("/Users/michael/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt05/R") #setwd("/Volumes/MichaelsHD/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt05") set.seed(23514) #1(a) #Daten einlesen data <- matrix(scan("./fieber.dat"),ncol=20,byrow=T) #t erzeugen t <- 1:100 #20 mal kleinen Datensatz anpassen #in T100 werden die Testgroessen gespeichert T100 <- numeric(20) #in p100 werden die p-Werte, die lm() berechnet gespeichert p100 <- numeric(20) for (i in 1:20) { fitH1 <- lm(data[,i] ~ 1+t) fitH0 <- lm(data[,i] ~ 1) # print(i) # print(summary(fitH1)) p100[i] <- summary(fitH1)$coeff[2,4] QR <- (summary(fitH1)$sigma)^2 * fitH1$df QRH0 <- (summary(fitH0)$sigma)^2 * fitH0$df T100[i] <- (QRH0 - QR)/QR*(100-2) } #95%-Quantil der F_1,98-Verteilung entsch100 <- qf(0.95,1,98) entsch100 #Wie oft wurde nicht abgelehnt sum(as.numeric(T100>entsch100)) #t100[4] wird nicht abgelehnt T100 T100[4] #p-Werte sum(as.numeric(p100<0.05)) p100[4] #summary(lm()) kommt zum selben Ergebnis #p-Wert von T100[4] (nicht gefordert) 1-pf(T100[4],1,98) #1(b) #selbst Daten erzeugen set.seed(23512) epsilons <- rnorm(2000,mean=36.8,sd=0.7) mydat <- matrix(epsilons,ncol=20) myT <- numeric(20) myp <- numeric(20) for (i in 1:20) { fitH1 <- lm(mydat[,i] ~ 1+t) fitH0 <- lm(mydat[,i] ~ 1) # print(i) # print(summary(fitH1)) QR <- (summary(fitH1)$sigma)^2 * fitH1$df QRH0 <- (summary(fitH0)$sigma)^2 * fitH0$df myT[i] <- (QRH0 - QR)/QR*(100-2) myp[i] <- summary(fitH1)$coeff[2,4] } #Wie oft wurde nicht abgelehnt sum(as.numeric(myT>entsch100)) #p-Werte sum(as.numeric(myp<0.05)) #sogar dreimal H_0 nicht abgelehnt, obwohl unter H_0 simuliert