simulate.gp {kergp}

Simulation of Paths from a gp Object

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 \widehat{\boldsymbol{\beta}} 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)