mleHomGP | R Documentation |
Gaussian process regression under homoskedastic noise based on maximum likelihood estimation of the hyperparameters. This function is enhanced to deal with replicated observations.
mleHomGP(
X,
Z,
lower = NULL,
upper = NULL,
known = NULL,
noiseControl = list(g_bounds = c(sqrt(.Machine$double.eps), 100)),
init = NULL,
covtype = c("Gaussian", "Matern5_2", "Matern3_2"),
maxit = 100,
eps = sqrt(.Machine$double.eps),
settings = list(return.Ki = TRUE, factr = 1e+07)
)
X |
matrix of all designs, one per row, or list with elements:
|
Z |
vector of all observations. If using a list with |
lower , upper |
optional bounds for the |
known |
optional list of known parameters, e.g., |
noiseControl |
list with element ,
|
init |
optional list specifying starting values for MLE optimization, with elements:
|
covtype |
covariance kernel type, either 'Gaussian', 'Matern5_2' or 'Matern3_2', see |
maxit |
maximum number of iteration for L-BFGS-B of |
eps |
jitter used in the inversion of the covariance matrix for numerical stability |
settings |
list with argument |
The global covariance matrix of the model is parameterized as nu_hat * (C + g * diag(1/mult)) = nu_hat * K
,
with C
the correlation matrix between unique designs, depending on the family of kernel used (see cov_gen
for available choices) and values of lengthscale parameters.
nu_hat
is the plugin estimator of the variance of the process.
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.
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,
nu_hat
: plugin estimator of the variance,
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 (not scaled by nu_hat
) (if return.Ki
is TRUE
in settings
)
time
: time to train the model, in seconds.
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.
predict.homGP
for predictions, update.homGP
for updating an existing model.
summary
and plot
functions are available as well.
mleHomTP
provide a Student-t equivalent.
##------------------------------------------------------------
## Example 1: Homoskedastic GP 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")
model <- mleHomGP(X = X, Z = Z, lower = 0.01, upper = 100)
## Display averaged observations
points(model$X0, model$Z0, pch = 20)
xgrid <- matrix(seq(0, 60, length.out = 301), ncol = 1)
predictions <- predict(x = xgrid, object = model)
## Display mean prediction
lines(xgrid, predictions$mean, col = 'red', lwd = 2)
## Display 95% confidence intervals
lines(xgrid, qnorm(0.05, predictions$mean, sqrt(predictions$sd2)), col = 2, lty = 2)
lines(xgrid, qnorm(0.95, predictions$mean, sqrt(predictions$sd2)), col = 2, lty = 2)
## Display 95% prediction intervals
lines(xgrid, qnorm(0.05, predictions$mean, sqrt(predictions$sd2 + predictions$nugs)),
col = 3, lty = 2)
lines(xgrid, qnorm(0.95, predictions$mean, sqrt(predictions$sd2 + predictions$nugs)),
col = 3, lty = 2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.