Fit.Kriging: Fit or Create Kriging Models

Description

This function provides a common interface to fit various kriging models from the data. OK (Ordinary.Kriging), UK (Universal.Kriging), LK (Limit.Kriging), RK (Rational.Kriging), and GRK (Generalized.Rational.Kriging) are supported in this package.

Usage

Fit.Kriging(
  X,
  y,
  interpolation = TRUE,
  fit = TRUE,
  model = "OK",
  model.parameters = list(),
  kernel = NULL,
  kernel.parameters = list(),
  nlopt.parameters = list()
)

Value

A Kriging Class Object.

Arguments

X: a matrix for input (feature)
y: a vector for output (target), only one-dimensional output is supported
interpolation: whether to interpolate, for noisy data please set interpolate=FALSE
fit: whether to fit the length scale parameters from data
model: choice of kriging model: OK, UK, LK, RK, GRK
model.parameters: a list of parameters required for the specific kriging model, e.g. basis functions for universal kriging
kernel: a kernel class object
kernel.parameters: a list of parameters required for the kernel, if no kernel class object is provided
nlopt.parameters: a list of parameters required for NLopt, including choice of optimization algorithm and maximum number of evaluation

Author

Chaofan Huang and V. Roshan Joseph

Details

Kriging gives the best linear unbiased predictor given the data. It can also be given a Bayesian interpretation by assuming a Gaussian process (GP) prior for the underlying function. Please see Santner et al. (2003), Rasmussen and Williams (2006), and Joseph (2024) for details.

For data from deterministic computer experiments, use interpolation=TRUE and will give an interpolator. For noisy data, use interpolation=FALSE, which will give an approximator of the underlying function. Currently, approximation is supported for OK (Ordinary.Kriging) and UK (Universal.Kriging).

The kernel choices are required and can be specified by (i) providing the kernel class object to kernel or (ii) specifying the kernel type and other parameters in kernel.parameters. Please see examples section for detail usages.

When the lengthscale / correlation parameters are known, simply provide them in (i) the kernel class object to kernel or (ii) in kernel.parameters, and set fit=FALSE. The kriging model will be fitted with the user provided parameters.

When the lengthscale / correlation parameters are unknown, they can be estimated via Maximum Likelihood method by setting fit=TRUE. The initial / lower bound / upper bound of the lengthscale parameters can be provided in kernel.parameters, otherwise a good initial and range would be estimated from the data. The optimization is performed via NLopt, a open-source library for nonlinear optimization. All gradient-free optimization methods in NLopt are supported and can be specified in nlopt.parameters. See nloptr::nloptr.print.options() for the list of available derivative-free algorithms (prefix with NLOPT_GN or NLOPT_LN). The maximum number of optimization steps can also be defined in nlopt.parameters. Please see examples section for detail usages.

References

Joseph, V. R. (2006). Limit kriging. Technometrics, 48(4), 458-466.

Joseph, V. R. (2024). Rational Kriging. Journal of the American Statistical Association.

Rasmussen, C. E. & Williams, C. K. (2006). Gaussian Processes for Machine Learning. The MIT Press.

Santner, T. J., Williams, B. J., Notz, W. I., & Williams, B. J. (2003). The design and analysis of computer experiments (Vol. 1). New York: Springer.

Examples

Run this code

# one dimensional example 
f <- function(x) {
  x <- 0.5 + 2*x
  y <- sin(10*pi*x)/(2*x) + (x-1)^4
  return (y)
}

set.seed(1234)
# train set
n <- 30
p <- 1
X <- matrix(runif(n),ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0,1,length=1001), ncol=p)

#############################################################################
################ Minimal Example for Fitting a Kriging Model ################
#############################################################################
# Ordinary Kriging
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="OK",
                       kernel.parameters=list(type="Gaussian"))
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

# Universal Kriging
basis.function <- function(x) {c(1,x[1],x[1]^2)}
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="UK",
                       model.parameters=list(basis.function=basis.function),
                       kernel.parameters=list(type="Gaussian"))
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

# Limit Kriging
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="LK",
                       kernel.parameters=list(type="Gaussian"))
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

# Rational Kriging
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="RK",
                       kernel.parameters=list(type="RQ",alpha=1))
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

# Generalized Rational Kriging
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="GRK",
                       kernel.parameters=list(type="RQ",alpha=1))
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

#############################################################################
################ Fitting a Kriging Model with Kernel Object #################
#############################################################################
kernel <- Gaussian.Kernel(rep(1,p))
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="OK", kernel=kernel)
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

#############################################################################
############### Creating a Kriging Model with Kernel Object #################
#############################################################################
# set fit = FALSE to create Kriging model with user provided kernel
# no optimization for the length scale parameters 
kernel <- Gaussian.Kernel(rep(1e-1,p))
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=FALSE, model="OK", kernel=kernel)
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3) 
Get.Kriging.Parameters(kriging)

#############################################################################
############ Fitting a Kriging Model with Range for Lengthscale #############
#############################################################################
kernel.parameters <- list(
    type = 'Gaussian',
    lengthscale = rep(1,p), # initial value
    lengthscale.lower.bound = rep(1e-3,p), # lower bound
    lengthscale.upper.bound = rep(1e3,p) # upper bound
) # if not provided, a good estimate would be computed from data
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="OK",
                       kernel.parameters=kernel.parameters)
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

#############################################################################
########### Fitting a Kriging Model with Different Optimization #############
#############################################################################
nlopt.parameters <- list(
    algorithm = 'NLOPT_GN_DIRECT', # optimization method
    maxeval = 250 # maximum number of evaluation
)
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="OK",
                       kernel.parameters=list(type="Gaussian"),
                       nlopt.parameters=nlopt.parameters)
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

# Global-Local Optimization
nlopt.parameters <- list(
    algorithm = 'NLOPT_GN_MLSL_LDS', # optimization method
    local.algorithm = 'NLOPT_LN_SBPLX', # local algorithm
    maxeval = 250 # maximum number of evaluation
)
kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="OK",
                       kernel.parameters=list(type="Gaussian"),
                       nlopt.parameters=nlopt.parameters)
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

#############################################################################
################# Fitting a Kriging Model from Noisy Data ###################
#############################################################################
y <- y + 0.1 * rnorm(length(y))
kriging <- Fit.Kriging(X, y, interpolation=FALSE, fit=TRUE, model="OK",
                       kernel.parameters=list(type="Gaussian"))
pred <- Predict.Kriging(kriging, newX)
plot(newX, f(newX), "l")
points(X, y, pch=16, col="blue")
lines(newX, pred$mean, col="red", lty=2)
lines(newX, pred$mean-2*pred$sd, col="red", lty=3)
lines(newX, pred$mean+2*pred$sd, col="red", lty=3)
Get.Kriging.Parameters(kriging)

Run the code above in your browser using DataLab