GPR - example 1
Simulating data from a GP with unidimensional input
We simulate 30 independent realisations from a zero-mean GP with a covariance function given by the sum of a liner kernel and a squared exponential kernel. Each observed curve has a sample size of 15 time points on [0, 1].
set.seed(123)
nrep <- 30
n <- 15
input <- seq(0, 1, length.out=n)
hp <- list('linear.a'=log(40), 'linear.i'=log(10),
'pow.ex.v'=log(5), 'pow.ex.w'=log(15),
'vv'=log(0.3))
Sigma <- cov.linear(hyper=hp, input=input) +
cov.pow.ex(hyper=hp, input=input, gamma=2) +
diag(exp(hp$vv), n, n)
Y <- t(mvrnorm(n=nrep, mu=rep(0,n), Sigma=Sigma))Estimation
Estimation of the GPR model can be carried out without using gradient:
set.seed(111)
fitNoGrad <- gpr(input=input, response=Y, Cov=c('linear','pow.ex'), gamma=2,
trace=4, nInitCandidates = 1, useGradient = F)
#>
#> --------- Initialising ----------
#> iter: -loglik: linear.a linear.i pow.ex.v pow.ex.w vv
#> 0: 3645.7590: 1.71278 4.17194 -2.38691 0.274907 -2.25353
#> 4: 939.61997: 3.32358 3.67538 -2.95628 -0.622218 0.392815
#> 8: 916.29300: 3.88924 2.30299 -2.79971 -0.321428 0.675116
#> 12: 827.16763: 3.94273 1.79195 -0.630069 3.24654 -0.621476
#> 16: 722.71406: 3.86892 1.97119 1.22332 2.87428 -1.18443
#> 20: 721.99682: 3.78325 2.01881 1.39491 2.76259 -1.20266
#> 24: 721.95187: 3.70787 2.04423 1.39566 2.76476 -1.20247
#> 28: 721.95061: 3.69770 2.05338 1.39932 2.76239 -1.20246
#> 32: 721.95041: 3.69609 2.05761 1.39636 2.76435 -1.20211
#>
#> optimization finished.If one wants to use gradient:
set.seed(111)
fit <- gpr(input=input, response=Y, Cov=c('linear','pow.ex'), gamma=2,
trace=4, nInitCandidates = 1, useGradient = T)
#>
#> --------- Initialising ----------
#> iter: -loglik: linear.a linear.i pow.ex.v pow.ex.w vv
#> 0: 3645.7590: 1.71278 4.17194 -2.38691 0.274907 -2.25353
#> 4: 932.96380: 3.37227 3.80452 -2.81914 -0.379572 0.527768
#> 8: 916.36097: 3.81598 2.03556 -2.68588 -0.119573 0.708268
#> 12: 836.83719: 4.28985 2.63705 0.224109 4.13627 0.193960
#> 16: 724.14027: 3.85944 2.40390 1.12431 2.96404 -1.16536
#> 20: 722.23742: 3.74330 2.29529 1.35200 2.78441 -1.17796
#> 24: 721.95351: 3.70114 2.05413 1.39781 2.75816 -1.19705
#> 28: 721.95041: 3.69634 2.05806 1.39680 2.76402 -1.20216
#>
#> optimization finished.Note the smaller number of iterations needed when the gradient analytical expressions are used in the optimisation.
We can see that the hyperparameter estimates are very accurate despite the fairly small sample size:
sapply(fit$hyper, exp)
#> linear.a linear.i pow.ex.v pow.ex.w vv
#> 40.2995 7.8307 4.0422 15.8635 0.3005The fitted model for the 10th realisation can be seen:
#> Predicted values not found, ploting fitted valuesPrediction
Predictions on a fine grid for the 10th realisation can be obtained as follows:
inputNew <- seq(0, 1, length.out = 1000)
pred1 <- gprPredict(train=fit, inputNew=inputNew, noiseFreePred=T)
plot(pred1, realisation=10)
If one wants to include noise variance in the predictions for the new time points:
