#==========================================================
# Einleitung in die Regressionsanalyse mittels Simulationen
#==========================================================

# Betrachtet wird ein SLR fuer n=11 feste x 0,1,...,10
x <- seq(0, 10)    # oder in Kurzform x <- 0:10
n <- length(x)     # hier n=11 mit Erwartungswerten
beta0 <- 3
beta1 <- 0.7
mu <- beta0 + beta1*x
sigma <- rep(2, n) # konstante Std-Abweichungen

# Zeichne das Mean-Modell mit den 2-sigma Grenzen
plot(x ,mu, type="l", ylim=c(-5,20),col="red",pch="*",ylab="mu+-2sigma")
lines(x, mu+2*sigma, col="red")
lines(x, mu-2*sigma, col="red")

# Simuliere einen Datensatz (mit n Elementen) unter diesen Annahmen
# Wir generieren dazu normalverteilte Responses
y <- rnorm(n, mean=mu, sd=sigma)
points(x, y)       # zeichne auch die Daten ein

# und berechnen die LS Schaetzer b_0, b_1 sowie MSE
lm(y ~ 1 + x)      # Kurzform: lm(y ~ x)
slr <- lm(y ~ x)   # speichere das Fit-Objekt unter "slr"
print(slr)         # und drucke "slr"

# mehr Information (Standardfehler, ...) mittels "summary(...)"
summary(slr)

# was so alles in "slr" und in "summary(slr)" ist, zeigt
attributes(slr); attributes(summary(slr))

abline(slr)        # zeichne auch die geschaetzte Regressionsgerade ein

# ein Freund erhebt auch Daten unter diesem Modell und erhaelt:
y2 <- rnorm(n, mean=mu, sd=sigma)
points(x, y2, col="blue")
abline(lm(y2 ~ x), col="blue") 

# die Stichprobe definiert die Parameterschaetzer b_0, b_1 !!
# nach oftmaligem Wiederholen dieses Experiments (S mal) erhaelt man
plot(x ,mu, type="l", ylim=c(-5,20),col="red",pch="*",ylab="mu+-2sigma")
lines(x, mu+2*sigma, col="red")
lines(x, mu-2*sigma, col="red")

S <- 10000
b0 <- b1 <- MSE <- rep(NA, S)
for(s in 1:S) 
  {
  y.sim <- rnorm(n, mean=mu, sd=sigma)
  slr.sim <- lm(y.sim ~ x)
  b0[s] <- slr.sim$coefficients[1] # oder slr.sim$coef["(Intercept)"]
  b1[s] <- slr.sim$coefficients[2] # oder slr.sim$coef["x"]
  MSE[s] <- summary(slr.sim)$sigma^2
  abline(slr.sim) 
  }

# Die Schaetzer b0, b1 sehen aus als waeren sie beide normalverteilt
hist(b0, freq=FALSE); summary(b0); sd(b0)
hist(b1, freq=FALSE); summary(b1); sd(b1)

# Der Schaetzer MSE für sieht aus als waere er chi^2-verteilt
hist(MSE, freq=FALSE); summary(MSE); sd(MSE)

# unsere Erst-Stichprobe ergab:
summary(slr)

# die theoretischen Egebnisse (uns immer unbekannt) waeren
Sxx <- sum((x-mean(x))^2); Sxx
var.b0 <- sigma[1]^2*sum(x^2)/(n*Sxx)  # Varianz von b_0
sqrt(var.b0)              # Standard Error von b_0
var.b1 <- sigma[1]^2/Sxx  # Varianz von b_1
sqrt(var.b1)              # Standard Error von b_1

# SSE/sigma^2 ~ chi^2_(n-2), wobei SSE = MSE*(n-2) 
# mit  Var(SSE/sigma^2) = 2(n-2)
# also Var(SSE) = 2(n-2)sigma^4
# oder Var(MSE) = Var(SSE/(n-2)) = 2sigma^4/(n-2)
var.MSE <- 2*sigma[1]^4/(n-2)  # Varianz von MSE
sqrt(var.MSE)                  # Standard Error von MSE


