design_experiment: Design Optimal Experiment using Gaussian Process Optimization
In GADGET: Gaussian Process Approximations for Designing Experiments

Description Usage Arguments Details Value References Examples

This function implements the method by Weaver, et al. (2016) that uses Gaussian process (GP) optimization to estimate an optimal design for a stochastic design criterion. This function can also optimize a deterministic design criterion as well. Validation of the fitted GP models is provided by the statistics described in Bastos and O'Hagan (2009). For sequential physical experiments or computer experiments with an expensive simulator, one can repeatedly use this function to compute the next step's optimal design point. For sequential computer experiments with inexpensive simulators, see sequential_experiment which will automatically run the simulator and continue the sequential design automatically.

design_experiment(
  criterion,
  lower_bound,
  upper_bound,
  stochastic = TRUE,
  init_budget = 10,
  optim_budget = 10,
  gp_options = list(formula = ~1, kernel = "matern5_2", optimizer = "gen"),
  genoud_options = list(pop.size = 1000, max.generations = 100, wait.generations = 10),
  diagnostics = 1,
  max_augment = 10,
  infile = NULL,
  outfile = NULL,
  verbose = TRUE,
  cluster = NULL
)

`criterion`	A function with vector input of length `d` (see details).
`lower_bound`	A vector of length `d`.
`upper_bound`	A vector of length `d`.
`stochastic`	Is the design criterion stochastic or deterministic (see details)?
`init_budget`	An integer defining the size of the initial training dataset and the size of the validation dataset for the GP model.
`optim_budget`	An integer defining the number of GP optimization iterations.
`gp_options`	A list specifying the type of GP model to fit (see `km`).
`genoud_options`	A list specifying the control options to optimizer (see `genoud`).
`diagnostics`	Type of GP diagnostics to perform before optimization occurs. There are currently three options: 0 (none), 1 (automatic) a simple Mahalanobis distance significance test, 2 (user inspected) execution is paused for visual inspection of pivoted-Cholesky residuals and QQ-plots.
`max_augment`	An integer defining the maximum number of design augmentations before terminating GP fitting.
`infile`	Saved evaluations of the DC from previous GADGET run
`outfile`	File to save all DC evaluations to while running
`verbose`	Print extra output during execution?
`cluster`	A `parallel` cluster object.

The design criterion (DC) is a stochastic or deterministic univariate function that measures the quality of a proposed design. GADGET assumes the design criterion must be minimized. For example, instead of maximizing the determinant of the Fisher-information matrix, GADGET would minimize the negative determinant of the Fisher-information matrix. If the DC is stochastic then the GP model is fit with a nugget effect and expected quantile improvement (EQI) is used to perform the GP optimization. The optimal design is taken to be the design that maximizes EQI on the final optimization iteration. If the DC is deterministic then the GP model is fit without a nugget effect and expected improvement (EI) is used to perform the optimization. The optimal design is taken to be the design with smallest observed DC over all evaluation of the DC.

The GADGET represents designs as a d-length vector. The user supplied DC function must translate this vector into the apporiate form for computing the DC. The upper_bound and lower_bound arguments define bounds of each element in the vectorized design.

A list containing the optimal design, diagnostic results, and final GP model.

Weaver, B. P., Williams, B. J., Anderson-Cook, C. M., Higdon, D. M. (2016). Computational enhancements to Bayesian design of experiments using Gaussian processes. Bayesian Analysis, 11(1), 191–213, <doi:10.1214/15-BA945>.

#--- Deterministic D-Optimal Design ---#
# simple linear regression model
# design = c(x1,x2,p) (two point design with weight) 

dc <- function(design) {
fisher_mat <- (1-design[3])*c(1,design[2]) %*% t(c(1,design[2]))
fisher_mat <- fisher_mat + design[3]*c(1,design[1]) %*% t(c(1,design[1]))
return(-log(det(fisher_mat)))
}
#optim_budget is small for demonstration purposes 
my_result <- design_experiment(dc,
                               c(0,0.5,0),
                               c(0.5,1,1),
                               FALSE,
                               optim_budget = 2) 

#- optimal design -#
print(my_result$experiment)
#- optimization report - #
print(my_result$optimization)
#- final gp model -#
print(my_result$gp)