gp: Gaussian process emulator construction
In dgpsi: Interface to 'dgpsi' for Deep and Linked Gaussian Process Emulations

View source: R/gp.R

gp	R Documentation

Gaussian process emulator construction

Description

This function builds and trains a GP emulator.

Usage

gp(
  X,
  Y,
  name = "sexp",
  lengthscale = rep(0.1, ncol(X)),
  bounds = NULL,
  prior = "ref",
  nugget_est = FALSE,
  nugget = ifelse(nugget_est, 0.01, 1e-08),
  scale_est = TRUE,
  scale = 1,
  training = TRUE,
  verb = TRUE,
  vecchia = FALSE,
  M = 25,
  ord = NULL,
  internal_input_idx = NULL,
  linked_idx = NULL,
  id = NULL
)

Arguments

`X`	a matrix where each row is an input data point and each column is an input dimension.
`Y`	a matrix with only one column and each row being an output data point.
`name`	kernel function to be used. Either `"sexp"` for squared exponential kernel or `"matern2.5"` for Matérn-2.5 kernel. Defaults to `"sexp"`.
`lengthscale`	initial values of lengthscales in the kernel function. It can be a single numeric value or a vector of length `ncol(X)`: if it is a single numeric value, it is assumed that kernel functions across input dimensions share the same lengthscale; if it is a vector, it is assumed that kernel functions across input dimensions have different lengthscales. Defaults to a vector of `0.1`.
`bounds`	the lower and upper bounds of lengthscales in the kernel function. It is a vector of length two where the first element is the lower bound and the second element is the upper bound. The bounds will be applied to all lengthscales in the kernel function. Defaults to `NULL` where no bounds are specified for the lengthscales.
`prior`	prior to be used for Maximum a Posterior for lengthscales and nugget of the GP: gamma prior (`"ga"`), inverse gamma prior (`"inv_ga"`), or jointly robust prior (`"ref"`). Defaults to `"ref"`. See the reference below for the jointly robust prior.
`nugget_est`	a bool indicating if the nugget term is to be estimated: `FALSE`: the nugget term is fixed to `nugget`. `TRUE`: the nugget term will be estimated. Defaults to `FALSE`.
`nugget`	the initial nugget value. If `nugget_est = FALSE`, the assigned value is fixed during the training. Set `nugget` to a small value (e.g., `1e-8`) and the corresponding bool in `nugget_est` to `FALSE` for deterministic computer models where the emulator should interpolate the training data points. Set `nugget` to a larger value and the corresponding bool in `nugget_est` to `TRUE` for stochastic emulation where the computer model outputs are assumed to follow a homogeneous Gaussian distribution. Defaults to `1e-8` if `nugget_est = FALSE` and `0.01` if `nugget_est = TRUE`.
`scale_est`	a bool indicating if the variance is to be estimated: `FALSE`: the variance is fixed to `scale`. `TRUE`: the variance term will be estimated. Defaults to `TRUE`.
`scale`	the initial variance value. If `scale_est = FALSE`, the assigned value is fixed during the training. Defaults to `1`.
`training`	a bool indicating if the initialized GP emulator will be trained. When set to `FALSE`, `gp()` returns an untrained GP emulator, to which one can apply `summary()` to inspect its specification or apply `predict()` to check its emulation performance before the training. Defaults to `TRUE`.
`verb`	a bool indicating if the trace information on GP emulator construction and training will be printed during function execution. Defaults to `TRUE`.
`vecchia`	a bool indicating whether to use Vecchia approximation for large-scale GP emulator construction and prediction. Defaults to `FALSE`. The Vecchia approximation implemented for the GP emulation largely follows Katzfuss et al. (2022). See reference below.
`M`	the size of the conditioning set for the Vecchia approximation in the GP emulator training. Defaults to `25`.
`ord`	an R function that returns the ordering of the input to the GP emulator for the Vecchia approximation. The function must satisfy the following basic rules: the first argument represents the input scaled by the lengthscales. the output of the function is a vector of indices that gives the ordering of the input to the GP emulator. If `ord = NULL`, the default random ordering is used. Defaults to `NULL`.
`internal_input_idx`	The column indices of `X` that are generated by the linked emulators in the preceding layers. Set `internal_input_idx = NULL` if the GP emulator is in the first layer of a system or all columns in `X` are generated by the linked emulators in the preceding layers. Defaults to `NULL`. The argument will be removed in the next release. To set up connections of emulators for linked emulations, please use the updated `lgp()` function instead.
`linked_idx`	Either a vector or a list of vectors: If `linked_idx` is a vector, it gives indices of columns in the pooled output matrix (formed by column-combined outputs of all emulators in the feeding layer) that feed into the GP emulator. The length of the vector shall equal to the length of `internal_input_idx` when `internal_input_idx` is not `NULL`. If the GP emulator is in the first layer of a linked emulator system, the vector gives the column indices of the global input (formed by column-combining all input matrices of emulators in the first layer) that the GP emulator will use. If the GP emulator is to be used in both the first and subsequent layers, one should initially set `linked_idx` to the appropriate values for the situation where the emulator is not in the first layer. Then, use the function `set_linked_idx()` to reset the linking information when the emulator is in the first layer. When the GP emulator is not in the first layer of a linked emulator system, `linked_idx` can be a list that gives the information on connections between the GP emulator and emulators in all preceding layers. The length of the list should equal to the number of layers before the GP emulator. Each element of the list is a vector that gives indices of columns in the pooled output matrix (formed by column-combined outputs of all emulators) in the corresponding layer that feed into the GP emulator. If the GP emulator has no connections to any emulator in a certain layer, set `NULL` in the corresponding position of the list. The order of input dimensions in `X[,internal_input_idx]` should be consistent with `linked_idx`. For example, a GP emulator in the second layer that is fed by the output dimension 1 and 3 of emulators in layer 1 should have `linked_idx = list( c(1,3) )`. In addition, the first and second columns of `X[,internal_input_idx]` should correspond to the output dimensions 1 and 3 from layer 1. Set `linked_idx = NULL` if the GP emulator will not be used for linked emulations. However, if this is no longer the case, one can use `set_linked_idx()` to add linking information to the GP emulator. Defaults to `NULL`. The argument will be removed in the next release. To set up connections of emulators for linked emulations, please use the updated `lgp()` function instead.
`id`	an ID to be assigned to the GP emulator. If an ID is not provided (i.e., `id = NULL`), a UUID (Universally Unique Identifier) will be automatically generated and assigned to the emulator. Default to `NULL`.

