R语言中xgbar(forecastxgb包)与forecast函数预测置信区间获取方法咨询

Great question! Let's break this down: the xgbar function from the forecastxgb package is built on XGBoost, a tree-based model that doesn't natively produce prediction intervals like traditional time series models (think ARIMA or exponential smoothing). That's why your current forecast() call only returns point predictions. But there are two solid ways to generate confidence intervals for your XGBoost time series forecast:

Method 1: Generate Intervals via Simulation #

You can simulate multiple prediction paths using your trained xgbar model, then calculate the quantiles of these simulations to get your confidence bounds. Here's how:

# First, train your model as you did before
model <- xgbar(train, seas_method = "fourier", maxlag=200, K = max(1, min(round(f/4 - 1), 10)))

# Generate 1000 simulated prediction paths for your forecast horizon
simulated_predictions <- predict(model, h = weeks_predicted, nsim = 1000)

# Calculate 95% confidence interval bounds (2.5th and 97.5th percentiles)
lower_bound <- apply(simulated_predictions, 1, quantile, probs = 0.025)
upper_bound <- apply(simulated_predictions, 1, quantile, probs = 0.975)

# Create a forecast object with the interval bounds included
fit_xgb_with_ci <- forecast(model, h = weeks_predicted)
fit_xgb_with_ci$lower <- matrix(lower_bound, ncol = 1)
fit_xgb_with_ci$upper <- matrix(upper_bound, ncol = 1)
colnames(fit_xgb_with_ci$lower) <- colnames(fit_xgb_with_ci$upper) <- "95%"

This method works by leveraging the model's ability to generate multiple plausible outcomes, then using the spread of these outcomes to estimate uncertainty. Adjust nsim (number of simulations) if you want more stable interval estimates (higher values = more stable, but slower computation).

Method 2: Train Quantile Regression Models #

Another approach is to directly train XGBoost models targeted at specific quantiles (e.g., 2.5%, 50%, 97.5%). This gives you explicit predictions for the bounds of your confidence interval:

# Train models for 2.5th, 50th (median/point forecast), and 97.5th quantiles
model_q025 <- xgbar(
  train, 
  seas_method = "fourier", 
  maxlag=200, 
  K = max(1, min(round(f/4 - 1), 10)),
  objective = "reg:quantileerror", 
  alpha = 0.025
)

model_q50 <- xgbar(
  train, 
  seas_method = "fourier", 
  maxlag=200, 
  K = max(1, min(round(f/4 - 1), 10))
) # Default objective is regression, which approximates the median

model_q975 <- xgbar(
  train, 
  seas_method = "fourier", 
  maxlag=200, 
  K = max(1, min(round(f/4 - 1), 10)),
  objective = "reg:quantileerror", 
  alpha = 0.975
)

# Generate predictions for each quantile
pred_q025 <- forecast(model_q025, h = weeks_predicted)$mean
pred_q50 <- forecast(model_q50, h = weeks_predicted)$mean
pred_q975 <- forecast(model_q975, h = weeks_predicted)$mean

# Build a forecast object with interval bounds
fit_xgb_with_ci <- forecast(model_q50, h = weeks_predicted)
fit_xgb_with_ci$lower <- matrix(pred_q025, ncol = 1)
fit_xgb_with_ci$upper <- matrix(pred_q975, ncol = 1)
colnames(fit_xgb_with_ci$lower) <- colnames(fit_xgb_with_ci$upper) <- "95%"

This method gives you more direct control over the quantiles you're predicting, but requires training multiple models (one per quantile of interest).

Once you've created fit_xgb_with_ci, you can use standard forecast package functions like plot(fit_xgb_with_ci) to visualize your predictions alongside the confidence intervals.

内容的提问来源于stack exchange，提问作者M_D