Package {rkriging}

ALM Sequential Design

Description

This function creates an Active Learning MacKay (ALM) sequential design object for emulation (space-filling based on predictive uncertainty).

Usage

ALM.Sequential(
  design.parameters,
  kriging = NULL,
  kriging.parameters = NULL,
  nlopt.parameters = list()
)

Arguments

Details

ALM selects the next design point as the one that maximizes the predictive standard deviation:

x_{next} = \arg\max_{x} \sigma(x),

where \sigma(x) is the kriging predictive standard deviation at x. This criterion is suitable for emulation tasks where the goal is to reduce overall prediction uncertainty.

The design space is specified via design.parameters, which must contain design.lower.bound and design.upper.bound vectors. During optimization, candidate points can be provided via design.candidates in Get.Next.Design. If the kriging object was fitted with scale=TRUE, design bounds, candidates, new design points, predictions, returned designs, and returned data are supplied and returned on the original input scale; scaling is handled internally for the sequential criterion.

Value

An ALM Sequential Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

MacKay, D. J. C. (1992). Information-based objective functions for active data selection. Neural computation, 4(4), 590-604.

See Also

Predict.Sequential, Get.Next.Design, Add.Next.Design.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- ALM.Sequential(design.parameters, kriging)

# sequential design loop with visualization
n.iterations <- 5
for (i in 1:n.iterations) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("ALM Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="ALM Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)

Add Next Design to Sequential Object

Description

This function adds a new design point and its corresponding output to the sequential design object and updates kriging (optionally refit the kriging hyperparameters).

Usage

Add.Next.Design(sequential, new.design, refit.kriging = TRUE)

Arguments

Details

The new design is provided as a list containing:

X

a matrix or vector of new design point(s); if a vector, it will be converted to a 1×p matrix

y

a vector of corresponding output value(s)

This function is typically used in a sequential design loop together with Get.Next.Design, allowing fine-grained control over each iteration. If the sequential object was created from a scaled kriging object, new.design$X should be provided on the original input scale and is scaled internally before updating the kriging model.

Value

The sequential object (modified in-place).

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

Get.Next.Design, Get.Sequential.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- EI.Sequential(design.parameters, kriging)

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("EI Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="EI Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)
best.design <- Get.Best.Design(sequential)
points(best.design$design, best.design$objective, pch=16, col="green")

EI Sequential Design

Description

This function creates an Expected Improvement (EI) sequential design object for optimization (finding the minimum of the objective function).

Usage

EI.Sequential(
  design.parameters,
  kriging = NULL,
  kriging.parameters = NULL,
  nlopt.parameters = list()
)

Arguments

Details

EI selects the next design point as the one that maximizes the expected improvement:

x_{next} = \arg\max_{x} \textrm{EI}(x) = \sigma(x)\left[u\Phi(u)+\phi(u)\right],

where u = (y_{min}-\mu(x))/\sigma(x), \mu(x) and \sigma(x) are the kriging predictive mean and standard deviation, and \Phi and \phi are the standard normal CDF and PDF. This is the most popular acquisition function for Bayesian optimization.

The design space is specified via design.parameters, which must contain design.lower.bound and design.upper.bound vectors. During optimization, candidate points can be provided via design.candidates in Get.Next.Design. If the kriging object was fitted with scale=TRUE, design bounds, candidates, new design points, predictions, returned designs, and returned data are supplied and returned on the original input scale; scaling is handled internally for the sequential criterion.

Value

An EI Sequential Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Jones, D. R., Schonlau, M., & Welch, W. J. (1998). Efficient global optimization of expensive black-box functions. Journal of Global optimization, 13(4), 455-492.

See Also

PI.Sequential, UCB.Sequential, UDF.Sequential, Get.Next.Design, Add.Next.Design, Get.Best.Design.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- EI.Sequential(design.parameters, kriging)

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("EI Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="EI Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)
best.design <- Get.Best.Design(sequential)
points(best.design$design, best.design$objective, pch=16, col="green")

Evaluate Kernel

Description

This function computes the kernel (correlation) matrix. Given the kernel class object and the input data X of size n, this function computes the corresponding n\times n kernel (correlation) matrix.

Usage

Evaluate.Kernel(kernel, X)

Arguments

Value

The kernel (correlation) matrix of X evaluated by the kernel.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

