-------------
#  Blatt 9  #
-------------

# das falsche Code-Fragment

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

# kann gefixt werden mit:

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


# probieren wir das aus: wir nehmen den Code 
# aus Verteilungen-der-Schaetzer.txt:


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

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, 
# N lineare Regressionen:

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)
}


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