fcalib: Objective function for performing large scale computer model calibration via optimization

Description

Defines an objective function for performing blackbox optimization towards solving a modularized calibration of of large computer model simulation to field data

Usage

fcalib(u, XU, Z, X, Y, da, d, g, uprior = NULL, 
  methods = rep("alc", 2), M = NULL, bias = TRUE, 
  omp.threads = 1, save.global = FALSE, verb=1)

Arguments

a vector of length ncol(XU) - ncol(X) containing a setting of the calibration parameter

a matrix or data.frame containing the full (large) design matrix of input locations to a computer simulator whose final ncol(XU) - ncol(X) columns are contains settings of a calibration or tuning paramete

a vector of responses/dependent values with length(Z) = ncol(XU) of computer model outputs at XU

a matrix or data.frame containing the full (large) design matrix of input locations

a vector values with length(Y) = ncol(X) containing the response from field data observations at X. A Y-vector with length(Y) = k*ncol(X), for positive integer k, can be supplied in w

for emulating Z at XU: a prior or initial setting for the (single/isotropic) lengthscale parameter in a Gaussian correlation function; a (default) NULL value triggers a sensible regularization (prior) and

for the discrepency between emulations Yhat at X, based on Z at XU, and the oputs Y observed at X. Otherwise, same description as da above.

for the nugget in the GP model for the discrepency between emulation Yhat at X, based on Z at XU, and the oputs Y observed at X: a prior or initial setting for the nugget

uprior

an optional function taking u arguments which returns a log prior density value for the calibration parameter.

methods

a sequence of local search methods to be deployed when emulating Z at XU via aGP; see aGP.seq for more details; provide met

a computer model simulation function taking two matrices as inputs, to be used in lieu of emulation; see aGP.seq for mode details

bias

a scalar logical indicating whether a GP discrepency or bias term should be estimated via discrep.est, or only a Gaussian (zero-mean) variance; see discrep.

omp.threads

a scalar positive integer indicating the number of threads to use for SMP parallel processing; see aGP for more details

save.global

an environment, e.g., .GlobalEnv if each evaluation of fcalib, say as called by a wrapper or optimization routine, should be saved. The variable used in that environment will be fcalib.save. Otherwise sav

verb

a positive integer specifying the verbosity level; verb = 0 is quiet, whereas a larger value causes each evaluation to be printed to the screen

Value

Returns a scalar measuring the negative log likelihood or posterior density of the calibration parameter u given the other inputs, for the purpose of optimization over u

Details

Gramacy, et al. (2014) defined an objective function which, when optimized, returns a setting of calibration parameters under a setup akin to the modularized calibration method of Liu, et al., (2009). The fcalib function return a log density (likelihood or posterior probability) value obtained by performing emulation at a set of inputs X augmented with a value of the calibration parameter, u. The emulator is trained on XU and Z, presumed to be very large relative relative to the size of the field data set X and Y, necesitating the use of approximate methods like aGP, via aGP.seq. The emulated values, call them Yhat are fed along with X and Y into the discrep.est function, whose whose likelihood or posterior calculation serves as a measure of merit for the value u, upon which the emulations are best.

The fcalib function is deterministic but, as Gramacy, et al. (2014) described, can result is a rugged objective surface for optimizing, meaning that conventional methods, like those in optim are unlikely to work well. They instead recommend using a blackbox derivative-free method, like NOMAD (Le Digabel, 2011). In our example below we use the implementation in the crs package, which provides an Rwrapper aroudn the underlying C library.

Note that while fcalib automates a call first to aGP.seq and then to discrep.est, it does not return enough information to complete, say, an out-of-sample prediction excercise like the one demonstrated in the discrep.est documentation. Therefore, after fcalib is used in an optimization to find the best setting of the calibration parameter, u, those functions must then be used in post-processing to complete a prediction ecercise. See demo("calib") or vignette("laGP") for more details

References

R.B. Gramacy, D. Bingham, JP. Holloway, M.J. Grosskopf, C.C. Kuranz, E. Rutter, M. Trantham, and P.R. Drake (2014). Calibrating a large computer experiment simulating radiative shock hydrodynamics. Preprint on arXiv:1410.3293 http://arxiv.org/abs/1410.3293

F. Liu, M. Bayarri, and J. Berger (2009). Modularization in Bayesian analysis, with emphasis on analysis of computer models. Bayesian Analysis, 4, 1, 119-150.

