# Fitting (Gradient-Enhanced) Kriging Models Estimation of a Kriging model with or without derivatives. UTF-8 ```r gekm(formula, data, deriv, covtype = c("matern5_2", "matern3_2", "gaussian"), theta = NULL, tol = NULL, optimizer = c("NMKB", "L-BFGS-B"), lower = NULL, upper = NULL, start = NULL, ncalls = 20, control = NULL, model = TRUE, x = FALSE, y = FALSE, dx = FALSE, dy = FALSE, ...) ## S3 method for class 'gekm' print(x, digits = 4L, ...) ``` ## Arguments - `formula`: a `formula` that defines the regression functions. Note that only formula expressions for which the derivations are contained in the derivatives table of `deriv` are supported in the gradient-enhanced Kriging model. In addition, formulas containing `I` also work for the gradient-enhanced Kriging model, although this function is not included in the derivatives table of `deriv`. See `derivModelMatrix` for some examples of the trend specification. - `data`: a `data.frame` with named columns of `n` training points of dimension `d`. Note, all variables contained in `data` are used for the construction of the correlation matrix without and with derivatives. - `deriv`: an optional `data.frame` with the derivatives, whose columns should be named like those of `data`. If not specified, a Kriging model without derivatives is estimated. - `covtype`: a `character` to specify the covariance structure to be used. One of `matern5_2`, `matern3_2` or `gaussian`. Default is `matern5_2`. - `theta`: a `numeric` vector of length `d` for the hyperparameters (optional). If not given, hyperparameters will be estimated via maximum likelihood. - `tol`: a tolerance for the conditional number of the correlation matrix, see `blockChol` for details. Default is `NULL`, i.e. no regularization is applied. - `optimizer`: an optional `character` that characterizes the optimization algorithms to be used for maximum likelihood estimation. See Details . - `lower`: an optional lower bound for the optimization of the correlation parameters. - `upper`: an optional upper bound for the optimization of the correlation parameters. - `start`: an optional `vector` of inital values for the optimization of the correlation parameters. - `ncalls`: an optional `integer` that defines the number of randomly selected initial values for the optimization. - `control`: a `list` of control parameters for the optimization routine. See `optim` or `nmkb`. - `model`: `logical`. Should the model frame be returned? Default is `TRUE`. - `x`: `logical`. Should the model matrix be returned? Default is `FALSE`. - `y`: `logical`. Should the response vector be returned? Default is `FALSE`. - `dx`: `logical`. Should the derivative of the model matrix be returned? Default is `FALSE`. - `dy`: `logical`. Should the derivatives of the response be returned? Default is `FALSE`. - `digits`: number of digits to be used for the `print` method. - ``...``: further arguments, currently not used. ## Details Parameter estimation is performed via maximum likelihood. The `optimizer` argument can be used to select one of the optimization algorithms `"NMBK"` or `"L-BFGS-B"`. In the case of the `"L-BFGS-B"`, analytical gradients of the concentrated log likelihood are used. For one-dimensional problems, `optimize` is called and the algorithm selected via `optimizer` is ignored. ## Returns `gekm` returns an object of `class` `"gekm"` whose underlying structure is a list containing at least the following components: - **coefficients**: the estimated regression coefficients. - **sigma**: the estimated process standard deviation. - **theta**: the (estimated) correlation parameters. - **covtype**: the name of the correlation function. - **chol**: (the components of) the upper triangular matrix of the Cholesky decomposition of the correlation matrix. - **optimizer**: the optimization algorithm. - **convergence**: the convergence code. - **message**: information from the optimizer. - **logLik**: the value of the negative concentrated log likelihood at the estimated parameters. - **derivatives**: `TRUE` if a gradient-enhanced Kriging model was adapted, otherwise `FALSE`. - **data**: the `data.frame` that was specified via the `data` argument. - **deriv**: if `derivatives = TRUE`, the `data.frame` with the derivatives that was specified via the `deriv` argument. - **call**: the matched call. - **terms**: the `terms` object used. - **model**: if requested (the default), the model frame used. - **x**: if requested, the model matrix. - **y**: if requested, the response vector. - **dx**: if requested, the derivatives of the model matrix. - **dx**: if requested, the vector of derivatives of the response. ## References Cressie, N. A. C. (1993). **Statistics for Spartial Data**. John Wiley & Sons. tools:::Rd_expr_doi("10.1002/9781119115151") . Koehler, J. and Owen, A. (1996). Computer Experiments. In Ghosh, S. and Rao, C. (eds.), **Design and Analysis of Experiments**, volume 13 of **Handbook of Statistics**, pp. 261–308. Elsevier Science. tools:::Rd_expr_doi("10.1016/S0169-7161(96)13011-X") . Krige, D. G. (1951). A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand. **Journal of the Southern African Institute of Mining and Metallurgy**, 52 (6):199--139. Laurent, L., Le Riche, R., Soulier, B., and Boucard, PA. (2019). An Overview of Gradient-Enhanced Metamodels with Applications. **Archives of Computational Methods in Engineering**, 26 (1):61--106. tools:::Rd_expr_doi("10.1007/s11831-017-9226-3") . Martin, J. D. and Simpson, T. W. (2005). Use of Kriging Models to Approximate Deterministic Computer Models. **AIAA Journal**, 43 (4):853--863. tools:::Rd_expr_doi("10.2514/1.8650") . Morris, M., Mitchell, T., and Ylvisaker, D. (1993). Bayesian Design and Analysis of Computer Experiments: Use of Derivatives in Surface Prediction. **Technometrics**, 35 (3):243--255. tools:::Rd_expr_doi("10.1080/00401706.1993.10485320") . Oakley, J. and O'Hagan, A. (2002). Bayesian Inference for the Uncertainty Distribution of Computer Model Outputs. **Biometrika**, 89 (4):769--784. tools:::Rd_expr_doi("10.1093/biomet/89.4.769") . O'Hagan, A., Kennedy, M. C., and Oakley, J. E. (1999). Uncertainty Analysis and Other Inference Tools for Complex Computer Codes. In **Bayesian Statistics 6**, Ed. J. M. Bernardo, J. O. Berger, A. P. Dawid and A .F. M. Smith, 503--524. Oxford University Press. O'Hagan, A. (2006). Bayesian Analysis of Computer Code Outputs: A Tutorial. **Reliability Engineering & System Safet**, 91 (10):1290--1300. tools:::Rd_expr_doi("10.1016/j.ress.2005.11.025") . Park, J.-S. and Beak, J. (2001). Efficient Computation of Maximum Likelihood Estimators in a Spatial Linear Model with Power Exponential Covariogram. **Computers & Geosciences**, 27 (1):1--7. tools:::Rd_expr_doi("10.1016/S0098-3004(00)00016-9") . Ranjan, P., Haynes, R. and Karsten, R. (2011). A Computationally Stable Approach to Gaussian Process Interpolation of Deterministic Computer Simulation Data. **Technometrics**, 53 :366--378. tools:::Rd_expr_doi("10.1198/TECH.2011.09141") . Rasmussen, C. E. and Williams, C. K. I. (2006). **Gaussian Processes for Machine Learning**. The MIT Press. [https://gaussianprocess.org/gpml/](https://gaussianprocess.org/gpml/). Ripley, B. D. (1981). **Spatial Statistics**. John Wiley & Sons. tools:::Rd_expr_doi("10.1002/0471725218") . Sacks, J., Welch, W. J., Mitchell, T. J., and Wynn, H. P. (1989). Design and Analysis of Computer Experiments. **Statistical Science**, 4 (4):409--423. tools:::Rd_expr_doi("10.1214/ss/1177012413") . Santner, T. J., Williams, B. J., and Notz, W. I. (2018). **The Design and Analysis of Computer Experiments**. 2nd edition. Springer-Verlag. Stein, M. L. (1999). **Interpolation of Spatial Data: Some Theory for Kriging**. Springer Series in Statistics. Springer-Verlag. tools:::Rd_expr_doi("10.1007/978-1-4612-1494-6") . ## Author(s) Carmen van Meegen ## See Also `predict.gekm` for prediction at new data points based on a model of class `"gekm"`. `simulate.gekm` for simulation of process paths conditional on a model of class `"gekm"`. ## Examples ```r ## 1-dimensional example: Oakley and O’Hagan (2002) # Define test function and its gradient f <- function(x) 5 + x + cos(x) fGrad <- function(x) 1 - sin(x) # Generate coordinates and calculate slopes x <- seq(-5, 5, length = 5) y <- f(x) dy <- fGrad(x) dat <- data.frame(x, y) deri <- data.frame(x = dy) # Fit Kriging model km.1d <- gekm(y ~ x, data = dat, covtype = "gaussian", theta = 1) km.1d # Fit Gradient-Enhanced Kriging model gekm.1d <- gekm(y ~ x, data = dat, deriv = deri, covtype = "gaussian", theta = 1) gekm.1d ## 2-dimensional example: Morris et al. (1993) # List of inputs with their distributions and their respective ranges inputs <- list("r_w" = list(dist = "norm", mean = 0.1, sd = 0.0161812, min = 0.05, max = 0.15), "r" = list(dist = "lnorm", meanlog = 7.71, sdlog = 1.0056, min = 100, max = 50000), "T_u" = list(dist = "unif", min = 63070, max = 115600), "H_u" = list(dist = "unif", min = 990, max = 1110), "T_l" = list(dist = "unif", min = 63.1, max = 116), "H_l" = list(dist = "unif", min = 700, max = 820), "L" = list(dist = "unif", min = 1120, max = 1680), # for a more nonlinear, nonadditive function, see Morris et al. (1993) "K_w" = list(dist = "unif", min = 1500, max = 15000)) # Generate design design <- data.frame("r_w" = c(0, 0.268, 1), "r" = rep(0, 3), "T_u" = rep(0, 3), "H_u" = rep(0, 3), "T_l" = rep(0, 3), "H_l" = rep(0, 3), "L" = rep(0, 3), "K_w" = c(0, 1, 0.268)) # Function to transform design onto input range transform <- function(x, data){ for(p in names(data)){ data[ , p] <- (x[[p]]$max - x[[p]]$min) * data[ , p] + x[[p]]$min } data } # Function to transform derivatives deriv.transform <- function(x, data){ for(p in colnames(data)){ data[ , p] <- data[ , p] * (x[[p]]$max - x[[p]]$min) } data } # Generate outcome and derivatives design.trans <- transform(inputs, design) design$y <- borehole(design.trans) deri.trans <- boreholeGrad(design.trans) deri <- data.frame(deriv.transform(inputs, deri.trans)) # Design and data cbind(design[ , c("r_w", "K_w", "y")], deri[ , c("r_w", "K_w")]) # Fit gradient-enhanced Kriging model with Gaussian correlation function mod <- gekm(y ~ 1, data = design[ , c("r_w", "K_w", "y")], deriv = deri[ , c("r_w", "K_w")], covtype = "gaussian") mod ## Compare results with Morris et al. (1993): # Estimated hyperparameters # in Morris et al. (1993): 0.429 and 0.467 1 / (2 * mod$theta^2) # Estimated intercept # in Morris et al. (1993): 69.15 coef(mod) # Estimated standard deviation # in Morris et al. (1993): 135.47 mod$sigma # Posterior mean and standard deviation at (0.5, 0.5) # in Morris et al. (1993): 69.4 and 2.7 predict(mod, data.frame("r_w" = 0.5, "K_w" = 0.5)) # Posterior mean and standard deviation at (1, 1) # in Morris et al. (1993): 230.0 and 19.2 predict(mod, data.frame("r_w" = 1, "K_w" = 1)) ## Graphical comparison: # Generate a 21 x 21 grid for prediction n_grid <- 21 x <- seq(0, 1, length.out = n_grid) grid <- expand.grid("r_w" = x, "K_w" = x) pred <- predict(mod, grid) # Compute ground truth on (transformed) grid newdata <- data.frame("r_w" = grid[ , "r_w"], "r" = 0, "T_u" = 0, "H_u" = 0, "T_l" = 0, "H_l" = 0, "L" = 0, "K_w" = grid[ , "K_w"]) newdata <- transform(inputs, newdata) truth <- borehole(newdata) # Contour plots of predicted and actual output par(mfrow = c(1, 2), oma = c(3.5, 3.5, 0, 0.2), mar = c(0, 0, 1.5, 0)) contour(x, x, matrix(pred$fit, nrow = n_grid, ncol = n_grid, byrow = TRUE), levels = c(seq(10, 50, 10), seq(100, 250, 50)), main = "Predicted output") points(design[ , c("K_w", "r_w")], pch = 16) contour(x, x, matrix(truth, nrow = n_grid, ncol = n_grid, byrow = TRUE), levels = c(seq(10, 50, 10), seq(100, 250, 50)), yaxt = "n", main = "Ground truth") points(design[ , c("K_w", "r_w")], pch = 16) mtext(side = 1, outer = TRUE, line = 2.5, "Normalized hydraulic conductivity of borehole") mtext(side = 2, outer = TRUE, line = 2.5, "Normalized radius of borehole") ```

gekm function