mleCRNGP: Gaussian process modeling with correlated noise
In hetGP: Heteroskedastic Gaussian Process Modeling and Design under Replication

mleCRNGP

R Documentation

Gaussian process modeling with correlated noise

Description

Gaussian process regression when seed (or trajectory) information is provided, based on maximum likelihood estimation of the hyperparameters. Trajectory handling involves observing all times for any given seed.

Usage

mleCRNGP(
  X,
  Z,
  T0 = NULL,
  stype = c("none", "XS"),
  lower = NULL,
  upper = NULL,
  known = NULL,
  noiseControl = list(g_bounds = c(sqrt(.Machine$double.eps) * 10, 100), rho_bounds =
    c(0.001, 0.9)),
  init = NULL,
  covtype = c("Gaussian", "Matern5_2", "Matern3_2"),
  maxit = 100,
  eps = sqrt(.Machine$double.eps),
  settings = list(return.Ki = TRUE, factr = 1e+07)
)

Arguments

`X`	matrix of all designs, one per row. The last column is assumed to contain the integer seed value.
`Z`	vector of all observations. If `ts` is provided, the `Z` is a matrix of size `nrow(X) x length(ts)`.
`T0`	optional vector of times (same for all `X`s)
`stype`	structural assumptions, options include: `none`: no structure, regular matrix inversion is used (only when no time is present); When time is present, the Kronecker structure is always used (the alternative is to provide times as an extra variable in `X`) Using the Kronecker structure becomes efficient when the product (nx x ns) x nt becomes large.
`lower, upper`	optional bounds for the `theta` parameter (see `cov_gen` for the exact parameterization). In the multivariate case, it is possible to give vectors for bounds (resp. scalars) for anisotropy (resp. isotropy)
`known`	optional list of known parameters, e.g., `beta0`, `theta`, `g` or `rho`.
`noiseControl`	list with element, `g_bounds`, vector providing minimal and maximal noise to signal ratio; `rho_bounds`, vector providing minimal and maximal correlation between seed values, in [0,1];
`init`	optional list specifying starting values for MLE optimization, with elements: `theta_init` initial value of the theta parameters to be optimized over (default to 10% of the range determined with `lower` and `upper`) `g_init` initial value of the nugget parameter to be optimized over. `rho_init` initial value of the seed correlation parameter.
`covtype`	covariance kernel type, either 'Gaussian', 'Matern5_2' or 'Matern3_2', see `cov_gen`
`maxit`	maximum number of iteration for L-BFGS-B of `optim`
`eps`	jitter used in the inversion of the covariance matrix for numerical stability
`settings`	list with argument `return.Ki`, to include the inverse covariance matrix in the object for further use (e.g., prediction). Arguments `factr` (default to 1e9) and `pgtol` are available to be passed to `control` for L-BFGS-B in `optim`.

Details

The global covariance matrix of the model is parameterized as nu_hat * (Cx + g Id) * Cs = nu_hat * K, with Cx the spatial correlation matrix between unique designs, depending on the family of kernel used (see cov_gen for available choices) and values of lengthscale parameters. Cs is the correlation matrix between seed values, equal to 1 if the seeds are equal, rho otherwise. nu_hat is the plugin estimator of the variance of the process.

Compared to mleHomGP, here the replications have a specific identifier, i.e., the seed.

Value

a list which is given the S3 class "CRNGP", with elements:

theta: maximum likelihood estimate of the lengthscale parameter(s),
g: maximum likelihood estimate of the nugget variance,
rho: maximum likelihood estimate of the seed correlation parameter,
trendtype: either "SK" if beta0 is given, else "OK"
beta0: estimated trend unless given in input,
nu_hat: plugin estimator of the variance,
ll: log-likelihood value,
X0, S0, T0: values for the spatial, seed and time designs
Z, eps, covtype, stype,: values given in input,
call: user call of the function
used_args: list with arguments provided in the call
nit_opt, msg: counts and msg returned by optim
Ki: inverse covariance matrix (not scaled by nu_hat) (if return.Ki is TRUE in settings)
Ct: if time is used, corresponding covariance matrix.
time: time to train the model, in seconds.

Note

This function is experimental at this time and could evolve in the future.

References

