Package 'gek'

Rosenbrock's Banana Function

Description

Rosenbrock's banana function is defined by

$f_{\rm banana}(x_1, ..., x_d) = \sum_{k = 1}^{d - 1} (100 (x_{k+1} - x_k^2)^2 + (x_k - 1)^2)$

with $x_k \in [-5, 10]$ for $k = 1, ..., d$ and $d \geq 2$.

Usage

banana(x)
bananaGrad(x)

Arguments

x	a numeric vector of length d or a numeric matrix with n rows and d columns, where d must be greater than 1.

Details

The gradient of Rosenbrock's banana function is

$\nabla f_{\rm banana}(x_1, ..., x_d) = \begin{pmatrix} -400 (x_2 - x_1)^2 x_1 + 2 (x_1 - 1) \\ 200 (x_2 - x_1)^2 - 400 x_2 (x_3 - x_2^2) + 2 (x_2 - 1) \\ \vdots \\ 200 (x_{d-1} - x_{d-2})^2 - 400 x_{d-1} (x_d - x_{d-1}^2) + 2 (x_{d-1} - 1) \\ 200 (x_d - x_{d -1}^2)\end{pmatrix}.$

Rosenbrock's banana function has one global minimum $f_{\rm banana}(x^{\star}) = 0$ at $x^{\star} = (1,\dots, 1)$.

Value

banana returns the function value of Rosenbrock's banana function at x.

bananaGrad returns the gradient of Rosenbrock's banana function at x.

Author(s)

Carmen van Meegen

References

Jamil, M. and Yang, X.-S. (2013). A Literature Survey of Benchmark Functions for Global Optimization Problems. International Journal of Mathematical Modelling and Numerical Optimisation, 4(2):150-–194. doi:10.1504/IJMMNO.2013.055204.

Rosenbrock, H. H. (1960). An Automatic Method for Finding the Greatest or least Value of a Function. The Computer Journal, 3(3):175–184. doi:10.1093/comjnl/3.3.175.

Surjanovic, S. and Bingham, D. (2013). Virtual Library of Simulation Experiments: Test Functions and Datasets. https://www.sfu.ca/~ssurjano/ (retrieved January 19, 2024).

Examples

# Contour plot of Rosenbrock's banana function 
n.grid <- 50
x1 <- seq(-2, 2, length.out = n.grid)
x2 <- seq(-1, 3, length.out = n.grid)
y <- outer(x1, x2, function(x1, x2) banana(cbind(x1, x2)))
contour(x1, x2, y, xaxs = "i", yaxs = "i", nlevels = 25, xlab = "x1", ylab = "x2")

# Perspective plot of Rosenbrock's banana function
col.pal <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
	"#FF7F00", "red", "#7F0000"))
colors <- col.pal(100)
y.facet.center <- (y[-1, -1] + y[-1, -n.grid] + y[-n.grid, -1] + y[-n.grid, -n.grid])/4
y.facet.range <- cut(y.facet.center, 100)
persp(x1, x2, y, phi = 30, theta = -315, expand = 0.75, ticktype = "detailed", 
	col = colors[y.facet.range])

---

Block Cholesky Decomposition

Description

blockChol calculates the block Cholesky decomposition of a partitioned matrix of the form

$A = \begin{pmatrix} K & R^{\rm T} \\ R & S \end{pmatrix}.$

Usage

blockChol(K, R = NULL, S = NULL, tol = NULL)

Arguments

K	a real symmetric positive-definite square submatrix.
R	an (optinal) rectangular submatrix.
S	an (optional) real symmetric positive-definite square submatrix.
tol	an (optional) numeric tolerance, see ‘Details’.

Details

To obtain the block Cholesky decomposition

$\begin{pmatrix} K & R^{\rm T} \\ R & S \end{pmatrix} = \begin{pmatrix} L^{\rm T} & 0 \\ Q^{\rm T} & M^{\rm T} \end{pmatrix} \begin{pmatrix} L & Q \\ 0 & M \end{pmatrix}$

the following steps are performed:

1. Calculate $K = L^{\rm T}L$ with upper triangular matrix $L$.

2. Solve $L^{\rm T}Q = R^{\rm T}$ via forward substitution.

3. Compute $N = S - Q^{\rm T}Q$ the Schur complement of the block $K$ of the matrix $A$.

4. Determine $N = M^{\rm T}M$ with upper triangular matrix $M$.

The upper triangular matrices $L$ and $M$ in step 1 and 4 are obtained by chol. Forward substitution in step 2 is carried out with backsolve and the option transpose = TRUE.

If tol is specified a regularization of the form $A_{\epsilon} = A + \epsilon I$ is conducted. Here, tol is the upper bound for the logarithmic condition number of $A_{\epsilon}$. Then

$\epsilon = \max\left\{ \frac{\lambda_{\max}(\kappa(A) - e^{\code{tol}})}{\kappa(A)(e^{\code{tol}} - 1)}, 0 \right\}$

is chosen as the minimal "nugget" that is added to the diagonal of $A$ to ensure $\log(\kappa(A_{\epsilon})) \leq$ tol.

Within gek this function is used to calculate the block Cholesky decomposition of the correlation matrix with derivatives. Here K is the Kriging correlation matrix. R is the matrix containing the first partial derivatives and S consists of the second partial derivatives of the correlation matrix K.

