#---------------------------------------------------
# Wir wollen verifizieren, dass die Schaetzer fuer #
# die Regressionskoeffizienten normalverteilt sind #
# und wir wollen die Formel fuer die Varianzen der #
# Schaetzer ueberpruefen:                          #
#--------------------------------------------------#


n = 11                               # Anzahl Daten
x = seq(from=-5,to=5,length=n)
sigma = 3

# die folgenden Werte sollen weiter 
# unten geschaetzt werden:
beta0 = 3
beta1 = 2
beta2 = 1 


eps = rnorm(n,mean=0,sd=sigma)
y = beta0 + beta1*x + beta2*(x^2) + eps 
plot(x,y)


# Jetzt: N Regressions-Experimente:

N = 10000

hatbeta0 = rep(0,N)
hatbeta1 = rep(0,N)
hatbeta2 = rep(0,N) 

sigsqML = rep(0,N)                # Maximum-Likelihood-Schaetzer
sigsqET = rep(0,N)                # erwartungstreuer Schaetzer, da n=11 
                                  # ziemlich klein ist, sollte man 
                                  # deutlichen Unterschied sehen

x0 = rep(1,n)
x1 = x
x2 = x^2

X12 = cbind(x1,x2)
X   = cbind(x0,x1,x2)

for(k in 1:N)
{
  eps = rnorm(n,mean=0,sd=sigma)
  y = beta0 + beta1*x + beta2*(x^2) + eps 

  res = lm( y ~ X12 )
  betas = res$coeff
  hatbeta0[k] = betas[1]
  hatbeta1[k] = betas[2]
  hatbeta2[k] = betas[3]

  # Regression-Fit:
  yhat = res$fit
  sigsqML[k] = sum((y-yhat)^2)/n
  sigsqET[k] = sum((y-yhat)^2)/(n-3)
}


hist(hatbeta0,breaks=20)
hist(hatbeta1,breaks=20)
hist(hatbeta2,breaks=20)


# Varianzen der Schaetzer gemaess 
# Formel aus der Vorlesung:

XTX = t(X) %*% X
XTXinv = solve(XTX)
V0 = sigma^2*XTXinv[1,1]
sd0 = sqrt(V0)
sd0
V1 = sigma^2*XTXinv[2,2]
sd1 = sqrt(V1)
sd1
V2 = sigma^2*XTXinv[3,3]
sd2 = sqrt(V2)
sd2


# Gauss'sche Glockenkurven:

curve( dnorm(x,mean=beta0,sd=sd0), from=-2, to=8 )

par(mfrow=c(1,3))                        # 3 Bilder in einem Fenster

hist(hatbeta0,breaks=20,prob=TRUE)
curve( dnorm(x,mean=beta0,sd=sd0) , add=TRUE , col="red" )

hist(hatbeta1,breaks=20,prob=TRUE)
curve( dnorm(x,mean=beta1,sd=sd1) , add=TRUE , col="red" )

hist(hatbeta2,breaks=20,prob=TRUE)
curve( dnorm(x,mean=beta2,sd=sd2) , add=TRUE , col="red" )


sigma^2
mean(sigsqML)
mean(sigsqET)