我正在尝试使用高斯过程回归(GPR)模型来预测河流中的每小时径流流量。我已经得到了很好的结果，我应用了脱：：内核实验室的列车()函数(感谢Kuhn！)。

由于不确定性概念是GPR的主要固有优势之一，我想知道是否有人可以帮助我访问与测试数据集的预测积分相关的结果。

我将摘录我一直使用的代码。由于我的真实数据非常庞大(老实说，我不知道怎么说才好)，我将以数据(空气质量)为例。此特定示例中的主要目标是使用空气质量$温度作为预测滞后变量来预测空气质量$臭氧。

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)

最后几行将生成以下图形：

下一步，也是最后一步，将提取基于每个预测点周围的高斯分布的预测间隔，以将其与最后一张图一起绘制。

例如，与仅使用kernLab：：gausspr Radial()或甚至tgp：：bgp()包相比，插入：：kernLab Train()设备返回更好的预测。我可以找到这两种情况的预测区间。

例如，要通过tgp：：bgp()获取预测间隔，可以输入：

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

因此，通过脱字符：：kernLab Train()，我希望可以找到所需的标准差，如

Model_train$...

或者可能，使用

Ozone_Pred$...

此外，在链接：https://stats.stackexchange.com/questions/414079/can-mad-median-absolute-deviation-or-mae-mean-absolute-error-be-used-to-calc处， 作者斯蒂芬·科拉萨解释说，我们可以通过MAE，甚至RMSE来估计预测区间。但我不明白这是否是我的意思，因为在这个例子中，我得到的MAE只是Obs x预测的臭氧数据之间的比较。

拜托，这个解决方案对我来说非常重要！我想我已经接近取得主要成果了，但我不知道该如何去尝试。 非常感谢，朋友们！

推荐答案

我真的不知道caret框架的工作原理，但手动获取具有高斯似然率的GP回归的预测间隔非常容易。

首先，我们只需要一个平方指数核的函数，也称为径向基函数核，这就是您所使用的。sf这里是比例因子(kernlab实现中未使用)，ell是长度比例，kernlab实现中称为sigma：

covSEiso <- function(x1, x2 = x1, sf = 1.0, ell = 1.0) {
    sf     <- sf^2
    ell    <- -0.5 * (1 / (ell^2))
    n      <- nrow(x1)
    m      <- nrow(x2)
    d      <- ncol(x1)
    result <- matrix(0, nrow = n, ncol = m)
    for ( j in 1:m ) {
        for ( i in 1:n ) {
            result[i, j] <- sf * exp(ell * sum((x1[i, ] - x2[j, ])^2))
        }
    }
    return(result)
}

我不确定您的代码对使用哪种长度比例表示了什么；下面我将使用长度比例25和比例因子50(通过GPML的超参数优化例程获得)。然后我们使用上面的covSEiso()函数得到相关的协方差，剩下的就是基本高斯恒等式的应用。我建议您参考Rasmussen and Williams (2006)的第二章(在线免费提供)。

data(airquality)
library(tidyverse)
library(magrittr)
df_lags <- airquality %>%
    mutate(Temp_lag1 = lag(n = 1L, Temp)) %>%
    na.omit()
ESM_train <- data.frame(df_lags[1:81, ])             # Training Data 75% dataset
ESM_test  <- data.frame(df_lags[82:nrow(df_lags), ]) # Testing  Data 25% dataset
## For convenience I'll define separately the training and test inputs
X <- ESM_train[ , c("Temp", "Temp_lag1")]
Xstar <- ESM_test[ , c("Temp", "Temp_lag1")]
## Get the kernel manually
K <- covSEiso(X, ell = 25, sf = 50)
## We also need covariance between the test cases
Kstar <- covSEiso(Xstar, X, ell = 25, sf = 50)
Ktest <- covSEiso(Xstar, ell = 25, sf = 50)
## Now the 95% credible region for the posterior is
predictive_mean <- Kstar %*% solve(K + diag(nrow(K))) %*% ESM_train$Ozone
predictive_var  <- Ktest - (Kstar %*% solve(K + diag(nrow(K))) %*% t(Kstar))
## Then for the prediction interval we only need to add the observation noise
z <- sqrt(diag(predictive_var)) + 25
interval_high <- predictive_mean + 2 * z
interval_low <- predictive_mean - 2 * z

然后我们可以检查预测区间

这一切都很容易通过我的gplmr包(available on GitHub)完成，如果您安装了Octave：

，它可以从R调用GPML

data(airquality)
library(tidyverse)
library(magrittr)
library(gpmlr)
df_lags <- airquality %>%
    mutate(Temp_lag1 = lag(n = 1L, Temp)) %>%
    na.omit()
ESM_train <- data.frame(df_lags[1:81, ])             # Training Data 75% dataset
ESM_test  <- data.frame(df_lags[82:nrow(df_lags), ]) # Testing  Data 25% dataset
X <-  as.matrix(ESM_train[ , c("Temp", "Temp_lag1")])
y <- ESM_train$Ozone
Xs <- as.matrix(ESM_test[ , c("Temp", "Temp_lag1")])
ys <- ESM_test$Ozone
hyp0 <- list(mean = numeric(), cov = c(0, 0), lik = 0)
hyp  <- set_hyperparameters(hyp0, "infExact", "meanZero", "covSEiso","likGauss",
                            X, y)
gp_res <- gp(hyp, "infExact", "meanZero", "covSEiso", "likGauss", X, y, Xs, ys)
predictive_mean <- gp_res$YMU
interval_high   <- gp_res$YMU + 2 * sqrt(gp_res$YS2)
interval_low    <- gp_res$YMU - 2 * sqrt(gp_res$YS2)

然后将预测绘制出来，如下所示：

plot(NULL, xlab = "", ylab = "", xaxt = "n", yaxt = "n",
     xlim = range(ESM_test$Temp), ylim = range(c(interval_high, interval_low)))
axis(1, tick = FALSE, line = -0.75)
axis(2, tick = FALSE, line = -0.75)
mtext("Temp", 1, 1.5)
mtext("Ozone", 2, 1.5)
idx <- order(ESM_test$Temp)
polygon(c(ESM_test$Temp[idx], rev(ESM_test$Temp[idx])),
        c(interval_high[idx], rev(interval_low[idx])),
        border = NA, col = "#80808080")
lines(ESM_test$Temp[idx], predictive_mean[idx])
points(ESM_test$Temp, ESM_test$Ozone, pch = 19)
plot(NULL, xlab = "", ylab = "", xaxt = "n", yaxt = "n",
     xlim = range(ESM_test$Temp_lag1), ylim = range(c(interval_high, interval_low)))
axis(1, tick = FALSE, line = -0.75)
axis(2, tick = FALSE, line = -0.75)
mtext("Temp_lag1", 1, 1.5)
mtext("Ozone", 2, 1.5)
idx <- order(ESM_test$Temp_lag1)
polygon(c(ESM_test$Temp_lag1[idx], rev(ESM_test$Temp_lag1[idx])),
        c(interval_high[idx], rev(interval_low[idx])),
        border = NA, col = "#80808080")
lines(ESM_test$Temp_lag1[idx], predictive_mean[idx])
points(ESM_test$Temp_lag1, ESM_test$Ozone, pch = 19)