Value

blockChol returns a list with the following components:

L	the upper triangular factor of the Cholesky decomposition of K.
Q	the solution of the triangular system t(L) %*% Q == t(R).
M	the upper triangular factor of the Cholesky decomposition of the Schur complement N.

If R or S are not specified, Q and M are set to NULL, i.e. only the Cholesky decomposition of K is calculated.

The attribute "eps" gives the minimum “nugget” that is added to the diagonal.

Warning

As in chol there is no check for symmetry.

Author(s)

Carmen van Meegen

References

Chen, J., Jin, Z., Shi, Q., Qiu, J., and Liu, W. (2013). Block Algorithm and Its Implementation for Cholesky Factorization.

Gustavson, F. G. and Jonsson, I. (2000). Minimal-storage high-performance Cholesky factorization via blocking and recursion. IBM Journal of Research and Development, 44(6):823–850. doi:10.1147/rd.446.0823.

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. doi:10.1198/TECH.2011.09141.

See Also

chol for the Cholesky decomposition.

backsolve for backward substitution.

blockCor for computing a correlation matrix with derivatives.

Examples

# Construct correlation matrix
x <- matrix(seq(0, 1, length = 5), ncol = 1)
res <- blockCor(x, theta = 1, covtype = "gaussian", derivatives = TRUE)
A <- cbind(rbind(res$K, res$R), rbind(t(res$R), res$S))

# Cholesky decomposition of correlation matix without derivatives
cholK <- blockChol(res$K)
cholK
cholK$L == chol(res$K)

# Cholesky decomposition of correlation matix with derivatives
cholA <- blockChol(res$K, res$R, res$S) 
cholA <- cbind(rbind(cholA$L, matrix(0, ncol(cholA$Q), nrow(cholA$Q))), 
	rbind(cholA$Q, cholA$M))
cholA
cholA == chol(A)

# Cholesky decomposition of correlation matix with derivatives with regularization
res <- blockCor(x, theta = 2, covtype = "gaussian", derivatives = TRUE)
A <- cbind(rbind(res$K, res$R), rbind(t(res$R), res$S))
try(blockChol(res$K, res$R, res$S))
blockChol(res$K, res$R, res$S, tol = 35)

---

Correlation Matrix without or with Derivatives

Description

Calculation of a correlation matrix with or without derivatives according to the specification of a correlation structure.

Usage

blockCor(x, ...)

## Default S3 method:
blockCor(x, theta, covtype = c("matern5_2", "matern3_2", "gaussian"), 
	derivatives = FALSE, ...)
## S3 method for class 'gekm'
blockCor(x, ...)

Arguments

x	a numeric matrix or an object of class gekm.
theta	numeric vector of length d for the hyperparameters.
covtype	character specifying the correlation function to be used. Must be one of "matern5_2", "matern3_2" or "gaussian". See ‘Details’.
derivatives	logical, if TRUE the first and second partial derivatives of the correlation matrix are calculated, otherwise not.
...	further arguments passed to or from other methods.

Details

The correlation matrix with derivatives is defined as a block matrix of the form

$\begin{pmatrix} K & R^{\rm T} \\ R & S \end{pmatrix}.$

Three correlation functions are currently implemented in blockCor:

• Matérn 5/2 kernel with covtype = "matern5_2":

  $\phi(x, x'; \theta) = \prod_{k = 1}^d \left(1 + \frac{\sqrt{5}|x_k - x_k'|}{\theta_k} + \frac{5 |x_k - x_k'|^2}{3 \theta_k^2}\right) \exp\left( -\frac{\sqrt{5}|x_k - x_k'|}{\theta_k}\right)$

• Matérn 3/2 kernel with covtype = "matern3_2":

  $\phi(x, x'; \theta) = \prod_{k = 1}^d \left(1 + \frac{\sqrt{3} |x_k - x_k'|}{\theta_k} \right) \exp\left( -\frac{\sqrt{3} |x_k - x_k'|}{\theta_k}\right)$

• Gaussian kernel with covtype = "gaussian":

  $\phi(x, x'; \theta) = \prod_{k = 1}^d \exp\left( -\frac{(x_k - x_k')^2}{2\theta_k^2}\right)$

Value

blockCor returns a list with the following components:

K	The correlation matrix without derivatives.
R	If derivatives = TRUE, the correlation matrix with first partial derivatives, otherwise NULL.
S	If derivatives = TRUE, the correlation matrix with second partial derivatives, otherwise NULL.

The components of the list can be combined in the following form to constructed the complete correlation matrix with derivatives: cbind(rbind(K, R), rbind(t(R), S)).

Author(s)

Carmen van Meegen

References

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. doi:10.1016/S0169-7161(96)13011-X.

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. The MIT Press. https://gaussianprocess.org/gpml/.

Stein, M. L. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer Series in Statistics. Springer-Verlag. doi:10.1007/978-1-4612-1494-6.

See Also

blockChol for block Cholesky decomposition.

tangents for drawing tangent lines.

Examples

# Some examples of correlation matrices:

