fit_BKP: Fit a Beta Kernel Process (BKP) Model

Description

Fits a Beta Kernel Process (BKP) model to binary or binomial response data using local kernel smoothing. The method constructs a flexible latent probability surface by updating Beta priors with kernel-weighted observations.

Usage

fit_BKP(
  X,
  y,
  m,
  Xbounds = NULL,
  prior = c("noninformative", "fixed", "adaptive"),
  r0 = 2,
  p0 = mean(y/m),
  kernel = c("gaussian", "matern52", "matern32"),
  loss = c("brier", "log_loss"),
  n_multi_start = NULL,
  theta = NULL
)

Value

A list of class "BKP" containing the fitted BKP model, including:

theta_opt: Optimized kernel hyperparameters (lengthscales).
kernel: Kernel function used, as a string.
loss: Loss function used for hyperparameter tuning.
loss_min: Minimum loss achieved during optimization, or NA if theta was user-specified.
X: Original input matrix (\(n \times d\)).
Xnorm: Normalized input matrix scaled to \([0,1]^d\).
Xbounds: Normalization bounds for each input dimension (\(d \times 2\)).
y: Observed success counts.
m: Observed binomial trial counts.
prior: Type of prior used.
r0: Prior precision parameter.
p0: Prior mean (for fixed priors).
alpha0: Prior Beta shape parameter \(\alpha_0(\mathbf{x})\).
beta0: Prior Beta shape parameter \(\beta_0(\mathbf{x})\).
alpha_n: Posterior shape parameter \(\alpha_n(\mathbf{x})\).
beta_n: Posterior shape parameter \(\beta_n(\mathbf{x})\).

Arguments

X: A numeric input matrix of size \(n \times d\), where each row corresponds to a covariate vector.
y: A numeric vector of observed successes (length n).
m: A numeric vector of total binomial trials (length n), corresponding to each y.
Xbounds: Optional \(d \times 2\) matrix specifying the lower and upper bounds of each input dimension. Used to normalize inputs to \([0,1]^d\). If NULL, inputs are assumed to be pre-normalized, and default bounds \([0,1]^d\) are applied.
prior: Type of prior: "noninformative" (default), "fixed", or "adaptive".
r0: Global prior precision (used when prior = "fixed" or "adaptive").
p0: Global prior mean (used when prior = "fixed"). Default is mean(y/m).
kernel: Kernel function for local weighting: "gaussian" (default), "matern52", or "matern32".
loss: Loss function for kernel hyperparameter tuning: "brier" (default) or "log_loss".
n_multi_start: Number of random initializations for multi-start optimization. Default is \(10 \times d\).
theta: Optional. A positive scalar or numeric vector of length d specifying kernel lengthscale parameters directly. If NULL (default), lengthscales are optimized using multi-start L-BFGS-B to minimize the specified loss.

References

Zhao J, Qing K, Xu J (2025). BKP: An R Package for Beta Kernel Process Modeling. arXiv. https://doi.org/10.48550/arXiv.2508.10447.

Examples

Run this code

#-------------------------- 1D Example ---------------------------
set.seed(123)

# Define true success probability function
true_pi_fun <- function(x) {
  (1 + exp(-x^2) * cos(10 * (1 - exp(-x)) / (1 + exp(-x)))) / 2
}

n <- 30
Xbounds <- matrix(c(-2,2), nrow=1)
X <- tgp::lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- sample(100, n, replace = TRUE)
y <- rbinom(n, size = m, prob = true_pi)

# Fit BKP model
model1 <- fit_BKP(X, y, m, Xbounds=Xbounds)
print(model1)


#-------------------------- 2D Example ---------------------------
set.seed(123)

# Define 2D latent function and probability transformation
true_pi_fun <- function(X) {
  if(is.null(nrow(X))) X <- matrix(X, nrow=1)
  m <- 8.6928
  s <- 2.4269
  x1 <- 4*X[,1]- 2
  x2 <- 4*X[,2]- 2
  a <- 1 + (x1 + x2 + 1)^2 *
    (19- 14*x1 + 3*x1^2- 14*x2 + 6*x1*x2 + 3*x2^2)
  b <- 30 + (2*x1- 3*x2)^2 *
    (18- 32*x1 + 12*x1^2 + 48*x2- 36*x1*x2 + 27*x2^2)
  f <- log(a*b)
  f <- (f- m)/s
  return(pnorm(f))  # Transform to probability
}

n <- 100
Xbounds <- matrix(c(0, 0, 1, 1), nrow = 2)
X <- tgp::lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- sample(100, n, replace = TRUE)
y <- rbinom(n, size = m, prob = true_pi)

# Fit BKP model
model2 <- fit_BKP(X, y, m, Xbounds=Xbounds)
print(model2)