#===========================================================================
# Zusatzmaterial zu Übungsblatt 3
# Statistik 4 für Nebenfachstudierende, SoSe 2017
#
# Datum: 11.5.2017
# Autor: Nora Fenske, Micha Schneider
#=============================================================================

load("decathlon.Rdata")

#=============================================================================
# Aufgabe 7: Tests im Einstichprobenfall
#=============================================================================

#---------------
# (a)  Kovarianzmatrix wird als bekannt vorausgesetzt
#---------------

# Reduziere Datensatz auf Stichprobe aus dem Jahr 2008
X2008<- decathlonO[decathlonO$year==2008,c("M100", "M400")]

head(X2008)

# Berechne Mittelwertsvektor 
mean2008 <- round(apply(X2008, MARGIN=2, mean), digits=2)
mean2008

# Definiere mu0 und Sigma0
mu0 <- c(11, 50)
Sigma0 <- matrix(c(2, 0, 0, 4), nrow=2)
Sigma0

n <- nrow(X2008)

# Berechne die Teststatistik
Tstat <- n*t(mean2008-mu0)%*%solve(Sigma0)%*%(mean2008-mu0)
Tstat

# Berechne 95%-Quantil der Chi-Quadrat-Verteilung mit 2 df
K <- qchisq(0.95, df=2)
K

# Vergleiche Tstat mit 95%-Quantil
Tstat > K

# Graphik:
plot(seq(0, 15, by = .1), dchisq(seq(0, 15, by = .1),2), type="l", main=expression(chi^2-Dichte))
abline(v=qchisq(0.95, df=2), lty = "dashed")

# Graphik(empirisch):
dat <- rchisq(10000,2)
hist(dat, prob=T, breaks=100)
abline(v=qchisq(0.95, df=2), lty = "dashed")

#---------------
# (b) Kovarianzmatrix wird als unbekannt vorausgesetzt
#---------------

# Berechne Kovarianzmatrix wie in der Angabe
S2008 <- round(cov(X2008), digits=2)
S2008

# Berechne Teststatistik
Tstatemp <- (n*(n-2))/(2*(n-1))*t(mean2008-mu0)%*%solve(S2008)%*%(mean2008-mu0)
Tstatemp

# Berechne 95%-Quantil der F-Verteilung mit 2 und 178 Freiheitsgraden
Kemp <- qf(0.95, df1=2, df2=278)
Kemp

# Vergleiche Teststatistik mit 95%-Quantil
Tstatemp  > Kemp

#=============================================================================
# Aufgabe 8: Test im Zweistichprobenfall - unabhängige Stichproben
#=============================================================================

# Reduziere Datensatz auf Stichprobe aus dem Jahr 1998
X1998<- decathlonO[decathlonO$year==1998,c("M100", "M400")]


# Berechne Mittelwertsvektor und Kovarianzmatrix wie in Angabe
mean1998 <- round(apply(X1998, MARGIN=2, mean), digits=2)
mean1998

S1998 <- round(cov(X1998), digits=2)
S1998

n1998 <- nrow(X1998)
n2008 <- nrow(X2008)


# Berechne gepoolte Kovarianzmatrix
Stot <- 1/(n1998+n2008-2)*((n1998-1)*S1998 + (n2008-1)*S2008)
Stot <- round(Stot, digits=2)
Stot 

# Berechne Teststatistik

# -> ungerundet
Tstat <- n1998*n2008/(n1998 + n2008) * t(mean1998 - mean2008)%*%solve(Stot)%*%(mean1998 - mean2008)
Tstat

# -> gerundet
StotInv <- round(solve(round(Stot, digits=2)), digits=2)

TstatR <- n1998*n2008/(n1998 + n2008) * t(mean1998 - mean2008)%*%StotInv%*%(mean1998 - mean2008)
TstatR


# Vorfaktor für Teststatistik 
(n1998 + n2008 - 2 - 1)/((n1998+n2008-2)*2)
T2 <- (n1998 + n2008 - 2 - 1)/((n1998+n2008-2)*2)*Tstat
T2

# Berechne Quantil der F-Verteilung
K2 <- qf(0.95, df1=2, df2=n1998+n2008-2-1)
K2

# Vergleiche Teststatistik mit 95%-Quantil
T2  > K2



#=============================================================================
# Aufgabe 9: Test im Zweistichprobenfall - verbundene Stichproben
#=============================================================================

# Extrahiere Daten von 4 Athleten
ind1 <- rownames(decathlon[decathlon$V15=="Bryan Ezra Clay" & decathlon$year%in%c(2007,2008),])
ind2 <- rownames(decathlon[decathlon$V15=="Dmitri Karpov" & decathlon$year%in%c(2007,2008),])
ind3 <- rownames(decathlon[decathlon$V15=="Andrei Krauchanka" & decathlon$year%in%c(2007,2008),])
ind4 <- rownames(decathlon[decathlon$V15=="Cory Roberts" & decathlon$year%in%c(2007,2008),])

a7dat <- decathlonO[rownames(decathlonO) %in% c(ind1, ind2, ind3, ind4),]
a7dat <- a7dat[c(1, 6, 2, 7, 3, 5, 4, 8),c("V16", "Disc", "Jav", "year")] 
a7dat

# Formatiere Daten vom long-Format ins wide-Format um
a7wide <- reshape(a7dat, v.names=c("Disc", "Jav"),
                  timevar="year", idvar="V16", direction="wide")
                  

# Berechne Datenmatrix mit Differenzen
Dmat <- cbind(a7wide$Disc.2008-a7wide$Disc.2007, a7wide$Jav.2008-a7wide$Jav.2007)

# Berechne Mittelwertsvektor und Kovarianzmatrix
dmean <- apply(Dmat, MARGIN=2, FUN=mean)
dmean

dmean <- round(dmean, digits=2)
dmean

Sd <- cov(Dmat)
Sd

SdInv <- round(solve(Sd), digits=6)
SdInv

# Berechne Teststatistik
Tstat <- (4-2)*4/((4-1)*2)*t(dmean)%*%SdInv%*%dmean
Tstat

# Berechne Quantil der F-Verteilung
K <- qf(0.95, df1=2, df2=2)
K

# Vergleiche Teststatistik mit 95%-Quantil
Tstat  > K