x <- matrix(1:4, ncol = 1)
blockCor(x, theta = 1)
blockCor(x, theta = 1, covtype = "gaussian")
blockCor(x, theta = 1, covtype = "gaussian", derivatives = TRUE)
blockCor(x, theta = 1, covtype = "matern3_2", derivatives = TRUE)

res <- blockCor(x, theta = 2, covtype = "gaussian", derivatives = TRUE)
cbind(rbind(res$K, res$R), rbind(t(res$R), res$S))


# Illustration of correlation functions and their derivatives:

x <- seq(0, 1, length = 100)
X <- matrix(x, ncol = 1)
gaussian <- blockCor(X, theta = 0.25, covtype = "gaussian", derivatives = TRUE)
matern5_2 <- blockCor(X, theta = 0.25, covtype = "matern5_2", derivatives = TRUE)
matern3_2 <- blockCor(X, theta = 0.25, covtype = "matern3_2", derivatives = TRUE)

# Correlation functions and first partial derivatives:
index <- c(10, 20, 40, 80)
par(mar = c(5.1, 5.1, 4.1, 2.1))

# Matern 3/2
plot(x, matern3_2$K[1, ], type = "l", xlab = expression(group("|", x - x*minute, "|")),
	ylab = expression(phi(x, x*minute, theta == 0.25)), lwd = 2)