Get.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

kernel <- Gaussian.Kernel(lengthscale)
Evaluate.Kernel(kernel, X)

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,
  nugget = NULL,
  model = "OK",
  model.parameters = list(),
  kernel = NULL,
  kernel.parameters = list(),
  nlopt.parameters = list(),
  scale = FALSE
)

Arguments

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), Joseph (2025), and Joseph (2026) 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. When scale=TRUE, the returned object stores the training input minima and maxima as scale_min and scale_max attributes.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Santner, T. J., Williams, B. J., & Notz, W. I. (2003). The Design and Analysis of Computer Experiments. Springer.

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

Joseph, V. R. (2025). Rational Kriging. Journal of the American Statistical Association. 120, 548-558.

Joseph, V. R. (2026). Experimental Design for Data Science and Engineering. Chapman & Hall/CRC, doi:10.1201/9781003661535.

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

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)

Gaussian Kernel

Description

This function specifies the Gaussian / Squared Exponential (SE) / Radial Basis Function (RBF) kernel.

Usage

Gaussian.Kernel(lengthscale)

Arguments

Details

The Gaussian kernel is given by

k(r)=\exp(-r^2/2),

where

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Gaussian Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- Gaussian.Kernel(lengthscale)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="Gaussian")
Evaluate.Kernel(kernel, X) 

Generalized Rational Kriging

Description

This functions fits the generalized rational kriging model to the data.

Usage

Generalized.Rational.Kriging(
  X,
  y,
  fit = TRUE,
  nugget = NULL,
  kernel = NULL,
  kernel.parameters = list(),
  nlopt.parameters = list()
)

Arguments

Details

Ordinary kriging and rational kriging can be obtained as special cases of generalized rational kriging. Please see Joseph (2025) for details. The Spectra library is used for fast computation of the first eigenvalues/vectors. Only interpolation is available. Noisy output is not supported for generalized rational 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 of Fit.Kriging for detail usages.

When the lengthscale / correlation parameters are unknown, all parameters including the constant mean 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 library, 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 of Fit.Kriging for detail usages.

Value

A Generalized Rational Kriging Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Joseph, V. R. (2025). Rational Kriging. Journal of the American Statistical Association. 120, 548-558.

Qiu, Y., Guennebaud, G., & Niesen, J. (2015). Spectra: C++ library for large scale eigenvalue problems.

See Also

Fit.Kriging, Predict.Kriging, Get.Kriging.Parameters.

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)

# approach 1
kriging <- Generalized.Rational.Kriging(X, y, fit=TRUE, 
                                        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)

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

Get Best Design from Optimization Sequential

Description

This function returns the best (minimum) objective value and its corresponding design from an optimization sequential design object (e.g., EI.Sequential, PI.Sequential, UCB.Sequential).

Usage

Get.Best.Design(sequential)

Arguments

Details

If the sequential object was created from a scaled kriging object, the returned best design is converted back to the original input scale.

Value

A list with components:

design

the design (input) corresponding to the best objective

objective

the best (minimum) objective value

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

EI.Sequential, PI.Sequential, UCB.Sequential.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# fit kriging
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- EI.Sequential(design.parameters, kriging)

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("EI Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="EI Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)
best.design <- Get.Best.Design(sequential)
points(best.design$design, best.design$objective, pch=16, col="green")

Get Data from Sequential

Description

This function retrieves the current training data (X and y) from a sequential design object.

Usage

Get.Data(sequential)

Arguments

Details

The sequential design object maintains the current training data internally. This function provides access to both the design matrix X and the corresponding output vector y. If the sequential object was created from a scaled kriging object, X is converted back to the original input scale.

Value

A list with components:

X

the design matrix (training inputs)

y

the output vector (training responses)

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

Get.Sequential.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- EI.Sequential(design.parameters, kriging)
Get.Data(sequential)

Kernel

Description

This function provides a common interface to specify various kernels. See arguments section for the available kernels in this package.

Usage

Get.Kernel(lengthscale, type, parameters = list())

Arguments

Value

A Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

Evaluate.Kernel, Gaussian.Kernel, RQ.Kernel, IM.Kernel, Matern12.Kernel, Matern32.Kernel, Matern52.Kernel Matern.Kernel, UDF.Kernel, MultiplicativeRQ.Kernel, MultiplicativeIM.Kernel, MultiplicativeMatern.Kernel, MultiplicativeUDF.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# Gaussian 
kernel <- Get.Kernel(lengthscale, type="Gaussian")
Evaluate.Kernel(kernel, X)