Details

See further examples and tutorials at https://mingdeyu.github.io/dgpsi-R/.

Value

An S3 class named gp that contains five slots:

id: A number or character string assigned through the id argument.
data: a list that contains two elements: X and Y which are the training input and output data respectively.
specs: a list that contains seven elements:
1. kernel: the type of the kernel function used. Either "sexp" for squared exponential kernel or "matern2.5" for Matérn-2.5 kernel.
2. lengthscales: a vector of lengthscales in the kernel function.
3. scale: the variance value in the kernel function.
4. nugget: the nugget value in the kernel function.
5. internal_dims: the column indices of X that correspond to the linked emulators in the preceding layers of a linked system. The slot will be removed in the next release.
6. external_dims: the column indices of X that correspond to global inputs to the linked system of emulators. It is shown as FALSE if internal_input_idx = NULL. The slot will be removed in the next release.
7. linked_idx: the value passed to argument linked_idx. It is shown as FALSE if the argument linked_idx is NULL. The slot will be removed in the next release.
8. vecchia: whether the Vecchia approximation is used for the GP emulator training.
9. M: the size of the conditioning set for the Vecchia approximation in the GP emulator training.
constructor_obj: a 'python' object that stores the information of the constructed GP emulator.
container_obj: a 'python' object that stores the information for the linked emulation.
emulator_obj: a 'python' object that stores the information for the predictions from the GP emulator.

The returned gp object can be used by

predict() for GP predictions.
validate() for LOO and OOS validations.
plot() for validation plots.
lgp() for linked (D)GP emulator constructions.
summary() to summarize the trained GP emulator.
write() to save the GP emulator to a .pkl file.
design() for sequential designs.
update() to update the GP emulator with new inputs and outputs.
alm(), mice(), and vigf() to locate next design points.

Note

Any R vector detected in X and Y will be treated as a column vector and automatically converted into a single-column R matrix. Thus, if X is a single data point with multiple dimensions, it must be given as a matrix.

References

Gu, M. (2019). Jointly robust prior for Gaussian stochastic process in emulation, calibration and variable selection. Bayesian Analysis, 14(3), 857-885.
Katzfuss, M., Guinness, J., & Lawrence, E. (2022). Scaled Vecchia approximation for fast computer-model emulation. SIAM/ASA Journal on Uncertainty Quantification, 10(2), 537-554.

Examples

## Not run: 
# load the package and the Python env
library(dgpsi)

# construct a step function
f <- function(x) {
   if (x < 0.5) return(-1)
   if (x >= 0.5) return(1)
  }

# generate training data
X <- seq(0, 1, length = 10)
Y <- sapply(X, f)

# training
m <- gp(X, Y)

# summarizing
summary(m)

# LOO cross validation
m <- validate(m)
plot(m)

# prediction
test_x <- seq(0, 1, length = 200)
m <- predict(m, x = test_x)

# OOS validation
validate_x <- sample(test_x, 10)
validate_y <- sapply(validate_x, f)
plot(m, validate_x, validate_y)

# write and read the constructed emulator
write(m, 'step_gp')
m <- read('step_gp')

## End(Not run)

dgpsi documentation built on April 3, 2025, 9:26 p.m.