tangents(x[index], matern3_2$K[index, 1], matern3_2$R[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern3_2$K[index, 1], pch = 16)

# Matern 5/2
lines(x, matern5_2$K[1, ], lwd = 2, col = 3)
tangents(x[index], matern5_2$K[index, 1], matern5_2$R[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern5_2$K[index, 1], pch = 16)

# Gaussian
lines(x, gaussian$K[1, ], lwd = 2, col = 4)
tangents(x[index], gaussian$K[index, 1], gaussian$R[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], gaussian$K[index, 1], pch = 16)

legend("topright", lty = 1, lwd = 2, col = c(1, 3, 4), bty = "n",
	legend = c("Matern 3/2", "Matern 5/2", "Gaussian"))


# First and second partial derivatives of correlation functions:
index <- c(5, 10, 20, 40, 80)
par(mar = c(5.1, 6.1, 4.1, 2.1))

# Gaussian
plot(x, matern3_2$R[1, ], type = "l", xlab = expression(group("|", x - x*minute, "|")),
	ylab = expression(frac(partialdiff * phi(x, x*minute, theta == 0.25), 
		partialdiff * x * minute)), lwd = 2)
tangents(x[index], matern3_2$R[1, index], matern3_2$S[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern3_2$R[1, index], pch = 16)

# Matern 5/2
lines(x, matern5_2$R[1, ], lwd = 2, col = 3)
tangents(x[index], matern5_2$R[1, index], matern5_2$S[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], matern5_2$R[1, index], pch = 16)

# Matern 3/2
lines(x, gaussian$R[1, ], lwd = 2, col = 4)
tangents(x[index], gaussian$R[1, index], gaussian$S[index, 1], 
	length = 0.15, lwd = 2, col = 2)
points(x[index], gaussian$R[1, index], pch = 16)

legend("topright", lty = 1, lwd = 2, col = c(1, 3, 4), bty = "n",
	legend = c("Matern 3/2", "Matern 5/2", "Gaussian"))

---

Borehole Function

Description

The borehole function is defined by

$f_{\rm borehole}(x) = \frac{2 \pi T_u (H_u - H_l)}{\log(r/r_w)\left(1 + \frac{2 L T_u}{\log(r / r_w) r_w^2 K_w} + \frac{T_u}{T_l}\right)}$

with $x = (r_w, r, T_u, H_u, T_l, H_l, L, K_w)$.

Usage

borehole(x)
boreholeGrad(x)

Arguments

x	a numeric vector of length 8 or a numeric matrix with n rows and 8 columns.

Details

The borehole function calculates the water flow rate $\rm [m^3/yr]$ through a borehole.

Input	Domain	Distribution	Description
$r_w$	$[0.05, 0.15]$	$\mathcal{N}(0.1, 0.0161812)$	radius of borehole in $\rm m$
$r$	$[100, 50\,000]$	$\mathcal{LN}(7.71, 1.0056)$	radius of influence in $\rm m$
$T_u$	$[63070, 115600]$	$\mathcal{U}(63070, 115600)$	transmissivity of upper aquifer in $\rm m^2/yr$
$H_u$	$[990, 1100]$	$\mathcal{U}(990, 1110)$	potentiometric head of upper aquifer in $\rm m$
$T_l$	$[63.1, 116]$	$\mathcal{U}(63.1, 116)$	transmissivity of lower aquifer in $\rm m^2/yr$
$H_l$	$[700, 820]$	$\mathcal{U}(700, 820)$	potentiometric head of lower aquifer in $\rm m$
$L$	$[1120, 1680]$	$\mathcal{U}(1120, 1680)$	length of borehole in $\rm m$
$K_w$	$[9855, 12045]$	$\mathcal{U}(9855, 12045)$	hydraulic conductivity of borehole in $\rm m/yr$

Note, $\mathcal{N}(\mu, \sigma)$ represents the normal distribution with expected value $\mu$ and standard deviation $\sigma$ and $\mathcal{LN}(\mu, \sigma)$ is the log-normal distribution with mean $\mu$ and standard deviation $\sigma$ of the logarithm. Further, $\mathcal{U}(a,b)$ denotes the continuous uniform distribution over the interval $[a,b]$.

Value

borehole returns the function value of borehole function at x.

boreholeGrad returns the gradient of borehole function at x.

Author(s)

Carmen van Meegen

References

Harper, W. V. and Gupta, S. K. (1983). Sensitivity/Uncertainty Analysis of a Borehole Scenario Comparing Latin Hypercube Sampling and Deterministic Sensitivity Approaches. BMI/ONWI-516, Office of Nuclear Waste Isolation, Battelle Memorial Institute, Columbus, OH.

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. doi:10.1080/00401706.1993.10485320.

See Also

gekm for another example.

Examples

# 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))

# Function for Monte Carlo simulation
samples <- function(x, N = 10^5){
	switch(x$dist,
		"norm" = rnorm(N, x$mean, x$sd),
		"lnorm" = rlnorm(N, x$meanlog, x$sdlog),
		"unif" = runif(N, x$min, x$max))
}

# Uncertainty distribution of the water flow rate
set.seed(1)
X <- sapply(inputs, samples)
y <- borehole(X)
hist(y, breaks = 50, xlab = expression(paste("Water flow rate ", group("[", m^3/yr, "]"))), 
	main = "", freq = FALSE)

---

Branin-Hoo Function

Description

The Branin-Hoo function is defined by

$f_{\rm branin}(x_1, x_2) = \left(x_2 - \frac{5.1}{4 \pi^2}x_1^2 + \frac{5}{\pi}x_1 - 6\right)^2 + 10 \left(1-\frac{1}{8\pi}\right)\cos(x_1) + 10$

with $x_1 \in [-5, 10]$ and $x_2 \in [0, 15]$.

Usage

branin(x)
braninGrad(x)

Arguments

x	a numeric vector of length 2 or a numeric matrix with n rows and 2 columns.

Details

The gradient of the Branin-Hoo function is

$\nabla f_{\rm branin}(x_1, x_2) = \begin{pmatrix} 2 \left(x_2 - \frac{5.1 x_1^2}{4 \pi^2} + \frac{5 x_1}{\pi} - 6\right) \left(-10.2 \frac{x_1}{4\pi^2} + \frac{5}{\pi}\right) - 10 \left(1 - \frac{1}{8\pi}\right) \sin(x_1) \\ 2 \left( x_2 - \frac{5.1 x_1^2}{4 \pi^2} + \frac{5 x_1}{\pi} - 6\right)\end{pmatrix}.$

The Branin-Hoo function has three global minima $f_{\rm branin}(x^{\star}) = 0.397887$ at $x^{\star} = (-\pi, 12.275)$, $x^{\star} = (\pi, 2.275)$ and $x^{\star} = (9.42478, 2.475)$.

Value

branin returns the function value of the Branin-Hoo function at x.

braninGrad returns the gradient of the Branin-Hoo function at x.

Author(s)

Carmen van Meegen

References

Branin, Jr., F. H. (1972). Widely Convergent Method of Finding Multiple Solutions of Simultaneous Nonlinear Equations. IBM Journal of Research and Development, 16(5):504–522.

Surjanovic, S. and Bingham, D. (2013). Virtual Library of Simulation Experiments: Test Functions and Datasets. https://www.sfu.ca/~ssurjano/ (retrieved January 19, 2024).

Examples

# Contour plot of Branin-Hoo function 
n.grid <- 50
x1 <- seq(-5, 10, length.out = n.grid)
x2 <- seq(0, 15, length.out = n.grid)
y <- outer(x1, x2, function(x1, x2) branin(cbind(x1, x2)))
contour(x1, x2, y, xaxs = "i", yaxs = "i", nlevels = 25, xlab = "x1", ylab = "x2")

# Perspective plot of Branin-Hoo function
col.pal <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
	"#FF7F00", "red", "#7F0000"))
colors <- col.pal(100)
y.facet.center <- (y[-1, -1] + y[-1, -n.grid] + y[-n.grid, -1] + y[-n.grid, -n.grid])/4
y.facet.range <- cut(y.facet.center, 100)
persp(x1, x2, y, phi = 30, theta = -315, expand = 0.75, ticktype = "detailed", 
	col = colors[y.facet.range])

---

Three-Hump Camel Function

Description

The three-hump camel function is defined by

$f_{\rm camel3}(x_1, x_2) = 2 x_1^2 - 1.05 x_1^4 + \frac{x_1^6}{6} + x_1 x_2 + x_2^2$

with $x_1, x_2 \in [-5, 5]$.

Usage

camel3(x)
camel3Grad(x)

Arguments

x	a numeric vector of length 2 or a numeric matrix with n rows and 2 columns.

Details

The gradient of the three-hump camel function is

$\nabla f_{\rm camel3}(x_1, x_2) = \begin{pmatrix} 4 x_1 - 4.2 x_1^3 + x_1^5 + x_2 \\ x_1 + 2 x_2 \end{pmatrix}.$

The three-hump camel function has one global minimum $f_{\rm camel3}(x^{\star}) = 0$ at $x^{\star} = (0, 0)$.

Value

camel3 returns the function value of the three-hump camel function at x.

camel3Grad returns the gradient of the three-hump camel function at x.

Author(s)

Carmen van Meegen

References

Jamil, M. and Yang, X.-S. (2013). A Literature Survey of Benchmark Functions for Global Optimization Problems. International Journal of Mathematical Modelling and Numerical Optimisation, 4(2):150-–194. doi:10.1504/IJMMNO.2013.055204.

Surjanovic, S. and Bingham, D. (2013). Virtual Library of Simulation Experiments: Test Functions and Datasets. https://www.sfu.ca/~ssurjano/ (retrieved January 19, 2024).

Examples

# Contour plot of three-hump camel function 
n.grid <- 50
x1 <- x2 <- seq(-2, 2, length.out = n.grid)
y <- outer(x1, x2, function(x1, x2) camel3(cbind(x1, x2)))
contour(x1, x2, y, xaxs = "i", yaxs = "i", nlevels = 25, xlab = "x1", ylab = "x2")

# Perspective plot of three-hump camel function
col.pal <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
	"#FF7F00", "red", "#7F0000"))
