------------------------------------------ # # # Vertauensintervalle fuer die # # Regressionskoeffizienten # # # ------------------------------------------ In der Vorlesung haben wir gezeigt: Die Test- groesse (hatbeta_j -beta_j) / sqrt( hatvar ) mit hatvar = hats^2 (X^T*X)^(-1)_{j,j} und hats^2 = 1/[n-(p+1)] * (y-X*hatbeta)^2 ist t_{n-(p+1)} - verteilt. Das koennen wir jetzt benutzen, um Vertrauensintervalle fuer die Schaetzungen der Regressionskoeffizienten zu bekommen. Wir muessen uns zunaechst daran erinnern, was Quantile sind: Da wir weiter unten wieder das Regressions- problem aus Verteilungen-der-Schaetzer.txt betrachten werden, betrachten wir als Beispiel t_{n-(p+1)} = t_{11-(2+1)} = t_8 - verteilte Zufallszahlen. Wir erzeugen N = 10000 Stueck davon: N = 10000 x = rt(N,df=8) plot(x) hist(x) hist(x,breaks=100,prob=TRUE) curve(dt(x,df=8),col="red",add=TRUE) summary(x) # In der Ausgabe von summary(x) gibt es die # Zahlen 1st Qu, das steht fuer "erstes Quartil", # Median und 3rd Qu fuer drittes Quartil. Den # Median koennte man auch als zweites Quartil # bezeichnen. Diese Zahlen haben die folgende # Bedeutung: wir ordnen die Zufallszahlen der # Groesse nach und schreiben das Resultat in den # Vektor sortx sortx = sort(x) head(sortx,50) # die ersten 50 Vektorelemente tail(sortx,50) # die letzten 50 Vektorelemente summary(x) summary(sortx) # ok, dieselben Zahlen wie in x, # nur eben geordnet plot(sortx) # das sieht nach einer deterministischen # Funktion aus, die bestimmen wir gleich # Ein q%-Quantil von x ist jetzt so eine Zu- # fallszahl x_q aus x, sodass q% aller Zahlen # in x kleiner sind als x_q und die restlichen # (100-q)% der Zahlen in x sind groesser als # x_q. Ein 25%-Quantil heisst Quartil und ein # 50%-Quantil heisst der Median, der Median liegt # also genau in der Mitte, die Haelfte der Zahlen # in x sind kleiner als der Median und die # andere Haelfte sind groesser der Median. # Ueberpruefen wir das: summary(x) sortx[N/2] # der Median - ok, nicht ganz.. sortx[N/2+1] # etwas zu gross.. (sortx[N/2]+sortx[N/2+1])/2 # das passt.. sortx[N/4] # erstes Quartil, bis auf obigen Effekt sortx[3*N/4] # drittes Quartil # Die Quantile koennen theoretisch berechnet werden: # Die Dichte der t_8-Verteilung ist gegeben durch # dt(x,df=8): curve(dt(x,df=8),from=-10,to=10,ylim=c(0,1)) # Das Integral ueber die Dichte, # # int_{-Infty}^x Dichte(y) dy # # ist gegeben durch die Funktion pt(x,df=8): curve(pt(x,df=8),col="blue",add=TRUE) # Wenn wir jetzt ein 25%-Quantil haben wollen, # muessen wir ein x_{q=25%} so bestimmen, dass # # int_{-Infty}^x_q Dichte(y) dy = 25% # # also: # # pt(x_q,df=8) = 25% (1) # # Wenn wir das x_q wissen wollen, brauchen wir # die Umkehrfunktion der pt()-Funktion, diese # Umkehrfunktion, auch Quantil-Funktion genannt, # ist ebenfalls in R eingebaut und hat die Syn- # tax qVerteilungsname, also hier qt: (1) ist # also aequivalent zu # # x_q = pt^{-1}(25%) = qt(25%,df=8) # # checken wir das: qt(0.25,df=8) summary(x) # ok, die x sind ja Zufallszahlen, # muss nicht exakt uebereinstimmen.. pt( qt(0.25,df=8) , df=8 ) # 0.25 qt(0.5,df=8) # nahe beim Median von x qt(0.75,df=8) # nahe beim 3.Quartil pt( qt(0.75,df=8) , df=8 ) # 0.75 curve(qt(x,df=8),from=0,to=1) plot(sortx) # das sieht sehr aehnlich aus, # versuchen wir das zu matchen: plot((1:N)/N,sortx) curve(qt(x,df=8),col="red",add=TRUE) # ok, das passt # Jetzt zurueck zum Regressions-Setting: Wir # haben: Die Testgroesse # # (hatbeta_j -beta_j) / sqrt( hatvar ) # # ist t_{n-(p+1)} - verteilt. Wenn man Ver- # trauensintervalle berechnen moechte, muss # man ein Vertrauenslevel alpha vorgeben, # etwa alpha = 90%: alpha = 0.9 # In 90% aller Faelle liegt eine t_8-verteilte # Zufallszahl x im Intervall [-x_q,x_q], wenn # # int_{-Infty}^{-x_q} Dichte(y)dy = 5% und # int_{-Infty}^{x_q} Dichte(y) dy = 95% # # Also: pt(-x_q,df=8) = 5% oder # pt(x_q,df=8) = 95% und damit # # x_q = pt^{-1}(95%,df=8) = qt(95%,df=8): xq = qt(0.95,df=8) xq minus_xq = qt(0.05,df=8) minus_xq # ok, ist gleich -xq # Wir haben also: In 90% aller Faelle ist # # (hatbeta_j -beta_j) / sqrt( hatvar ) in [-xq,xq] # # <=> hatbeta_j -beta_j in [-xq*sqrt(hatvar),xq*sqrt(hatvar)] # # <=> beta_j in (2) # [ hatbeta_j - xq*sqrt(hatvar) , hatbeta_j + xq*sqrt(hatvar) ] # # Die Gleichung (2) sagt also folgendes: Wenn wir # das Intervall auf der rechten Seite von (2) mit # xq = qt(0.95,df=8) berechnen, dann liegt in 90% # aller Faelle das wirkliche beta_j in diesem # Intervall. Ueberpruefen wir das mit einer geeigneten # Simulation: Wir betrachte dasselbe Regressions- # problem wie in Verteilungen-der-Schaetzer.txt : # 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 # folgendes brauchen wir fuer die Vertrauens- # intervalle: X0 = cbind(rep(1,n),X) # Regressoren mit Konstante XTXinv = solve(t(X0)%*%X0) # die Matrix (X^T*X)^{-1} XTXinv00 = XTXinv[1,1] XTXinv11 = XTXinv[2,2] XTXinv22 = XTXinv[3,3] xq = qt(0.95,df=8) # Vertrauenslevel alpha = 90% # obere und untere Intervallgrenzen: beta0up = rep(0,N) beta1up = rep(0,N) beta2up = rep(0,N) beta0down = rep(0,N) beta1down = rep(0,N) beta2down = rep(0,N) 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) # ist das hat(s^2) # jetzt die Vertrauensintervalle: beta0up[k] = hatbeta0[k] + xq*sqrt(hats[k]*XTXinv00) beta1up[k] = hatbeta1[k] + xq*sqrt(hats[k]*XTXinv11) beta2up[k] = hatbeta2[k] + xq*sqrt(hats[k]*XTXinv22) beta0down[k] = hatbeta0[k] - xq*sqrt(hats[k]*XTXinv00) beta1down[k] = hatbeta1[k] - xq*sqrt(hats[k]*XTXinv11) beta2down[k] = hatbeta2[k] - xq*sqrt(hats[k]*XTXinv22) } # Schauen wir uns die Intervallgrenzen mal an: plot(rep(beta0,N),col="red",ylim=c(-12,0)) # wirkliches beta points(beta0up) points(beta0down,col="blue") # ok, sieht nicht unplausibel aus.. # jetzt etwas quantitativer: test0 = ifelse( beta0beta0down , 1 , 0 ) sum(test0)/N # ok, sehr nahe bei alpha = 90%, all good.. # jetzt dasselbe fuer beta1 und beta2: # fuer beta1: plot(rep(beta1,N),col="red",ylim=c(-2,0)) # wirkliches beta points(beta1up) points(beta1down,col="blue") # ok, sieht nicht unplausibel aus.. test1 = ifelse( beta1beta1down , 1 , 0 ) sum(test1)/N # ok, sehr nahe bei alpha = 90% # fuer beta2: plot(rep(beta2,N),col="red",ylim=c(0,1)) # wirkliches beta points(beta2up) points(beta2down,col="blue") # ok, sieht nicht unplausibel aus.. test2 = ifelse( beta2beta2down , 1 , 0 ) sum(test2)/N # ok, sehr nahe bei alpha = 90%