# R/optim.tgp.R In tgp: Bayesian Treed Gaussian Process Models

#### Documented in optim.ptgpfoptim.step.tgptgp.cands

```#*******************************************************************************
#
# Bayesian Regression and Adaptive Sampling with Gaussian Process Trees
# Copyright (C) 2005, University of California
#
# 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 (rbgramacy@ams.ucsc.edu)
#
#*******************************************************************************

## optim.ptgpf:
##
## find the minima of the MAP predictive (kriging) surface
## encoded in the tgp object, starting at the specified spot
## restrected to the provided rectangle with the specified
## method  -- eventually we should be calclating and using GP
## derivative information

optim.ptgpf <- function(start, rect, tgp.obj,
"CG", "SANN", "optimize"))
{
## check the method argument
method <- match.arg(method)

## ptgpf:
##
## predict at x for the MAP tgp object, to be used by optim
## for finding the minimum of the MAP kriging surface

ptgpf <- function(x, tgp.obj, rect=NULL)
{
## only need to check rectangle when 2-d or more
if(!is.null(rect)) for(i in nrow(rect))
if(x[i] < rect[i,1] || x[i] > rect[i,2]) return(Inf)

## necessary b/c check.matrix doesn't know correct ncol
if(!is.null(rect)) x <- matrix(x, ncol=nrow(rect))

## run predict
out <- predict(tgp.obj, XX=x, pred.n=FALSE)
return(as.vector(out\$ZZ.km))
}

## optimize is for 1-d data only
if(method == "optimize") {
if(nrow(rect) != 1) ## check if optimize method is appropriate
stop("method=\"optimize\" only valid for 1-d functions")
opt <- optimize(ptgpf, interval=rect[1,], tgp.obj=tgp.obj)
return(list(par=opt\$minimum, value=opt\$objective, convergence=1))
}

## otherwise use optim in some way
if(method == "L-BFGS-B") {  ## use the boundary informatoin in rect
opt <- optim(par=start, ptgpf, method=method, tgp.obj=tgp.obj,
rect=rect, lower=rect[,1], upper=rect[,2])

} else { ## otheriwise, apply a method without boundaries
opt <- optim(par=start, ptgpf, method=method, tgp.obj=tgp.obj,
rect=rect)
}

## return
return(opt)
}

## tgp.cands:
##
## create NN candidate locations (XX) either via Latin Hypercube
## sample (LHS), or sequential treed D-optimal design (based on an
## initial LHS

tgp.cands <- function(tgp.obj, NN, cands=c("lhs", "tdopt"), rect, X, verb=0)
{
## check the cands argument
cands <- match.arg(cands)

## return a latin hypercibe sample
if(cands == "lhs") return(lhs(NN, rect))

## return a sequential treed D-optimal sample from initial LHS cands
Xcand <- lhs(10*NN, rect)
if(is.null(tgp.obj)) XX <- dopt.gp(NN, X=X, Xcand, verb=verb)\$XX
else XX <- tgp.design(NN, Xcand, tgp.obj, verb=verb)
XX <- matrix(XX, ncol=ncol(X))
return(XX)
}

## optim.tgp:

## execute one step in a search for the global optimum (minimum) of a
## noisy function (f) bounded in rect with starting (X,Z) data provided:
## fit a tgp model and predict creating NN+{1,2} candidates and
## select the one (or two) which give have the highest expected improv
## statistic.  NN of the candidates come from cands (lhs or tdopt),
## plus one which is the location of the minima found (e.g.,) via calling
## optim (with particular method) on the MAP btgpm predictive surface
## (passed in with prev).  When as != "none" an additional candidate
## is also selected, which has the highest expected alc or alm statistic
## The new X (which may be 1-3 rows) are returned

optim.step.tgp <-
function(f, rect, model=btgp, prev=NULL, X=NULL, Z=NULL,
NN=20*length(rect), improv=c(1,5), cands=c("lhs", "tdopt"),
method=c("L-BFGS-B", "Nelder-Mead", "BFGS", "CG", "SANN", "optimize"),
...)
{
## lhs should verify that the rect makes sense
rect <- matrix(rect, ncol=2)

## XX a predictive grid
XX <- tgp.cands(prev\$obj, NN, cands, rect, X)

## add optim results in as a predictive location
XX <- rbind(XX, as.numeric(prev\$progress[1,1:nrow(rect)]))
Xboth <- rbind(X,XX)

## fit a tgp model
out <- model(X=X, Z=Z, XX=XX, improv=improv, ...)

## find the predicted minimum
m <- which.min(c(out\$Zp.mean, out\$ZZ.mean))
Xm <- Xboth[m,]

## find the optimum with kriging, and record in opt
opt <- optim.ptgpf(Xm, rect, out, method)
opt <- data.frame(matrix(c(opt\$par, opt\$value), nrow=1))
names(opt) <- c(paste("x", 1:nrow(rect), sep=""), "z")

## X & from tgp-improv
ir <- out\$improv[,2]
Ximprov <- matrix(XX[ir <= improv[2],], nrow=sum(ir <= improv[2]))

## assemble the info about the current minimum, and return
as <- data.frame(improv=max(out\$improv[,1]))
r <- list(X=Ximprov, progress=cbind(opt, as), obj=out)
return(r)
}
```

## Try the tgp package in your browser

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

tgp documentation built on Sept. 7, 2020, 5:10 p.m.