colors <- col.pal(100)
y.facet.center <- (y[-1, -1] + y[-1, -n.grid] + y[-n.grid, -1] + y[-n.grid, -n.grid])/4
y.facet.range <- cut(y.facet.center, 100)
persp(x1, x2, y, phi = 30, theta = -315, expand = 0.75, ticktype = "detailed", 
	col = colors[y.facet.range])

---

Six-Hump Camel Function

Description

The six-hump camel function is defined by

$f_{\rm camel6}(x_1, x_2) = \left(4 - 2.1 x_1^2 + \frac{x_1^4}{3}\right) x_1^2 + x_1 x_2 + (-4 + 4 x_2^2) x_2^2$

with $x_1 \in [-3, 3]$ and $x_2 \in [-2, 2]$.

Usage

camel6(x)
camel6Grad(x)

Arguments

x	a numeric vector of length 2 or a numeric matrix with n rows and 2 columns.

Details

The gradient of the six-hump camel function is

$\nabla f_{\rm camel6}(x_1, x_2) = \begin{pmatrix} 8 x_1 - 8.4 x_1^3 + 2 x_1^5 + x_2 \\ x_1 - 8 x_2 + 16 x_2^3 \end{pmatrix}.$

The six-hump camel function has two global minima $f_{\rm camel6}(x^{\star}) = -1.031628$ at $x^{\star} = (0.0898, -0.7126)$ and $x^{\star} = (-0.0898, 0.7126)$.

Value

camel6 returns the function value of the six-hump camel function function at x.

camel6Grad returns the gradient of the six-hump camel function function at x.

Author(s)

Carmen van Meegen

References

Jamil, M. and Yang, X.-S. (2013). A Literature Survey of Benchmark Functions for Global Optimization Problems. International Journal of Mathematical Modelling and Numerical Optimisation, 4(2):150-–194. doi:10.1504/IJMMNO.2013.055204.

Surjanovic, S. and Bingham, D. (2013). Virtual Library of Simulation Experiments: Test Functions and Datasets. https://www.sfu.ca/~ssurjano/ (retrieved January 19, 2024).

Examples

# Contour plot of six-hump camel function 
n.grid <- 50
x1 <- seq(-2, 2, length.out = n.grid)
x2 <- seq(-1, 1, length.out = n.grid)
y <- outer(x1, x2, function(x1, x2) camel6(cbind(x1, x2)))
contour(x1, x2, y, xaxs = "i", yaxs = "i", nlevels = 25, xlab = "x1", ylab = "x2")

# Perspective plot of six-hump camel function
col.pal <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
	"#FF7F00", "red", "#7F0000"))
colors <- col.pal(100)
y.facet.center <- (y[-1, -1] + y[-1, -n.grid] + y[-n.grid, -1] + y[-n.grid, -n.grid])/4
y.facet.range <- cut(y.facet.center, 100)
persp(x1, x2, y, phi = 30, theta = -315, expand = 0.75, ticktype = "detailed", 
	col = colors[y.facet.range])

---

Bent Cigar Function

Description

The Bent cigar function is defined by

$f_{\rm cigar}(x_1, ..., x_d) = x_1^2 + 10^6 \sum_{k = 2}^{d} x_k^2$

with $x_k \in [-100, 100]$ for $k = 1, ..., d$.

Usage

cigar(x)
cigarGrad(x)

Arguments

x	a numeric vector of length 2 or a numeric matrix with n rows and 2 columns.

Details

The gradient of the bent cigar function is

$\nabla f_{\rm cigar}(x_1, ..., x_d) = \begin{pmatrix} 2x_1 \\ 20^6x_2\\ \vdots \\ 20^6x_d\end{pmatrix}.$

The bent cigar function has one global minimum $f_{\rm cigar}(x^{\star}) = 0$ at $x^{\star} = (1,\dots, 1)$.

Value

cigar returns the function value of the bent cigar function at x.

cigarGrad returns the gradient of the bent cigar function at x.

Author(s)

Carmen van Meegen

References

Plevris, V. and Solorzano, G. (2022). A Collection of 30 Multidimensional Functions for Global Optimization Benchmarking. Data, 7(4):46. doi:10.3390/data7040046.

See Also

