Loading co2 data
In this example, we use a real dataset and apply a customisedcovariance kernel.
For visualisation purposes, we will use only the first five years ofthe sample.
y<- co2[1:60]n<-length(y)input<-5*(1:n)/n
We want to fit a GPR model using the observed data and then makepredictions at a\(1000\) equallyspaced time points:
inputNew<-5*seq(1/n,1,length.out=1000)# input test for prediction
Using exponential covariance kernel
We will first fit a GPR model with a linear trend line for the meanfunction by settingmeanModel='t', and a exponentialcovariance function by settingCov='pow.ex' andgamma=1.
set.seed(789)fit1<-gpr(input=input,response=y,Cov='pow.ex',gamma=1,meanModel='t',nInitCandidates=50,trace=2)
sapply(fit1$hyper, exp)#> pow.ex.v pow.ex.w vv#> 3.432e+00 2.444e+00 2.794e-09
Since the noise variancevv estimate is extremely low,this fitted model basically interpolates the datapoints:
plot(fit1,main="exponential kernel - fitted model")
See below the large uncertainty in predictions for test inputpoints:
pred1<-gprPredict(train = fit1,inputNew=inputNew)plot(pred1,main="exponential kernel - predictions")
These results suggest that the exponential kernel is likely not to bethe best choice. A kernel which takes into account a periodic patternmay be wanted.
Using a customised covariance kernel
To capture the periodic pattern, we use the following covariancefunction:\[\begin{equation}k(t,t')=v \exp(-w (t-t')^2 - u \sin^2(\pi (t-t'))), \quadv,w,u > 0\label{cus-cov}\end{equation}\] (see Williams, C. K., & Rasmussen, C. E.(2006). “Gaussian Processes for Machine Learning”, the MIT press).
This is not a standard covariance kernel and is not included in thepackage. Therefore, we will define it manually.
Defining a customised covariance kernel
The first step is to define the covariance kernel itself. The name ofthe kernel must be constituted by the prefixcov. and 6letters ( e.g.,cov.custom).
Here the C++ functiondistMat is used to calculatedistances\(\exp(w)(t-t')^2\). ThesimilardistMatSq function is used whent andt' are identical. See package documentation for moredetails.
cov.custom<-function(hyper, input,inputNew=NULL){ hyper<-lapply(hyper, exp)if(is.null(inputNew)){ inputNew<-as.matrix(input) }else{ inputNew<-as.matrix(inputNew) } input<-as.matrix(input) A1<-distMat(input=input,inputNew=inputNew,A=as.matrix(hyper$custom.w),power=2) sept<-outer(input[,1], inputNew[,1],"-") A2<- hyper$custom.u*(sin(pi*sept))^2 customCov<- hyper$custom.v*exp(-A1-A2)return(customCov)}
Defining the first derivatives
The second step is to define the first derivative of thelog-likelihood with respect to each of the hyperparameters of thecovariance kernel.
The name of the functions should look likeDloglik.custom.w, where the middle affix is the customdefined name, and the last letter (i.e. w) names the vectorof the hyperparameters.
For details about derivatives, see Shi, J. Q., and Choi, T. (2011),“Gaussian Process Regression Analysis for Functional Data”, CRCPress.
Dloglik.custom.w<-function(hyper, input, AlphaQ){ A1<-distMatSq(input=input,A=as.matrix(exp(hyper$custom.w)),power=2) Dcov<--cov.custom(hyper, input)*A1 out<-0.5*sum(diag(AlphaQ%*%Dcov))return(out)}Dloglik.custom.u<-function(hyper, input, AlphaQ){ sept<-outer(input[,1], input[,1],"-") A2<-exp(hyper$custom.u)*(sin(pi*sept))^2 Dcov<--cov.custom(hyper, input)*A2 out<-0.5*sum(diag(AlphaQ%*%Dcov))return(out)}Dloglik.custom.v<-function(hyper, input, AlphaQ){ Dcov<-cov.custom(hyper, input) out<-0.5*sum(diag(AlphaQ%*%Dcov))return(out)}
Defining the second derivatives
The third step is to define the second derivatives of thelog-likelihood.
D2custom.w<-function(hyper, input, inv.Q, Alpha.Q){ Cov<-cov.custom(hyper, input) A1<-distMatSq(input=input,A=as.matrix(exp(hyper$custom.w)),power=2) D1cov<--Cov*A1 D2cov<--(D1cov+ Cov)*A1 D2c.w<-D2(D1cov, D2cov, inv.Q, Alpha.Q)return(D2c.w)}D2custom.u<-function(hyper, input, inv.Q, Alpha.Q){ Cov<-cov.custom(hyper, input) sept<-outer(input[,1], input[,1],"-") A2<-exp(hyper$custom.u)*(sin(pi*sept))^2 D1cov<-- Cov*A2 D2cov<--(D1cov+ Cov)*A2 D2c.u<-D2(D1cov, D2cov, inv.Q, Alpha.Q)return(D2c.u)}D2custom.v<-function(hyper, input, inv.Q, Alpha.Q){ out<-cov.custom(hyper, input)return(out)}
Defining a diagonal matrix including the signal variance
Finally, we define a function to calculate a diagonal matrixcontanining the variance of the customised covariance function. This isused for making predictions at new input points.
diag.custom<-function(hyper, input){ Qstar<-rep(exp(hyper$custom.v),nrow(input))return(Qstar)}
Fitting the GPR model with the customised kernel and makingpredictions
fit2<-gpr(input=input,response=y,Cov='custom',NewHyper=c('custom.w','custom.u','custom.v'),gamma=2,meanModel='t',nInitCandidates=50,trace=2,useGradient = F)
sapply(fit2$hyper, exp)#> custom.u custom.v custom.w vv#> 0.4299023 15.4127823 0.0009464 0.0677307
Note that now the noise variancevv estimate is nolonger so small. Both the fitted model and the predictions seem muchmore reasonable for this dataset:
plot(fit2,main="customised kernel - fitted model")
#> Predicted values not found, ploting fitted values
pred2<-gprPredict(train = fit2,inputNew=inputNew)plot(pred2,main="customised kernel - predictions")