Xi Chen, Bruce E Ankenman, and Barry L Nelson. The effects of common random numbers on stochastic kriging metamodels. ACM Transactions on Modeling and Computer Simulation (TOMACS), 22(2):1-20, 2012.

Michael Pearce, Matthias Poloczek, and Juergen Branke. Bayesian simulation optimization with common random numbers. In 2019 Winter Simulation Conference (WSC), pages 3492-3503. IEEE, 2019.

A Fadikar, M Binois, N Collier, A Stevens, KB Toh, J Ozik. Trajectory-oriented optimization of stochastic epidemiological models. arXiv preprint arXiv:2305.03926

Examples

##------------------------------------------------------------
## Example 1: CRN GP modeling on 1d sims
##------------------------------------------------------------
#' set.seed(42)
nx <- 50
ns <- 5
x <- matrix(seq(0,1, length.out = nx), nx)
s <- matrix(seq(1, ns, length.out = ns))
g <- 1e-3
theta <- 0.01
KX <- cov_gen(x, theta = theta)
rho <- 0.3
KS <- matrix(rho, ns, ns)
diag(KS) <- 1
YY <- MASS::mvrnorm(n = 1, mu = rep(0, nx*ns), Sigma = kronecker(KX, KS) + g * diag(nx*ns))
YYmat <- matrix(YY, ns, nx)
matplot(x, t(YYmat), pch = 1, type = "b", lty = 3)

Xgrid <- as.matrix(expand.grid(s, x))
Xgrid <- cbind(Xgrid[,2], Xgrid[,1])
ids <- sample(1:nrow(Xgrid), 20)
X0 <- Xgrid[ids,]
Y0 <-  YY[ids]
points(X0[,1], Y0, pch = 20, col = 1 + ((X0[,2] - 1) %% 6))

model <- mleCRNGP(X0, Y0, known = list(theta = 0.01, g = 1e-3, rho = 0.3))

preds <- predict(model, x = Xgrid, xprime = Xgrid)
matlines(x, t(matrix(preds$mean, ns, nx)), lty = 1)
# prediction on new seed (i.e., average prediction)
xs1 <- cbind(x, ns+1)
predsm <- predict(model, x = xs1)
lines(x, predsm$mean, col = "orange", lwd = 3)
lines(x, predsm$mean + 2 * sqrt(predsm$sd2), col = "orange", lwd = 2, lty = 3)
lines(x, predsm$mean - 2 * sqrt(predsm$sd2), col = "orange", lwd = 2, lty = 3)

# Conditional realizations
sims <- MASS::mvrnorm(n = 1, mu = preds$mean, Sigma = 1/2 * (preds$cov + t(preds$cov)))
plot(Xgrid[,1], sims, col = 1 + ((Xgrid[,2] - 1) %% 6))
points(X0[,1], Y0, pch = 20, col = 1 + ((X0[,2] - 1) %% 6))
## Not run: 
##------------------------------------------------------------
## Example 2: Homoskedastic GP modeling on 2d sims
##------------------------------------------------------------
set.seed(2)
nx <- 31
ns <- 5
d <- 2
x <- as.matrix(expand.grid(seq(0,1, length.out = nx), seq(0,1, length.out = nx)))
s <- matrix(seq(1, ns, length.out = ns))
Xgrid <- as.matrix(expand.grid(seq(1, ns, length.out = ns), seq(0,1, length.out = nx), 
                               seq(0,1, length.out = nx)))
Xgrid <- Xgrid[,c(2, 3, 1)]
g <- 1e-3
theta <- c(0.02, 0.05)
KX <- cov_gen(x, theta = theta)
rho <- 0.33
KS <- matrix(rho, ns, ns)
diag(KS) <- 1
YY <- MASS::mvrnorm(n = 1, mu = rep(0, nx*nx*ns), Sigma = kronecker(KX, KS) + g * diag(nx*nx*ns))
YYmat <- matrix(YY, ns, nx*nx)
filled.contour(matrix(YYmat[1,], nx))
filled.contour(matrix(YYmat[2,], nx))

ids <- sample(1:nrow(Xgrid), 80)
X0 <- Xgrid[ids,]
Y0 <-  YY[ids]

