Metodą przedstawioną poniżej jest metoda opisana w sekcji 6.3.3 Davidson i Hinckley (1997),
Metody ładowania początkowego i ich zastosowanie . Dzięki Glen_b i jego komentarzowi tutaj . Biorąc pod uwagę, że było kilka pytań na temat Cross Validated na ten temat, pomyślałem, że warto to napisać.
Model regresji liniowej to:
Yja= Xjaβ+ ϵja
Mamy dane, , których używamy do oszacowania jako:
beta beta OLSi = 1 , 2 , … , Nβ
β^OLS= ( X′X)- 1X′Y
Teraz chcemy przewidzieć, jaki będzie dla nowego punktu danych, biorąc pod uwagę, że znamy dla niegoTo jest problem z prognozowaniem. Nazwijmy nowy (który znamy) i nowy (który chcielibyśmy przewidzieć), . Zwykła prognoza (jeśli zakładamy, że są iid i nieskorelowane z ) to:
X X X N + 1 Y Y N + 1 ϵ i X Y p N + 1YXXXN.+ 1YYN.+ 1ϵjaX
YpN.+ 1= XN.+ 1β^OLS
Błąd prognozy spowodowany tą prognozą to:
mipN.+ 1= YN.+ 1- YpN.+ 1
Możemy ponownie napisać to równanie, takie jak:
YN.+ 1= YpN.+ 1+ epN.+ 1
Teraz już obliczyliśmy. Tak więc, jeśli chcemy związać w przedziale, powiedzmy, w 90% przypadków, wystarczy oszacować konsekwentnie i percentyle / kwantyle , nazwij je , a przedział przewidywania będzie . Y N + 1 5 t h 95 t h e p N + 1 e 5 , e 95 [ Y p N + 1 + e 5 , Y p N + 1 + e 95 ]YpN.+ 1YN.+ 15t godz95t godzmipN.+ 1mi5, e95[ YpN.+ 1+ e5, YpN.+ 1+ e95]
Jak oszacować kwantyle / percentyle ? Cóż, możemy napisać:
e p N + 1mipN.+ 1
mipN.+ 1= YN.+ 1- YpN.+ 1= XN.+ 1β+ ϵN.+ 1- XN.+ 1β^OLS= XN.+ 1( β- β^OLS) + ϵN.+ 1
Strategia polega na próbkowaniu (w pewnym sensie bootstrapu) wiele razy z a następnie obliczaniu percentyli w zwykły sposób. Może więc spróbujemy 10 000 razy z , a następnie oszacujemy percentyle i jako i najmniejszych członków próbka. e p N + 1 5 t h 95 t h 500 t h 9 , 500 t hmipN.+ 1mipN.+ 15t godz95t godz500t godz9 , 500t godz
Aby narysować , możemy załadować błędy (przypadki też byłyby w porządku, ale zakładamy, że iid błędy mimo to). Tak więc przy każdej replikacji bootstrap rysujesz razy z zamianą z reszt skorygowanych o wariancję (patrz następny akapit), aby uzyskać , a następnie utworzyć nowy , a następnie uruchom OLS na nowym zestawie danych, aby uzyskać tej repliki . Nareszcie losowanie tej repliki na to N ε * i Y * i = X i β OLS + ε * i ( Y * , X ) β * r X N + 1 ( β - β OLS ) X N + 1 ( β OLS - β * rXN.+ 1( β- β^OLS)N.ϵ∗jaY∗ja= Xjaβ^OLS+ ϵ∗ja( Y∗, X)β∗rXN.+ 1( β- β^OLS)XN.+ 1( β^OLS- β∗r)
Biorąc pod uwagę, że zakładamy, że iid , naturalnym sposobem próbkowania z części równania jest użycie resztek z regresji, . Resztki mają różne i ogólnie zbyt małe wariancje, więc będziemy chcieli próbować z , reszty z korektą wariancji, gdzie i jest dźwignią obserwacji .ϵϵN.+ 1{ e∗1, e∗2), … , E∗N.}{ s1- s¯¯¯, s2)- s¯¯¯, … , SN.- s¯¯¯}sja= e∗ja/ ( 1 - godzja)------√hjaja
I wreszcie algorytm tworzenia 90% przedziału predykcji dla , biorąc pod uwagę, że to to:YN.+ 1XXN.+ 1
- Dokonaj prognozy .YpN.+ 1= XN.+ 1β^OLS
- Dokonaj reszt skorygowanych o wariancję, , gdzie .{ s1- s¯¯¯, s2)- s¯¯¯, … , SN.- s¯¯¯}sja= eja/ (√1 - godzja)
- Dla replikacji :
r = 1 , 2 , … , R
- Narysuj razy na skorygowanych resztkach, aby zrobić resztki ładowania początkowego
N.{ ϵ∗1, ϵ∗2), … , Ε∗N.}
- Wygeneruj bootstrapY∗= Xβ^OLS+ ϵ∗
- Oblicz estymator OLS bootstrap dla tej replikacji,
β∗r= ( X′X)- 1X′Y∗
- Uzyskaj resztki ładowania początkowego z tej replikacji,mi∗r= Y∗- Xβ∗r
- Oblicz resztki skorygowane o wariancję bootstrap z tej replikacji,s∗- s∗¯¯¯¯¯
- Narysuj jedną z reszt skorygowanych o wariancję bootstrapu z tej replikacji,ϵ∗N.+ 1 , r
- Oblicz losowanie tej replikacji na
,mipN.+ 1mip ∗r= XN.+ 1( β^OLS- β∗r) + ϵ∗N.+ 1 , r
- Znajdź i percentyli ,5t godz95t godzmipN.+ 1mi5, e95
- 90% przedział predykcji dla to
.YN.+ 1[ YpN.+ 1+ e5, YpN.+ 1+ e95]
Oto R
kod:
# This script gives an example of the procedure to construct a prediction interval
# for a linear regression model using a bootstrap method. The method is the one
# described in Section 6.3.3 of Davidson and Hinckley (1997),
# _Bootstrap Methods and Their Application_.
#rm(list=ls())
set.seed(12344321)
library(MASS)
library(Hmisc)
# Generate bivariate regression data
x <- runif(n=100,min=0,max=100)
y <- 1 + x + (rexp(n=100,rate=0.25)-4)
my.reg <- lm(y~x)
summary(my.reg)
# Predict y for x=78:
y.p <- coef(my.reg)["(Intercept)"] + coef(my.reg)["x"]*78
y.p
# Create adjusted residuals
leverage <- influence(my.reg)$hat
my.s.resid <- residuals(my.reg)/sqrt(1-leverage)
my.s.resid <- my.s.resid - mean(my.s.resid)
reg <- my.reg
s <- my.s.resid
the.replication <- function(reg,s,x_Np1=0){
# Make bootstrap residuals
ep.star <- sample(s,size=length(reg$residuals),replace=TRUE)
# Make bootstrap Y
y.star <- fitted(reg)+ep.star
# Do bootstrap regression
x <- model.frame(reg)[,2]
bs.reg <- lm(y.star~x)
# Create bootstrapped adjusted residuals
bs.lev <- influence(bs.reg)$hat
bs.s <- residuals(bs.reg)/sqrt(1-bs.lev)
bs.s <- bs.s - mean(bs.s)
# Calculate draw on prediction error
xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"]
xb.xb <- xb.xb + (coef(my.reg)["x"] - coef(bs.reg)["x"])*x_Np1
return(unname(xb.xb + sample(bs.s,size=1)))
}
# Do bootstrap with 10,000 replications
ep.draws <- replicate(n=10000,the.replication(reg=my.reg,s=my.s.resid,x_Np1=78))
# Create prediction interval
y.p+quantile(ep.draws,probs=c(0.05,0.95))
# prediction interval using normal assumption
predict(my.reg,newdata=data.frame(x=78),interval="prediction",level=0.90)
# Quick and dirty Monte Carlo to see which prediction interval is better
# That is, what are the 5th and 95th percentiles of Y_{N+1}
#
# To do it properly, I guess we would want to do the whole procedure above
# 10,000 times and then see what percentage of the time each prediction
# interval covered Y_{N+1}
y.np1 <- 1 + 78 + (rexp(n=10000,rate=0.25)-4)
quantile(y.np1,probs=c(0.05,0.95))