| Title: | Bayesian Treed Gaussian Process Models |
|---|---|
| Description: | Bayesian nonstationary, semiparametric nonlinear regression and design by treed Gaussian processes (GPs) with jumps to the limiting linear model (LLM). Special cases also implemented include Bayesian linear models, CART, treed linear models, stationary separable and isotropic GPs, and GP single-index models. Provides 1-d and 2-d plotting functions (with projection and slice capabilities) and tree drawing, designed for visualization of tgp-class output. Sensitivity analysis and multi-resolution models are supported. Sequential experimental design and adaptive sampling functions are also provided, including ALM, ALC, and expected improvement. The latter supports derivative-free optimization of noisy black-box functions. For details and tutorials, see Gramacy (2007) <doi:10.18637/jss.v019.i09> and Gramacy & Taddy (2010) <doi:10.18637/jss.v033.i06>. |
| Authors: | Robert B. Gramacy [aut, cre], Matt A. Taddy [aut] |
| Maintainer: | Robert B. Gramacy <[email protected]> |
| License: | LGPL |
| Version: | 2.4-23 |
| Built: | 2024-10-04 06:41:25 UTC |
| Source: | CRAN |
A Bayesian nonstationary nonparametric regression and design package implementing an array of models of varying flexibility and complexity.
This package implements Bayesian nonstationary, semiparametric nonlinear regression with “treed Gaussian process models” with jumps to the limiting linear model (LLM). The package contains functions which facilitate inference for seven regression models of varying complexity using Markov chain Monte Carlo (MCMC): linear model, CART (Classification and Regression Tree), treed linear model, Gaussian process (GP), GP with jumps to the LLM, GP single-index models, treed GPs, treed GP LLMs, and treed GP single-index models. R provides an interface to the C/C++ backbone, and a serves as mechanism for graphically visualizing the results of inference and posterior predictive surfaces under the models. A Bayesian Monte Carlo based sensitivity analysis is implemented, and multi-resolution models are also supported. Sequential experimental design and adaptive sampling functions are also provided, including ALM, ALC, and expected improvement. The latter supports derivative-free optimization of noisy black-box functions.
For a fuller overview including a complete list of functions, demos and
vignettes, please use help(package="tgp").
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. (See Chapter 9.) https://bobby.gramacy.com/surrogates/
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/ doi:10.18637/jss.v033.i06
Gramacy, R. B., Lee, H. K. H. (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), pp. 1119-1130. Also available as ArXiv article 0710.4536 https://arxiv.org/abs/0710.4536
Robert B. Gramacy, Heng Lian (2011). Gaussian process single-index models as emulators for computer experiments. Available as ArXiv article 1009.4241 https://arxiv.org/abs/1009.4241
Gramacy, R. B., Lee, H. K. H. (2006). Adaptive design of supercomputer experiments. Available as UCSC Technical Report ams2006-02.
Gramacy, R.B., Samworth, R.J., and King, R. (2007) Importance Tempering. ArXiV article 0707.4242 https://arxiv.org/abs/0707.4242
Gray, G.A., Martinez-Canales, M., Taddy, M.A., Lee, H.K.H., and Gramacy, R.B. (2007) Enhancing Parallel Pattern Search Optimization with a Gaussian Process Oracle, SAND2006-7946C, Proceedings of the NECDC
https://bobby.gramacy.com/r_packages/tgp/
The seven functions described below implement Bayesian regression models of varying complexity: linear model, linear CART, Gaussian process (GP), GP with jumps to the limiting linear model (LLM), treed GP, and treed GP LLM.
blm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", BTE = c(1000, 4000, 3), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, trace = FALSE, verb = 1, ...) btlm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", tree = c(0.5, 2), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, trace = FALSE, verb = 1, ...) bcart(X, Z, XX = NULL, bprior = "bflat", tree = c(0.5, 2), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv=FALSE, sens.p = NULL, trace = FALSE, verb = 1, ...) bgp(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", BTE = c(1000, 4000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...) bgpllm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", gamma=c(10,0.2,0.7), BTE = c(1000, 4000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...) btgp(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", tree = c(0.5, 2), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, linburn = FALSE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...) btgpllm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", tree = c(0.5, 2), gamma=c(10,0.2,0.7), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, linburn = FALSE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...)blm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", BTE = c(1000, 4000, 3), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, trace = FALSE, verb = 1, ...) btlm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", tree = c(0.5, 2), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, trace = FALSE, verb = 1, ...) bcart(X, Z, XX = NULL, bprior = "bflat", tree = c(0.5, 2), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv=FALSE, sens.p = NULL, trace = FALSE, verb = 1, ...) bgp(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", BTE = c(1000, 4000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...) bgpllm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", gamma=c(10,0.2,0.7), BTE = c(1000, 4000, 2), R = 1, m0r1 = TRUE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...) btgp(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", tree = c(0.5, 2), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, linburn = FALSE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...) btgpllm(X, Z, XX = NULL, meanfn = "linear", bprior = "bflat", corr = "expsep", tree = c(0.5, 2), gamma=c(10,0.2,0.7), BTE = c(2000, 7000, 2), R = 1, m0r1 = TRUE, linburn = FALSE, itemps = NULL, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, nu = 1.5, trace = FALSE, verb = 1, ...)
Each of the above functions takes some subset of the following arguments...
X |
|
Z |
Vector of output responses |
XX |
Optional |
meanfn |
A choice of mean function for the process. When
where
|
bprior |
Linear (beta) prior, default is |
tree |
a 2-vector containing the tree process prior parameterization
automatically giving zero probability to trees
with partitions containing less than |
gamma |
Limiting linear model parameters
|
corr |
Gaussian process correlation model. Choose between the isotropic
power exponential family ( |
BTE |
3-vector of Monte-carlo parameters (B)urn in, (T)otal, and (E)very. Predictive samples are saved every E MCMC rounds starting at round B, stopping at T. |
R |
Number of repeats or restarts of |
m0r1 |
If |
linburn |
If |
itemps |
Importance tempering (IT) inverse temperature ladder,
or powers to improve mixing. See |
pred.n |
|
krige |
|
zcov |
If |
Ds2x |
|
improv |
|
sens.p |
Either |
nu |
“beta” functionality: fixed smoothness parameter for
the Matern correlation function; |
trace |
|
verb |
Level of verbosity of R-console print statements: from 0 (none); 1 (default) which shows the “progress meter”; 2 includes an echo of initialization parameters; up to 3 and 4 (max) with more info about successful tree operations |
... |
These ellipses arguments are interpreted as augmentations to the prior specification generated by
You may use these to specify a custom setting of any of default
parameters in the output list |
The functions and their arguments can be categorized by whether or not they use treed partitioning (T), GP models, and jumps to the LLM (or LM)
| blm | LM | Linear Model |
| btlm | T, LM | Treed Linear Model |
| bcart | T | Treed Constant Model |
| bgp | GP | GP Regression |
| bgpllm | GP, LLM | GP with jumps to the LLM |
| btgp | T, GP | treed GP Regression |
| btgpllm | T, GP, LLM | treed GP with jumps to the LLM |
Each function implements a special case of the generic function
tgp which is an interface to C/C++ code for treed Gaussian process
modeling of varying parameterization. Documentation for tgp
has been declared redundant, and has subsequently been removed. To see
how the b* functions use tgp simply examine the
function. In the latest version, with the addition of the ellipses
“...” argument, there is nothing that can be done
with the direct tgp function that cannot also be done with a
b* function
Only functions in the T (tree) category take the tree argument;
GP category functions take the corr argument; and LLM category
functions take the gamma argument. Non-tree class functions omit
the parts output, see below
bcart is the same as btlm except that only the
intercept term in the LM is estimated; the others are zero, thereby
implementing a Bayesian version of the original CART model
The sens.p argument contains a vector of parameters for
sensitivity analysis. It should be NULL unless created by the
sens function. Refer to help(sens) for details.
If itemps =! NULL then importance tempering (IT) is performed
to get better mixing. After each restart (when R > 1) the
observation counts are used to update the pseudo-prior. Stochastic
approximation is performed in the first burn-in rounds (for B-T
rounds, not B) when c0 and n0 are positive.
Every subsequent burn-in after the first restart is for B
rounds in order to settle-in after using the observation counts. See
default.itemps for more details and an example
Please see vignette("tgp") for a detailed illustration
bgp returns an object of class "tgp".
The function plot.tgp
can be used to help visualize results.
An object of class "tgp" is a list containing at least the
following components... The parts output is unique to the T
(tree) category functions. Tree viewing is supported by
tgp.trees
X |
Input argument: |
n |
Number of rows in |
d |
Number of cols in |
Z |
Vector of output responses |
XX |
Input argument: |
nn |
Number of rows in |
BTE |
Input argument: Monte-carlo parameters |
R |
Input argument: restarts |
linburn |
Input argument: initialize MCMC with linear CART |
params |
|
dparams |
Double-representation of model input parameters used by the C-code |
itemps |
|
Zp.mean |
Vector of mean predictive estimates at |
Zp.q1 |
Vector of 5% predictive quantiles at |
Zp.q2 |
Vector of 95% predictive quantiles at |
Zp.q |
Vector of quantile norms |
Zp.s2 |
If input |
Zp.km |
Vector of (expected) kriging means at |
Zp.vark |
Vector of posterior variance for kriging surface (no additive noise) at |
Zp.ks2 |
Vector of (expected) predictive kriging variances at |
ZZ.mean |
Vector of mean predictive estimates at |
ZZ.q1 |
Vector of 5% predictive quantiles at |
ZZ.q2 |
Vector of 95% predictive quantiles at |
ZZ.q |
Vector of quantile norms |
ZZ.s2 |
If input |
ZpZZ.s2 |
If input |
ZZ.km |
Vector of (expected) kriging means at |
ZZ.vark |
Vector of posterior variance for kriging surface (no additive noise) at |
ZZ.ks2 |
Vector of (expected) predictive kriging variances at |
Ds2x |
If argument |
improv |
If argument |
response |
Name of response |
parts |
Internal representation of the regions depicted by partitions of the maximum a' posteriori (MAP) tree |
trees |
|
trace |
If |
ess |
Importance tempering effective sample size (ESS). If
Otherwise the ESS will be lower due to a non-zero coefficient of variation of the calculated importance tempering weights |
sens |
See |
Inputs X, XX, Z containing NaN, NA, or Inf are
discarded with non-fatal warnings
Upon execution, MCMC reports are made every 1,000 rounds to indicate progress
Stationary (non-treed) processes on larger inputs (e.g., X,Z)
of size greater than 500, *might* be slow in execution, especially on
older machines. Once the C code starts executing, it can be interrupted
in the usual way: either via Ctrl-C (Unix-alikes) or pressing the Stop
button in the R-GUI. When this happens, interrupt messages will
indicate which required cleanup measures completed before returning
control to R.
Whereas most of the tgp models will work reasonably well with
little or no change to the default prior specification, GP's with the
"mrexpsep" correlation imply a very specific relationship between
fine and coarse data, and a careful prior specification is usually
required.
The ranks provided in the second column of the improv field
of a tgp object are based on the expectation of a multivariate
improvement that may or may not be raised to a positive integer power.
They can thus differ significantly from a simple ranking of the first
column of expected univariate improvement values.
Regarding trace=TRUE: Samples from the posterior will be
collected for all parameters in the model. GP parameters are collected
with reference to the locations in XX, resulting
nn=nrow{XX} traces of d,g,s2,tau2, etc. Therefore, it
is recommended that nn is chosen to be a small, representative,
set of input locations. Besides GP parameters, traces are saved for
the tree partitions, areas under the LLM, log posterior (as a function
of tree height), and samples from the posterior predictive
distributions. Note that since some traces are stored in
files, multiple tgp/R sessions in the same working
directory can clobber the trace files of other sessions
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. (See Chapter 9.) https://bobby.gramacy.com/surrogates/
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/ doi:10.18637/jss.v033.i06
Gramacy, R. B., Lee, H. K. H. (2007). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), pp. 1119-1130. Also available as ArXiv article 0710.4536 https://arxiv.org/abs/0710.4536
Gramacy, R. B. and Lee, K.H. (2008). Gaussian Processes and Limiting Linear Models. Computational Statistics and Data Analysis, 53, pp. 123-136. Also available as ArXiv article 0804.4685 https://arxiv.org/abs/0804.4685
Gramacy, R. B., Lee, H. K. H. (2009). Adaptive design and analysis of supercomputer experiments. Technometrics, 51(2), pp. 130-145. Also avaliable on ArXiv article 0805.4359 https://arxiv.org/abs/0805.4359
Robert B. Gramacy, Heng Lian (2011). Gaussian process single-index models as emulators for computer experiments. Available as ArXiv article 1009.4241 https://arxiv.org/abs/1009.4241
Chipman, H., George, E., & McCulloch, R. (1998). Bayesian CART model search (with discussion). Journal of the American Statistical Association, 93, 935–960.
Chipman, H., George, E., & McCulloch, R. (2002). Bayesian treed models. Machine Learning, 48, 303–324.
M. Schonlau and Jones, D.R. and Welch, W.J. (1998). Global versus local search in constrained optimization of computer models. In "New Developments and applications in experimental design", IMS Lecture Notes - Monograph Series 34. 11–25.
https://bobby.gramacy.com/r_packages/tgp/
plot.tgp, tgp.trees,
predict.tgp, sens, default.itemps
## ## Many of the examples below illustrate the above ## function(s) on random data. Thus it can be fun ## (and informative) to run them several times. ## # # simple linear response # # input and predictive data X <- seq(0,1,length=50) XX <- seq(0,1,length=99) Z <- 1 + 2*X + rnorm(length(X),sd=0.25) out <- blm(X=X, Z=Z, XX=XX) # try Linear Model plot(out) # plot the surface # # 1-d Example # # construct some 1-d nonstationary data X <- seq(0,20,length=100) XX <- seq(0,20,length=99) Z <- (sin(pi*X/5) + 0.2*cos(4*pi*X/5)) * (X <= 9.6) lin <- X>9.6; Z[lin] <- -1 + X[lin]/10 Z <- Z + rnorm(length(Z), sd=0.1) out <- btlm(X=X, Z=Z, XX=XX) # try Linear CART plot(out) # plot the surface tgp.trees(out) # plot the MAP trees out <- btgp(X=X, Z=Z, XX=XX) # use a treed GP plot(out) # plot the surface tgp.trees(out) # plot the MAP trees # # 2-d example # (using the isotropic correlation function) # # construct some 2-d nonstationary data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z XX <- exp2d.data$XX # try a GP out <- bgp(X=X, Z=Z, XX=XX, corr="exp") plot(out) # plot the surface # try a treed GP LLM out <- btgpllm(X=X, Z=Z, XX=XX, corr="exp") plot(out) # plot the surface tgp.trees(out) # plot the MAP trees # # Motorcycle Accident Data # # get the data require(MASS) # try a GP out <- bgp(X=mcycle[,1], Z=mcycle[,2]) plot(out) # plot the surface # try a treed GP LLM # best to use the "b0" beta linear prior to capture common # common linear process throughout all regions (using the # ellipses "...") out <- btgpllm(X=mcycle[,1], Z=mcycle[,2], bprior="b0") plot(out) # plot the surface tgp.trees(out) # plot the MAP trees## ## Many of the examples below illustrate the above ## function(s) on random data. Thus it can be fun ## (and informative) to run them several times. ## # # simple linear response # # input and predictive data X <- seq(0,1,length=50) XX <- seq(0,1,length=99) Z <- 1 + 2*X + rnorm(length(X),sd=0.25) out <- blm(X=X, Z=Z, XX=XX) # try Linear Model plot(out) # plot the surface # # 1-d Example # # construct some 1-d nonstationary data X <- seq(0,20,length=100) XX <- seq(0,20,length=99) Z <- (sin(pi*X/5) + 0.2*cos(4*pi*X/5)) * (X <= 9.6) lin <- X>9.6; Z[lin] <- -1 + X[lin]/10 Z <- Z + rnorm(length(Z), sd=0.1) out <- btlm(X=X, Z=Z, XX=XX) # try Linear CART plot(out) # plot the surface tgp.trees(out) # plot the MAP trees out <- btgp(X=X, Z=Z, XX=XX) # use a treed GP plot(out) # plot the surface tgp.trees(out) # plot the MAP trees # # 2-d example # (using the isotropic correlation function) # # construct some 2-d nonstationary data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z XX <- exp2d.data$XX # try a GP out <- bgp(X=X, Z=Z, XX=XX, corr="exp") plot(out) # plot the surface # try a treed GP LLM out <- btgpllm(X=X, Z=Z, XX=XX, corr="exp") plot(out) # plot the surface tgp.trees(out) # plot the MAP trees # # Motorcycle Accident Data # # get the data require(MASS) # try a GP out <- bgp(X=mcycle[,1], Z=mcycle[,2]) plot(out) # plot the surface # try a treed GP LLM # best to use the "b0" beta linear prior to capture common # common linear process throughout all regions (using the # ellipses "...") out <- btgpllm(X=mcycle[,1], Z=mcycle[,2], bprior="b0") plot(out) # plot the surface tgp.trees(out) # plot the MAP trees
Parameterized by the minimum desired inverse temperature, this
function generates a ladder of inverse temperatures k[1:m]
starting at k[1] = 1, with m steps down to the final
temperature k[m] = k.min progressing sigmoidally,
harmonically or geometrically.
The output is in a format convenient for the b* functions
in the tgp package (e.g. btgp), including
stochastic approximation parameters and
for tuning the uniform pseudo-prior output by this function
default.itemps(m = 40, type = c("geometric", "harmonic","sigmoidal"), k.min = 0.1, c0n0 = c(100, 1000), lambda = c("opt", "naive", "st"))default.itemps(m = 40, type = c("geometric", "harmonic","sigmoidal"), k.min = 0.1, c0n0 = c(100, 1000), lambda = c("opt", "naive", "st"))
m |
Number of temperatures in the ladder; |
type |
Choose from amongst two common defaults for simulated tempering and Metropolis-coupled MCMC, i.e., geometric (default) or harmonic, or a sigmoidal ladder (default) that concentrates more inverse temperatures near 1 |
k.min |
Minimum inverse temperature desired |
c0n0 |
Stochastic approximation parameters used to tune
the simulated tempering pseudo-prior ( |
lambda |
Method for combining the importance samplers at each
temperature. Optimal combination (
Setting |
The geometric and harmonic inverse temperature ladders are usually defined
by an index and a parameter
. The geometric ladder is defined by
and the harmonic ladder by
Alternatively, specifying the minimum temperature
in the ladder can be used to
uniquely determine . E.g., for the geometric
ladder
and for the harmonic
In a similar spirit, the sigmoidal ladder is specified by first
situating indices so that
and
under
The remaining are spaced evenly
between and to fill out the ladder
.
For more details, see the Importance tempering paper cited
below and a full demonstration in vignette("tgp2")
The return value is a list which is compatible with the input argument
itemps to the b* functions (e.g. btgp),
containing the following entries:
c0n0 |
A copy of the |
k |
The generated inverse temperature ladder; a vector
with |
pk |
A vector with |
lambda |
IT method, as specified by the input argument |
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R.B., Samworth, R.J., and King, R. (2010) Importance Tempering. ArXiV article 0707.4242 Statistics and Computing, 20(1), pp. 1-7; https://arxiv.org/abs/0707.4242.
For stochastic approximation and simulated tempering (ST):
Geyer, C.~and Thompson, E.~(1995). Annealing Markov chain Monte Carlo with applications to ancestral inference. Journal of the American Statistical Association, 90, 909–920.
For the geometric temperature ladder:
Neal, R.M.~(2001) Annealed importance sampling. Statistics and Computing, 11, 125–129
Justifying geometric and harmonic defaults:
Liu, J.S.~(1002) Monte Carlo Strategies in Scientific Computing. New York: Springer. Chapter 10 (pages 213 & 233)
https://bobby.gramacy.com/r_packages/tgp/
## comparing the different ladders geo <- default.itemps(type="geometric") har <- default.itemps(type="harmonic") sig <- default.itemps(type="sigmoidal") par(mfrow=c(2,1)) matplot(cbind(geo$k, har$k, sig$k), pch=21:23, main="inv-temp ladders", xlab="indx", ylab="itemp") legend("topright", pch=21:23, c("geometric","harmonic","sigmoidal"), col=1:3) matplot(log(cbind(sig$k, geo$k, har$k)), pch=21:23, main="log(inv-temp) ladders", xlab="indx", ylab="itemp") ## Not run: ## using Importance Tempering (IT) to improve mixing ## on the motorcycle accident dataset library(MASS) out.it <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,22000,2), R=3, itemps=default.itemps(), bprior="b0", trace=TRUE, pred.n=FALSE) ## compare to regular tgp w/o IT out.reg <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,22000,2), R=3, bprior="b0", trace=TRUE, pred.n=FALSE) ## compare the heights explored by the three chains: ## REG, combining all temperatures, and IT p <- out.it$trace$post L <- length(p$height) hw <- suppressWarnings(sample(p$height, L, prob=p$wlambda, replace=TRUE)) b <- hist2bar(cbind(out.reg$trace$post$height, p$height, hw)) par(mfrow=c(1,1)) barplot(b, beside=TRUE, xlab="tree height", ylab="counts", col=1:3, main="tree heights encountered") legend("topright", c("reg MCMC", "All Temps", "IT"), fill=1:3) ## End(Not run)## comparing the different ladders geo <- default.itemps(type="geometric") har <- default.itemps(type="harmonic") sig <- default.itemps(type="sigmoidal") par(mfrow=c(2,1)) matplot(cbind(geo$k, har$k, sig$k), pch=21:23, main="inv-temp ladders", xlab="indx", ylab="itemp") legend("topright", pch=21:23, c("geometric","harmonic","sigmoidal"), col=1:3) matplot(log(cbind(sig$k, geo$k, har$k)), pch=21:23, main="log(inv-temp) ladders", xlab="indx", ylab="itemp") ## Not run: ## using Importance Tempering (IT) to improve mixing ## on the motorcycle accident dataset library(MASS) out.it <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,22000,2), R=3, itemps=default.itemps(), bprior="b0", trace=TRUE, pred.n=FALSE) ## compare to regular tgp w/o IT out.reg <- btgpllm(X=mcycle[,1], Z=mcycle[,2], BTE=c(2000,22000,2), R=3, bprior="b0", trace=TRUE, pred.n=FALSE) ## compare the heights explored by the three chains: ## REG, combining all temperatures, and IT p <- out.it$trace$post L <- length(p$height) hw <- suppressWarnings(sample(p$height, L, prob=p$wlambda, replace=TRUE)) b <- hist2bar(cbind(out.reg$trace$post$height, p$height, hw)) par(mfrow=c(1,1)) barplot(b, beside=TRUE, xlab="tree height", ylab="counts", col=1:3, main="tree heights encountered") legend("topright", c("reg MCMC", "All Temps", "IT"), fill=1:3) ## End(Not run)
Create sequential D-Optimal design for a stationary Gaussian process model of fixed parameterization by subsampling from a list of candidates
dopt.gp(nn, X=NULL, Xcand, iter=5000, verb=0)dopt.gp(nn, X=NULL, Xcand, iter=5000, verb=0)
nn |
Number of new points in the design. Must
be less than or equal to the number of candidates contained in
|
X |
|
Xcand |
|
iter |
number of iterations of stochastic accent algorithm,
default |
verb |
positive integer indicating after how many rounds of
stochastic approximation to print each progress statement;
default |
Design is based on a stationary Gaussian process model with stationary isotropic
exponential correlation function with parameterization fixed as a function
of the dimension of the inputs. The algorithm implemented is a simple stochastic
ascent which maximizes det(K)– the covariance matrix constructed
with locations X and a subset of Xcand of size nn.
The selected design is locally optimal
The output is a list which contains the inputs to, and outputs of, the C code
used to find the optimal design. The chosen design locations can be
accessed as list members XX or equivalently Xcand[fi,].
X |
Input argument: |
nn |
Input argument: number new points in the design |
Xcand |
Input argument: |
ncand |
Number of rows in |
fi |
Vector of length |
XX |
|
Inputs X, Xcand containing NaN, NA, Inf are discarded with non-fatal
warnings. If nn > dim(Xcand)[1] then a non-fatal warning is displayed
and execution commences with nn = dim(Xcand)[1]
In the current version there is no progress indicator. You will have to be patient. Creating D-optimal designs is no speedy task
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. (See Chapter 6.) https://bobby.gramacy.com/surrogates/
Chaloner, K. and Verdinelli, I. (1995). Bayesian experimental design: A review. Statist. Sci., 10, (pp. 273–304).
# # 2-d Exponential data # (This example is based on random data. # It might be fun to run it a few times) # # get the data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z Xcand <- exp2d.data$XX # find a treed sequential D-Optimal design # with 10 more points dgp <- dopt.gp(10, X, Xcand) # plot the d-optimally chosen locations # Contrast with locations chosen via # the tgp.design function plot(X, pch=19, xlim=c(-2,6), ylim=c(-2,6)) points(dgp$XX)# # 2-d Exponential data # (This example is based on random data. # It might be fun to run it a few times) # # get the data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z Xcand <- exp2d.data$XX # find a treed sequential D-Optimal design # with 10 more points dgp <- dopt.gp(10, X, Xcand) # plot the d-optimally chosen locations # Contrast with locations chosen via # the tgp.design function plot(X, pch=19, xlim=c(-2,6), ylim=c(-2,6)) points(dgp$XX)
A 2-dimensional data set that can be used to validate non-stationary models.
data(exp2d)data(exp2d)
A data frame with 441 observations on the following 4 variables.
X1Numeric vector describing the first dimension of X inputs
X2Numeric vector describing the second dimension of X inputs
ZNumeric vector describing the response Z(X)+N(0,sd=0.001)
ZtrueNumeric vector describing the true response Z(X),
without noise
The true response is evaluated as
Zero-mean normal noise
with sd=0.001 has been added to the true response
This data is used in the examples of the functions listed below in
the “See Also” section via the exp2d.rand function
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. https://bobby.gramacy.com/surrogates/
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/. doi:10.18637/jss.v033.i06
Gramacy, R. B., Lee, H. K. H. (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), pp. 1119-1130. Also available as ArXiv article 0710.4536 https://arxiv.org/abs/0710.4536
https://bobby.gramacy.com/r_packages/tgp/
exp2d.rand, exp2d.Z,
btgp, and other b* functions
A Random subsample of data(exp2d), or
Latin Hypercube sampled data evaluated with exp2d.Z
exp2d.rand(n1 = 50, n2 = 30, lh = NULL, dopt = 1)exp2d.rand(n1 = 50, n2 = 30, lh = NULL, dopt = 1)
n1 |
Number of samples from the first, interesting, quadrant |
n2 |
Number of samples from the other three, uninteresting, quadrants |
lh |
If |
dopt |
If |
When is.null(lh), data is subsampled without replacement from
data(exp2d). Of the n1 + n2 <= 441
input/response pairs X,Z, there are n1 are taken from the
first quadrant, i.e., where the response is interesting,
and the remaining n2 are taken from the other three
quadrants. The remaining 441 - (n1 + n2) are treated as
predictive locations
Otherwise, when !is.null(lh), Latin Hypercube Sampling
(lhs) is used
If dopt >= 2 then n1*dopt LH candidates are used
for to get a D-optimal subsample of size n1 from the
first (interesting) quadrant. Similarly n2*dopt in the
rest of the un-interesting region.
A total of lh*dopt candidates will be used for sequential D-optimal
subsampling for predictive locations XX in all four
quadrants assuming the already-sampled X locations will
be in the design.
In all three cases, the response is evaluated as
thus creating the outputs Ztrue and ZZtrue.
Zero-mean normal noise with sd=0.001 is added to the
responses Z and ZZ
Output is a list with entries:
X |
2-d |
Z |
Numeric vector describing the responses (with noise) at the
|
Ztrue |
Numeric vector describing the true responses (without
noise) at the |
XX |
2-d |
ZZ |
Numeric vector describing the responses (with noise) at
the |
ZZtrue |
Numeric vector describing the responses (without
noise) at the |
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/ doi:10.18637/jss.v033.i06
Gramacy, R. B., Lee, H. K. H. (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), pp. 1119-1130. Also available as ArXiv article 0710.4536 https://arxiv.org/abs/0710.4536
https://bobby.gramacy.com/r_packages/tgp/
lhs, exp2d, exp2d.Z,
btgp, and other b* functions
## randomly subsampled data ## ------------------------ eds <- exp2d.rand() # higher span = 0.5 required because the data is sparse # and was generated randomly eds.g <- interp.loess(eds$X[,1], eds$X[,2], eds$Z, span=0.5) # perspective plot, and plot of the input (X & XX) locations par(mfrow=c(1,2), bty="n") persp(eds.g, main="loess surface", theta=-30, phi=20, xlab="X[,1]", ylab="X[,2]", zlab="Z") plot(eds$X, main="Randomly Subsampled Inputs") points(eds$XX, pch=19, cex=0.5) ## Latin Hypercube sampled data ## ---------------------------- edlh <- exp2d.rand(lh=c(20, 15, 10, 5)) # higher span = 0.5 required because the data is sparse # and was generated randomly edlh.g <- interp.loess(edlh$X[,1], edlh$X[,2], edlh$Z, span=0.5) # perspective plot, and plot of the input (X & XX) locations par(mfrow=c(1,2), bty="n") persp(edlh.g, main="loess surface", theta=-30, phi=20, xlab="X[,1]", ylab="X[,2]", zlab="Z") plot(edlh$X, main="Latin Hypercube Sampled Inputs") points(edlh$XX, pch=19, cex=0.5) # show the quadrants abline(h=2, col=2, lty=2, lwd=2) abline(v=2, col=2, lty=2, lwd=2) ## Not run: ## D-optimal subsample with a factor of 10 (more) candidates ## --------------------------------------------------------- edlhd <- exp2d.rand(lh=c(20, 15, 10, 5), dopt=10) # higher span = 0.5 required because the data is sparse # and was generated randomly edlhd.g <- interp.loess(edlhd$X[,1], edlhd$X[,2], edlhd$Z, span=0.5) # perspective plot, and plot of the input (X & XX) locations par(mfrow=c(1,2), bty="n") persp(edlhd.g, main="loess surface", theta=-30, phi=20, xlab="X[,1]", ylab="X[,2]", zlab="Z") plot(edlhd$X, main="D-optimally Sampled Inputs") points(edlhd$XX, pch=19, cex=0.5) # show the quadrants abline(h=2, col=2, lty=2, lwd=2) abline(v=2, col=2, lty=2, lwd=2) ## End(Not run)## randomly subsampled data ## ------------------------ eds <- exp2d.rand() # higher span = 0.5 required because the data is sparse # and was generated randomly eds.g <- interp.loess(eds$X[,1], eds$X[,2], eds$Z, span=0.5) # perspective plot, and plot of the input (X & XX) locations par(mfrow=c(1,2), bty="n") persp(eds.g, main="loess surface", theta=-30, phi=20, xlab="X[,1]", ylab="X[,2]", zlab="Z") plot(eds$X, main="Randomly Subsampled Inputs") points(eds$XX, pch=19, cex=0.5) ## Latin Hypercube sampled data ## ---------------------------- edlh <- exp2d.rand(lh=c(20, 15, 10, 5)) # higher span = 0.5 required because the data is sparse # and was generated randomly edlh.g <- interp.loess(edlh$X[,1], edlh$X[,2], edlh$Z, span=0.5) # perspective plot, and plot of the input (X & XX) locations par(mfrow=c(1,2), bty="n") persp(edlh.g, main="loess surface", theta=-30, phi=20, xlab="X[,1]", ylab="X[,2]", zlab="Z") plot(edlh$X, main="Latin Hypercube Sampled Inputs") points(edlh$XX, pch=19, cex=0.5) # show the quadrants abline(h=2, col=2, lty=2, lwd=2) abline(v=2, col=2, lty=2, lwd=2) ## Not run: ## D-optimal subsample with a factor of 10 (more) candidates ## --------------------------------------------------------- edlhd <- exp2d.rand(lh=c(20, 15, 10, 5), dopt=10) # higher span = 0.5 required because the data is sparse # and was generated randomly edlhd.g <- interp.loess(edlhd$X[,1], edlhd$X[,2], edlhd$Z, span=0.5) # perspective plot, and plot of the input (X & XX) locations par(mfrow=c(1,2), bty="n") persp(edlhd.g, main="loess surface", theta=-30, phi=20, xlab="X[,1]", ylab="X[,2]", zlab="Z") plot(edlhd$X, main="D-optimally Sampled Inputs") points(edlhd$XX, pch=19, cex=0.5) # show the quadrants abline(h=2, col=2, lty=2, lwd=2) abline(v=2, col=2, lty=2, lwd=2) ## End(Not run)
Evaluate the functional (mean) response for the 2-d
exponential data (truth) at the X inputs, and randomly
sample noisy Z–values having normal error with standard
deviation provided.
exp2d.Z(X, sd=0.001)exp2d.Z(X, sd=0.001)
X |
Must be a |
sd |
Standard deviation of iid normal noise added to the responses |
The response is evaluated as
thus creating the outputs Z and Ztrue.
Zero-mean normal noise with sd=0.001 is added to the
responses Z and ZZ
Output is a data.frame with columns:
Z |
Numeric vector describing the responses (with noise) at the
|
Ztrue |
Numeric vector describing the true responses (without
noise) at the |
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. https://bobby.gramacy.com/surrogates/
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/ doi:10.18637/jss.v033.i06
Gramacy, R. B., Lee, H. K. H. (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), pp. 1119-1130. Also available as ArXiv article 0710.4536 https://arxiv.org/abs/0710.4536
https://bobby.gramacy.com/r_packages/tgp/
N <- 20 x <- seq(-2,6,length=N) X <- expand.grid(x, x) Zdata <- exp2d.Z(X) persp(x,x,matrix(Zdata$Ztrue, nrow=N), theta=-30, phi=20, main="Z true", xlab="x1", ylab="x2", zlab="Ztrue")N <- 20 x <- seq(-2,6,length=N) X <- expand.grid(x, x) Zdata <- exp2d.Z(X) persp(x,x,matrix(Zdata$Ztrue, nrow=N), theta=-30, phi=20, main="Z true", xlab="x1", ylab="x2", zlab="Ztrue")
Generate X and Y values from the 10-dim “first” Friedman data set used to validate the Multivariate Adaptive Regression Splines (MARS) model, and a variation involving boolean indicators. This test function has three non-linear and interacting variables, along with two linear, and five which are irrelevant. The version with indicators has parts of the response turned on based on the setting of the indicators
friedman.1.data(n = 100) fried.bool(n = 100)friedman.1.data(n = 100) fried.bool(n = 100)
n |
Number of samples desired |
In the original formulation, as implemented by friedman.1.data
the function has 10-dim inputs X are drawn from Unif(0,1), and responses
are where
and
The variation fried.bool uses indicators
. The function also has 10-dim
inputs X with columns distributed as Unif(0,1) and responses
are where
and
where
The indicator I is coded in binary in the output data frame as:
c(0,0,0) for I=1,
c(0,0,1) for I=2,
c(0,1,0) for I=3, and
c(1,0,0) for I=4.
Output is a data.frame with columns
X.1, ..., X.10
|
describing the 10-d randomly sampled inputs |
I.1, ..., I.3
|
boolean version of the indicators provided only
for |
Y |
sample responses (with N(0,1) noise) |
Ytrue |
true responses (without noise) |
An example using the original version of the data
(friedman.1.data) is contained in the first package vignette:
vignette("tgp"). The boolean version fried.bool
is used in second vignette vignette("tgp2")
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/. doi:10.18637/jss.v033.i06
Friedman, J. H. (1991). Multivariate adaptive regression splines. “Annals of Statistics”, 19, No. 1, 1–67.
Gramacy, R. B., Lee, H. K. H. (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), pp. 1119-1130. Also available as ArXiv article 0710.4536 https://arxiv.org/abs/0710.4536
Chipman, H., George, E., & McCulloch, R. (2002). Bayesian treed models. Machine Learning, 48, 303–324.
https://bobby.gramacy.com/r_packages/tgp/
bgpllm, btlm,
blm, bgp, btgpllm, bgp
Use the loess function to interpolate the
two-dimensional x, y, and z data onto a uniform grid. The output
produced is an object directly usable by the plotting functions
persp, image,
and contour, etc.
This function is designed as an alternative to the
interp functions from the akima
library.
interp.loess(x, y, z, gridlen = c(40,40), span = 0.1, ...)interp.loess(x, y, z, gridlen = c(40,40), span = 0.1, ...)
x |
Vector of |
y |
Vector of |
z |
Vector of |
gridlen |
Size of the interpolated grid to be produced in x and y.
The default of |
span |
Kernel span argument to the |
... |
Further arguments to be passed to the
|
Uses expand.grid function to produce a uniform
grid of size gridlen with domain equal to the rectangle implied
by X and Y. Then, a loess a smoother
is fit to the data Z = f(X,Y). Finally,
predict.loess is used to predict onto the grid.
The output is a list compatible with the 2-d plotting functions
persp, image,
and contour, etc.
The list contains...
x |
Vector of with |
y |
Vector of with |
z |
|
As mentioned above, the default span = 0.1 parameter is
significantly smaller that the default loess setting.
This asserts a tacit assumption that
the input is densely packed and that the noise in z's is small.
Such should be the case when the data are output from a tgp regression –
this function was designed specifically for this situation.
For data that is random or sparse, simply choose higher setting,
e.g., the default loess setting of span =
0.75, or a more intermediate setting of span = 0.5 as in the example below
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
https://bobby.gramacy.com/r_packages/tgp/
interp, loess,
persp, image, contour
# random data ed <- exp2d.rand() # higher span = 0.5 required because the data is sparse # and was generated randomly ed.g <- interp.loess(ed$X[,1], ed$X[,2], ed$Z, span=0.5) # perspective plot persp(ed.g)# random data ed <- exp2d.rand() # higher span = 0.5 required because the data is sparse # and was generated randomly ed.g <- interp.loess(ed$X[,1], ed$X[,2], ed$Z, span=0.5) # perspective plot persp(ed.g)
Functions for making barplots summarizing the progress of
importance tempering. The itemps.barplot function can be
used to make a histogram of the inverse temperatures visited
in the trans-temporal Markov chain. The hist2bar function
is useful for making a histogram of integer-valued samples
(e.g., tree heights) encountered in one or several Markov chains
itemps.barplot(obj, main = NULL, xlab = "itemps", ylab = "counts", plot.it = TRUE, ...) hist2bar(x)itemps.barplot(obj, main = NULL, xlab = "itemps", ylab = "counts", plot.it = TRUE, ...) hist2bar(x)
obj |
|
main |
Main plot label to be augmented by |
xlab |
Label for the x-axis |
ylab |
Label for the y-axis |
plot.it |
whether to plot the |
... |
other arguments passed to |
x |
|
itemps.barplot specifically works with the $trace
field of a "tgp"-class object. An error will be produced
if this field is NULL, i.e., if the b* function used
the create the object was not run with the argument trace=TRUE
The hist2bar function can be used on any integer (or discrete)
valued matrix. The columns are interpreted as different realizations
of similar processes for comparison with one another via a histogram.
The histogram is obtained with the barplot command used
with the argument beside=TRUE. See the examples section of
default.itemps
Both functions return a data.frame that can be used
within the barplot function with argument
beside=TRUE
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R.B., Samworth, R.J., and King, R. (2007) Importance Tempering. ArXiv article 0707.4242 https://arxiv.org/abs/0707.4242
https://bobby.gramacy.com/r_packages/tgp/
default.itemps, vignette(tgp2),
barplot
Draw a (random) Latin Hypercube (LH) sample of size n from in
the region outlined by the provided rectangle
lhs(n, rect, shape=NULL, mode=NULL)lhs(n, rect, shape=NULL, mode=NULL)
n |
Size of the LH sample |
rect |
Rectangle describing the domain from which the LH sample
is to be taken. The rectangle should be a |
shape |
Optional vector of shape parameters for the Beta distribution.
Vector of length equal to the dimension of the domain, with elements > 1.
If this is specified, the LH sample is proportional to a joint pdf formed by
independent Beta distributions in each dimension of the domain,
scaled and shifted to have support defined by |
mode |
Optional vector of mode values for the Beta distribution.
Vector of length equal to the dimension of the domain, with elements within
the support defined by |
The output is a matrix with n rows and
nrow(rect) columns. Each of the n rows represents
a sample point.
The domain bounds specified by the rows of rect can
be specified backwards with no change in effect.
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. (See Chapter 4.) https://bobby.gramacy.com/surrogates/
McKay, M. D., W. J. Conover and R. J. Beckman. (1979). A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code, Technometrics 21: (pp. 239–245).
# get and plot a 2-d LH design s1 <- lhs(10, rbind(c(-2,3), c(0.5, 0.8))) plot(s1) # plot a grid to show that there is one sample # in each grid location abline(v=seq(-2,3,length=11), lty=2, col=3) abline(h=seq(0.5,0.8,length=11), lty=2, col=3)# get and plot a 2-d LH design s1 <- lhs(10, rbind(c(-2,3), c(0.5, 0.8))) plot(s1) # plot a grid to show that there is one sample # in each grid location abline(v=seq(-2,3,length=11), lty=2, col=3) abline(h=seq(0.5,0.8,length=11), lty=2, col=3)
Plot the maximum a' posteriori (MAP) tree from a "tgp"-class
object, or add one on top of an existing plot. Like plot.tgp,
projections and slices of trees can be plotted as specified
mapT(out, proj = NULL, slice = NULL, add = FALSE, lwd = 2, ...)mapT(out, proj = NULL, slice = NULL, add = FALSE, lwd = 2, ...)
out |
|
proj |
1-or-2-Vector describing the dimensions to be shown in a
projection. The argument is ignored for 1-d data, i.e., if |
slice |
|
add |
Specify whether the to add partitions to an existing plot
( |
lwd |
Plotting argument specifying the width of the lines used to depict the partitions |
... |
Additional arguments to |
The only output of this function is a beautiful region-representation of the MAP tree.
For examples, see vignette("tgp") and the examples provided
in the documentation for the tgp.design function
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
https://bobby.gramacy.com/r_packages/tgp/
plot.tgp, tgp.trees,
tgp.design, vignette("tgp")
Optimize (minimize) a noisy black-box function (i.e., a function which
may not be differentiable, and may not execute deterministically).
A b* tgp model is used as a surrogate for
adaptive sampling via improvement (and other) statistics. Note that
this function is intended as a skeleton to be tailored as required
for a particular application
optim.step.tgp(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"), ...) optim.ptgpf(start, rect, tgp.obj, method=c("L-BFGS-B", "Nelder-Mead", "BFGS", "CG", "SANN", "optimize"))optim.step.tgp(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"), ...) optim.ptgpf(start, rect, tgp.obj, method=c("L-BFGS-B", "Nelder-Mead", "BFGS", "CG", "SANN", "optimize"))
f |
A function to be optimized, having only one free argument |
rect |
|
model |
The |
prev |
The output from a previous call to |
X |
|
Z |
Vector of current output responses |
NN |
Number of candidate locations ( |
improv |
Indicates the |
cands |
The type of candidates ( |
method |
A method from |
... |
Further arguments to the |
start |
A starting value for optimization of the MAP predictive
(kriging) surface of a |
tgp.obj |
A |
optim.step.tgp executes one step in a search for the global
optimum (minimum) of a noisy function (Z~f(X)) in a bounded
rectangle (rect). The procedure essentially fits a tgp
model and samples from
the posterior distribution of improvement
statistics at NN+1 candidates locations.
NN of the candidates come from
cands, i.e., "lhs" or "tdopt", plus one which
is the location of the minima found in a previous run via
prev by using optim (with a particular
method or optimize instead) on the MAP
model predictive surface using the "tgp"-class object
contained therein.
The improv[2] with the the highest expected improvement are
recommended for adding into the design on output.
optim.ptgpf is the subroutine used by
optim.step.tgp to find optimize on the MAP (surrogate)
predictive surface for the "tgp"-class object contained in
prev.
Please see vignette("tgp2") for a detailed illustration
The list return has the following components.
X |
A |
progress |
A one-row |
obj |
the |
The ellipses (...) argument is used differently here, as
compared to optim, and optimize. It
allows further arguments to be passed to the b* model
function, whereas for optim it would describe
further (static) arguments to the function f to be optimized.
If static arguments need to be set for f, then we recommend
setting defaults via the formals of f
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. (See Chapter 7.) https://bobby.gramacy.com/surrogates/
Matthew Taddy, Herbert K.H. Lee, Genetha A. Gray, and Joshua D. Griffin. (2009) Bayesian guided pattern search for robust local optimization. Technometrics, 51(4), pp. 389-401
https://bobby.gramacy.com/r_packages/tgp/
btgp, etc., optim,
optimize, tgp.design,
predict.tgp, dopt.gp
## optimize the simple exponential function f <- function(x) { exp2d.Z(x)$Z } ## create the initial design with D-optimal candidates rect <- rbind(c(-2,6), c(-2,6)) Xcand <- lhs(500, rect) X <- dopt.gp(50, X=NULL, Xcand)$XX Z <- f(X) ## do 10 rounds of adaptive sampling out <- progress <- NULL for(i in 1:10) { ## get recommendations for the next point to sample out <- optim.step.tgp(f, X=X, Z=Z, rect=rect, prev=out) ## add in the inputs, and newly sampled outputs X <- rbind(X, out$X) Z <- c(Z, f(out$X)) ## keep track of progress and best optimum progress <- rbind(progress, out$progress) print(progress[i,]) } ## plot the progress so far par(mfrow=c(2,2)) plot(out$obj, layout="surf") plot(out$obj, layout="as", as="improv") matplot(progress[,1:nrow(rect)], main="optim results", xlab="rounds", ylab="x[,1:2]", type="l", lwd=2) plot(log(progress$improv), type="l", main="max log improv", xlab="rounds", ylab="max log(improv)")## optimize the simple exponential function f <- function(x) { exp2d.Z(x)$Z } ## create the initial design with D-optimal candidates rect <- rbind(c(-2,6), c(-2,6)) Xcand <- lhs(500, rect) X <- dopt.gp(50, X=NULL, Xcand)$XX Z <- f(X) ## do 10 rounds of adaptive sampling out <- progress <- NULL for(i in 1:10) { ## get recommendations for the next point to sample out <- optim.step.tgp(f, X=X, Z=Z, rect=rect, prev=out) ## add in the inputs, and newly sampled outputs X <- rbind(X, out$X) Z <- c(Z, f(out$X)) ## keep track of progress and best optimum progress <- rbind(progress, out$progress) print(progress[i,]) } ## plot the progress so far par(mfrow=c(2,2)) plot(out$obj, layout="surf") plot(out$obj, layout="as", as="improv") matplot(progress[,1:nrow(rect)], main="optim results", xlab="rounds", ylab="x[,1:2]", type="l", lwd=2) plot(log(progress$improv), type="l", main="max log improv", xlab="rounds", ylab="max log(improv)")
Partition data according to the maximum a' posteriori (MAP)
tree contained in a "tgp"-class object.
partition(X, out)partition(X, out)
X |
|
out |
|
Output is a list of data.frames populated with the inputs
X contained in each region of the partition of the MAP tree
in the "tgp"-class object out
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
https://bobby.gramacy.com/r_packages/tgp/
# # 2-d Exponential data # (This example is based on random data. # It might be fun to run it a few times) # # get the data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z Xcand <- exp2d.data$XX # fit treed GP LLM model to data w/o prediction # basically just to get MAP tree (and plot it) out <- btgpllm(X=X, Z=Z, pred.n=FALSE, BTE=c(2000,3000,2)) tgp.trees(out) # find a treed sequential D-Optimal design # with 10 more points Xcand.parts <- partition(Xcand, out)# # 2-d Exponential data # (This example is based on random data. # It might be fun to run it a few times) # # get the data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z Xcand <- exp2d.data$XX # fit treed GP LLM model to data w/o prediction # basically just to get MAP tree (and plot it) out <- btgpllm(X=X, Z=Z, pred.n=FALSE, BTE=c(2000,3000,2)) tgp.trees(out) # find a treed sequential D-Optimal design # with 10 more points Xcand.parts <- partition(Xcand, out)
A generic function for plotting of "tgp"-class objects.
1-d posterior mean and error plots, 2-d posterior mean and
error image and perspective plots, and 3+-dimensional mean and error
image and perspective plots are supported via projection
and slicing.
## S3 method for class 'tgp' plot(x, pparts = TRUE, proj = NULL, slice = NULL, map = NULL, as = NULL, center = "mean", layout = "both", main = NULL, xlab = NULL, ylab = NULL, zlab = NULL, pc = "pc", gridlen = c(40,40), span = 0.1, pXX = TRUE, legendloc = "topright", maineff = TRUE, mrlayout="both", rankmax = 20, ...)## S3 method for class 'tgp' plot(x, pparts = TRUE, proj = NULL, slice = NULL, map = NULL, as = NULL, center = "mean", layout = "both", main = NULL, xlab = NULL, ylab = NULL, zlab = NULL, pc = "pc", gridlen = c(40,40), span = 0.1, pXX = TRUE, legendloc = "topright", maineff = TRUE, mrlayout="both", rankmax = 20, ...)
x |
|
pparts |
If |
proj |
1-or-2-Vector describing the dimensions to be shown in a
projection. The argument is ignored for 1-d data, i.e., if |
slice |
|
map |
Optional 2-d map (longitude and latitude) from maps to be shown on top of image plots |
center |
Default |
as |
Optional string indicator for plotting of adaptive sampling
statistics: specifying |
layout |
Specify whether to plot the mean predictive surface
( |
main |
Optional |
xlab |
Optional |
ylab |
Optional |
zlab |
Optional |
pc |
Selects perspective-posterior mean and image-error plots
( |
(only valid for 2-d plots)
gridlen |
Number of regular grid points for 2-d slices and
projections in x and y. The default of |
span |
Span for |
pXX |
scalar logical indicating if |
legendloc |
Location of the |
maineff |
Format for the plots of sensitivity analyses produced
with |
mrlayout |
The plot layout for double resolution
tgp objects with |
rankmax |
When |
... |
Extra arguments to 1-d ( |
The only output of this function is beautiful plots
This plotting function is provided with the intention that it
will be used as an aid in the visualization of "tgp"-class
objects. Users are encouraged to use the source code for
this function in order to develop custom plotting functions.
1-d projections for 3-d or higher data are also available
by specifying a 1-d projection argument (e.g. proj=c(1)
for projecting onto the first input variable).
For examples, see vignette("tgp") and the
help files of those functions in "See Also", below
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
https://bobby.gramacy.com/r_packages/tgp/
plot, bgpllm, btlm,
blm, bgp, btgpllm,
predict.tgp,
tgp.trees, mapT, loess, sens
This generic prediction method was designed to obtain samples
from the posterior predictive distribution after the b*
functions have finished. Samples, or kriging mean and variance
estimates, can be obtained from the MAP model encoded in the
"tgp"-class object, or this parameterization can be used
as a jumping-off point in obtaining further samples from the
joint posterior and posterior predictive distributions
## S3 method for class 'tgp' predict(object, XX = NULL, BTE = c(0, 1, 1), R = 1, MAP = TRUE, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, trace = FALSE, verb = 0, ...)## S3 method for class 'tgp' predict(object, XX = NULL, BTE = c(0, 1, 1), R = 1, MAP = TRUE, pred.n = TRUE, krige = TRUE, zcov = FALSE, Ds2x = FALSE, improv = FALSE, sens.p = NULL, trace = FALSE, verb = 0, ...)
object |
|
XX |
Optional |
BTE |
3-vector of Monte-carlo parameters (B)urn in, (T)otal, and
(E)very. Predictive samples are saved every E MCMC rounds starting
at round B, stopping at T. The default |
R |
Number of repeats or restarts of |
MAP |
When |
pred.n |
|
krige |
|
zcov |
If |
Ds2x |
|
improv |
|
sens.p |
Either |
trace |
|
verb |
Level of verbosity of R-console print statements: from 0 (default: none); 1 which shows the “progress meter”; 2 includes an echo of initialization parameters; up to 3 and 4 (max) with more info about successful tree operations |
... |
Ellipses are not used in the current version
of |
While this function was designed with prediction in mind, it is
actually far more general. It allows a continuation of
MCMC sampling where the b* function left off (when
MAP=FALSE) with a possibly new set of predictive locations
XX. The intended use of this function is to obtain quick
kriging-style predictions for a previously-fit MAP estimate
(contained in a "tgp"-class object)
on a new set of predictive locations XX. However,
it can also be used simply to extend the search for an MAP model
when MAP=FALSE, pred.n=FALSE, and XX=NULL
The output is the same, or a subset of, the output produced
by the b* functions, for example see btgp
It is important to note that this function is not a replacement
for supplying XX to the b* functions, which is the only
way to get fully Bayesian samples from the posterior prediction
at new inputs. It is only intended as a post-analysis (diagnostic)
tool.
Inputs XX containing NaN, NA, or Inf are
discarded with non-fatal warnings. Upon execution, MCMC reports are
made every 1,000 rounds to indicate progress.
If XXs are provided which fall outside the range of X
inputs provided to the original b* function, then those will
not be extrapolated properly, due to the way that bounding rectangles
are defined in the original run. For a workaround, supply
out$Xsplit <- rbind(X, XX) before running predict on
out.
See note for btgp or another b* function
regarding the handling and appropriate specification of traces.
The "tgp" class output produced by predict.tgp can
also be used as input to predict.tgp, as well as others (e.g.,
plot.tgp.
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
https://bobby.gramacy.com/r_packages/tgp/
predict, blm, btlm,
bgp, btgp, bgpllm,
btgpllm, plot.tgp
## revisit the Motorcycle data require(MASS) ## fit a btgpllm without predictive sampling (for speed) out <- btgpllm(X=mcycle[,1], Z=mcycle[,2], bprior="b0", pred.n=FALSE) ## nothing to plot here because there is no predictive data ## save the "tgp" class output object for use later and save(out, file="out.Rsave") ## then remove it (for illustrative purposes) out <- NULL ## (now imagine emailing the out.Rsave file to a friend who ## then performs the following in order to use your fitted ## tgp model on his/her own predictive locations) ## load in the "tgp" class object we just saved load("out.Rsave") ## new predictive locations XX <- seq(2.4, 56.7, length=200) ## now obtain kriging estimates from the MAP model out.kp <- predict(out, XX=XX, pred.n=FALSE) plot(out.kp, center="km", as="ks2") ## actually obtain predictive samples from the MAP out.p <- predict(out, XX=XX, pred.n=FALSE, BTE=c(0,1000,1)) plot(out.p) ## use the MAP as a jumping-off point for more sampling out2 <- predict(out, XX, pred.n=FALSE, BTE=c(0,2000,2), MAP=FALSE, verb=1) plot(out2) ## (generally you would not want to remove the file) unlink("out.Rsave")## revisit the Motorcycle data require(MASS) ## fit a btgpllm without predictive sampling (for speed) out <- btgpllm(X=mcycle[,1], Z=mcycle[,2], bprior="b0", pred.n=FALSE) ## nothing to plot here because there is no predictive data ## save the "tgp" class output object for use later and save(out, file="out.Rsave") ## then remove it (for illustrative purposes) out <- NULL ## (now imagine emailing the out.Rsave file to a friend who ## then performs the following in order to use your fitted ## tgp model on his/her own predictive locations) ## load in the "tgp" class object we just saved load("out.Rsave") ## new predictive locations XX <- seq(2.4, 56.7, length=200) ## now obtain kriging estimates from the MAP model out.kp <- predict(out, XX=XX, pred.n=FALSE) plot(out.kp, center="km", as="ks2") ## actually obtain predictive samples from the MAP out.p <- predict(out, XX=XX, pred.n=FALSE, BTE=c(0,1000,1)) plot(out.p) ## use the MAP as a jumping-off point for more sampling out2 <- predict(out, XX, pred.n=FALSE, BTE=c(0,2000,2), MAP=FALSE, verb=1) plot(out2) ## (generally you would not want to remove the file) unlink("out.Rsave")
Fully Bayesian Monte Carlo sensitivity analysis scheme, based upon any of the regression models contained in the tgp package. Random Latin hypercube samples are drawn at each MCMC iteration in order to estimate main effects as well as 1st order and total sensitivity indices.
sens(X, Z, nn.lhs, model = btgp, ngrid = 100, span = 0.3, BTE = c(3000,8000,10), rect = NULL, shape = NULL, mode = NULL, ...)sens(X, Z, nn.lhs, model = btgp, ngrid = 100, span = 0.3, BTE = c(3000,8000,10), rect = NULL, shape = NULL, mode = NULL, ...)
X |
|
Z |
Vector of output responses |
nn.lhs |
Size of each Latin hypercube drawn for
use in the Monte Carlo integration scheme. Total number of locations
for prediction is |
model |
Either the regression model used for prediction, or
|
ngrid |
The number of grid points in each input dimension upon which main effects will be estimated. |
span |
Smoothing parameter for main effects integration:
the fraction of |
BTE |
3-vector of Monte-Carlo parameters (B)urn in, (T)otal, and (E)very. Predictive samples are saved every E MCMC rounds starting at round B, stopping at T |
rect |
Rectangle describing the domain of the uncertainty
distribution with respect to which the sensitivity is to be
determined. This defines the domain from which the LH sample
is to be taken. The rectangle should be a |
shape |
Optional vector of shape parameters for the Beta
distribution. Vector of length equal to the dimension of the domain,
with elements > 1. If specified, the uncertainty distribution
(i.e. the LH sample) is proportional to a joint pdf formed by
independent Beta distributions in each dimension of the domain,
scaled and shifted to have support defined by |
mode |
Optional vector of mode values for the Beta
uncertainty distribution. Vector of length equal to the dimension
of the domain, with elements within the support defined by
|
... |
Extra arguments to the tgp |
Saltelli (2002) describes a Latin Hypercube sampling based method for estimation of the 'Sobal' sensitivity indices:
1st Order for input ,
where is the -th input.
Total Effect for input ,
where is all inputs except for the -th.
All moments are with respect to the appropriate marginals of the
uncertainty distribution – that is, the probability
distribution on the inputs with respect to which sensitivity is being
investigated.
Under this approach, the integrals involved are approximated through
averages over properly chosen samples based on two LH samples
proportional to U. If nn.lhs is the sample size for the
Monte Carlo estimate, this scheme requires nn.lhs*(ncol(X)+2)
function evaluations.
The sens function implements the method for unknown functions
, through prediction via one of the tgp regression
models conditional on an observed set of X locations.
At each MCMC iteration of the tgp model fitting,
the nn.lhs*(ncol(X)+2) locations are drawn randomly from the
LHS scheme and realizations of the sensitivity indices are
calculated. Thus we obtain a posterior sample of the indices,
incorporating variability from both the Monte Carlo estimation and
uncertainty about the function output. Since a subset of the
predictive locations are actually an LHS proportional to the
uncertainty distribution, we can also estimate the main effects
through simple non-parametric regression (a moving average).
Please see vignette("tgp2") for a detailed illustration
The output is a "tgp"-class object. The details for btgp
contain a complete description of this output. The list entry that
is relevance to sensitivity analysis is sens, which itself has entries:
par |
This contains a |
Xgrid |
A |
ZZ.mean |
A |
ZZ.q1 |
A |
ZZ.q2 |
A |
S |
A |
T |
A |
The quality of sensitivity analysis is dependent on the size of
the LH samples used for integral approximation; as with any Monte
Carlo integration scheme, the sample size (nn.lhs) must
increase with the dimensionality of the problem. The total
sensitivity indices are forced non-negative,
and if negative values occur it is necessary to increase
nn.lhs. The plot.tgp function replaces negative values with zero
for illustration.
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
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.D. Morris, A. Kottas, M. Taddy, R. Furfaro, and B. Ganapol. (2009) A statistical framework for the sensitivity analysis of radiative transfer models. IEEE Transactions on Geoscience and Remote Sensing, to appear.
Saltelli, A. (2002) Making best use of model evaluations to compute sensitivity indices. Computer Physics Communications, 145, 280-297.
btgp, plot.tgp,
predict.tgp, lhs
# Take a look at the air quality in New York: Sensitivity of # ozone levels with respect to solar radiation, wind, and # temperature. See help(airquality) for details. X <- airquality[,2:4] Z <- airquality$Ozone # Uncertainty distribution is the default: uniform over range(X) # There is missing data, which is removed automatically by tgp # range(X). s <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, model=btgp, ngrid=100, span=0.3, BTE=c(5000,10000,10))) # plot the results plot(s, layout="sens", ylab="Ozone", main="main effects") # plot only the sensitivity indices plot(s, layout="sens", ylab="Ozone", maineff=FALSE) # plot only the main effects, side by side plot(s, layout="sens", ylab="Ozone", main="", maineff=t(1:3)) # build a 'sens.p' parameter vector for a data-dependent # informative uncertainty distribution. For each variable, # the input distribution will be a scaled Beta with shape=2, # and mode equal to the data mean rect <- t(apply(X, 2, range, na.rm=TRUE)) mode <- apply(X , 2, mean, na.rm=TRUE) shape <- rep(2,3) # plot a sample from the marginal uncertainty distribution. Udraw <- lhs(300, rect=rect, mode=mode, shape=shape) par(mfrow=c(1,3)) for(i in 1:3) hist(Udraw[,i], breaks=15,xlab=names(X)[i]) # build sens.p with the 'sens' function. sens.p <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, model=NULL, ngrid=100, rect=rect, shape=shape, mode=mode)) # Use predict.tgp to quickly analyze with respect to this new # uncertainty distribution without re-running the MCMC, then # plot the results. s.new <- predict(s, BTE=c(1,1000,1), sens.p=sens.p, verb=1) plot(s.new, layout="sens", ylab="Ozone", main="main effects")# Take a look at the air quality in New York: Sensitivity of # ozone levels with respect to solar radiation, wind, and # temperature. See help(airquality) for details. X <- airquality[,2:4] Z <- airquality$Ozone # Uncertainty distribution is the default: uniform over range(X) # There is missing data, which is removed automatically by tgp # range(X). s <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, model=btgp, ngrid=100, span=0.3, BTE=c(5000,10000,10))) # plot the results plot(s, layout="sens", ylab="Ozone", main="main effects") # plot only the sensitivity indices plot(s, layout="sens", ylab="Ozone", maineff=FALSE) # plot only the main effects, side by side plot(s, layout="sens", ylab="Ozone", main="", maineff=t(1:3)) # build a 'sens.p' parameter vector for a data-dependent # informative uncertainty distribution. For each variable, # the input distribution will be a scaled Beta with shape=2, # and mode equal to the data mean rect <- t(apply(X, 2, range, na.rm=TRUE)) mode <- apply(X , 2, mean, na.rm=TRUE) shape <- rep(2,3) # plot a sample from the marginal uncertainty distribution. Udraw <- lhs(300, rect=rect, mode=mode, shape=shape) par(mfrow=c(1,3)) for(i in 1:3) hist(Udraw[,i], breaks=15,xlab=names(X)[i]) # build sens.p with the 'sens' function. sens.p <- suppressWarnings(sens(X=X, Z=Z, nn.lhs=300, model=NULL, ngrid=100, rect=rect, shape=shape, mode=mode)) # Use predict.tgp to quickly analyze with respect to this new # uncertainty distribution without re-running the MCMC, then # plot the results. s.new <- predict(s, BTE=c(1,1000,1), sens.p=sens.p, verb=1) plot(s.new, layout="sens", ylab="Ozone", main="main effects")
Construct a default list of parameters to the b*
functions– the interfaces to treed Gaussian process
modeling
tgp.default.params(d, meanfn = c("linear", "constant"), corr = c("expsep", "exp", "mrexpsep", "matern", "sim", "twovar"), splitmin = 1, basemax = d, ...)tgp.default.params(d, meanfn = c("linear", "constant"), corr = c("expsep", "exp", "mrexpsep", "matern", "sim", "twovar"), splitmin = 1, basemax = d, ...)
d |
number of input dimensions |
meanfn |
A choice of mean function for the process. When
where
|
corr |
Gaussian process correlation model. Choose between the isotropic
power exponential family ( |
splitmin |
Indicates which column of the inputs |
basemax |
Indicates which column of the inputs |
... |
These ellipses arguments are interpreted as augmentations to the prior specification. You may use these to specify a custom setting of any of default parameters in the output list detailed below |
The output is the following list of params...
col |
dimension of regression coefficients |
meanfn |
copied from the inputs |
corr |
copied from the inputs |
bprior |
Linear (beta) prior, default is |
beta |
|
tree |
with zero probability given to trees
with partitions containing less than |
s2.p |
|
tau2.p |
|
d.p |
c(1.0,20.0,10.0,10.0) Mixture of gamma prior parameter (initial values)
for the range parameter(s) |
nug.p |
|
gamma |
|
d.lam |
|
nug.lam |
|
s2.lam |
|
tau2.lam |
|
delta.p |
|
nugf.p |
|
dp.sim |
|
Please refer to the examples for the functions in
"See Also" below, vignette("tgp") and vignette(tgp2)
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/ doi:10.18637/jss.v033.i06
Gramacy, R. B., Lee, H. K. H. (2008). Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103(483), pp. 1119-1130. Also available as ArXiv article 0710.4536 https://arxiv.org/abs/0710.4536
Robert B. Gramacy, Heng Lian (2011). Gaussian process single-index models as emulators for computer experiments. Available as ArXiv article 1009.4241 https://arxiv.org/abs/1009.4241
https://bobby.gramacy.com/r_packages/tgp/
blm, btlm, bgp,
btgp, bgpllm, btgpllm
Based on the maximum a' posteriori (MAP)
treed partition extracted from a "tgp"-class object,
calculate independent sequential treed D-Optimal designs in each of the regions.
tgp.design(howmany, Xcand, out, iter = 5000, verb = 0)tgp.design(howmany, Xcand, out, iter = 5000, verb = 0)
howmany |
Number of new points in the design. Must
be less than the number of candidates contained in
|
Xcand |
|
out |
|
iter |
number of iterations of stochastic accent algorithm,
default |
verb |
positive integer indicating after how many rounds of
stochastic approximation in |
This function partitions Xcand and out$X based on
the MAP tree (obtained on "tgp"-class out with
partition) and calls
dopt.gp in order to obtain a D-optimal design under
independent stationary Gaussian processes models defined in each
region. The aim is to obtain a design where new points from Xcand
are spaced out relative to themselves, and relative to
the existing locations (out$X) in the region.
The number of new points from each region of the partition is
proportional to the number of candidates Xcand in the region.
Output is a list of data.frames containing XX design
points for each region of the MAP tree in out
Input Xcand containing NaN, NA, Inf are discarded with non-fatal
warnings
D-Optimal computation in each region is preceded by a print statement
indicated the number of new locations to be chosen and the number of candidates
in the region. Other than that, there are no other indicators of progress.
You will have to be patient.
Creating treed sequential D-optimal designs is no speedy task. At least it
faster than the non-treed version (see dopt.gp).
The example below is also part of vignette("tgp").
Please see vignette("tgp2") for a similar example based on
optimization using the optim.step.tgp
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
Gramacy, R. B. (2020) Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Boca Raton, Florida: Chapman Hall/CRC. (See Chapter 9.) https://bobby.gramacy.com/surrogates/
Gramacy, R. B. (2007). tgp: An R Package for Bayesian Nonstationary, Semiparametric Nonlinear Regression and Design by Treed Gaussian Process Models. Journal of Statistical Software, 19(9). https://www.jstatsoft.org/v19/i09 doi:10.18637/jss.v019.i09
Robert B. Gramacy, Matthew Taddy (2010). Categorical Inputs, Sensitivity Analysis, Optimization and Importance Tempering with tgp Version 2, an R Package for Treed Gaussian Process Models. Journal of Statistical Software, 33(6), 1–48. https://www.jstatsoft.org/v33/i06/ doi:10.18637/jss.v033.i06
Gramacy, R. B., Lee, H. K. H. (2006). Adaptive design and analysis of supercomputer experiments. Technometrics, 51(2), pp. 130-145. Also avaliable on ArXiv article 0805.4359 https://arxiv.org/abs/0805.4359
Gramacy, R. B., Lee, H. K. H., & Macready, W. (2004). Parameter space exploration with Gaussian process trees. ICML (pp. 353–360). Omnipress & ACM Digital Library.
https://bobby.gramacy.com/r_packages/tgp/
bgpllm, btlm, blm,
bgp, btgpllm, plot.tgp,
dopt.gp, lhs,
partition, optim.step.tgp
# # 2-d Exponential data # (This example is based on random data. # It might be fun to run it a few times) # # get the data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z Xcand <- exp2d.data$XX # fit treed GP LLM model to data w/o prediction # basically just to get MAP tree (and plot it) out <- btgpllm(X=X, Z=Z, pred.n=FALSE, corr="exp") tgp.trees(out) # find a treed sequential D-Optimal design # with 10 more points. It is interesting to # contrast this design with one obtained via # the dopt.gp function XX <- tgp.design(10, Xcand, out) # now fit the model again in order to assess # the predictive surface at those new design points dout <- btgpllm(X=X, Z=Z, XX=XX, corr="exp") plot(dout)# # 2-d Exponential data # (This example is based on random data. # It might be fun to run it a few times) # # get the data exp2d.data <- exp2d.rand() X <- exp2d.data$X; Z <- exp2d.data$Z Xcand <- exp2d.data$XX # fit treed GP LLM model to data w/o prediction # basically just to get MAP tree (and plot it) out <- btgpllm(X=X, Z=Z, pred.n=FALSE, corr="exp") tgp.trees(out) # find a treed sequential D-Optimal design # with 10 more points. It is interesting to # contrast this design with one obtained via # the dopt.gp function XX <- tgp.design(10, Xcand, out) # now fit the model again in order to assess # the predictive surface at those new design points dout <- btgpllm(X=X, Z=Z, XX=XX, corr="exp") plot(dout)
Plot the maximum a' posteriori (MAP) tree as a function of tree height, and show the log posterior probabilities for comparison.
tgp.trees(out, heights = NULL, main = NULL, ...)tgp.trees(out, heights = NULL, main = NULL, ...)
out |
|
heights |
Index vector of length less than |
main |
Optional character string to add to the main title of the plot |
... |
Extra arguments to the |
The maximum a' posteriori (MAP) tree encountered at each height
(in the MCMC chain) is plotted, and the log posterior probabilities
are shown for comparison. The text at the branches in the tree show
the splitting variable and value. The text at the leaves show the
number of input data points (X and Z) that fall
into the region(s) along with an estimate of the variability therein.
The only output of this function is beautiful tree diagrams.
Plotting trees that the maptree library is installed, which itself requires that the combinat library also be installed.
See vignette("tgp") and the examples sections of the functions
under “See Also”, below
Robert B. Gramacy, [email protected], and Matt Taddy, [email protected]
https://bobby.gramacy.com/r_packages/tgp/
bgpllm, btlm, blm,
bgp, btgpllm,
plot.tgp, mapT, vignette("tgp")