## 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 correlation parameters
# 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
sigma(mod)
# Predicted 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))
# Predicted 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, sd.fit = FALSE)
# 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, 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")
Run the code above in your browser using DataLab