# Rational Quadratic (RQ)
kernel <- Get.Kernel(lengthscale, type="RQ", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

# Inverse Multiquadric (IM)
kernel <- Get.Kernel(lengthscale, type="IM", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

# Matern(1/2)
kernel <- Get.Kernel(lengthscale, type="Matern12")
Evaluate.Kernel(kernel, X) 

# Matern(3/2)
kernel <- Get.Kernel(lengthscale, type="Matern32")
Evaluate.Kernel(kernel, X) 

# Matern(5/2)
kernel <- Get.Kernel(lengthscale, type="Matern52")
Evaluate.Kernel(kernel, X) 

# Generalized Matern
kernel <- Get.Kernel(lengthscale, type="Matern", parameters=list(nu=2.01))
Evaluate.Kernel(kernel, X) 

# User Defined Function (UDF) Kernel
kernel.function <- function(sqdist) {return (exp(-sqrt(sqdist)))} 
kernel <- Get.Kernel(lengthscale, type="UDF", 
                     parameters=list(kernel.function=kernel.function))
Evaluate.Kernel(kernel, X) 

# Multiplicative Rational Quadratic (RQ)
kernel <- Get.Kernel(lengthscale, type="MultiplicativeRQ", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

# Multiplicative Inverse Multiquadric (IM)
kernel <- Get.Kernel(lengthscale, type="MultiplicativeIM", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

# Multiplicative Generalized Matern
kernel <- Get.Kernel(lengthscale, type="MultiplicativeMatern", parameters=list(nu=2.01))
Evaluate.Kernel(kernel, X)

# Multiplicative User Defined Function (UDF)
kernel.function <- function(sqdist) {return (exp(-sqrt(sqdist)))} 
kernel <- Get.Kernel(lengthscale, type="MultiplicativeUDF", 
                     parameters=list(kernel.function=kernel.function))
Evaluate.Kernel(kernel, X) 

Get Kriging Parameters

Description

This function can be used for extracting the estimates of the kriging parameters.

Usage

Get.Kriging.Parameters(kriging)

Arguments

Value

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

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

Joseph, V. R. (2025). Rational Kriging. Journal of the American Statistical Association. 120, 548-558.

Joseph, V. R. (2026). Experimental Design for Data Science and Engineering. Chapman & Hall/CRC, doi:10.1201/9781003661535.

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

See Also

Fit.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)

kriging <- Fit.Kriging(X, y, interpolation=TRUE, fit=TRUE, model="OK",
                       kernel.parameters=list(type="Gaussian"))
Get.Kriging.Parameters(kriging)

Get Next Design from Sequential

Description

This function selects the next design point using the sequential criterion.

Usage

Get.Next.Design(sequential, design.candidates = NULL)

Arguments

Details

If design.candidates is provided, the next design is selected from the candidate set by evaluating the criterion at each candidate point. Otherwise, NLopt optimization is used within the design space bounds specified in design.parameters. If the sequential object was created from a scaled kriging object, design.candidates should be provided on the original input scale and the returned next design is also on the original input scale.

Value

A vector representing the next design point.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

Get.Sequential, Add.Next.Design.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# fit kriging
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- EI.Sequential(design.parameters, kriging)

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("EI Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="EI Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)
best.design <- Get.Best.Design(sequential)
points(best.design$design, best.design$objective, pch=16, col="green")

# alternatively can pass a set of candidate points for finding the next design
design.candidates <- matrix(seq(0,1,length=1001), ncol=p)
Get.Next.Design(sequential, design.candidates=design.candidates)

Get Sequential Design

Description

This function provides a common interface to create various sequential design objects. ALM (ALM.Sequential), EI (EI.Sequential), PI (PI.Sequential), UCB (UCB.Sequential), and UDF (UDF.Sequential) are supported.

Usage

Get.Sequential(
  design.parameters,
  kriging = NULL,
  kriging.parameters = NULL,
  method = "ALM",
  method.parameters = list(),
  nlopt.parameters = list()
)

Arguments

Details

Sequential design methods select the next design point iteratively based on an acquisition criterion. For emulation tasks (space-filling), use method="ALM" which maximizes predictive uncertainty. For optimization tasks (finding the minimum), use method="EI", "PI", or "UCB". A custom criterion can be specified via method="UDF".

The design space is specified via design.parameters, which must contain design.lower.bound and design.upper.bound vectors. During optimization, candidate points can be provided via design.candidates in Get.Next.Design. If the kriging object was fitted with scale=TRUE, design bounds, candidates, new design points, predictions, returned designs, and returned data are supplied and returned on the original input scale; scaling is handled internally for the sequential criterion. For method="UDF", the user-defined criterion function receives the current training X on the original input scale.

The kriging model can be provided either as:

An existing kriging object

e.g., from Ordinary.Kriging or Fit.Kriging

Kriging parameters

as a list to pass to Fit.Kriging, must contain at least X and y

Value

A Sequential Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

ALM.Sequential, EI.Sequential, PI.Sequential, UCB.Sequential, UDF.Sequential, Predict.Sequential, Get.Next.Design, Add.Next.Design.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- Get.Sequential(design.parameters, kriging, method="EI")

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("EI Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="EI Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)
best.design <- Get.Best.Design(sequential)
points(best.design$design, best.design$objective, pch=16, col="green")

# alternatively can create sequential by passing kriging parameters
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
kriging.parameters <- list(X=X, y=y, interpolation=TRUE, fit=TRUE,
                           kernel.parameters=list(type="Gaussian"))
sequential <- Get.Sequential(design.parameters, kriging.parameters=kriging.parameters, method="EI")

Inverse Multiquadric (IM) Kernel

Description

This function specifies the Inverse Multiquadric (IM) kernel.

Usage

IM.Kernel(lengthscale, alpha = 1)

Arguments

Details

The Inverse Multiquadric (IM) kernel is given by

k(r;\alpha)=\left(1+\frac{r^2}{2}\right)^{-\alpha},

where \alpha is a positive shape parameter and

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

An Inverse Multiquadric (IM) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

Hardy, R. L. & Göpfert, W. M. (1975). Least squares prediction of gravity anomalies, geoidal undulations, and deflections of the vertical with multiquadric harmonic functions. Geophysical Research Letters, 2(10), 423–426. doi:10.1029/GL002i010p00423

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

See Also

MultiplicativeIM.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- IM.Kernel(lengthscale, alpha=1)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="IM", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

Limit Kriging

Description

This functions fits the limit kriging model to the data.

Usage

Limit.Kriging(
  X,
  y,
  fit = TRUE,
  nugget = NULL,
  kernel = NULL,
  kernel.parameters = list(),
  nlopt.parameters = list()
)

Arguments

Details

Limit kriging avoids the mean reversion issue of ordinary kriging. Please see Joseph (2006) for details. Only interpolation is available. Noisy output is not supported for limit 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 of Fit.Kriging for detail usages.

When the lengthscale / correlation parameters are unknown, all parameters including the constant mean 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 of Fit.Kriging for detail usages.

Value

A Limit Kriging Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

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

See Also

Fit.Kriging, Predict.Kriging, Get.Kriging.Parameters.

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)

# approach 1
kriging <- Limit.Kriging(X, y, fit=TRUE, 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)

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

Generalized Matern Kernel

Description

This function specifies the (Generalized) Matern kernel with any smoothness parameter \nu.

Usage

Matern.Kernel(lengthscale, nu = 2.01)

Arguments

Details

The Generalized Matern kernel is given by

k(r;\nu)=\frac{2^{1-\nu}}{\Gamma(\nu)}(\sqrt{2\nu}r)^{\nu}K_{\nu}(\sqrt{2\nu}r),

where \nu is the smoothness parameter, K_{\nu} is the modified Bessel function, \Gamma is the gamma function, and

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s. As \nu\to\infty, it converges to the Gaussian.Kernel.

Value

A Generalized Matern Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

Matern12.Kernel, Matern32.Kernel, Matern52.Kernel, MultiplicativeMatern.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- Matern.Kernel(lengthscale, nu=2.01)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="Matern", parameters=list(nu=2.01))
Evaluate.Kernel(kernel, X) 

Matern(1/2) Kernel

Description

This function specifies the Matern kernel with smoothness parameter \nu=1/2. It is also known as the Exponential kernel.

Usage

Matern12.Kernel(lengthscale)

Arguments

Details

The Matern(1/2) kernel is given by

k(r)=\exp(-r),

where

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Matern(1/2) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

Matern.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- Matern12.Kernel(lengthscale)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="Matern12")
Evaluate.Kernel(kernel, X) 