tangents for drawing tangent lines.

Examples

# 1-dimensional Cigar function with tangents
curve(cigar(x), from = -5, to = 5, n = 200)
x <- seq(-4.5, 4.5, length = 5)
y <- cigar(x)
dy <- cigarGrad(x)
tangents(x, y, dy, length = 2, lwd = 2, col = "red")
points(x, y, pch = 16)

# Contour plot of Cigar function 
n.grid <- 50
x1 <- x2 <- seq(-100, 100, length.out = n.grid)
y <- outer(x1, x2, function(x1, x2) cigar(cbind(x1, x2)))
contour(x1, x2, y, xaxs = "i", yaxs = "i", nlevels = 25, xlab = "x1", ylab = "x2")

# Perspective plot of Cigar function
col.pal <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
	"#FF7F00", "red", "#7F0000"))
colors <- col.pal(100)
y.facet.center <- (y[-1, -1] + y[-1, -n.grid] + y[-n.grid, -1] + y[-n.grid, -n.grid])/4
y.facet.range <- cut(y.facet.center, 100)
persp(x1, x2, y, phi = 30, theta = -315, expand = 0.75, ticktype = "detailed", 
	col = colors[y.facet.range])

---

Derivatives of Model Matrix

Description

Determine the derivatives of a model matrix.

Usage

derivModelMatrix(object, ...)

## Default S3 method:
derivModelMatrix(object, data, ...)
## S3 method for class 'gekm'
derivModelMatrix(object, ...)
## S3 method for class 'gekm'
model.matrix(object, ...)

Arguments

object	an object of an appropriate class. For the default method, a formula defining the regression functions.
data	a data.frame with named columns.
...	further arguments, yet not used.

Details

derivModelMatrix makes use of the function deriv. Accordingly, the calculation of derivatives is only possible for functions that are contained in the derivatives table of deriv.

Note, in contrast to model.matrix, factors are not supported.

Value

The derivatives of the model (or design) matrix.

As in model.matrix there is an attribute "assign".

Author(s)

Carmen van Meegen

See Also

deriv for more details on supported arithmetic operators and functions.

model.matrix for construction of a design (or model) matrix.

Examples

## Several examples for the derivatives of a model matrix

dat <- data.frame(x1 = seq(-2, 2, length.out = 5))

model.matrix(~ 1, dat)
derivModelMatrix(~ 1, dat)

model.matrix(~ ., dat)
derivModelMatrix(~ ., dat)

model.matrix(~ . - 1, dat)
derivModelMatrix(~ . - 1, dat)

model.matrix(~ sin(x1) + I(x1^2), dat)
derivModelMatrix(~ sin(x1) + I(x1^2), dat)

dat <- cbind(dat, x2 = seq(1, 5, length.out = 5))

model.matrix(~ 1, dat)
derivModelMatrix(~ 1, dat)

model.matrix(~ .^2, dat)
derivModelMatrix(~ .^2, dat)

model.matrix(~ log(x2), dat)
derivModelMatrix(~ log(x2), dat)

model.matrix(~ x1:x2, dat)
derivModelMatrix(~ x1:x2, dat)

model.matrix(~ I(x1^2) * I(x2^3), dat)
derivModelMatrix(~ I(x1^2) * I(x2^3), dat)

model.matrix(~ sin(x1) + cos(x2) + atan(x1 * x2), dat)
derivModelMatrix(~ sin(x1) + cos(x2) + atan(x1 * x2), dat)

---

Fitting (Gradient-Enhanced) Kriging Models

Description

Estimation of a Kriging model with or without derivatives.

Usage

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.

Value

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.

Author(s)

Carmen van Meegen

References

Cressie, N. A. C. (1993). Statistics for Spartial Data. John Wiley & Sons. 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. 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. 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. 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. 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. 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. 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. 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. 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/.

Ripley, B. D. (1981). Spatial Statistics. John Wiley & Sons. 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. 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. doi:10.1007/978-1-4612-1494-6.

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

## 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")

---

Griewank Function

Description

Griewank function is defined by

$f_{\rm griewank}(x_1, ..., x_d) = \sum_{k = 1}^{d} \frac{x_k^2}{4000} - \prod_{k = 1}^d \cos\left(\frac{x_k}{\sqrt{k}}\right) + 1$

with $x_k \in [-600, 600]$ for $k = 1, ..., d$.

Usage

griewank(x)
griewankGrad(x)

Arguments

x	a numeric vector or a numeric matrix with n rows and d columns. If a vector is passed, the 1-dimensional version of the Griewank function is calculated.

Details

The gradient of Griewank function is

$\nabla f_{\rm griewank}(x_1, ..., x_d) = \begin{pmatrix} \frac{x_1}{2000} + \frac{1}{\sqrt{1}} \sin\left( \frac{x_1}{\sqrt{1}}\right) \prod_{k = 2}^d \cos\left(\frac{x_k}{\sqrt{k}}\right) \\ \vdots \\ \frac{x_d}{2000} + \frac{1}{\sqrt{d}} \sin\left( \frac{x_d}{\sqrt{d}}\right) \prod_{k = 1}^{d-1} \cos\left(\frac{x_k}{\sqrt{k}}\right) \end{pmatrix}.$

