mleHomTP {hetGP}

Student-T process modeling with homoskedastic noise

Description

Student-t process regression under homoskedastic noise based on maximum likelihood estimation of the hyperparameters. This function is enhanced to deal with replicated observations.

Usage

mleHomTP(
  X,
  Z,
  lower = NULL,
  upper = NULL,
  known = list(beta0 = 0),
  noiseControl = list(g_bounds = c(sqrt(.Machine$double.eps), 10000), nu_bounds = c(2 +
    0.001, 30), sigma2_bounds = c(sqrt(.Machine$double.eps), 10000)),
  init = list(nu = 3),
  covtype = c("Gaussian", "Matern5_2", "Matern3_2"),
  maxit = 100,
  eps = sqrt(.Machine$double.eps),
  settings = list(return.Ki = TRUE, factr = 1e+09)
)

Arguments

X

matrix of all designs, one per row, or list with elements:

• X0 matrix of unique design locations, one point per row

• Z0 vector of averaged observations, of length nrow(X0)

• mult number of replicates at designs in X0, of length nrow(X0)

Z

vector of all observations. If using a list with X, Z has to be ordered with respect to X0, and of length sum(mult)

lower, upper

bounds for the theta parameter (see cov_gen for the exact parameterization). In the multivariate case, it is possible to give vectors for bounds (resp. scalars) for anisotropy (resp. isotropy)

known

optional list of known parameters, e.g., beta0 (default to 0), theta, g, sigma2 or nu

noiseControl

list with element,

• g_bound, vector providing minimal and maximal noise variance

• sigma2_bounds, vector providing minimal and maximal signal variance

• nu_bounds, vector providing minimal and maximal values for the degrees of freedom. The mininal value has to be stricly greater than 2. If the mle optimization gives a large value, e.g., 30, considering a GP with mleHomGP may be better.

init

list specifying starting values for MLE optimization, with elements:

• theta_init initial value of the theta parameters to be optimized over (default to 10% of the range determined with lower and upper)

• g_init initial value of the nugget parameter to be optimized over (based on the variance at replicates if there are any, else 10% of the variance)

• sigma2 initial value of the variance paramter (default to 1)

• nu initial value of the degrees of freedom parameter (default to 3)

covtype

covariance kernel type, either 'Gaussian', 'Matern5_2' or 'Matern3_2', see cov_gen

maxit

maximum number of iteration for L-BFGS-B of optim

eps

jitter used in the inversion of the covariance matrix for numerical stability

settings

list with argument return.Ki, to include the inverse covariance matrix in the object for further use (e.g., prediction). Arguments factr (default to 1e9) and pgtol are available to be passed to control for L-BFGS-B in optim.

Details

The global covariance matrix of the model is parameterized as K = sigma2 * C + g * diag(1/mult), with C the correlation matrix between unique designs, depending on the family of kernel used (see cov_gen for available choices).

It is generally recommended to use find_reps to pre-process the data, to rescale the inputs to the unit cube and to normalize the outputs.

Value

a list which is given the S3 class "homGP", with elements:

• theta: maximum likelihood estimate of the lengthscale parameter(s),

• g: maximum likelihood estimate of the nugget variance,

• trendtype: either "SK" if beta0 is given, else "OK"

• beta0: estimated trend unless given in input,

• sigma2: maximum likelihood estimate of the scale variance,

• nu2: maximum likelihood estimate of the degrees of freedom parameter,

• ll: log-likelihood value,

• X0, Z0, Z, mult, eps, covtype: values given in input,

• call: user call of the function

• used_args: list with arguments provided in the call

• nit_opt, msg: counts and msg returned by optim

• Ki, inverse covariance matrix (if return.Ki is TRUE in settings)

• time: time to train the model, in seconds.

References

M. Binois, Robert B. Gramacy, M. Ludkovski (2018), Practical heteroskedastic Gaussian process modeling for large simulation experiments, Journal of Computational and Graphical Statistics, 27(4), 808–821.
Preprint available on arXiv:1611.05902.

A. Shah, A. Wilson, Z. Ghahramani (2014), Student-t processes as alternatives to Gaussian processes, Artificial Intelligence and Statistics, 877–885.

M. Chung, M. Binois, RB Gramacy, DJ Moquin, AP Smith, AM Smith (2019). Parameter and Uncertainty Estimation for Dynamical Systems Using Surrogate Stochastic Processes. SIAM Journal on Scientific Computing, 41(4), 2212-2238.
Preprint available on arXiv:1802.00852.

See Also

predict.homTP for predictions. summary and plot functions are available as well.

Examples

##------------------------------------------------------------
## Example 1: Homoskedastic Student-t modeling on the motorcycle data
##------------------------------------------------------------
set.seed(32)

## motorcycle data
library(MASS)
X <- matrix(mcycle$times, ncol = 1)
Z <- mcycle$accel
plot(X, Z, ylim = c(-160, 90), ylab = 'acceleration', xlab = "time")

noiseControl = list(g_bounds = c(1e-3, 1e4))
model <- mleHomTP(X = X, Z = Z, lower = 0.01, upper = 100, noiseControl = noiseControl)
summary(model)
  
## Display averaged observations
points(model$X0, model$Z0, pch = 20) 
xgrid <- matrix(seq(0, 60, length.out = 301), ncol = 1) 
preds <- predict(x = xgrid, object =  model)

## Display mean prediction
lines(xgrid, preds$mean, col = 'red', lwd = 2)
## Display 95% confidence intervals
lines(xgrid, preds$mean + sqrt(preds$sd2) * qt(0.05, df = model$nu + nrow(X)), col = 2, lty = 2)
lines(xgrid, preds$mean + sqrt(preds$sd2) * qt(0.95, df = model$nu + nrow(X)), col = 2, lty = 2)
## Display 95% prediction intervals
lines(xgrid, preds$mean + sqrt(preds$sd2 + preds$nugs) * qt(0.05, df = model$nu + nrow(X)), 
  col = 3, lty = 2)
lines(xgrid, preds$mean + sqrt(preds$sd2 + preds$nugs) * qt(0.95, df = model$nu + nrow(X)), 
  col = 3, lty = 2)