------------------------------------------
#                                        #
#     Die Verteilungen der Schaetzer     #  
#        fuer beta und sigma^2           #
#                                        #
------------------------------------------


In der Vorlesung haben wir gezeigt: 

1) hatbeta_j ist normalverteilt mit Mittelwert 
   beta_j und Varianz sigma^2 * (X^T*X)^(-1)_{j,j}

2) hats^2 ist sigma^2/(n-p-1) mal einer 
   chi^2-Verteilung

Dabei ist X = (x0,x1,...,xp) die Matrix der 
Regressoren, p+1 Stueck wenn man die Konstante 
mit dazu zaehlt, und n ist der Stichproben-
Umfang.

Wir wollen die Resultate (1) und (2) durch 
eine Simulation verifizieren:

Nehmen wir etwa an, das wir folgenden exakten 
Zusammenhang zwischen den y und den xj's haben: 

  y = beta0 + beta1*x + beta2*x^2 + eps       (3)

mit (mean=0,sd=sigma)-normalverteilten epsilons, 
wir haben dann also p+1 = 3 Regressoren. Die x's 
waehlen wir etwa aequidistant aus dem Intervall 
[-5,5], und zwar n Stueck, das ist dann der 
Stichprobenumfang. Also:


n = 11       # das ist ein ziemlich kleiner Stichprobenumfang, 
             # also man wird deutlich die Varianzen von 
             # hatbeta und hats^2 sehen koennen. 

x = seq( from=-5, to=5, length=n )
x

# wir probieren mal folgende betas:
beta0 = -6
beta1 = -1 
beta2 = 0.5

sigma = 2
eps = rnorm(n,mean=0,sd=sigma)

y = beta0 + beta1*x + beta2*x^2 + eps
plot(x,y)    


# Wir machen jetzt N Zufallsexperimente, ein 
# Zufallsexperiment besteht aus dem Generieren 
# der y's gemaess (3) und dem Durchfuehren einer
# linearen Regression. Die geschaetzen Werte 
# fuer die betas und das sigma^2 schreiben wir 
# jeweils in Vektoren der Laenge N:

N = 10000
hatbeta0 = rep(0,N)
hatbeta1 = rep(0,N)
hatbeta2 = rep(0,N)
hats     = rep(0,N)        # soll das hat(s^2) sein

X = cbind(x,x^2)           # Matrix der Regressoren, ohne die 
                           # Konstante, die wird automatisch 
                           # lm()-Befehl mitgenommen

for(k in 1:N)
{
  eps = rnorm(n,mean=0,sd=sigma)
  y = beta0 + beta1*x + beta2*x^2 + eps
  # die lineare Regression:
  res = lm(y ~ X)
  betas = res$coeff
  hatbeta0[k] = betas[1]
  hatbeta1[k] = betas[2]
  hatbeta2[k] = betas[3]
  hats[k] = 1/(n-3) * sum(res$residuals^2)
}


# bevor wir uns die Verteilungen anschauen, 
# berechnen wir noch die theoretischen Varianzen 
# der hatbeta:

X0 = cbind(rep(1,n),X)     # Regressoren mit Konstante
X0
covmatrix = sigma^2* solve(t(X0)%*%X0)
covmatrix
var0 = covmatrix[1,1]
var1 = covmatrix[2,2]
var2 = covmatrix[3,3]
sd0 = sqrt(var0)
sd1 = sqrt(var1)
sd2 = sqrt(var2)
sd0
sd1
sd2

# wir machen 4 Histogramme, 3 fuer die betas 
# und eins fuer das hat(s^2), packen wir die 
# etwa in ein 2x2 plot-array:

par(mfrow=c(2,2))          # sets up plot-array, 
                           # "make frame by row"

info1 = "Verteilung hat(beta0)"
hist(hatbeta0,breaks=50,prob=TRUE,main=info1)
curve( dnorm(x,mean=beta0,sd=sd0), col="red", add=TRUE)

info2 = "Verteilung hat(beta1)"
hist(hatbeta1,breaks=50,prob=TRUE,main=info2)
curve( dnorm(x,mean=beta1,sd=sd1), col="red", add=TRUE)

info3 = "Verteilung hat(beta2)"
hist(hatbeta2,breaks=50,prob=TRUE,main=info3)
curve( dnorm(x,mean=beta2,sd=sd2), col="red", add=TRUE)

info4 = "Verteilung   (n-3) * hat(s^2) / sigma^2"
hist((n-3)*hats/sigma^2,breaks=50,prob=TRUE,main=info4)
curve( dchisq(x,df=n-3), col="red", add=TRUE)