# Tests

x<-rnorm(10,3,1); x; 
mean(x); sd(x)

# Gauss-Test auf mu>=4 gegen mu<4, bekannte Varianz 1
n<-length(x)
(T<-sqrt(n)*(mean(x)-4)/1)	

(c<-qnorm(0.05))  		# -1.644854

if (T<=c)  print("Ablehnung der Nullhypothese, da T<= krit. Wert.") else print("keine Ablehnung, da T>krit. Wert.")

# -------------------------------------------------

# t-test bei unbekannter Varianz
n<-length(x)

(T<-sqrt(n)*(mean(x)-4)/sd(x))	# unbekannte Varinz 2.65535, T=-1.372221

(c<-qt(0.05,df=n-1))		# -2.131847

if (T<=c)  print("Ablehnung, da T<= krit. Wert.") else print("keine Ablehnung, da T>krit. Wert.")

# --------------------------------------------------

# Nochmal t-test mit t.test
s<-t.test(x,mu=4,alternative="less")
s
names(s)

alpha<-0.05
if (s$p.value<alpha)  print("Ablehnung, da p-Wert<=alpha") else print("keine Ablehnung, da p-Wert>alpha")

# --------------------------------------------------

x<-rnorm(10,0,2)
y<-rnorm(6,1,2)
n<-length(x)
m<-length(y)

# t-Test mu1<=mu2  gegen  mu1>mu2
mean(x); mean(y)

# Schaetzer fuer Varianz
v<- ( sum((x-mean(x))^2) + sum((y-mean(y))^2) ) / (n+m-2); v

# Teststatistik
T<-(mean(x)-mean(y)) / (sqrt(v *(1/n+1/m))); T

# Kritischer Wert = Quantil t_18 - Verteilung
c<-qt(0.975,n+m-2); c

# Testentscheidung
if (T>=c)  print("Ablehnung der Nullhypothese, da T>= krit. Wert.") else print("keine Ablehnung, da T<krit. Wert.")

# -----------------------------------------------------

# Auswertung mit t.test
s<-t.test(x,y,alternative="greater",var.equal = T); s

# Testentscheidung
alpha<-0.025
if (s$p.value<alpha)  print("Ablehnung, da p-Wert<=alpha") else print("keine Ablehnung, da p-Wert>alpha")

# -----------------------------------------------------------

# Rank-Test

x<-rnorm(10,0,2)
y<-rnorm(6,1,1)
n<-length(x)
m<-length(y)

x; rank(x)
data.frame(x,rank(x));
data.frame(y,rank(y))

# Rang von y in der gemeinsamen Stichprobe
ry<-rank(c(y,x))[1:m]; ry
sum(ry)

sum(rank(c(x,y))[1:n])

# Wilcoxon-Test X stoch <= Y gegen X stoch > Y
# (entspricht F_X>=F_Y  gegen F_X<F_Y)
wilcox.test(x,y,alternative="g")

# Wert der Teststatistik W: Rangsumme von x-n(n+1)/2
sum(rank(c(x,y))[1:n]) - n*(n+1)/2


# ----------------------------------------------------------------
# ----------------------------------------------------------------


# Regression mit R

# Daten erzeugen und Plotten
n<-100; epsi<-rnorm(n)/2; x<-sort(rnorm(n)) 
y<-2+0.5*x+epsi
plot(x,y)

# Modell
mo<-lm(y ~ x)
mo
abline(mo)

# Residuen plotten
plot(x,mo$residuals)
abline(h=0)

# weitere Angaben
summary(mo)
(cor(x,y))^2

# Nun zu Fuss
# Design-Matrix mit 
model.matrix(mo)
# oder
X<-cbind(rep(1,n),x)

# KQ-Schaetzer (X^T*X)^{-1}*X^T*y 
beta<-solve(t(X) %*% X) %*% t(X) %*% y
beta

# RSS/(n-p)
v<-sum(mo$residuals^2)/mo$df.residual

# Berechnung der t-Statistiken zum Test beta_i=0 gegen <>0
beta/sqrt(diag( solve(t(X) %*% X) ))/sqrt(v)

--------------------------------------------------------------

# Polynomiale Regression mit R

# Daten erzeugen und Plotten
n<-100; epsi<-rnorm(n)/2; x<-sort(rnorm(n)); z<-x^2
y<-2+0.5*x-z+epsi

# Lineare Regression passt hier nicht
# par(mfrow=c(2,1))
plot(x,y)
mo<-lm(y~x)
abline(mo,col="red")

plot(x,mo$residuals)
abline(h=0)

# Quadratische Regression  
plot(x,y)
mo<-lm(y ~ x+z)
# oder auch
mo<-lm(formula = y ~ x + I(x^2))
mo
lines(x,fitted(mo),col="red")


# Residuen plotten
plot(mo$residuals)
abline(h=0)

# weitere Angaben
summary(mo)
(cor(x,y))^2
(cor(mo$fitted.values,y))^2  # Bestimmtheitsmaß

# Design-Matrix
X<-model.matrix(mo)
# oder
X<-cbind(rep(1,n),x,z)

# KQ-Schaetzer (X^T*X)^{-1}*X^T*y 
beta<-solve(t(X) %*% X) %*% t(X) %*% y
beta

------------------------------------------