# R/calib.R In laGP: Local Approximate Gaussian Process Regression

#### Documented in aGP.seqdiscrep.estfcalib

```#*******************************************************************************
#
# Local Approximate Gaussian Process Regression
# Copyright (C) 2013, The University of Chicago
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# version 2.1 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
#
# Questions? Contact Robert B. Gramacy (rbg@vt.edu)
#
#*******************************************************************************

## approxGP.seq:
##
## fit a sequence of approx GPs, following the method progression specified
## -- assumes 1/1000 nugget

aGP.seq <- function(X, Z, XX, d, methods=rep("alc", 2), M=NULL, ncalib=0,
## aGP=aGP,
...)
{
## base case for not emulating, M function must be defined
if(length(methods) == 1 && methods == FALSE) {
nc <- ncol(X)
return(data.frame(mean=M(XX[,1:(nc-ncalib)], XX[,(nc-ncalib+1):nc]), var=0))
}

M.hat <- NULL
for(m in methods) {
M.hat <- aGP(X=X, Z=Z, XX=XX, method=m, d=d, ...)
if(length(d\$start) == 1) d\$start <- rep(d\$start, nrow(XX))
nna <- is.finite(M.hat\$mle\$d)
d\$start[nna] <- M.hat\$mle\$d[nna]
}
return(M.hat)
}

## discrep.est:
##
## switch function for estimating a GP discrepency, or simply
## using zero-mean IID and estimating variance (bias=FALSE)

discrep.est <- function(X, Y, Yhat, d, g, bias=TRUE, clean=TRUE)
{
## WHEN REPS > 1 COULD USE DUPLICATES TO GET INITIAL NUGGET ESTIMATE

## remove columns from X which have fewer than three
lun <- apply(X, 2, function(x) { length(unique(x)) })
lun2 <- which(lun <= 2)
if(length(lun2) > 0) X <- X[,-lun2]

## calculate discrepency: possibly dealing with replicates
reps <- length(Y)/nrow(X)
if(reps != 1) {
X <- matrix(rep(t(X), reps), byrow=TRUE, ncol=ncol(X))
Ym <- Y - rep(Yhat, reps)
} else
Ym <- Y - Yhat

## estimate discrepency or not
if(bias) {
gp <- newGP(X, Ym, d=d\$start, g=g\$start, dK=TRUE)
cmle <- jmleGP(gp, drange=c(d\$min, d\$max), grange=c(g\$min, g\$max),
dab=d\$ab, gab=g\$ab)
} else { ## case for no discrepency estimated - just variance
gp <- newGP(X, Ym, d=0, g=0)
cmle <- list(d=0, g=0, dits=0, gits=0)
}

## calculate the log likelihood
cmle\$ll <- llikGP(gp, dab=d\$ab, gab=g\$ab)

## clean up
if(clean) deleteGP(gp)
else cmle\$gp <- gp

## done
return(cmle)
}

## fcalib:
##
## function evaluation by an optimizer in order to find u.hat
## for calibration -- returns the negative log likelihood

fcalib <- function(u, XU, Z, X, Y, da, d, g, uprior=NULL,
methods=rep("alc", 2), M=NULL, bias=TRUE,
{
## fit an approxGP emulator to the computer model data
Xu <- cbind(as.matrix(X), matrix(rep(as.numeric(u), nrow(X)),
ncol=length(u), byrow=TRUE))

## loop over methods
M.hat <- aGP.seq(X=XU, Z=Z, XX=Xu, d=da, methods=methods, M=M,

## find mle parameters for discrepency model
cmle <- discrep.est(X, Y, M.hat\$mean, d, g, bias, TRUE)

## possibly add in prior for u
if(!is.null(uprior)) cmle\$ll <- cmle\$ll + uprior(u)

## progress meter
if(verb > 0)
cat("u=(", paste(u, collapse=", "), "), its=", cmle\$dits+cmle\$gits,
" (d,g)=(", cmle\$d, ",", cmle\$g, "), llik=", cmle\$ll, "\n", sep="")

## possibly save to a global variable
if(!is.logical(save.global)) {
if(exists("fcalib.save", envir=save.global)) {
assign("fcalib.save", rbind(get("fcalib.save", save.global),
c(u, cmle\$ll)), save.global)
} else {
assign("fcalib.save", data.frame(u=matrix(u, nrow=1), ll=cmle\$ll),
save.global)
}
} else if(is.logical(save.global) && save.global == TRUE)
stop("when not FALSE, save.global must be an environment")

## done
return(-cmle\$ll)
}
```

## Try the laGP package in your browser

Any scripts or data that you put into this service are public.

laGP documentation built on March 31, 2023, 9:46 p.m.