Matern(3/2) Kernel

Description

This function specifies the Matern kernel with smoothness parameter \nu=3/2.

Usage

Matern32.Kernel(lengthscale)

Arguments

Details

The Matern(3/2) kernel is given by

k(r)=(1+\sqrt{3}r)\exp(-\sqrt{3}r),

where

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Matern(3/2) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

Matern.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- Matern32.Kernel(lengthscale)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="Matern32")
Evaluate.Kernel(kernel, X) 

Matern(5/2) Kernel

Description

This function specifies the Matern kernel with smoothness parameter \nu=5/2.

Usage

Matern52.Kernel(lengthscale)

Arguments

Details

The Matern(5/2) kernel is given by

k(r)=(1+\sqrt{5}r+5r^2/3)\exp(-\sqrt{5}r),

where

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Matern(5/2) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

Matern.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- Matern52.Kernel(lengthscale)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="Matern52")
Evaluate.Kernel(kernel, X) 

Multiplicative Inverse Multiquadric (IM) Kernel

Description

This function specifies the Multiplicative Inverse Multiquadric (IM) kernel.

Usage

MultiplicativeIM.Kernel(lengthscale, alpha = 1)

Arguments

Details

The Multiplicative Inverse Multiquadric (IM) kernel is given by

k(r;\alpha)=\prod_{i=1}^{p}\left(1+\frac{r_{i}^2}{2}\right)^{-\alpha},

