¿Los "intervalos de incertidumbre" del Profeta son intervalos de confianza o intervalos de predicción?

Question

De forma predeterminada, Prophet devolverá intervalos de incertidumbre para el pronóstico yhat.

Lamentablemente, la documentación sobre esos "intervalos de incertidumbre" es extremadamente vaga, y no ayuda que los estadísticos tengan definiciones precisas tanto de intervalos de confianza como de intervalos de predicción (¡hay una gran diferencia!), pero no de "intervalos de incertidumbre". De hecho, la documentación podría interpretarse como refiriéndose a cualquiera de los dos, aunque una interpretación de "intervalo de predicción" me parece un poco más probable.

¿Alguien puede decirnos (por ejemplo, investigando en el código) si estos son intervalos de confianza o de predicción?

Answer 1

Al revisar el código como sugieres, parece que son intervalos de predicción. Específicamente, el modelo se ajusta muestreando un posterior con Stan (línea 1249-1266 en prophet.R) y el intervalo se basa en extracciones de la distribución predictiva posterior. Ver línea 1542 del archivo prophet.R, que reproduzco a continuación.

#' Intervalos de incertidumbre de Prophet para yhat y tendencia
#'
#' @param m Objeto Prophet.
#' @param df Dataframe de predicción.
#'
#' @return Dataframe con intervalos de incertidumbre.
#'
#' @keywords internal
predict_uncertainty <- function(m, df) {
  sim.values <- sample_posterior_predictive(m, df)
  # Agregar estimaciones de incertidumbre
  lower.p <- (1 - m$interval.width)/2
  upper.p <- (1 + m$interval.width)/2

  intervals <- cbind(
    t(apply(t(sim.values$yhat), 2, stats::quantile, c(lower.p, upper.p),
        na.rm = TRUE)),
t(apply(t(sim.values$trend), 2, stats::quantile, c(lower.p, upper.p),
            na.rm = TRUE))
  )

  colnames(intervals) <- paste(rep(c('yhat', 'trend'), each=2),
                               c('lower', 'upper'), sep = "_")

  return(dplyr::as_tibble(intervals))
}

Ver también la función sample_model que es llamada desde sample_posterior_predictive, la cual muestra exactamente cómo se calculan los yhats en línea 1573:

sample_model <- function(m, df, seasonal.features, iteration, s_a, s_m) {
  trend <- sample_predictive_trend(m, df, iteration)

  beta <- m$params$beta[iteration,]
  Xb_a = as.matrix(seasonal.features) %*% (beta * s_a) * m$y.scale
  Xb_m = as.matrix(seasonal.features) %*% (beta * s_m)

  sigma <- m$params$sigma_obs[iteration]
  noise <- stats::rnorm(nrow(df), mean = 0, sd = sigma) * m$y.scale

  return(list("yhat" = trend * (1 + Xb_m) + Xb_a + noise,
              "trend" = trend))
}