Interwał przewidywania ładowania początkowego

Czy jest dostępna technika ładowania początkowego do obliczania przedziałów predykcji dla prognoz punktowych uzyskanych np. Z regresji liniowej lub innej metody regresji (k-najbliższy sąsiad, drzewa regresji itp.)?

Jakoś wydaje mi się, że czasami proponowanym sposobem, aby po prostu wyrzucić prognozę punktową (patrz np. Przedziały predykcji dla regresji kNN ), nie zapewnia przedziału prognoz, ale przedział ufności.

Przykład w R.

# STEP 1: GENERATE DATA

set.seed(34345)

n <- 100 
x <- runif(n)
y <- 1 + 0.2*x + rnorm(n)
data <- data.frame(x, y)


# STEP 2: COMPUTE CLASSIC 95%-PREDICTION INTERVAL
fit <- lm(y ~ x)
plot(fit) # not shown but looks fine with respect to all relevant aspects

# Classic prediction interval based on standard error of forecast
predict(fit, list(x = 0.1), interval = "p")
# -0.6588168 3.093755

# Classic confidence interval based on standard error of estimation
predict(fit, list(x = 0.1), interval = "c")
# 0.893388 1.54155


# STEP 3: NOW BY BOOTSTRAP
B <- 1000
pred <- numeric(B)
for (i in 1:B) {
  boot <- sample(n, n, replace = TRUE)
  fit.b <- lm(y ~ x, data = data[boot,])
  pred[i] <- predict(fit.b, list(x = 0.1))
}
quantile(pred, c(0.025, 0.975))
# 0.8699302 1.5399179

Oczywiście 95% podstawowy interwał ładowania początkowego pasuje do 95% przedziału ufności, a nie 95% przedziału predykcji. Więc moje pytanie: jak to zrobić poprawnie?

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:

$$\begin{align} Y_i &= X_i\beta+\epsilon_i \end{align}$$

Mamy dane, , których używamy do oszacowania jako: beta beta$i=1,2,\ldots,N$$\beta$

$$\begin{align} \hat{\beta}_{\text{OLS}} &= \left( X'X \right)^{-1}X'Y \end{align}$$

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 + $Y$$X$$X$$X_{N+1}$$Y$$Y_{N+1}$$\epsilon_i$$X$

$$\begin{align} Y^p_{N+1} &= X_{N+1}\hat{\beta}_{\text{OLS}} \end{align}$$

Błąd prognozy spowodowany tą prognozą to:

$$\begin{align} e^p_{N+1} &= Y_{N+1}-Y^p_{N+1} \end{align}$$

Możemy ponownie napisać to równanie, takie jak:

$$\begin{align} Y_{N+1} &= Y^p_{N+1} + e^p_{N+1} \end{align}$$

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$Y^p_{N+1}$$Y_{N+1}$$5^{th}$$95^{th}$$e^p_{N+1}$$e^5,e^{95}$$\left[Y^p_{N+1}+e^5,Y^p_{N+1}+e^{95} \right]$

Jak oszacować kwantyle / percentyle ? Cóż, możemy napisać: e p N + $e^p_{N+1}$

$$\begin{align} e^p_{N+1} &= Y_{N+1}-Y^p_{N+1}\\ &= X_{N+1}\beta + \epsilon_{N+1} - X_{N+1}\hat{\beta}_{\text{OLS}}\\ &= X_{N+1}\left( \beta-\hat{\beta}_{\text{OLS}} \right) + \epsilon_{N+1} \end{align}$$

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 $e^p_{N+1}$$e^p_{N+1}$$5^{th}$$95^{th}$$500^{th}$$9,500^{th}$

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 - β * $X_{N+1}\left( \beta-\hat{\beta}_{\text{OLS}} \right)$$N$$\epsilon^*_i$$Y^*_i=X_i\hat{\beta}_{\text{OLS}}+\epsilon^*_i$$\left(Y^*,X \right)$$\beta^*_r$$X_{N+1}\left( \beta-\hat{\beta}_{\text{OLS}} \right)$$X_{N+1}\left( \hat{\beta}_{\text{OLS}}-\beta^*_r \right)$

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 .$\epsilon$$\epsilon_{N+1}$$\left\{ e^*_1,e^*_2,\ldots,e^*_N \right\}$$\left\{ s_1-\overline{s},s_2-\overline{s},\ldots,s_N-\overline{s} \right\}$$s_i=e^*_i/\sqrt{(1-h_i)}$$h_i$$i$

I wreszcie algorytm tworzenia 90% przedziału predykcji dla , biorąc pod uwagę, że to to:$Y_{N+1}$$X$$X_{N+1}$

1. Dokonaj prognozy .$Y^p_{N+1}=X_{N+1}\hat{\beta}_{\text{OLS}}$
2. Dokonaj reszt skorygowanych o wariancję, , gdzie .$\left\{ s_1-\overline{s},s_2-\overline{s},\ldots,s_N-\overline{s}\right\}$$s_i=e_i/\sqrt(1-h_{i})$
3. Dla replikacji : $r=1,2,\ldots,R$
   • Narysuj razy na skorygowanych resztkach, aby zrobić resztki ładowania początkowego $N$$\left\{\epsilon^*_1,\epsilon^*_2,\ldots,\epsilon^*_N \right\}$
   • Wygeneruj bootstrap$Y^*=X\hat{\beta}_{\text{OLS}}+\epsilon^*$
   • Oblicz estymator OLS bootstrap dla tej replikacji, $\beta^*_r=\left( X'X \right)^{-1}X'Y^*$
   • Uzyskaj resztki ładowania początkowego z tej replikacji,$e^*_r=Y^*-X\beta^*_r$
   • Oblicz resztki skorygowane o wariancję bootstrap z tej replikacji,$s^*-\overline{s^*}$
   • Narysuj jedną z reszt skorygowanych o wariancję bootstrapu z tej replikacji,$\epsilon^*_{N+1,r}$
   • Oblicz losowanie tej replikacji na ,$e^p_{N+1}$$e^{p*}_r=X_{N+1}\left( \hat{\beta}_{\text{OLS}}-\beta^*_r \right)+\epsilon^*_{N+1,r}$
4. Znajdź i percentyli ,$5^{th}$$95^{th}$$e^p_{N+1}$$e^5,e^{95}$
5. 90% przedział predykcji dla to .$Y_{N+1}$$\left[Y^p_{N+1}+e^5,Y^p_{N+1}+e^{95} \right]$

Oto Rkod:

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