S. Le Digabel (2011). Algorithm 909: NOMAD: Nonlinear Optimization with the MADS algorithm. ACM Transactions on Mathematical Software, 37, 4, 44:1-44:15.

J.S. Racine, Z. and Nie (2012). crs: Categorical regression splines. Rpackage version 0.15-18.

Examples

Run this code

## the example here illustrates how fcalib combines aGP.seq and 
## discrep.est functions, duplicating the example in the discrep.est
## documentation file.  it is comprised of snippets from demo("calib"), 
## which contains code from the Calibration Section of vignette("laGP")

## Here we generate calibration data using a true calibration
## parameter, u, and then evaluate log posterior probabities; 
## the discrep.est documentation repeats this with by first calling
## aGP.seq and then discrep.est.  The answers should be idential, however
## note that a call first to example("fcalib") and then 
## example("discrep.est") will generate two random data sets, causing
## the results not to match

## begin data-generation code identical to aGP.seq, discrep.est, fcalib
## example sections and demo("calib")

## M: computer model test functon used in Goh et al, 2013 (Technometrics)
## an elaboration of one from Bastos and O'Hagan, 2009 (Technometrics) 
M <- function(x,u) 
  {
    x <- as.matrix(x)
    u <- as.matrix(u)
    out <- (1-exp(-1/(2*x[,2]))) 
    out <- out * (1000*u[,1]*x[,1]^3+1900*x[,1]^2+2092*x[,1]+60) 
    out <- out / (100*u[,2]*x[,1]^3+500*x[,1]^2+4*x[,1]+20)  
    return(out)
  }
  
## bias: discrepancy function from Goh et al, 2013 
bias <- function(x) 
  {
    x<-as.matrix(x)   
    out<- 2*(10*x[,1]^2+4*x[,2]^2) / (50*x[,1]*x[,2]+10)
    return(out)
  }

## beta.prior: marginal beta prior for u, 
## defaults to a mode at 1/2 in hypercube
beta.prior <- function(u, a=2, b=2, log=TRUE)
{
  if(length(a) == 1) a <- rep(a, length(u))
  else if(length(a) != length(u)) stop("length(a) must be 1 or length(u)")
  if(length(b) == 1) b <- rep(b, length(u))
  else if(length(b) != length(u)) stop("length(b) must be 1 or length(u)")
  if(log) return(sum(dbeta(u, a, b, log=TRUE)))
  else return(prod(dbeta(u, a, b, log=FALSE)))
}

## tgp for LHS sampling
library(tgp)
rect <- matrix(rep(0:1, 4), ncol=2, byrow=2)

## training inputs
ny <- 50; 
X <- lhs(ny, rect[1:2,])    ## computer model train

## true (but unknown) setting of the calibration parameter
## for the computer model
u <- c(0.2, 0.1)
Zu <- M(X, matrix(u, nrow=1)) 

## field data response, biased and replicated
sd <- 0.5
## Y <- computer output + bias + noise
reps <- 2 ## example from paper uses reps <- 10
Y <- rep(Zu,reps) + rep(bias(X),reps) + rnorm(reps*length(Zu), sd=sd) 
## variations: remove the bias or change 2 to 1 to remove replicates

## computer model design
nz <- 10000
XU <- lhs(nz, rect)
nth <- 1 ## number of threads to use in emulation, demo uses 8

## augment with physical model design points 
## with various u settings
XU2 <- matrix(NA, nrow=10*ny, ncol=4)
for(i in 1:10) {
  I <- ((i-1)*ny+1):(ny*i)
  XU2[I,1:2] <- X
}
XU2[,3:4] <- lhs(10*ny, rect[3:4,])
XU <- rbind(XU, XU2)

## evaluate the computer model
Z <- M(XU[,1:2], XU[,3:4])

## flag indicating if estimating bias/discrepency or not
bias.est <- TRUE
## two passes of ALC with MLE calculations for aGP.seq
methods <- rep("alcray", 2) ## demo uses rep("alc", 2)

## set up priors
da <- d <- darg(NULL, XU)
g <- garg(list(mle=TRUE), Y) 

## end identical data generation code

## now calculate log posterior for true calibration parameter 
## value (u).  You could repeat this for an estimate value 
## from demo("calib"), for example u.hat <- c(0.8236673, 0.1406989)

fcalib(u, XU, Z, X, Y, da, d, g, beta.prior, methods, M, bias.est, nth)

Run the code above in your browser using DataLab