Interval prediksi bootstrap

Apakah ada teknik bootstrap yang tersedia untuk menghitung interval prediksi untuk prediksi titik yang diperoleh misalnya dari regresi linier atau metode regresi lainnya (k-tetangga terdekat, pohon regresi dll)?

Entah bagaimana saya merasa bahwa cara yang kadang-kadang diusulkan untuk hanya menghapus prediksi titik (lihat misalnya interval prediksi untuk regresi kNN ) tidak memberikan interval prediksi tetapi interval kepercayaan.

Contoh dalam 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

Jelas, interval bootstrap dasar 95% cocok dengan interval kepercayaan 95%, bukan interval prediksi 95%. Jadi pertanyaan saya: Bagaimana cara melakukannya dengan benar?

Metode yang dijelaskan di bawah ini adalah yang dijelaskan dalam Bagian 6.3.3 dari Davidson dan Hinckley (1997), Metode Bootstrap dan Penerapannya . Terima kasih untuk Glen_b dan komentarnya di sini . Mengingat ada beberapa pertanyaan tentang Cross yang divalidasi pada topik ini, saya pikir itu layak ditulis.

Model regresi linier adalah:

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

Kami memiliki data, , yang kita gunakan untuk memperkirakan β sebagai: β$i=1,2,\ldots,N$$\beta$

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

Sekarang, kami ingin memprediksi apa yang akan menjadi untuk titik data baru, mengingat bahwa kami tahu X untuk itu. Ini adalah masalah prediksi. Mari panggilan baru X (yang kita tahu) X N + 1 dan baru Y (yang kita ingin memprediksi), Y N + 1 . Prediksi yang biasa (jika kita mengasumsikan bahwa ϵ i adalah iid dan tidak berkorelasi dengan X ) adalah: 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}$$

Kesalahan perkiraan yang dibuat oleh prediksi ini adalah:

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

Kita dapat menulis ulang persamaan ini seperti:

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

Sekarang, kita sudah dihitung. Jadi, jika kita ingin terikat Y N + 1 di interval, katakanlah, 90% dari waktu, semua yang perlu kita lakukan adalah perkiraan konsisten 5 t h dan 95 t h persentil / quantiles dari e p N + 1 , panggilan mereka e 5 , e 95 , dan interval prediksi akan [ Y p N + 1 + e 5 , Y p $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]$

Bagaimana cara memperkirakan kuantil / persentil dari ? Kita bisa menulis: 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}$$

Strateginya adalah dengan mengambil sampel (dengan cara bootstrap) berkali-kali dari dan kemudian menghitung persentil dengan cara biasa. Jadi, mungkin kita akan sampel 10.000 kali dari e p N + 1 , dan kemudian perkirakan persentil 5 t jam dan 95 t jam sebagai 500 t jam dan 9 , 500 t jam anggota terkecil dari sampel.$e^p_{N+1}$$e^p_{N+1}$$5^{th}$$95^{th}$$500^{th}$$9,500^{th}$

Untuk menggambar pada , kita dapat bootstrap kesalahan (kasus akan baik-baik saja, juga, tapi kita mengasumsikan kesalahan iid pula). Jadi, pada setiap replikasi bootstrap, Anda menggambar N kali dengan penggantian dari residual variance-disesuaikan (lihat butir berikutnya) untuk mendapatkan ε * i , kemudian membuat yang baru Y * i = X i β OLS + ε * i , kemudian jalankan OLS pada dataset baru, ( Y ∗ , X $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)$untuk mendapatkan replikasi ini . Akhirnya, hasil imbang ini replikasi pada X N + 1 ( β - β OLS ) adalah X N + 1 ( β OLS - β * r$\beta^*_r$$X_{N+1}\left( \beta-\hat{\beta}_{\text{OLS}} \right)$$X_{N+1}\left( \hat{\beta}_{\text{OLS}}-\beta^*_r \right)$

Mengingat kita mengasumsikan iid , cara alami untuk mengambil sampel dari bagian ϵ N + 1 dari persamaan adalah dengan menggunakan residu yang kita miliki dari regresi, { e ∗ 1 , e ∗ 2 , ... , e ∗ N } . Residual memiliki varian yang berbeda dan umumnya terlalu kecil, jadi kami ingin mengambil sampel dari { s 1 - ¯ s , s 2 - ¯ s , ... , s N - ¯ s$\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\}$, varians-residuals yang dikoreksi, di mana danhiadalah leverage dari pengamatani.$s_i=e^*_i/\sqrt{(1-h_i)}$$h_i$$i$

Dan, akhirnya, algoritma untuk membuat interval prediksi 90% untuk , mengingat bahwa X adalah X N + 1 adalah:$Y_{N+1}$$X$$X_{N+1}$

1. Membuat prediksi .$Y^p_{N+1}=X_{N+1}\hat{\beta}_{\text{OLS}}$
2. Buat residu yang disesuaikan dengan varians, , di mana s i = e i / $\left\{ s_1-\overline{s},s_2-\overline{s},\ldots,s_N-\overline{s}\right\}$.$s_i=e_i/\sqrt(1-h_{i})$
3. Untuk replikasi : $r=1,2,\ldots,R$
   • Gambar kali pada residu yang disesuaikan untuk membuat resid bootstrap { ϵ ∗ 1 , ϵ ∗ 2 , … , ϵ ∗ N$N$$\left\{\epsilon^*_1,\epsilon^*_2,\ldots,\epsilon^*_N \right\}$
   • Menghasilkan bootstrap $Y^*=X\hat{\beta}_{\text{OLS}}+\epsilon^*$
   • Hitung estimator bootstrap OLS untuk replikasi ini, $\beta^*_r=\left( X'X \right)^{-1}X'Y^*$
   • Dapatkan sisa bootstrap dari replikasi ini, $e^*_r=Y^*-X\beta^*_r$
   • Hitung residu yang disesuaikan dengan varian bootstrap dari replikasi ini, $s^*-\overline{s^*}$
   • Menarik salah satu bootstrap varians-disesuaikan residual dari replikasi ini, $\epsilon^*_{N+1,r}$
   • $e^p_{N+1}$$e^{p*}_r=X_{N+1}\left( \hat{\beta}_{\text{OLS}}-\beta^*_r \right)+\epsilon^*_{N+1,r}$
4. $5^{th}$$95^{th}$$e^p_{N+1}$$e^5,e^{95}$
5. $Y_{N+1}$$\left[Y^p_{N+1}+e^5,Y^p_{N+1}+e^{95} \right]$.

Here is R code:

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