where \alpha is a positive shape parameter and

r_{i}(x,x^{\prime})=\sqrt{\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the dimension-wise euclidean distances between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Multiplicative Inverse Multiquadric (IM) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

Song, D. & Joseph, V. R. (2026). Efficient Active Learning Strategies for Computer Experiments. Technometrics, 68(1), 65–78. doi:10.1080/00401706.2025.2546366

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

See Also

IM.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- MultiplicativeIM.Kernel(lengthscale, alpha=1)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="MultiplicativeIM", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

Multiplicative Generalized Matern Kernel

Description

This function specifies the Multiplicative Generalized Matern kernel.

Usage

MultiplicativeMatern.Kernel(lengthscale, nu = 2.01)

Arguments

Details

The Multiplicative Generalized Matern kernel is given by

k(r;\nu)=\prod_{i=1}^{p}\frac{2^{1-\nu}}{\Gamma(\nu)}(\sqrt{2\nu}r_{i})^{\nu}K_{\nu}(\sqrt{2\nu}r_{i}),

where \nu is the smoothness parameter, K_{\nu} is the modified Bessel function, \Gamma is the gamma function, and

r_{i}(x,x^{\prime})=\sqrt{\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the dimension-wise euclidean distances between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Multiplicative Generalized Matern Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

Matern.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- MultiplicativeMatern.Kernel(lengthscale, nu=2.01)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="MultiplicativeMatern", parameters=list(nu=2.01))
Evaluate.Kernel(kernel, X) 

Multiplicative Rational Quadratic (RQ) Kernel

Description

This function specifies the Multiplicative Rational Quadratic (RQ) kernel.

Usage

MultiplicativeRQ.Kernel(lengthscale, alpha = 1)

Arguments

Details

The Multiplicative Rational Quadratic (RQ) kernel is given by

k(r;\alpha)=\prod_{i=1}^{p}\left(1+\frac{r_{i}^2}{2\alpha}\right)^{-\alpha},

where \alpha is the scale mixture parameter and

r_{i}(x,x^{\prime})=\sqrt{\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the dimension-wise euclidean distances between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Multiplicative Rational Quadratic (RQ) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Song, D. and Joseph, V. R. (2026). "Efficient Active Learning Strategies for Computer Experiments." Technometrics, 68, 65-78.

See Also

RQ.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- MultiplicativeRQ.Kernel(lengthscale, alpha=1)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="MultiplicativeRQ", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

Multiplicative User Defined Function (UDF) Kernel

Description

This function specifies the Multiplicative kernel with the user defined R function.

Usage

MultiplicativeUDF.Kernel(lengthscale, kernel.function)

Arguments

Details

The Multiplicative User Defined Function (UDF) kernel is given by

k(r)=\prod_{i=1}^{p}f(r_{i}^2),

where f is the user defined kernel function that takes r_{i}^2 as input, where

r_{i}(x,x^{\prime})=\sqrt{\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the dimension-wise euclidean distances between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A Multiplicative User Defined Function (UDF) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

UDF.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

kernel.function <- function(sqdist) {return (exp(-sqrt(sqdist)))} 

# approach 1
kernel <- MultiplicativeUDF.Kernel(lengthscale, kernel.function=kernel.function)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="MultiplicativeUDF", 
                     parameters=list(kernel.function=kernel.function))
Evaluate.Kernel(kernel, X) 

Ordinary Kriging

Description

This functions fits the ordinary kriging model to the data.

Usage

Ordinary.Kriging(
  X,
  y,
  interpolation = TRUE,
  fit = TRUE,
  nugget = NULL,
  kernel = NULL,
  kernel.parameters = list(),
  nlopt.parameters = list()
)

Arguments

Details

Ordinary kriging assumes a constant mean. Please see Joseph (2026) 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.

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 of Fit.Kriging for detail usages.

When the lengthscale / correlation parameters are unknown, all parameters including the constant mean 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 of Fit.Kriging for detail usages.

Value

A Ordinary Kriging Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Joseph, V. R. (2026). Experimental Design for Data Science and Engineering. Chapman & Hall/CRC, doi:10.1201/9781003661535.

See Also

Fit.Kriging, Predict.Kriging, Get.Kriging.Parameters.

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)

# approach 1
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE, 
                            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)

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

