Fit.Kriging {rkriging} | R Documentation |
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()
)
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 |
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 |
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.
Value
A Kriging Class Object.
Author(s)
Chaofan Huang and V. Roshan Joseph
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.
See Also
Predict.Kriging, Get.Kriging.Parameters, Get.Kernel, Ordinary.Kriging, Universal.Kriging, Limit.Kriging, Rational.Kriging, Generalized.Rational.Kriging.
Examples
# 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)