bhetGP {bhetGP} | R Documentation |
MCMC sampling for Heteroskedastic GP
Description
Conducts MCMC sampling of hyperparameters and latent noise
process llam
for a hetGP. Separate length scale
parameters theta_lam
and theta_y
govern the correlation
strength of the hidden layer and outer layer respectively.
lam
layer may have a non-zero nugget g
which governs
noise for the latent noise layer. tau2_y
and tau2_lam
control the amplitude of the mean and noise process respectively.
In Matern covariance, v
governs smoothness.
Usage
bhetGP(
x = NULL,
y = NULL,
reps_list = NULL,
nmcmc = 1000,
sep = TRUE,
inits = NULL,
priors = NULL,
reps = TRUE,
cov = c("exp2", "matern", "ARD matern"),
v = 2.5,
stratergy = c("default", "flat"),
vecchia = FALSE,
m = min(25, length(y) - 1),
ordering = NULL,
verb = TRUE,
omp_cores = 4
)
Arguments
x |
vector or matrix of input locations |
y |
vector of response values |
reps_list |
list object from hetGP::find_reps |
nmcmc |
number of MCMC iterations |
sep |
logical indicating whether to fit isotropic GP ( |
inits |
set initial values for hyparameters:
|
priors |
hyperparameters for priors and proposals (see details) |
reps |
logical; if |
cov |
covariance kernel, either Matern, ARD Matern
or squared exponential ( |
v |
Matern smoothness parameter (only used if |
stratergy |
choose initialization stratergy; "default" uses hetGP for
|
vecchia |
logical indicating whether to use Vecchia approximation |
m |
size of Vecchia conditioning sets (only used if
|
ordering |
optional ordering for Vecchia approximation, must correspond
to rows of |
verb |
logical indicating whether to print progress |
omp_cores |
logical; if |
Details
Maps inputs x
to mean response y
and noise levels
llam
. Conducts sampling of the latent noise process using Elliptical
Slice sampling. Utilizes Metropolis Hastings sampling of the length
scale and nugget parameters with proposals and priors controlled by
priors
. g
for the noise process is set to a specific
value, and by default, is not estimated. When vecchia = TRUE
,
all calculations leverage the Vecchia approximation with
specified conditioning set size m
. tau2_y
is always
inferred from likelihood; tau2_lam
is inferred by default but
may be pre-specified and fixed.
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.
theta_star <- runif(1, l * theta_t / u, u * theta_t / l)
,
with defaults l = 1
and u = 2
provided in priors
.
To adjust these, set priors = list(l = new_l, u = new_u)
.
Priors on g
, theta_y
, and theta_lam
follow Gamma
distributions with shape parameters (alpha
) and rate parameters
(beta
) controlled within the priors
list object.
Defaults are
-
priors$alpha$theta_lam <- 1.5
-
priors$beta$theta_lam <- 4
-
priors$alpha$theta_y <- 1.5
-
priors$beta$theta_y <- 4
-
priors$alpha$g <- 1.5
-
priors$beta$g <- 4
tau2_y
and tau2_lam
are not sampled; rather directly inferred
under conjugate Inverse Gamma prior with shape (alpha
) and scale
parameters (beta
) within the priors
list object
-
priors$a$tau2_y <- 10
-
priors$a$tau2_y <- 4
-
priors$a$tau2_lam <- 10
-
priors$a$tau2_lam <- 4
These priors are designed for x
scaled to
[0, 1] and y
having mean mean_y
. These may be
adjusted using the priors
input.
Initial values for theta_y
, theta_lam
, llam
may be
specified by the user. If no initial values are specified, stratergy
will determine the initialization method. stratergy = "default"
leverages mleHetGP for initial values of hyper-parameters if
vecchia = FALSE
and Scaled-Vecchia with Stochastic Kriging (Sk-Vec)
hybrid approach if vecchia = TRUE
.
For SK-Vec hybrid approach, scaled Vecchia code from https://github.com/katzfuss-group/scaledVecchia/blob/master/vecchia_scaled.R is used to fit two GPs using the Vecchia approximation. The first for (x, y) pairs, which result in estimated residual sums of squares based on predicted y values. Another GP on (x, s) to obtain latent noise estimates which are smoothed. A script is leveraged internally within this package that fits this method.
Optionally, choose stratergy = "flat" which which will start at
uninformative initial values; llam
= log(var(y) * 0.1) or
specify initial values.
The output object of class bhetgp
or bhetgp_vec
is designed for
use with trim
, predict
, and plot
.
Value
a list of the S3 class bhetgp
or bhetgp_vec
with elements:
-
x
: copy of input matrix -
y
: copy of response mean at inputs (x) -
Ylist
: list of all responses observed per location (x) -
A
: number of replicates at each location -
nmcmc
: number of MCMC iterations -
priors
: copy of proposal/priors -
settings
: setting for predictions usingbhetgp
orbhetgp_vec
object -
theta_y
: vector of MCMC samples fortheta_y
(length scale of mean process) -
theta_lam
: matrix of MCMC samples fortheta_lam
(length scale of latent noise process) -
llam_samples
: matrix of ESS samples forlog lambda
(latent noise process samples) -
g
: vector of MCMC samples forg
if infered -
tau2
: vector of MAP estimates fortau2
(scale parameter of mean process) -
tau2_lam
: vector of MAP estimates fortau2_lam
(scale parameter of latent noise process) -
llik_y
: vector of MVN log likelihood of Y for reach Gibbs iteration -
llik_lam
: vector of MVN log likelihood ofllam
i.e. the latent noise process for reach Gibbs iteration -
x_approx
: conditioning set, NN and ordering forvecchia = TRUE
-
m
: copy of size of conditioning set forvecchia = TRUE
-
time
: computation time in seconds
References
Binois, Mickael, Robert B. Gramacy, and Mike Ludkovski. "Practical heteroscedastic Gaussian process modeling for large simulation experiments." Journal of Computational and Graphical Statistics 27.4 (2018): 808-821.
Katzfuss, Matthias, Joseph Guinness, and Earl Lawrence. "Scaled Vecchia approximation for fast computer-model emulation." SIAM/ASA Journal on Uncertainty Quantification 10.2 (2022): 537-554.
Sauer, Annie Elizabeth. "Deep Gaussian process surrogates for computer experiments." (2023).
Examples
# 1D function with 1D noise
# Truth
fx <- function(x){
result <- (6 * x - 2)^2* sin(12 * x - 4)
}
# Noise
rx <- function(x){
result <- (1.1 + sin(2 * pi * x))^2
return(result)
}
# Training data
r <- 10 # replicates
xn <- seq(0, 1, length = 25)
x <- rep(xn, r)
rn <- drop(rx(x))
noise <- as.numeric(t(mvtnorm::rmvnorm(r, sigma = diag(rn, length(xn)))))
f <- fx(x)
y <- f + noise
# Testing data
xx <- seq(0, 1, length = 100)
yy <- fx(xx)
#---------------------------------------------------------------------------
# Example 1: Full model, no Vecchia
#---------------------------------------------------------------------------
# Fitting a bhetGP model using all the data
fit <- bhetGP(x, y, nmcmc = 100, verb = FALSE)
# Trimming the object to remove burn in and thin samples
fit <- trim(fit, 50, 10)
# Predition using the bhetGP object (indepedent predictions)
fit <- predict(fit, xx, cores = 2)
# Visualizing the mean predictive surface.
# Can run plot(fit, trace = TRUE) to view trace plots
plot(fit)
#---------------------------------------------------------------------------
# Example 2: Vecchia approximated model
#---------------------------------------------------------------------------
# Fitting a bhetGP model with vecchia approximation. Two cores for OpenMP
fit <- bhetGP(x, y, nmcmc = 100, vecchia = TRUE, m = 5, omp_cores = 2, verb = FALSE)
# Trimming the object to remove burn in and thin samples
fit <- trim(fit, 50, 10)
# Predition using the bhetGP_vec object with joint predictions (lite = FALSE)
# Two cores for OpenMP, default setting (omp_cores = 2). No SNOW
fit <- predict(fit, xx, lite = FALSE, vecchia = TRUE)
# Visualizing the mean predictive surface
plot(fit)
#---------------------------------------------------------------------------
# Example 3: Vecchia inference, non-vecchia predictions
#---------------------------------------------------------------------------
# Fitting a bhetGP model with vecchia approximation. Two cores for OpenMP
fit <- bhetGP(x, y, nmcmc = 200, vecchia = TRUE, m = 5, omp_cores = 2)
# Trimming the object to remove burn in and thin samples
fit <- trim(fit, 100, 10)
# Predition using the bhetGP object with joint predictions (lite = FALSE)
# Two cores for OpenMP which is default setting (omp_cores = 2)
# Two cores for SNOW (cores = 2)
fit <- predict(fit, xx, vecchia = FALSE, cores = 2, lite = FALSE)
# Visualizing the mean predictive surface
plot(fit)