时间序列预测的高斯过程回归不确定性估计

问题描述

我有一个 2043 小时的时间序列。我的目标是通过高斯过程回归 (GPR) 模型训练和测试预测模型，包括其不确定性估计。我在应用 R 包 caret + kernlab::train( ) 联合模型的准时预测方面取得了很好的结果。我如何使用这些 R 包来估计其预测区间？

我的 .csv 数据可通过三个 DropBox 链接获得：

Training data

Testing data

Complete Time series data.frame (Time,train and test data)

目的是使用“ST3”和“ST4_lag1”回归量预测“ST4_lead1”。它们是我们在 data.frames 上的列。

从现在开始，我将展示我的代码。第一个窗口代码对我们俩来说都是通用的。您可以完全运行它。

# Forecasting dataseries: caret + kernlab packages
################ Load R packages
library(readr)
library(tidyverse)
library(caret)
library(kernlab)
library(hydroGOF)
library(lubridate)


################ Load data
setwd("C:/...") # Set the folder you saved the downloaded .csv data

data_train_emd <- read.csv("data_train_emd.csv",sep = ';')
data_test_emd <- read.csv("data_test_emd.csv",sep = ';')
final_df_emd <- read.csv("final_df_emd.csv",sep = ';')


################ GaussprRadial model Preprocess
set.seed(111)
cvCtrl <- trainControl(
  method = "repeatedcv",# Cross-correlation preprocess
  repeats = 1,number = 3,allowParallel = TRUE,verboseIter = TRUE,savePredictions = "final")


################ Set gaussprRadial grid for tuning parameters
grid_gaussprRadial <- expand.grid(sigma = c(1*1e-3,5*1e-3))


################ Implement Gaussian Process
################ Train
set.seed(111)
attach(data_train_emd)
system.time(Fitting_model <- train(ST4_lead1 ~ ST3 + ST4_lag1,# Predict "ST4_lead1" based on regressors "ST3" and "ST4_lag1"
                              method ="gaussprRadial",# Radial Basis kernel function
                              trControl = cvCtrl,data = data_train_emd,metric = "MAE",# Using "MAE" as fitting parameters,since I'm focusing on low data
                              preProcess = c("center","scale"),# Centering and scaling data Preprocess
                              tuneGrid = grid_gaussprRadial,maxit = 1000,linout = 1))                          # 1 for regression model

#The training process lasts 14 seconds on my machine.



################ Test
attach(final_df_emd)
forecasted <- data_test_emd %>%                                           
  mutate(Qpred = predict(Fitting_model,newdata = data_test_emd),Time = final_df_emd %>% filter(Type == "Test") %>% pull(Time))   


################ Tidy training results
################ View metrics for training process
results_gaussprRadial <- Fitting_model$resample %>% mutate(Model = "")


colnames(results_gaussprRadial) = c("RMSE","R2","MAE","Sigma","Model")
results_gaussprRadial %>% 
  select(RMSE,R2,MAE,Model) %>% 
  gather(a,b,-Model) %>% 
  ggplot(aes(Model,color = Model,fill = Model)) + 
  geom_Boxplot(alpha = 0.3,show.legend = FALSE) + 
  facet_wrap(~ a,scales = "free") + 
  labs(x = "Train simulation performances",y = NULL)


QObs_Tr = Fitting_model$pred$obs # Observed data to be trained
Qsim = Fitting_model$pred$pred   # Simulated data in the training process
Residual_Tr = QObs_Tr - Qsim 
train_obs_calc <- data.frame(QObs_Tr,Qsim,Residual_Tr)

ggplot(data = train_obs_calc,aes(x = QObs_Tr,y = Qsim)) +
  geom_point() +
  geom_abline(intercept = 0,slope = 1,color = "black") +
  labs(title = "Error Scatterplot",x = "Q(m³/s)",y = "Q(m³/s)")

我们已经完成了常见的培训和测试流程，并获得了测试结果。

在第二个窗口代码中，我们将从预测对象中提取一些测试结果。顺便说一下，我希望 forecasted 会给我们带来所需的预测间隔，我需要帮助。

################ Tidy testing results
################ Set time period
N.hours <- as.numeric(ymd_h("2016-10-23 10") - ymd_h("2016-10-02 05"))*24 # Set date limits
Time = ymd_h("2016-10-02 05") + hours(0:N.hours)


################ View testing process results
QObs_Test = forecasted$ST4_lead1 # Observed data to be tested
QPred = forecasted$Qpred   # Simulated data in the testing process
Residual_Test = QObs_Test - QPred 
test_obs_calc <- data.frame(Time,QObs_Test,QPred,Residual_Test)

ggplot(data = test_obs_calc,aes(x = Time,y = QObs_Test)) +
  geom_line() +
  geom_point(aes(y = QPred),color = "coral") +
  labs(title = "Testing Results",x = "",y = "Q(m³/s)")

它为我们带来了测试结果图，如下所示。实线为观测数据，红点为预测估计值。

最后，在第三个窗口代码中，我们只计算回归的指标：

################ Metrics performance
################ Train
gof(QObs_Tr,Qsim)

################ Test
gof(QObs_Test,QPred)

如我们所见，对于训练过程，R² = 0.76 和 MAE = 1.06；而对于测试对应物，R² = 0.80 和 MAE = 1.36。这些对我们来说已经足够了。

在 GPR 中，预测 invervals 基于每个预测点周围的高斯分布，如 Rasmussen and Williams (2006) 提供的那样。我曾通过其他三种方式尝试 GPR：仅使用 kernlab::train( )、tgp::bgp( ) 和 gplmr::train( )。 duckmayr 慷慨地建议了最后一个。对于所有这些，我找到了预测区间。然而，“caret + kernlab::train( )”（如我的代码所示）的联合使用返回了更好的准时预测。我不知道为什么。

正如我在第一个和第二个窗口代码之间提到的，我可能期望从“预测”对象中，我们可以从每个预测点中提取平均数据，然后加（或减）其标准偏差的 2 倍，允许 y = 平均值 +- 2*sd。但是，我不知道如何提取这个平均值，甚至每个预测点的标准差。

简而言之，如何使用 caret + kernlab 包的联合使用来获得高斯过程回归的预测区间？

解决方法

您是否尝试输入“forecasted$”并查看这些平均值和标准差是否出现？

gaussian-process kernlab machine-learning r r r-caret

时间序列预测的高斯过程回归不确定性估计

问题描述

解决方法

相关问答