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, sep = FALSE, 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 
sep 
logical indicating whether to use separable ( 
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"
.
NOTE on OpenMP: The Vecchia implementation relies on OpenMP parallelization for efficient computation. This function will produce a warning message if the package was installed without OpenMP (this is the default for CRAN packages installed on Apple machines). To set up OpenMP parallelization, download the package source code and install using the gcc/g++ compiler.
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. "Vecchiaapproximated Deep Gaussian
Processes for Computer Experiments." preprint on arXiv:2204.02904
Gramacy, RB. Surrogates: Gaussian Process Modeling, Design, and
Optimization for the Applied Sciences. Chapman Hall, 2020.
# Examples of realworld implementations are available at: # https://bitbucket.org/gramacylab/deepgpex/ # 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 = 1e6) 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.