Griewank function has one global minimum $f_{\rm griewank}(x^{\star}) = 0$ at $x^{\star} = (0,\dots, 0)$.

Value

griewank returns the function value of Griewank function at x.

griewankGrad returns the gradient of Griewank function at x.

Author(s)

Carmen van Meegen

References

Plevris, V. and Solorzano, G. (2022). A Collection of 30 Multidimensional Functions for Global Optimization Benchmarking. Data, 7(4):46. doi:10.3390/data7040046.

Surjanovic, S. and Bingham, D. (2013). Virtual Library of Simulation Experiments: Test Functions and Datasets. https://www.sfu.ca/~ssurjano/ (retrieved January 19, 2024).

See Also

tangents for drawing tangent lines.

Examples

# 1-dimensional Griewank function with tangents
curve(griewank(x), from = -5, to = 5, n = 200)
x <- seq(-5, 5, length = 5)
y <- griewank(x)
dy <- griewankGrad(x)
tangents(x, y, dy, length = 2, lwd = 2, col = "red")
points(x, y, pch = 16)

# Contour plot of Griewank function 
n.grid <- 50
x1 <- x2 <- seq(-5, 5, length.out = n.grid)
y <- outer(x1, x2, function(x1, x2) griewank(cbind(x1, x2)))
contour(x1, x2, y, xaxs = "i", yaxs = "i", nlevels = 25, xlab = "x1", ylab = "x2")

# Perspective plot of Griewank function
col.pal <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
	"#FF7F00", "red", "#7F0000"))
colors <- col.pal(100)
y.facet.center <- (y[-1, -1] + y[-1, -n.grid] + y[-n.grid, -1] + y[-n.grid, -n.grid])/4
y.facet.range <- cut(y.facet.center, 100)
persp(x1, x2, y, phi = 30, theta = -315, expand = 0.75, ticktype = "detailed", 
	col = colors[y.facet.range])

---

Himmelblau's Function

Description

Himmelblau's function is defined by

$f_{\rm himmelblau}(x_1, x_2) = (x_1^2 + x_2 - 11)^2 + (x_1 + x_2^2 -7)^2$

with $x_1, x_2 \in [-5, 5]$.

Usage

himmelblau(x)
himmelblauGrad(x)

Arguments

x	a numeric vector of length 2 or a numeric matrix with n rows and 2 columns.

Details

The gradient of Himmelblau's function is

$\nabla f_{\rm himmelblau}(x_1, x_2) = \begin{pmatrix} 4 x_1 (x_1^2 + x_2 - 11) + 2 (x_1 + x_2^2 - 7) \\ 2 (x_1^2 + x_2 - 11) + 4 x_2 (x_1 + x_2^2 - 7) \end{pmatrix}.$

Himmelblau's function has four global minima $f_{\rm himmelblau}(x^{\star}) = 0$ at $x^{\star} = (3, 2)$, $x^{\star} = (-2.805118, 3.131312)$, $x^{\star} = (-3.779310, -3.283186)$ and $x^{\star} = (3.584428, -1.848126)$.

Value

himmelblau returns the function value of Himmelblau's function at x.

himmelblauGrad returns the gradient of Himmelblau's function at x.

Author(s)

Carmen van Meegen

References

Himmelblau, D. (1972). Applied Nonlinear Programming. McGraw-Hill. ISBN 0-07-028921-2.

Jamil, M. and Yang, X.-S. (2013). A Literature Survey of Benchmark Functions for Global Optimization Problems. International Journal of Mathematical Modelling and Numerical Optimisation, 4(2):150-–194. doi:10.1504/IJMMNO.2013.055204.

Surjanovic, S. and Bingham, D. (2013). Virtual Library of Simulation Experiments: Test Functions and Datasets. https://www.sfu.ca/~ssurjano/ (retrieved January 19, 2024).

Examples

# Contour plot of Himmelblau's function 
n.grid <- 50
x1 <- x2 <- seq(-5, 5, length.out = n.grid)
y <- outer(x1, x2, function(x1, x2) himmelblau(cbind(x1, x2)))
contour(x1, x2, y, xaxs = "i", yaxs = "i", nlevels = 25, xlab = "x1", ylab = "x2")

# Perspective plot of Himmelblau's function
col.pal <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
	"#FF7F00", "red", "#7F0000"))
colors <- col.pal(100)
y.facet.center <- (y[-1, -1] + y[-1, -n.grid] + y[-n.grid, -1] + y[-n.grid, -n.grid])/4
y.facet.range <- cut(y.facet.center, 100)
persp(x1, x2, y, phi = 30, theta = -315, expand = 0.75, ticktype = "detailed", 
	col = colors[y.facet.range])

---

Predict Method for (Gradient-Enhanced) Kriging Model Fits

Description

Predicted values and standard deviations based on a gekm model.

Usage

## S3 method for class 'gekm'
predict(object, newdata, ...)

Arguments

object	an object of class "gekm".
newdata	a data.frame containing the points where to perform predictions.
...	further arguments, currently not used.

Value

fit	predicted mean computed at newdata.
sd.fit	predicted standard deviation of predicted mean.

Author(s)

Carmen van Meegen

References

