fit_one_layer | R Documentation |
Conducts MCMC sampling of hyperparameters for a one layer
GP. Length scale parameter theta
governs
the strength of the correlation and nugget parameter g
governs noise. In Matern covariance, v
governs smoothness.
fit_one_layer( x, y, nmcmc = 10000, verb = TRUE, g_0 = 0.01, theta_0 = 0.1, true_g = NULL, settings = NULL, cov = c("matern", "exp2"), v = 2.5, vecchia = FALSE, m = min(25, length(y) - 1) )
x |
vector or matrix of input locations |
y |
vector of response values |
nmcmc |
number of MCMC iterations |
verb |
logical indicating whether to print iteration progress |
g_0 |
initial value for |
theta_0 |
initial value for |
true_g |
if true nugget is known it may be specified here (set to a small value to make fit deterministic). Note - values that are too small may cause numerical issues in matrix inversions. |
settings |
hyperparameters for proposals and priors (see details) |
cov |
covariance kernel, either Matern or squared exponential
( |
v |
Matern smoothness parameter (only used if |
vecchia |
logical indicating whether to use Vecchia approximation |
m |
size of Vecchia conditioning sets (only used if
|
Utilizes Metropolis Hastings sampling of the length scale and
nugget parameters with proposals and priors controlled by
settings
. When true_g
is set to a specific value, the
nugget is not estimated. When vecchia = TRUE
, all calculations
leverage the Vecchia approximation with specified conditioning set size
m
. Vecchia approximation is only implemented for
cov = "matern"
.
Proposals for g
and theta
follow a uniform sliding window
scheme, e.g.
g_star <- runif(1, l * g_t / u, u * g_t / l)
,
with defaults l = 1
and u = 2
provided in settings
.
To adjust these, set settings = list(l = new_l, u = new_u)
.
Priors on g
and theta
follow Gamma distributions with
shape parameters (alpha
) and rate parameters (beta
)
controlled within the settings
list object. Defaults are
settings$alpha$g <- 1.5
settings$beta$g <- 3.9
settings$alpha$theta <- 1.5
settings$beta$theta <- 3.9 / 1.5
These priors are designed for x
scaled
to [0, 1] and y
scaled to have mean 0 and variance 1. These may
be adjusted using the settings
input.
The output object of class gp
is designed for use with
continue
, trim
, and predict
.
a list of the S3 class gp
or gpvec
with elements:
x
: copy of input matrix
y
: copy of response vector
nmcmc
: number of MCMC iterations
settings
: copy of proposal/prior settings
v
: copy of Matern smoothness parameter (v = 999
indicates cov = "exp2"
)
g
: vector of MCMC samples for g
theta
: vector of MCMC samples for theta
tau2
: vector of MLE estimates for tau2
(scale parameter)
time
: computation time in seconds
Sauer, A, RB Gramacy, and D Higdon. 2020. "Active Learning for Deep Gaussian
Process Surrogates." Technometrics, to appear; arXiv:2012.08015.
Sauer, A, A Cooper, and RB Gramacy. 2022. "Vecchia-approximated Deep Gaussian
Processes for Computer Experiments." pre-print on arXiv:2204.02904
Gramacy, RB. Surrogates: Gaussian Process Modeling, Design, and
Optimization for the Applied Sciences. Chapman Hall, 2020.
# Examples of real-world implementations are available at: # https://bitbucket.org/gramacylab/deepgp-ex/ # G function (https://www.sfu.ca/~ssurjano/gfunc.html) f <- function(xx, a = (c(1:length(xx)) - 1) / 2) { new1 <- abs(4 * xx - 2) + a new2 <- 1 + a prod <- prod(new1 / new2) return((prod - 1) / 0.86) } # Training data d <- 1 n <- 20 x <- matrix(runif(n * d), ncol = d) y <- apply(x, 1, f) # Testing data n_test <- 100 xx <- matrix(runif(n_test * d), ncol = d) yy <- apply(xx, 1, f) plot(xx[order(xx)], yy[order(xx)], type = "l") points(x, y, col = 2) # Example 1: full model (nugget fixed) fit <- fit_one_layer(x, y, nmcmc = 2000, true_g = 1e-6) plot(fit) fit <- trim(fit, 1000, 2) fit <- predict(fit, xx, cores = 1) plot(fit) # Example 2: full model (nugget estimated, EI calculated) fit <- fit_one_layer(x, y, nmcmc = 2000) plot(fit) fit <- trim(fit, 1000, 2) fit <- predict(fit, xx, cores = 1, EI = TRUE) plot(fit) par(new = TRUE) # overlay EI plot(xx[order(xx)], fit$EI[order(xx)], type = 'l', lty = 2, axes = FALSE, xlab = '', ylab = '') # Example 3: Vecchia approximated model fit <- fit_one_layer(x, y, nmcmc = 2000, vecchia = TRUE, m = 10) plot(fit) fit <- trim(fit, 1000, 2) fit <- predict(fit, xx, cores = 1) plot(fit)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.