# Simulation of Paths from a `gp` Object Simulation of paths from a `gp` object. ```r ## 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. ## Returns 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 $\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. ## Examples ```r 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) ``` ## Author(s) Yves Deville

simulate.gp function