PI Sequential Design

Description

This function creates a Probability of Improvement (PI) sequential design object for optimization (finding the minimum of the objective function).

Usage

PI.Sequential(
  design.parameters,
  kriging = NULL,
  kriging.parameters = NULL,
  nlopt.parameters = list()
)

Arguments

Details

PI selects the next design point as the one that maximizes the probability of improvement:

x_{next} = \arg\max_{x} \textrm{PI}(x) = \Phi(u),

where u = (y_{min}-\mu(x))/\sigma(x), \mu(x) and \sigma(x) are the kriging predictive mean and standard deviation, and \Phi is the standard normal CDF.

The design space is specified via design.parameters, which must contain design.lower.bound and design.upper.bound vectors. During optimization, candidate points can be provided via design.candidates in Get.Next.Design. If the kriging object was fitted with scale=TRUE, design bounds, candidates, new design points, predictions, returned designs, and returned data are supplied and returned on the original input scale; scaling is handled internally for the sequential criterion.

Value

A PI Sequential Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Kushner, H. J. (1964). A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise. Journal of Basic Engineering, 86(1), 97-106.

See Also

EI.Sequential, UCB.Sequential, UDF.Sequential, Get.Next.Design, Add.Next.Design, Get.Best.Design.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- PI.Sequential(design.parameters, kriging)

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("PI Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization and best design
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="PI Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)
best.design <- Get.Best.Design(sequential)
points(best.design$design, best.design$objective, pch=16, col="green")

Kriging Prediction

Description

This function gives prediction and uncertainty quantification of the kriging model on a new input.

Usage

Predict.Kriging(kriging, X, Neumann = FALSE)

Arguments

Details

If kriging was fitted with scale=TRUE in Fit.Kriging, X is scaled to [0, 1] using the training minima and maxima stored in the fitted object before prediction.

Value

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

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

Joseph, V. R. (2025). Rational Kriging. Journal of the American Statistical Association. 120, 548-558.

Joseph, V. R. (2026). Experimental Design for Data Science and Engineering. Chapman & Hall/CRC, doi:10.1201/9781003661535.

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

See Also

Fit.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)

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)

Predict from Sequential Design

Description

This function computes kriging predictions from the sequential design object at new input locations.

Usage

Predict.Sequential(sequential, newX)

Arguments

Details

If the sequential object was created from a scaled kriging object, newX should be provided on the original input scale and is scaled internally before prediction.

Value

A list with components:

mean

the predictive mean

sd

the predictive standard deviation

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

Get.Sequential.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- EI.Sequential(design.parameters, kriging)
# get prediction and visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="EI Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)

