deepgp (version 0.2.0)

fit_one_layer: MCMC sampling for one layer GP

Description

Conducts MCMC sampling of hyperparameters for a one layer GP. Covariance structure is based on inverse exponentiated squared euclidean distance with length scale parameter "theta" governing the strength of the correlation and nugget parameter "g" governing noise.

Usage

fit_one_layer(
  x,
  y,
  nmcmc = 10000,
  trace = TRUE,
  g_0 = 0.01,
  theta_0 = 0.5,
  true_g = NULL,
  settings = list(l = 1, u = 2, alpha = list(g = 1.5, theta = 1.5), beta = list(g =
    3.9, theta = 3.9/1.5))
)

Arguments

x

vector or matrix of input locations

y

vector of response values

nmcmc

number of MCMC iterations

trace

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)

settings

hyperparameters for proposals and priors on g and theta

Value

a list of the S3 class "gp" with elements:

• x: copy of input matrix

• y: copy of response vector

• nmcmc: number of MCMC iterations

• settings: copy of proposal/prior settings

• g: vector of MCMC samples for g

• theta: vector of MCMC samples for theta

• time: computation time in seconds

Details

Utilizes Metropolis Hastings sampling of the length scale and nugget parameters with proposals and priors controlled by settings. 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. Priors on g and theta follow Gamma distributions with shape parameter (alpha) and rate parameter (beta) provided in settings. These priors are designed for "x" scaled to [0,1] and "y" scaled to have mean zero and variance 1.

The output object of class "gp" is designed for use with continue, trim, and predict.

References

Sauer, A, RB Gramacy, and D Higdon. 2020. "Active Learning for Deep Gaussian Process Surrogates." arXiv:2012.08015. Gramacy, RB. Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences. Chapman Hall, 2020.

Examples

# NOT RUN {
# Toy example (runs in less than 5 seconds) -----------------------------------
# This example uses a small number of MCMC iterations in order to run quickly
# More iterations are required to get appropriate fits
# Function defaults are recommended (see additional example below)

f <- function(x) {
  if (x <= 0.4) return(-1)
  if (x >= 0.6) return(1)
  if (x > 0.4 & x < 0.6) return(10*(x-0.5))
}
x <- seq(0.05, 0.95, length = 7)
y <- sapply(x, f)
x_new <- seq(0, 1, length = 100)

# Fit model and calculate EI
fit <- fit_one_layer(x, y, nmcmc = 500)
fit <- trim(fit, 400)
fit <- predict(fit, x_new, lite = TRUE, store_all = TRUE)
ei <- EI(fit)

# }
# NOT RUN {
# One Layer and EI ------------------------------------------------------------

f <- function(x) {
  sin(5 * pi * x) / (2 * x) + (x - 1) ^ 4
}
  
# Training data
x <- seq(0.5, 2, length = 30)
y <- f(x) + rnorm(30, 0, 0.01)
  
# Testing data
xx <- seq(0.5, 2, length = 100)
yy <- f(xx)
  
# Standardize inputs and outputs
xx <- (xx - min(x)) / (max(x) - min(x))
x <- (x - min(x)) / (max(x) - min(x))
yy <- (yy - mean(y)) / sd(y)
y <- (y - mean(y)) / sd(y)
  
# Conduct MCMC
fit <- fit_one_layer(x, y, nmcmc = 10000)
plot(fit) # investigate trace plots
fit <- trim(fit, 8000, 2)
  
# Predict and calculate EI
fit <- predict(fit, xx, lite = TRUE, store_all = TRUE)
ei <- EI(fit)
  
# Visualize Fit
plot(fit)
par(new = TRUE) # overlay EI
plot(xx, ei$value, type = 'l', lty = 2, axes = FALSE, xlab = '', ylab = '')

# Select next design point
x_new <- xx[which.max(ei$value)]

# Evaluate fit
rmse(yy, fit$mean) # lower is better
# }
# NOT RUN {
# }