fit_one_layer {deepgp}

R Documentation

MCMC sampling for one layer GP

Description

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.

Usage

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)
)

Arguments

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 g
theta_0	initial value for theta
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 ("exp2")
v	Matern smoothness parameter (only used if cov = "matern")
vecchia	logical indicating whether to use Vecchia approximation
m	size of Vecchia conditioning sets (only used if vecchia = TRUE)

Details

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.

Value

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

References

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

# 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)

---