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 "hetGP"-class or "homGP"-class object as output from one of the main fitting functions, e.g., mleHetGP and mleHomGP. TP variants are not yet supported.

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 Xnew.

Znew

optional averaged observations at Xnew. If not provided, only the MSPE is returned

multnew

number of replicates at each design in Xnew

covreturn

boolean to return the predicte covariance matrix

forceSym

boolean to enforce the return of a symmetric predictive covariance matrix (default to TRUE)

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)

[Package hetGP version 1.1.8 Index]