Rational Quadratic (RQ) Kernel

Description

This function specifies the Rational Quadratic (RQ) kernel.

Usage

RQ.Kernel(lengthscale, alpha = 1)

Arguments

Details

The Rational Quadratic (RQ) kernel is given by

k(r;\alpha)=\left(1+\frac{r^2}{2\alpha}\right)^{-\alpha},

where \alpha is the scale mixture parameter and

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2}

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s. As \alpha\to\infty, it converges to the Gaussian.Kernel.

Value

A Rational Quadratic (RQ) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

MultiplicativeRQ.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

# approach 1
kernel <- RQ.Kernel(lengthscale, alpha=1)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="RQ", parameters=list(alpha=1))
Evaluate.Kernel(kernel, X) 

Rational Kriging

Description

This functions fits the rational kriging model to the data.

Usage

Rational.Kriging(
  X,
  y,
  fit = TRUE,
  nugget = NULL,
  kernel = NULL,
  kernel.parameters = list(),
  nlopt.parameters = list()
)

Arguments

Details

Rational kriging gives a rational predictor in terms of the inputs, which gives a more stable estimate of the mean. Please see Joseph (2025) for details. Only interpolation is available. Noisy output is not supported for rational 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 of Fit.Kriging for detail usages.

When the lengthscale / correlation parameters are unknown, all parameters including the constant mean 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 of Fit.Kriging for detail usages.

Value

A Rational Kriging Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Joseph, V. R. (2025). Rational Kriging. Journal of the American Statistical Association. 120, 548-558.

See Also

Fit.Kriging, Predict.Kriging, Get.Kriging.Parameters.

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)

# approach 1
kriging <- Rational.Kriging(X, y, fit=TRUE, 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)

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

UCB Sequential Design

Description

This function creates an Upper Confidence Bound (UCB) sequential design object for optimization (finding the minimum of the objective function).

Usage

UCB.Sequential(
  design.parameters,
  kriging = NULL,
  kriging.parameters = NULL,
  lambda = 2,
  nlopt.parameters = list()
)

Arguments

Details

UCB selects the next design point as the one that minimizes the lower confidence bound:

x_{next} = \arg\min_{x} \mu(x) - \lambda\sigma(x),

where \mu(x) and \sigma(x) are the kriging predictive mean and standard deviation, and \lambda controls the exploration-exploitation tradeoff. Larger \lambda encourages more exploration.

The design space is specified via design.parameters, which must contain design.lower.bound and design.upper.bound vectors. During optimization, candidate points can be provided via design.candidates in Get.Next.Design. If the kriging object was fitted with scale=TRUE, design bounds, candidates, new design points, predictions, returned designs, and returned data are supplied and returned on the original input scale; scaling is handled internally for the sequential criterion.

Value

A UCB Sequential Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Srinivas, N., Krause, A., Kakade, S. M., & Seeger, M. (2009). Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995.

See Also

