fit_one_layer {deepgp} | 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)