#=========================================================
# Konfidenzintervalle fuer unbekannte Parameter:
# Eine Einfuehrung in das Konzept mittels Simulation
#=========================================================
# (1-alpha)*100% Prozent KIV(beta_1)
# ===> b1 +- t(1-alpha/2; n-2)*sqrt(MSE/Sxx)
alpha <- 0.05
t <- qt(1-alpha/2, n-2) # t=2.262 (z=1.960)

u <- o <- rep(NA, S)
for (s in 1:S){
  u[s] <- b1[s] - t*sqrt(MSE[s]/Sxx)
  o[s] <- b1[s] + t*sqrt(MSE[s]/Sxx)
}
wrong.o <- sum(o < beta1); wrong.o
wrong.u <- sum(beta1 < u); wrong.u
alpha.mc <- (wrong.o + wrong.u)/S; alpha.mc

R <- 200
plot(seq(1:R), b1[1:R],
     ylim=c(min(u[1:R]), max(o[1:R])),xlab="repetition",ylab="KIV(beta1)")
for (s in 1:R)
  {
  col <- 1*(u[s]<beta1 & o[s]>beta1) + 2*(o[s]<beta1) + 3*(u[s]>beta1)
  l <- 1*(col == 1) + 4*(col > 1) # line width
  lines(c(s, s), c(u[s], o[s]), col=col, lwd=l)
}
abline(h=beta1, lwd=2)


#=========================================================
# Hypothesentests fuer unbekannte Parameter:
# Eine Einfuehrung in das Konzept mittels Simulation
#=========================================================
# Teste H0: beta1 = 0 gegen Ha: beta1 > 0 mit t* = b1/(sqrt(MSE/Sxx))

tstar <- rep(NA, S)
for (s in 1:S){
  tstar[s] <- b1[s]/sqrt(MSE[s]/Sxx)
}

R <- S
plot(seq(1:R), tstar[1:R], xlab="repetition", ylab="t*")

t <- qt(1-alpha, n-2) # t=1.833
abline(h=t, col="red", lwd=2)

# Vergleiche die hypothetische Verteilung von t* unter H0 mit der
# durch die simulierten t* Werte geschaetzten:
plot((-30:100)/10,dt((-30:100)/10,n-2),type="l",lwd=2,xlab="t*",ylab="Dichte")
lines(density(tstar), col="blue", lwd=2)
abline(v=t, col="red", lwd=2)

# die richtige Enscheidung trifft man mit Wahrscheinlichkeit (~95%)
sum(tstar < t)/S

#=========================================================
# Jetzt testen wir H0: beta1 = 0.7 gegen Ha: beta1 \= 0.7
tstar1 <- rep(NA, S)
for (s in 1:S){
  tstar1[s] <- (b1[s] - 0.7)/sqrt(MSE[s]/Sxx)
}

plot(seq(1:R), tstar1[1:R], xlab="repetition", ylab="t*")
t <- qt(1-alpha/2, n-2)
abline(h=t, col="red", lwd=2); abline(h=-t, col="red", lwd=2)

# Vergleiche die hypothetische Verteilung von t* unter H0 mit der
# durch die simulierten t* Werte geschaetzten:
plot((-40:40)/10,dt((-40:40)/10,n-2),type="l",lwd=2,xlab="t*",ylab="Dichte")
lines(density(tstar1), col="blue", lwd=2)
abline(v=t, col="red", lwd=2); abline(v=-t, col="red", lwd=2)

# die richtige Enscheidung trifft man mit Wahrscheinlichkeit (~95%)
sum(abs(tstar1) < t)/S

#=========================================================
#=========================================================
# vergleiche Normal mit t Verteilungen
plot((-30:30)/10, dnorm((-30:30)/10))
lines((-30:30)/10, dt((-30:30)/10, 2), lty=2, lwd=2)
lines((-30:30)/10, dt((-30:30)/10, 10), lty=3, lwd=2)
lines((-30:30)/10, dt((-30:30)/10, 50), lty=4, lwd=2)

# vergleiche chi^2 Verteilungen
plot((0:500)/10, dchisq((0:500)/10, 5))
lines((0:500)/10, dchisq((0:500)/10, 10), lty=2, lwd=2)
lines((0:500)/10, dchisq((0:500)/10, 20), lty=2, lwd=2)
lines((0:500)/10, dchisq((0:500)/10, 40), lty=2, lwd=2)