GPR - example 2

library(GPFDA)
require(MASS)
# packages required for visualisation:
require(interp)
require(fields)

Simulating data from a GP with 2-dimensional input

We simulate \(10\) independent realisations (surfaces) from a zero-mean GP with a Matern \((\nu=3/2)\) covariance function. Each observed surface has a sample size of \(30 \times 30 = 900\) points on \([0,1]^2\).

set.seed(123)
nrep <- 10
n1 <- 30
n2 <- 30
n <- n1*n2
input1 <- seq(0,1,len=n1)
input2 <- seq(0,1,len=n2)
input <- as.matrix(expand.grid(input1=input1, input2=input2))
hp <- list('matern.v'=log(2),'matern.w'=c(log(20), log(25)),'vv'=log(0.2))
nu <- 1.5
Sigma <- cov.matern(hyper = hp, input = input, nu = nu) + diag(exp(hp$vv), n)
Y <- t(mvrnorm(n=nrep, mu=rep(0,n), Sigma=Sigma))

We now split the dataset into training and test sets, leaving about 80% of the observations for the test set.

idx <- expand.grid(1:n1, 1:n2)
n1test <- floor(n1*0.8)
n2test <- floor(n2*0.8)
idx1 <- sort(sample(1:n1, n1test))
idx2 <- sort(sample(1:n2, n2test))
whichTest <- idx[,1]%in%idx1 & idx[,2]%in%idx2

inputTest <- input[whichTest, ]
Ytest <- Y[whichTest, ]
inputTrain <- input[!whichTest, ]
Ytrain <- Y[!whichTest, ]

Estimation

Estimation of the GPR model is done by:

fit <- gpr(input=inputTrain, response=Ytrain, Cov='matern', trace=4, useGradient=T,
            iter.max=50, nu=nu, nInitCandidates=50)
#> 
#>  --------- Initialising ---------- 
#> iter:  -loglik:     matern.v     matern.w1     matern.w2     vv     
#>   0:     5377.7042:  4.36887  8.17782 0.587245 -0.792392
#>   4:     3531.8449: -0.235568  5.97922  5.31298 -2.41557
#>   8:     3075.2097: 0.698879  2.73899  3.33005 -1.78165
#>  12:     3045.5697: 0.687056  2.99280  3.12893 -1.63410
#>  16:     3044.4250: 0.567100  3.12758  3.26257 -1.62071
#>  20:     3044.4095: 0.557519  3.14699  3.27350 -1.62145
#>  24:     3044.4089: 0.554286  3.14900  3.27584 -1.62167
#> 
#>      optimization finished.

The hyperparameter estimates are:

sapply(fit$hyper, exp)
#> $matern.v
#> [1] 1.741
#> 
#> $matern.w
#> [1] 23.31 26.47
#> 
#> $vv
#> [1] 0.1976

Prediction

Predictions for the test set can then be found:

pred <- gprPredict(train=fit, inputNew=inputTest, noiseFreePred=T)

zlim <- range(c(pred$pred.mean, Ytest))

plotImage(response = Ytest, input = inputTest, realisation = 1, 
            n1 = n1test, n2 = n2test,
            zlim = zlim, main = "observed")


plotImage(response = pred$pred.mean, input = inputTest, realisation = 1, 
            n1 = n1test, n2 = n2test,
            zlim = zlim, main = "prediction")