Cressie, N. A. C. (1993). Statistics for Spartial Data. John Wiley & Sons. 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. 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. 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. 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. 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. 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. 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. 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. 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/.

Ripley, B. D. (1981). Spatial Statistics. John Wiley & Sons. 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. 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. doi:10.1007/978-1-4612-1494-6.

See Also

gekm for fitting a (gradient-enhanced) Kriging model.

tangents for drawing tangent lines.

Examples

## 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)

# Fit gradient-enhanced Kriging model
gekm.1d <- gekm(y ~ x, data = dat, deriv = deri, covtype = "gaussian", theta = 1)

# Generate new data for prediction
newdat <- data.frame(x = seq(-6, 6, length = 100))

# Prediction for both models
pred.km.1d <- predict(km.1d, newdat)
pred.gekm.1d  <- predict(gekm.1d, newdat)

# Draw curves
curve(f(x), from = -6, to = 6, lwd = 2)
lines(newdat$x, pred.km.1d$fit, type = "l", col = rgb(0.5, 0.5, 0.5), lwd = 2, lty = 2)
lines(newdat$x, pred.gekm.1d$fit, type = "l", col = rgb(0, 0.5, 1), lwd = 2, lty = 2)

# Add pointswise confidence intervals
lower.km.1d <- pred.km.1d$fit - qnorm(0.975) * pred.km.1d$sd
upper.km.1d <- pred.km.1d$fit + qnorm(0.975) * pred.km.1d$sd
lower.gekm.1d <- pred.gekm.1d$fit - qnorm(0.975) * pred.gekm.1d$sd
upper.gekm.1d <- pred.gekm.1d$fit + qnorm(0.975) * pred.gekm.1d$sd
 
polygon(c(newdat$x, rev(newdat$x)),	c(lower.km.1d, rev(upper.km.1d)),
	col = rgb(0.5, 0.5, 0.5, alpha = 0.1), border = NA)
polygon(c(newdat$x, rev(newdat$x)),	c(lower.gekm.1d, rev(upper.gekm.1d)),
	col = rgb(0, 0.5, 1, alpha = 0.1), border = NA)

# Add tangent lines and observations
tangents(x, y, dy, lwd = 2, col = "red")
points(x, y, pch = 16)

# Add legend
legend("topleft", lty = c(1, 2, 2), col = c("black", rgb(0.5, 0.5, 0.5), rgb(0, 0.5, 1)),
	legend = c("Test function", "Kriging", "GEK"), lwd = 2, bty = "n")


## 2-dimensional example: Branin-Hoo function

# Generate a grid for training
n <- 5
x1 <- seq(-5, 10, length = n)
x2 <- seq(0, 15, length = n)
x <- expand.grid(x1 = x1, x2 = x2)
y <- branin(x)
dy <- braninGrad(x)
dat <- data.frame(x, y)
deri <- data.frame(dy)

# Fit Kriging model
km.2d <- gekm(y ~ ., data = dat)
km.2d

# Fit gradient-enhanced Kriging model
gekm.2d <- gekm(y ~ ., data = dat, deriv = deri)
gekm.2d

# Generate new data for prediction
n.grid <- 50
x1.grid <- seq(-5, 10, length = n.grid)
x2.grid <- seq(0, 15, length = n.grid)
newdat <- expand.grid(x1 = x1.grid, x2 = x2.grid)

# Prediction for both models and actual outcome
pred.km.2d <- predict(km.2d, newdat)
pred.gekm.2d  <- predict(gekm.2d, newdat)
truth <- outer(x1.grid, x2.grid, function(x1, x2) branin(cbind(x1, x2)))

# Contour plots of predicted and actual output
par(mfrow = c(1, 3), oma = c(3.5, 3.5, 0, 0.2), mar = c(0, 0, 1.5, 0))
contour(x1.grid, x2.grid, truth, nlevels = 30, 
	levels = seq(0, 300, 10), main = "Branin-Hoo")
points(x, pch = 16)
contour(x1.grid, x2.grid, matrix(pred.km.2d$fit, nrow = n.grid, ncol = n.grid),
	levels = seq(0, 300, 10), main = "Kriging", yaxt = "n")
points(x, pch = 16)
contour(x1.grid, x2.grid, matrix(pred.gekm.2d$fit, nrow = n.grid, ncol = n.grid),
	levels = seq(0, 300, 10), main = "GEK", yaxt = "n")
points(x, pch = 16)
mtext(side = 1, outer = TRUE, line = 2.5, expression(x[1]))
mtext(side = 2, outer = TRUE, line = 2, expression(x[2]))

# Contour plots of predicted variance
par(mfrow = c(1, 2), oma = c(3.5, 3.5, 0, 0.2), mar = c(0, 0, 1.5, 0))
contour(x1.grid, x2.grid, matrix(pred.km.2d$sd.fit^2, nrow = n.grid, ncol = n.grid),
	main = "Kriging variance")
points(x, pch = 16)
contour(x1.grid, x2.grid, matrix(pred.gekm.2d$sd.fit^2, nrow = n.grid, ncol = n.grid),
	main = "GEK variance", yaxt = "n")
points(x, pch = 16)
mtext(side = 1, outer = TRUE, line = 2.5, expression(x[1]))
mtext(side = 2, outer = TRUE, line = 2, expression(x[2]))

---