The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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:     5368.8149:  4.36887  8.17782 0.587245 -0.792392
#>   4:     3899.6108: 0.403812  7.29654  3.05398 -1.81042
#>   8:     3128.8965: -0.0888745  4.21497  4.13713 -2.00594
#>  12:     3043.2620: 0.342497  3.45996  3.47264 -1.71053
#>  16:     3039.2156: 0.499630  3.22883  3.36662 -1.59966
#>  20:     3038.2633: 0.526918  3.20605  3.31268 -1.63734
#>  24:     3038.2241: 0.550658  3.17911  3.28203 -1.63260
#> 
#>      optimization finished.

The hyperparameter estimates are:

sapply(fit$hyper, exp)
#> $matern.v
#> [1] 1.734
#> 
#> $matern.w
#> [1] 24.03 26.63
#> 
#> $vv
#> [1] 0.1954

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")

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.