simulate.gp: Simulation of Paths from a 'gp' Object

In kergp: Gaussian Process Laboratory

Description

Simulation of paths from a gp object.

Usage

## S3 method for class 'gp'
simulate(object, nsim = 1L, seed = NULL,
         newdata = NULL,
         cond = TRUE,
         trendKnown = FALSE,
         newVarNoise = NULL,
         nuggetSim = 1e-8,
         checkNames = TRUE,
         output = c("list", "matrix"),
         label = "y", unit = "",
         ...)

Arguments

object	An object with class "gp".
nsim	Number of paths wanted.
seed	Not used yet.
newdata	A data frame containing the inputs values used for simulation as well as the required trend covariates, if any. This is similar to the newdata formal in predict.gp.
cond	Logical. Should the simulations be conditional on the observations used in the object or not?
trendKnown	Logical. If TRUE the vector of trend coefficients will be regarded as known so all simulated paths share the same trend. When FALSE, the trend must have been estimated so that its estimation covariance is known. Then each path will have a different trend.
newVarNoise	Variance of the noise for the "new" simulated observations. For the default NULL, the noise variance found in object is used. Note that if a very small positive value is used, each simulated path is the sum of the trend the smooth GP part and an interval containing say 95% of the simulated responses can be regarded as a confidence interval rather than a prediction interval.
nuggetSim	Small positive number ("nugget") added to the diagonal of conditional covariance matrices before computing a Cholesky decomposition, for numerical lack of positive-definiteness. This may happen when the covariance kernel is not (either theoretically or numerically) positive definite.
checkNames	Logical. It TRUE the colnames of X and the input names of the covariance in object as given by inputNames(object) must be identical sets.
output	The type of output wanted. A simple matrix as in standard simulation methods may be quite poor, since interesting intermediate results are then lost.
label, unit	A label and unit that will be copied into the output object when output is "list".
...	Further arguments to be passed to the simulate method of the "covAll" class.

Value

A matrix with the simulated paths as its columns or a more complete list with more results. This list which is given the S3 class "simulate.gp" has the following elements.

• X, F, y Inputs, trend covariates and response.

• XNew, FNew New inputs, new trend covariates.

• sim Matrix of simulated paths.

• trend Matrix of simulated trends.

• trendKnown, noise, newVarNoise Values of the formals.

• Call The call.

Note

When betaKnown is FALSE, the trend and the smooth GP parts of a simulation are usually correlated, and their sum will show less dispersion than each of the two components. The covariance of the vector β can be regarded as the posterior distribution corresponding to a non-informative prior, the distribution from which a new path is drawn being the predictive distribution.

Author(s)

Yves Deville

Examples

set.seed(314159)
n <- 40
x <- sort(runif(n))
y <- 2 + 4 * x  + 2 * x^2 + 3 * sin(6 * pi * x ) + 1.0 * rnorm(n)
df <- data.frame(x = x, y = y)

##-------------------------------------------------------------------------
## use a Matern 3/2 covariance. With model #2, the trend is mispecified,
## so the smooth GP part of the model captures a part of the trend.
##-------------------------------------------------------------------------

myKern <- k1Matern3_2
inputNames(myKern) <- "x"
mygp <- list()
mygp[[1]] <- gp(formula = y ~ x + I(x^2) + sin(6 * pi * x), data = df, 
                parCovLower = c(0.01, 0.01), parCovUpper = c(10, 100),
                cov = myKern, estim = TRUE, noise = TRUE)
mygp[[2]] <- gp(formula = y ~ sin(6 * pi * x), data = df, 
                parCovLower = c(0.01, 0.01), parCovUpper = c(10, 100),
                cov = myKern, estim = TRUE, noise = TRUE)

##-------------------------------------------------------------------------
## New data
##-------------------------------------------------------------------------

nNew <- 150
xNew <- seq(from = -0.2, to= 1.2, length.out = nNew)
dfNew <- data.frame(x = xNew)

opar <- par(mfrow = c(2L, 2L))

nsim <- 40
for (i in 1:2) {

    ##--------------------------------------------------------------------
    ## beta known or not, conditional
    ##--------------------------------------------------------------------

    simTU <- simulate(object = mygp[[i]], newdata = dfNew,  nsim = nsim,
                      trendKnown = FALSE)
    plot(simTU, main = "trend unknown, conditional")

    simTK <- simulate(object = mygp[[i]], newdata = dfNew, nsim = nsim,
                      trendKnown = TRUE)
    plot(simTK, main = "trend known, conditional")

    ##--------------------------------------------------------------------
    ## The same but UNconditional
    ##--------------------------------------------------------------------

    simTU <- simulate(object = mygp[[i]], newdata = dfNew,  nsim = nsim,
                     trendKnown = FALSE, cond = FALSE)
    plot(simTU, main = "trend unknown, unconditional")
    simTK <- simulate(object = mygp[[i]], newdata = dfNew, nsim = nsim,
                      trendKnown = TRUE, cond = FALSE)
    plot(simTK, main = "trend known, unconditional")
}

par(opar)

Example output

Loading required package: Rcpp
Loading required package: testthat
Loading required package: nloptr
Warning messages:
1: executing %dopar% sequentially: no parallel backend registered 
2: In nloptr.add.default.options(opts.user = opts, x0 = x0, num_constraints_ineq = num_constraints_ineq,  :
  No termination criterium specified, using default (relative x-tolerance = 1e-04)
Warning message:
In nloptr.add.default.options(opts.user = opts, x0 = x0, num_constraints_ineq = num_constraints_ineq,  :
  No termination criterium specified, using default (relative x-tolerance = 1e-04)
Loading required package: MASS

kergp documentation built on March 18, 2021, 5:06 p.m.