update_pred {hetGP} | R Documentation |
Prediction update with new designs and observations
Description
Fast updated prediction formulas (fixed hyperparameters)
Usage
update_pred(
model,
Xnew,
x = NULL,
Znew = NULL,
multnew = rep(1, nrow(Xnew)),
covreturn = TRUE,
forceSym = TRUE
)
Arguments
model |
a fitted |
Xnew |
matrix of new design location(s) |
x |
optional matrix of designs locations to predict at (one point per row). If not provided, the prediction is at |
Znew |
optional averaged observations at |
multnew |
number of replicates at each design in |
covreturn |
boolean to return the predicte covariance matrix |
forceSym |
boolean to enforce the return of a symmetric predictive covariance matrix (default to |
Details
This is alpha functionality at this time. Note that the nu_hat
parameter is not updated.
References
Gramacy, R. B. Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences
CRC Press, 2020
Chevalier, C., Ginsbourger, D., & Emery, X. (2013). Corrected kriging update formulae for batch-sequential data assimilation. In Mathematics of Planet Earth: Proceedings of the 15th Annual Conference of the International Association for Mathematical Geosciences (pp. 119-122). Berlin, Heidelberg: Springer Berlin Heidelberg.
Lyu, X., Binois, M. & Ludkovski, M. (2021). Evaluating Gaussian Process Metamodels and Sequential Designs for Noisy Level Set Estimation. Statistics and Computing, 31(4), 43. arXiv:1807.06712.
Examples
## Not run:
##------------------------------------------------------------
## Illustration of the update procedure
##------------------------------------------------------------
set.seed(1)
nvar <- 2
## test function
obj_fun <- function(x) return(sin(2 * x[1] * pi) + cos(2 * x[2] * pi))
noise_fun <- function(x, a = 0.1) return(a * sqrt(sum(x)))
ftest <- function(x) return(obj_fun(x) + rnorm(1, sd = noise_fun(x)))
n.grid <- 51
xgrid <- seq(0, 1, length.out = n.grid)
Xgrid <- as.matrix(expand.grid(xgrid, xgrid))
Yref <- apply(Xgrid, 1, obj_fun)
contour(xgrid, xgrid, matrix(Yref, n.grid))
## Unique (randomly chosen) design locations
n <- 25
Xu <- matrix(runif(n * 2), n)
## Select replication sites randomly
X <- Xu[sample(1:n, 20*n, replace = TRUE),]
## obtain training data response at design locations X
Z <- apply(X, 1, ftest)
## Formating of data for model creation (find replicated observations)
prdata <- find_reps(X, Z)
## Model fitting
model <- mleHetGP(X = list(X0 = prdata$X0, Z0 = prdata$Z0, mult = prdata$mult), Z = prdata$Z,
lower = rep(0.01, nvar), upper = rep(5, nvar), covtype = "Matern5_2")
predictions <- predict(x = Xgrid, object = model)
## Suppose that a new point is added, with some replicates
nnew <- 5
# i) a new unique design
Xnew1 <- matrix(runif(2), 1)
Znew1 <- apply(Xnew1[rep(1, nnew),, drop = FALSE], 1, ftest)
prednew1 <- update_pred(model = model, x = Xgrid, Xnew = Xnew1, Znew = mean(Znew1), multnew = nnew)
prednew1_sd2_only <- update_pred(model = model, x = Xgrid, Xnew = Xnew1, multnew = nnew)
# Verification
new_nug1 <- log(predict(x = Xnew1, model)$nugs/model$nu_hat)
model_v1 <- mleHetGP(X = list(X0 = rbind(prdata$X0, Xnew1),
Z0 = c(prdata$Z0, mean(Znew1)),mult = c(prdata$mult, nnew)),
Z = c(prdata$Z, Znew1),
lower = rep(0.01, nvar), upper = rep(5, nvar),
known = list(theta = model$theta, k_theta_g = model$k_theta_g,
theta_g = model$theta_g, g = model$g,
Delta = c(model$Delta, new_nug1)),
covtype = "Matern5_2")
pred_v1 <- predict(x = Xgrid, object = model_v1)
print(max(abs(pred_v1$mean - prednew1$mean)))
print(max(abs(pred_v1$sd2 - prednew1$sd2)))
print(max(abs(predictions$nugs - prednew1$nugs)))
## End(Not run)