GADGET (version 0.2.0)

gp_residuals: Gaussian Process Residuals

Description

This function computes standardized and pivoted-Cholesky residuals of a Gaussian process (GP) model on a validation data set. Mahalanobis distance and Mahalanobis p-value are calculated. These statistics provide evidence of lack-of-fit in the GP model. The residuals can be plotted against predicted values as well as QQ-plots to check the normality assumption.

Usage

gp_residuals(design, response, model, plot = TRUE, type = "SK")

Arguments

design

A matrix of n rows and d columns.

response

A column vector of length n.

model

A GP model of class km (see km-class).

plot

Plot residuals and QQ-plots (with outliers are highlighted)?

type

Kriging type: Simple Kriging "SK" or Universal Kriging "UK".

Value

A list including the Mahalanobis distance (MD), MD F-statistic, MD p-value, pivoted-Cholesky residuals, and standardized residuals.

References

Bastos, L. S., & O'Hagan, A. (2009). Diagnostics for gaussian process emulators. Technometrics, 51(4), 425<U+2013>438, <doi:10.1198/TECH.2009.08019>.

Examples

# NOT RUN {
#--- Simple iid Normal Example ---#
#model assumputions hold
set.seed(123)
# training data
x           <- matrix(runif(20,-1.5,1.5),ncol=1)
y           <- matrix(rnorm(20),ncol = 1)
my_model    <- DiceKriging::km(formula=~1,
                               design=x,
                               response=y,
                               covtype='matern5_2',
                               optim.method='BFGS',
                               nugget.estim=TRUE)
# validation data
v_x         <- matrix(runif(25,-1,1),ncol=1)
v_y         <- matrix(rnorm(25),ncol = 1)
diagnostics <-gp_residuals(design = v_x, response = v_y,my_model)

#--- Bastos and O'Hagan (2009) Two-input Toy Model ---#
# needs more than 20 training points
set.seed(123)
# training data
x   <- lhs::randomLHS(20,2)
y   <- space_eval(x,bo09_toy)
# validation data
v_x <- lhs::randomLHS(25,2)
v_y <- space_eval(v_x,bo09_toy)
my_model    <- DiceKriging::km(formula=~1,
                               design=x,
                               response=y,
                               covtype='matern5_2',
                               optim.method='BFGS',
                               nugget.estim=TRUE)
diagnostics <- gp_residuals(v_x,v_y,my_model)
# }