Package {RQK}

Fit a rational quadratic kriging model

Description

Fits the rational quadratic kriging (RQK) model using compiled C++ routines for factor-weight estimation, leave-one-out tuning, and conformal prediction.

The implementation is based on RcppArmadillo. It is recommended to use a high performance BLAS library for optimized speed.

Usage

RQK(
  X,
  y,
  theta = NULL,
  theta_fit = TRUE,
  noise = FALSE,
  nugget = NULL,
  q = NULL,
  theta_bound = NULL,
  scale = FALSE
)

Arguments

Details

The model uses a rational quadratic kernel on a p-dimensional input space:

R(\mathbf{x}_i, \mathbf{x}_j) = \left( 1 + \frac{\|\mathbf{x}_i - \mathbf{x}_j\|_w^2}{\theta^2} \right)^{-q} + \mathrm{nugget}\,\mathbb{I}\{\mathbf{x}_i = \mathbf{x}_j\}

and estimates anisotropic factor weights over the p coordinates. The returned fit stores both user-facing quantities such as theta, q, and w, and cached matrices used by predict.RQK() to evaluate point predictions and conformal prediction intervals. When scale = TRUE, fitting is performed on the scaled inputs and the returned X_train cache is also on the scaled input scale.

Value

An object of class "RQK" with elements including:

• w: estimated nonnegative factor weights summing to one.

• theta: fitted kernel length scale.

• nugget: fixed or estimated nugget term.

• q: kernel shape parameter.

• LOO_pred: leave-one-out point predictions for the training responses.

• scores: training conformity scores used by predict.RQK() to construct interval types 1–4.

• X_train, y_train: training inputs and responses.

• scale: logical flag indicating whether the inputs were scaled.

• scale_min, scale_max: training minima and maxima for each input variable, present when scale = TRUE.

• additional cached matrices used internally by predict.RQK().

References

Song, D., and Joseph, V. R. (2026). Rational Quadratic Kriging. arXiv preprint arXiv:2610.18077. doi:10.48550/arXiv.2610.18077

Examples

borehole <- function(xx) {
  lower <- c(0.05, 100, 63070, 990, 63.1, 700, 1120, 9855)
  upper <- c(0.15, 50000, 115600, 1110, 116, 820, 1680, 12045)
  xx <- xx * (upper - lower) + lower
  rw <- xx[1]
  r <- xx[2]
  Tu <- xx[3]
  Hu <- xx[4]
  Tl <- xx[5]
  Hl <- xx[6]
  L <- xx[7]
  Kw <- xx[8]
  frac1 <- 2 * pi * Tu * (Hu - Hl)
  frac2a <- 2 * L * Tu / (log(r / rw) * rw^2 * Kw)
  frac2b <- Tu / Tl
  frac2 <- log(r / rw) * (1 + frac2a + frac2b)
  frac1 / frac2
}

set.seed(1)
X_train <- matrix(runif(200 * 8), ncol = 8)
y_train <- apply(X_train, 1, borehole)

fit <- RQK(X_train, y_train, theta = 0.2, theta_fit = FALSE)
fit$theta
fit$q
head(fit$w)

X_test <- matrix(runif(1000 * 8), ncol = 8)
pred <- predict(fit, X_test)
head(data.frame(mean = pred$mean, lower = pred$lower, upper = pred$upper))

Obtain an initial length scale and factor weights

Description

Computes a starting point for RQK fitting, including the kernel length scale and the anisotropic factor weights over the p input coordinates.

Usage

initial_lengthscale(
  X,
  y,
  theta = NULL,
  theta_fit = TRUE,
  noise = FALSE,
  q = NULL
)

Arguments

Value

A list with elements:

• theta: initial kernel length scale on the original scale.

• logtheta: the same value on the log scale.

• q: kernel shape parameter used for the initialization.

• u: unconstrained optimization variables for the factor weights.

• w: factor weights on the simplex.

• par: full optimization vector returned by the initializer.

• logtheta_lb, logtheta_ub: log-scale bounds used for theta.

References

Song, D., and Joseph, V. R. (2026). Rational Quadratic Kriging. arXiv preprint arXiv:2610.18077. doi:10.48550/arXiv.2610.18077

Examples

borehole <- function(xx) {
  lower <- c(0.05, 100, 63070, 990, 63.1, 700, 1120, 9855)
  upper <- c(0.15, 50000, 115600, 1110, 116, 820, 1680, 12045)
  xx <- xx * (upper - lower) + lower
  rw <- xx[1]
  r <- xx[2]
  Tu <- xx[3]
  Hu <- xx[4]
  Tl <- xx[5]
  Hl <- xx[6]
  L <- xx[7]
  Kw <- xx[8]
  frac1 <- 2 * pi * Tu * (Hu - Hl)
  frac2a <- 2 * L * Tu / (log(r / rw) * rw^2 * Kw)
  frac2b <- Tu / Tl
  frac2 <- log(r / rw) * (1 + frac2a + frac2b)
  frac1 / frac2
}

set.seed(1)
X_train <- matrix(runif(200 * 8), ncol = 8)
y_train <- apply(X_train, 1, borehole)

init <- initial_lengthscale(X_train, y_train, noise = FALSE)
init$theta
head(init$w)

Predict from a fitted RQK model

Description

Computes point predictions and, unless disabled, conformal prediction intervals from an object returned by RQK().

Usage

## S3 method for class 'RQK'
predict(object, newdata, alpha = 0.05, interval_type = 4, ...)

Arguments

Details

If object was fitted with scale = TRUE, newdata is scaled to ⁠[0, 1]⁠ using the training minima and maxima saved in the fitted object before predictions are computed.

Value

A list with:

• mean: point predictions. If interval_type = "none", this is the only returned component. Otherwise the list also contains:

• lower, upper: numeric prediction bounds for interval types 1, 3, and 4. For type 2 they are NA_real_ because the prediction set may be disjoint.

• intervals: a list with one element per prediction point. Each element is a two-column matrix with columns lower and upper. Type 2 may contain multiple rows, whereas types 1, 3, and 4 contain a single row.

References

Song, D., and Joseph, V. R. (2026). Rational Quadratic Kriging. arXiv preprint arXiv:2610.18077. doi:10.48550/arXiv.2610.18077

Examples

borehole <- function(xx) {
  lower <- c(0.05, 100, 63070, 990, 63.1, 700, 1120, 9855)
  upper <- c(0.15, 50000, 115600, 1110, 116, 820, 1680, 12045)
  xx <- xx * (upper - lower) + lower
  rw <- xx[1]
  r <- xx[2]
  Tu <- xx[3]
  Hu <- xx[4]
  Tl <- xx[5]
  Hl <- xx[6]
  L <- xx[7]
  Kw <- xx[8]
  frac1 <- 2 * pi * Tu * (Hu - Hl)
  frac2a <- 2 * L * Tu / (log(r / rw) * rw^2 * Kw)
  frac2b <- Tu / Tl
  frac2 <- log(r / rw) * (1 + frac2a + frac2b)
  frac1 / frac2
}

set.seed(1)
X_train <- matrix(runif(200 * 8), ncol = 8)
y_train <- apply(X_train, 1, borehole)
fit <- RQK(X_train, y_train, theta = 0.2, theta_fit = FALSE)

X_test <- matrix(runif(1000 * 8), ncol = 8)
pred <- predict(fit, X_test, alpha = 0.1)
head(pred$mean)