Cómo encontrar un intervalo de predicción GBM

Estoy trabajando con modelos de GBM usando el paquete caret y buscando encontrar un método para resolver los intervalos de predicción de mis datos pronosticados. He buscado mucho, pero solo se me ocurren algunas ideas para encontrar intervalos de predicción para Random Forest. ¡Cualquier código de ayuda / R sería muy apreciado!

EDITAR: Como se señala en los comentarios a continuación, esto proporciona los intervalos de confianza para las predicciones y no estrictamente los intervalos de predicción . Estaba un poco desencadenante contento con mi respuesta y debería haber pensado un poco más en esto.

Siéntase libre de ignorar esta respuesta o intente construir sobre el código para obtener los intervalos de predicción.

He usado el programa de arranque simple para crear intervalos de predicción varias veces, pero puede haber otras (mejores) formas.

Considere los oildatos en el caretpaquete y suponga que queremos generar dependencias parciales e intervalos del 95% para el efecto de Stearic en Palmitic. A continuación se muestra solo un ejemplo simple, pero puede jugar con él para satisfacer sus necesidades. Asegúrese de que el gbmpaquete esté actualizado para permitir el grid.pointsargumento enplot.gbm

library(caret)
data(oil)
#train the gbm using just the defaults.
tr <- train(Palmitic ~ ., method = "gbm" ,data = fattyAcids, verbose = FALSE)

#Points to be used for prediction. Use the quartiles here just for illustration
x.pt <- quantile(fattyAcids$Stearic, c(0.25, 0.5, 0.75))

#Generate the predictions, or in this case, the partial dependencies at the selected points. Substitute plot() for predict() to get predictions
p <- plot(tr$finalModel, "Stearic", grid.levels = x.pt, return.grid = TRUE)

#Bootstrap the process to get prediction intervals
library(boot)

bootfun <- function(data, indices) {
  data <- data[indices,]

  #As before, just the defaults in this example. Palmitic is the first variable, hence data[,1]
  tr <- train(data[,-1], data[,1], method = "gbm", verbose=FALSE)

  # ... other steps, e.g. using the oneSE rule etc ...
  #Return partial dependencies (or predictions)

  plot(tr$finalModel, "Stearic", grid.levels = x.pt, return.grid = TRUE)$y
  #or predict(tr$finalModel, data = ...)
}

#Perform the bootstrap, this can be very time consuming. Just 99 replicates here but we usually want to do more, e.g. 500. Consider using the parallel option
b <- boot(data = fattyAcids, statistic = bootfun, R = 99)

#Get the 95% intervals from the boot object as the 2.5th and 97.5th percentiles
lims <- t(apply(b$t, 2, FUN = function(x) quantile(x, c(0.025, 0.975))))

Esta es una forma de hacerlo que, al menos, intenta dar cuenta de las incertidumbres que surgen al ajustar el gbm. Un enfoque similar se ha utilizado en http://onlinelibrary.wiley.com/doi/10.2193/2006-503/abstract

A veces, la estimación puntual está fuera del intervalo, pero la modificación de la cuadrícula de ajuste (es decir, el aumento del número de árboles y / o la profundidad) generalmente resuelve eso.

¡Espero que esto ayude!