fcalib | R Documentation |
Defines an objective function for performing blackbox optimization towards solving a modularized calibration of large computer model simulation to field data
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)
u |
a vector of length |
XU |
a |
Z |
a vector of responses/dependent values with |
X |
a |
Y |
a vector of values with |
da |
for emulating |
d |
for the discrepancy between emulations |
g |
for the nugget in the GP model for the discrepancy between emulation
|
uprior |
an optional function taking |
methods |
a sequence of local search methods to be deployed when emulating
|
M |
a computer model “simulation” function taking two matrices as
inputs, to be used in lieu of emulation; see |
bias |
a scalar logical indicating whether a GP discrepancy or bias term should
be estimated via |
omp.threads |
a scalar positive integer indicating the number
of threads to use for SMP parallel processing; see |
save.global |
an environment, e.g., |
verb |
a non-negative integer specifying the verbosity level; |
Gramacy, et al. (2015) 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 returns 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 to the size of the field data set X
and Y
,
necessitating 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
likelihood or posterior calculation serves as a measure of merit for
the value u
.
The fcalib
function is deterministic but, as Gramacy, et al. (2015)
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 R wrapper around 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 exercise 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 exercise. See demo("calib")
or vignette("laGP")
for more details
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
Note that in principle a separable correlation function could be used
(e.g, via newGPsep
and mleGPsep
),
however this is not implemented at this time
Robert B. Gramacy rbg@vt.edu
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. (See Chapter 8.) https://bobby.gramacy.com/surrogates/
R.B. Gramacy (2016). laGP: Large-Scale Spatial Modeling via
Local Approximate Gaussian Processes in R., Journal of Statistical
Software, 72(1), 1-46; \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v072.i01")}
or see vignette("laGP")
R.B. Gramacy, D. Bingham, JP. Holloway, M.J. Grosskopf, C.C. Kuranz, E. Rutter, M. Trantham, and P.R. Drake (2015). Calibrating a large computer experiment simulating radiative shock hydrodynamics. Annals of Applied Statistics, 9(3) 1141-1168; preprint on arXiv:1410.3293 https://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. R package version 0.15-18.
vignette("laGP")
,
jmleGP
, newGP
, aGP.seq
, discrep.est
,
snomadr
## 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 probabilities;
## the discrep.est documentation repeats this with by first calling
## aGP.seq and then discrep.est. The answers should be identical, 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 function 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/discrepancy 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.