## Uncomment below for For 3D visualisation
# library(rgl)
# plot3d(Xgrid[,1], Xgrid[,2], YY, col = 1 + (Xgrid[,3] - 1) %% 6)
# points3d(X0[,1], X0[,2], Y0, size = 10, col = 1 + ((X0[,3] - 1) %% 6))

model <- mleCRNGP(X0, Y0, know = list(beta0 = 0))

preds <- predict(model, x = Xgrid, xprime = Xgrid)
# surface3d(unique(Xgrid[1:nx^2,1]),unique(Xgrid[,2]), matrix(YY[Xgrid[,3]==1], nx), 
#   front = "lines", back = "lines")
# aspect3d(1, 1, 1)
# surface3d(unique(Xgrid[1:nx^2,1]),unique(Xgrid[,2]), matrix(preds$mean[Xgrid[,3]==1], nx), 
#   front = "lines", back = "lines", col = "red")
plot(preds$mean, YY)

# prediction on new seed (i.e., average prediction)
xs1 <- cbind(x, ns+1)
predsm <- predict(model, x = xs1)
# surface3d(unique(x[,1]), unique(x[,2]), matrix(predsm$mean, nx), col = "orange", 
#   front = "lines", back = "lines")

# Conditional realizations
sims <- MASS::mvrnorm(n = 1, mu = preds$mean, Sigma = 1/2 * (preds$cov + t(preds$cov)))
# plot3d(X0[,1], X0[,2], Y0, size = 10, col = 1 + ((X0[,3] - 1) %% 6))
# surface3d(unique(x[,1]), unique(x[,2]), matrix(sims[Xgrid[,3] == 1], nx), col = 1, 
#   front = "lines", back = "lines")
# surface3d(unique(x[,1]), unique(x[,2]), matrix(sims[Xgrid[,3] == 2], nx), col = 2, 
#   front = "lines", back = "lines")

# Faster alternative for conditional realizations 
# (note: here the design points are part of the simulation points)
Xgrid0 <- unique(Xgrid[, -(d + 1), drop = FALSE])
sims2 <- simul(object = model,Xgrid = Xgrid, ids = ids, nsim = 5, check = TRUE) 

##------------------------------------------------------------
## Example 3: Homoskedastic GP modeling on 1d trajectories (with time)
##------------------------------------------------------------
set.seed(42)
nx <- 11
nt <- 9
ns <- 7
x <- matrix(sort(seq(0,1, length.out = nx)), nx)
s <- matrix(sort(seq(1, ns, length.out = ns)))
t <- matrix(sort(seq(0, 1, length.out = nt)), nt)
covtype <- "Matern5_2"
g <- 1e-3
theta <- c(0.3, 0.5)
KX <- cov_gen(x, theta = theta[1], type = covtype)
KT <- cov_gen(t, theta = theta[2], type = covtype)
rho <- 0.3
KS <- matrix(rho, ns, ns)
diag(KS) <- 1
XST <- as.matrix(expand.grid(x, s, t))

Kmc <- kronecker(chol(KT), kronecker(chol(KS), chol(KX)))
YY <- t(Kmc) %*% rnorm(nrow(Kmc))

ninit <- 50
XS <- as.matrix(expand.grid(x, s))
ids <- sort(sample(1:nrow(XS), ninit))
XST0 <- cbind(XS[ids[rep(1:ninit, each = nt)],], rep(t[,1], times = ninit))
X0 <- XST[which(duplicated(rbind(XST, XST0), fromLast = TRUE)),]
Y0 <-  YY[which(duplicated(rbind(XST, XST0), fromLast = TRUE))]

# tmp <- hetGP:::find_reps(X = X0[,-3], Y0)
model <- mleCRNGP(X = XS[ids,], T0=t, Z = matrix(Y0, ncol = nt), covtype = covtype)

preds <- predict(model, x = XS, xprime = XS)

# compare with regular CRN GP
mref <- mleCRNGP(X = X0[, c(1, 3, 2)], Z = Y0, covtype = covtype)
pref <- predict(mref, x = XST[, c(1, 3, 2)], xprime = XST[, c(1, 3, 2)])

print(model$time) # Use Kronecker structure for time
print(mref$time)

plot(as.vector(preds$mean), YY)
plot(pref$mean, YY) 


## End(Not run)

hetGP documentation built on Oct. 3, 2023, 1:07 a.m.