问如何通过插入核包提取高斯过程回归的预测区间？

EN

rm(list = ls())
data(airquality)

airquality = na.omit(as.data.frame(airquality)); str(airquality)

library(tidyverse)
library(magrittr)

airquality$Ozone %>% plot(type = 'l')
lines(airquality$Temp, col = 2)
legend("topleft", legend = c("Ozone", "Temperature"),
   col=c(1, 2), lty = 1:1, cex = 0.7, text.font = 4, inset = 0.01, 
   box.lty=0, lwd = 1)

attach(airquality)
df_lags <- airquality %>%
  mutate(Temp_lag1 = lag(n = 1L, Temp)) %>%
  na.omit()

ESM_train = data.frame(df_lags[1:81, ])            # Training Observed 75% dataset
ESM_test = data.frame(df_lags[82:nrow(df_lags), ]) # Testing Observed 25% dataset

grid_gaussprRadial = expand.grid(.sigma = c(0.001, 0.01, 0.05, 0.1, 0.5, 1, 2)) # Sigma parameters searching for GPR

# TRAIN MODEL ############################
# Tuning set
library(caret)
set.seed(111)
cvCtrl <- trainControl(
  method ="repeatedcv",
  repeats = 1,
  number = 20,
  allowParallel = TRUE,
  verboseIter = TRUE,
  savePredictions = "final")

# Train (aprox. 4 seconds time-simulation)
attach(ESM_train)
set.seed(111)
system.time(Model_train <- caret::train(Ozone ~  Temp + Temp_lag1,
                                        trControl = cvCtrl,
                                        data = ESM_train,
                                        metric = "MAE", # Using MAE since I intend minimum values are my focus 
                                        preProcess = c("center", "scale"),
                                        method = "gaussprRadial", # Setting RBF kernel function
                                        tuneGrid = grid_gaussprRadial,
                                        maxit = 1000,
                                        linout = 1)) # Regression type

plot(Model_train)
Model_train
ESM_results_train <- Model_train$resample %>% mutate(Model = "") # K-fold Training measures

# Select the interested TRAIN data and arrange them as dataframe
Ozone_Obs_Tr = Model_train$pred$obs
Ozone_sim = Model_train$pred$pred
Resid = Ozone_Obs_Tr - Ozone_sim
train_results = data.frame(Ozone_Obs_Tr,
                           Ozone_sim,
                           Resid)

# Plot Obs x Simulated train results
library(ggplot2)
ggplot(data = train_results, aes(x = Ozone_Obs_Tr, y = Ozone_sim)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "black")


# TEST MODEL ############################
# From "ESM_test" dataframe, we predict ESM Ozone time series, adding it in "ESM_forecasted" dataframe
ESM_forecasted = ESM_test %>%                                              
  mutate(Ozone_Pred = predict(Model_train, newdata = ESM_test, variance.model = TRUE))
str(ESM_forecasted)

# Select the interested TEST data and arrange them as a dataframe
Ozone_Obs = ESM_forecasted$Ozone
Ozone_Pred = ESM_forecasted$Ozone_Pred

# Plot Obs x Predicted TEST results
ggplot(data = ESM_forecasted, aes(x = Ozone_Obs, y = Ozone_Pred)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "black")


# Model performance #####
library(hydroGOF)
gof_TR = gof(Ozone_sim, Ozone_Obs_Tr)
gof_TEST = gof(Ozone_Pred,Ozone_Obs)
Performances = data.frame(
                          Train = gof_TR,
                          Test = gof_TEST
                          ); Performances
# Plot the TEST prediction
attach(ESM_forecasted)
plot(Ozone_Obs, type = "l", xlab = "", ylab = "", ylim = range(Ozone_Obs, Ozone_Pred))
lines(Ozone_Pred , col = "coral2", lty = 2, lwd = 2)
legend("top", legend = c("Ozone Obs Test", "Ozone Pred Test"),
       col=c(1, "coral2"), lty = 1:2, cex = 0.7, text.font = 4, inset = 0.01, box.lty=0, lwd = 2)

    Upper_Bound <- Ozone_Pred$ZZ.q2 #Ozone_Pred - 2 * sigma^2 
    Lower_Bound <- Ozone_Pred$ZZ.q1 #Ozone_Pred + 2 * sigma^2

Model_train$...