EI.Sequential, PI.Sequential, UDF.Sequential, Get.Next.Design, Add.Next.Design, Get.Best.Design.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))
sequential <- UCB.Sequential(design.parameters, kriging, lambda=2)

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("UCB Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization and best design
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="UCB Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)
best.design <- Get.Best.Design(sequential)
points(best.design$design, best.design$objective, pch=16, col="green")

User Defined Function (UDF) Kernel

Description

This function specifies a kernel with the user defined R function.

Usage

UDF.Kernel(lengthscale, kernel.function)

Arguments

Details

The User Defined Function (UDF) kernel is given by

k(r) = f(r^2)

where f is the user defined kernel function that takes r^2 as input, where

r(x,x^{\prime})=\sqrt{\sum_{i=1}^{p}\left(\frac{x_{i}-x_{i}^{\prime}}{l_{i}}\right)^2},

is the euclidean distance between x and x^{\prime} weighted by the length scale parameters l_{i}'s.

Value

A User Defined Function (UDF) Kernel Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Duvenaud, D. (2014). The kernel cookbook: Advice on covariance functions.

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

See Also

MultiplicativeUDF.Kernel, Get.Kernel, Evaluate.Kernel.

Examples

n <- 5
p <- 3
X <- matrix(rnorm(n*p), ncol=p)
lengthscale <- c(1:p)

kernel.function <- function(sqdist) {return (exp(-sqrt(sqdist)))} 

# approach 1
kernel <- UDF.Kernel(lengthscale, kernel.function=kernel.function)
Evaluate.Kernel(kernel, X)

# approach 2
kernel <- Get.Kernel(lengthscale, type="UDF", 
                     parameters=list(kernel.function=kernel.function))
Evaluate.Kernel(kernel, X) 

UDF Sequential Design

Description

This function creates a User Defined Function (UDF) sequential design object with a custom acquisition function.

Usage

UDF.Sequential(
  design.parameters,
  kriging = NULL,
  kriging.parameters = NULL,
  criterion.function,
  criterion.parameters = list(),
  nlopt.parameters = list()
)

Arguments

Details

The UDF sequential design allows users to specify a custom criterion function for selecting the next design point. The criterion function receives the predictive mean, predictive standard deviation, the current training data, and any additional user-defined parameters. If the kriging object was fitted with scale=TRUE, the criterion function receives the current training X on the original input scale.

The design space is specified via design.parameters, which must contain design.lower.bound and design.upper.bound vectors. During optimization, candidate points can be provided via design.candidates in Get.Next.Design. If the kriging object was fitted with scale=TRUE, design bounds, candidates, new design points, predictions, returned designs, and returned data are supplied and returned on the original input scale; scaling is handled internally for the sequential criterion.

Value

A UDF Sequential Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

See Also

ALM.Sequential, EI.Sequential, PI.Sequential, UCB.Sequential, Predict.Sequential, Get.Next.Design, Add.Next.Design.

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)
n <- 10
p <- 1
X <- matrix(runif(n), ncol=p)
y <- apply(X, 1, f)
newX <- matrix(seq(0, 1, length=1001), ncol=p)

# create sequential
kriging <- Ordinary.Kriging(X, y, interpolation=TRUE, fit=TRUE,
                            kernel.parameters=list(type="Gaussian"))
design.parameters <- list(design.lower.bound=c(0), design.upper.bound=c(1))

# define custom criterion function (e.g., ALM-like using predictive SD)
my.criterion <- function(mean, sd, X, y, clist) { return(sd) }
sequential <- UDF.Sequential(design.parameters, kriging,
                             criterion.function=my.criterion,
                             criterion.parameters=list())

# sequential design loop with visualization
for (i in 1:5) {
  # visualize current kriging
  pred <- Predict.Sequential(sequential, newX)
  plot(newX, f(newX), "l", main=sprintf("UDF Iteration %d", i))
  points(Get.Data(sequential)$X, Get.Data(sequential)$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 next design and evaluate
  Xn <- Get.Next.Design(sequential)
  new.design <- list(X=Xn, y=f(Xn))
  points(new.design$X, new.design$y, pch=16, col="green")
  # add to kriging and update
  sequential <- Add.Next.Design(sequential, new.design)
}

# final visualization
pred <- Predict.Sequential(sequential, newX)
plot(newX, f(newX), "l", main="UDF Final")
points(Get.Data(sequential)$X, Get.Data(sequential)$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)

Universal Kriging

Description

This functions fits the universal kriging model to the data.

Usage

Universal.Kriging(
  X,
  y,
  basis.function,
  interpolation = TRUE,
  fit = TRUE,
  nugget = NULL,
  kernel = NULL,
  kernel.parameters = list(),
  nlopt.parameters = list()
)

Arguments

Details

Universal kriging permits a more general function of mean, which can be specified using basis.function. Please see Joseph (2006) 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.

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 of Fit.Kriging for detail usages.

When the lengthscale / correlation parameters are unknown, all parameters including the constant mean 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 of Fit.Kriging for detail usages.

Value

A Universal Kriging Class Object.

Author(s)

Chaofan Huang, Difan Song, and V. Roshan Joseph

References

Joseph, V. R. (2026). Experimental Design for Data Science and Engineering. Chapman & Hall/CRC, doi:10.1201/9781003661535.

See Also

Fit.Kriging, Predict.Kriging, Get.Kriging.Parameters.

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)

basis.function <- function(x) {c(1,x[1],x[1]^2)}

# approach 1
kriging <- Universal.Kriging(X, y, basis.function=basis.function,
                             interpolation=TRUE, fit=TRUE, 
                             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)

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