## 2-dimensional example
# Generate coordinates and calculate slopes
x1 <- seq(-1.75, 1.75, length = 3)
x2 <- seq(-0.75, 0.75, length = 3)
X <- expand.grid(x1 = x1, x2 = x2)
y <- camel6(X)
dy <- camel6Grad(X)
dat <- data.frame(X, y)
deri <- data.frame(dy)
# Fit (gradient-enhanced) Kriging model
km.2d <- gekm(y ~ 1, data = dat, covtype = "gaussian", optimizer = "L-BFGS-B")
gekm.2d <- gekm(y ~ 1, data = dat, deriv = deri, covtype = "gaussian", optimizer = "L-BFGS-B")
# Compute negative 'concentrated' log-likelihood values
n.grid <- 20
theta1.grid <- seq(0.5, 4, length = n.grid)
theta2.grid <- seq(0.5, 2, length = n.grid)
params <- expand.grid(theta1 = theta1.grid, theta2 = theta2.grid)
logLik.km.2d <- apply(params, 1, logLikFun, km.2d)
logLik.gekm.2d <- apply(params, 1, logLikFun, gekm.2d)
logLikGrad.km.2d <- t(apply(params, 1, logLikGrad, km.2d))
logLikGrad.gekm.2d <- t(apply(params, 1, logLikGrad, gekm.2d))
# Plot negative 'concentrated' log-likelihood
par(mfrow = c(1, 2), oma = c(3.6, 3.5, 1.5, 0.2), mar = c(0, 0, 1.5, 0))
contour(theta1.grid, theta2.grid, matrix(logLik.km.2d, nrow = n.grid, ncol = n.grid),
nlevels = 50, main = "Kriging")
vectorfield(params, logLikGrad.km.2d, col = 4, lwd = 2, length = 0.1)
points(km.2d$theta[1], km.2d$theta[2], col = "red", pch = 16)
contour(theta1.grid, theta2.grid, matrix(logLik.gekm.2d, nrow = n.grid, ncol = n.grid),
nlevels = 50, main = "GEK", yaxt = "n")
points(gekm.2d$theta[1], gekm.2d$theta[2], col = "red", pch = 16)
vectorfield(params, logLikGrad.gekm.2d, col = 4, lwd = 2, length = 0.1)
title(main = "Negative 'concentrated' log-likelihood", outer = TRUE)
mtext(side = 1, outer = TRUE, line = 2.5, expression(theta[1]))
mtext(side = 2, outer = TRUE, line = 2.5, expression(theta[2]))
Run the code above in your browser using DataLab