crit_optim: Criterion optimization

In hetGP: Heteroskedastic Gaussian Process Modeling and Design under Replication

Criterion optimization

Description

Search for the best value of available criterion, possibly using a h-steps lookahead strategy to favor designs with replication

Usage

crit_optim(
  model,
  crit,
  ...,
  h = 2,
  Xcand = NULL,
  control = list(multi.start = 10, maxit = 100),
  seed = NULL,
  ncores = 1
)

Arguments

model	homGP or hetGP model
crit	considered criterion, one of "crit_cSUR", "crit_EI", "crit_ICU", "crit_MCU" and "crit_tMSE". Note that crit_IMSPE has its dedicated method, see IMSPE_optim.
...	additional parameters of the criterion
h	horizon for multi-step ahead framework. The decision is made between: sequential crit search starting by a new design (optimized first) then adding h replicates sequential crit searches starting by 1 to h replicates before adding a new point Use h = 0 for the myopic criterion, i.e., not looking ahead.
Xcand	optional discrete set of candidates (otherwise a maximin LHS is used to initialize continuous search)
control	list in case Xcand == NULL, with elements multi.start, to perform a multi-start optimization based on optim, with maxit iterations each. Also, tol_dist defines the minimum distance to an existing design for a new point to be added, otherwise the closest existing design is chosen. In a similar fashion, tol_dist is the minimum relative change of crit for adding a new design.
seed	optional seed for the generation of LHS designs with maximinSA_LHS
ncores	number of CPU available (> 1 mean parallel TRUE), see mclapply

Details

When looking ahead, the kriging believer heuristic is used, meaning that the non-observed value is replaced by the mean prediction in the update.

Value

list with elements:

• par: best first design,

• value: criterion h-steps ahead starting from adding par,

• path: list of elements list(par, value, new) at each step h

References

M. Binois, J. Huang, R. B. Gramacy, M. Ludkovski (2019), Replication or exploration? Sequential design for stochastic simulation experiments, Technometrics, 61(1), 7-23.
Preprint available on arXiv:1710.03206.

Examples

###############################################################################
## Bi-variate example (myopic version)
###############################################################################

nvar <- 2 

set.seed(42)
ftest <- function(x, coef = 0.1) return(sin(2*pi*sum(x)) + rnorm(1, sd = coef))

n <- 25 # must be a square
xgrid0 <- seq(0.1, 0.9, length.out = sqrt(n))
designs <- as.matrix(expand.grid(xgrid0, xgrid0))
X <- designs[rep(1:n, sample(1:10, n, replace = TRUE)),]
Z <- apply(X, 1, ftest)

model <- mleHomGP(X, Z, lower = rep(0.1, nvar), upper = rep(1, nvar))

ngrid <- 51
xgrid <- seq(0,1, length.out = ngrid)
Xgrid <- as.matrix(expand.grid(xgrid, xgrid))

preds <- predict(x = Xgrid, object =  model)

## Initial plots
contour(x = xgrid,  y = xgrid, z = matrix(preds$mean, ngrid),
        main = "Predicted mean", nlevels = 20)
points(model$X0, col = 'blue', pch = 20)

crit <- "crit_EI"
crit_grid <- apply(Xgrid, 1, crit, model = model)
filled.contour(x = xgrid, y = xgrid, matrix(crit_grid, ngrid),
               nlevels = 20, color.palette = terrain.colors, 
               main = "Initial criterion landscape",
plot.axes = {axis(1); axis(2); points(model$X0, pch = 20)})

## Sequential crit search
nsteps <- 1 # Increase for better results

for(i in 1:nsteps){
  res <- crit_optim(model, crit = crit, h = 0, control = list(multi.start = 50, maxit = 30))
  newX <- res$par
  newZ <- ftest(newX)
  model <- update(object = model, Xnew = newX, Znew = newZ)
}

## Final plots
contour(x = xgrid,  y = xgrid, z = matrix(preds$mean, ngrid),
        main = "Predicted mean", nlevels = 20)
points(model$X0, col = 'blue', pch = 20)

crit_grid <- apply(Xgrid, 1, crit, model = model)
filled.contour(x = xgrid, y = xgrid, matrix(crit_grid, ngrid),
               nlevels = 20, color.palette = terrain.colors, 
               main = "Final criterion landscape",
plot.axes = {axis(1); axis(2); points(model$X0, pch = 20)})

###############################################################################
## Bi-variate example (look-ahead version)
###############################################################################
## Not run:   
nvar <- 2 

set.seed(42)
ftest <- function(x, coef = 0.1) return(sin(2*pi*sum(x)) + rnorm(1, sd = coef))

n <- 25 # must be a square
xgrid0 <- seq(0.1, 0.9, length.out = sqrt(n))
designs <- as.matrix(expand.grid(xgrid0, xgrid0))
X <- designs[rep(1:n, sample(1:10, n, replace = TRUE)),]
Z <- apply(X, 1, ftest)

model <- mleHomGP(X, Z, lower = rep(0.1, nvar), upper = rep(1, nvar))

ngrid <- 51
xgrid <- seq(0,1, length.out = ngrid)
Xgrid <- as.matrix(expand.grid(xgrid, xgrid))

nsteps <- 5 # Increase for more steps
crit <- "crit_EI"

# To use parallel computation (turn off on Windows)
library(parallel)
parallel <- FALSE #TRUE #
if(parallel) ncores <- detectCores() else ncores <- 1

for(i in 1:nsteps){
  res <- crit_optim(model, h = 3, crit = crit, ncores = ncores,
                    control = list(multi.start = 100, maxit = 50))
  
  # If a replicate is selected
  if(!res$path[[1]]$new) print("Add replicate")
  
  newX <- res$par
  newZ <- ftest(newX)
  model <- update(object = model, Xnew = newX, Znew = newZ)
  
  ## Plots 
  preds <- predict(x = Xgrid, object =  model)
  contour(x = xgrid,  y = xgrid, z = matrix(preds$mean, ngrid),
          main = "Predicted mean", nlevels = 20)
  points(model$X0, col = 'blue', pch = 20)
  points(newX, col = "red", pch = 20)
  
  crit_grid <- apply(Xgrid, 1, crit, model = model)
  filled.contour(x = xgrid, y = xgrid, matrix(crit_grid, ngrid),
                 nlevels = 20, color.palette = terrain.colors,
  plot.axes = {axis(1); axis(2); points(model$X0, pch = 20)})
}

## End(Not run)

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