Beta Kernel Process Modeling

Implements the Beta Kernel Process (BKP) for nonparametric modeling of covariate-dependent binomial probabilities, and the Dirichlet Kernel Process (DKP) for categorical or multinomial response data. Scalable global-local approximations are provided through TwinBKP and TwinDKP, using twinning-selected global subsets and local nearest-neighbour updates. Functions are included for model fitting, predictive inference with uncertainty quantification, posterior simulation, and visualization in one- and two-dimensional input spaces. Gaussian, Matern 5/2, Matern 3/2, and Wendland kernels are supported, with hyperparameters selected by multi-start derivative-free optimization. For more details, see Zhao, Qing, and Xu (2025) .

BKP: Beta Kernel Process Modeling

Total downloads

Model covariate-dependent binomial and multinomial probabilities directly using flexible kernels and closed-form Beta or Dirichlet updates—without introducing latent Gaussian variables or relying on MCMC.

BKP provides probability-scale kernel models for binary and aggregated binomial responses. Kernel-weighted conjugate updating yields closed-form, pointwise Beta posterior summaries, including posterior means, variances, credible intervals, quantiles, and simulations.

The package also implements the Dirichlet Kernel Process (DKP) for categorical and multinomial responses. For larger datasets, TwinBKP and TwinDKP combine twinning-selected global subsets with location-specific nearest-neighbour updates.

Software paper · Reproducibility materials · CRAN · Issue tracker

BKP posterior summaries for a two-dimensional binomial probability surface.

Why BKP?

Closed-form inference: obtain pointwise Beta or Dirichlet posterior summaries without MCMC or numerical posterior approximation.
Direct probability-scale modeling: model covariate-dependent probabilities without introducing latent Gaussian variables or a link function.
Multiple response types: handle binary, binomial, categorical, and multinomial observations within a common framework.
Flexible kernels: choose Gaussian, Matérn 5/2, Matérn 3/2, or compactly supported Wendland kernels with isotropic or anisotropic length scales.
Automatic hyperparameter tuning: select kernel length scales by leave-one-out cross-validation using the Brier score or log-loss.
Scalable approximations: use TwinBKP or TwinDKP for twinning-based global-local modeling of larger datasets.

Standard S3 interfaces are provided for prediction, simulation, visualization, model summaries, fitted values, posterior quantiles, and parameter extraction.

Which model should I use?

Response type	Full model	Scalable approximation	Shepard ESS calibration
Binary or binomial	`fit_BKP()`	`fit_TwinBKP()`	`fit_BKP()` only
Categorical or multinomial	`fit_DKP()`	`fit_TwinDKP()`	`fit_DKP()` only

Installation

Install the stable release from CRAN:

install.packages("BKP")

Install the development version from GitHub:

# install.packages("pak")
pak::pak("Jiangyan-Zhao/BKP")

The GitHub development version may contain changes that have not yet been released on CRAN.

Quick start

For binomial data:

X is an $n \times d$ matrix of covariates;
y contains the observed success counts;
m contains the corresponding numbers of trials;
Xbounds gives the lower and upper bounds of each covariate.

library(BKP)

fit <- fit_BKP(
  X = X,
  y = y,
  m = m,
  Xbounds = Xbounds
)

summary(fit)
plot(fit, engine = "ggplot")

# Xnew must have the same number of columns as X.
pred <- predict(
  fit,
  Xnew = Xnew
)

pred

When theta is omitted, fit_BKP() selects the kernel length scale by LOOCV. Supply a positive theta to skip optimization and fit the model using a fixed length scale.

Complete reproducible one-dimensional example

library(BKP)

set.seed(123)

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 <- matrix(sort(runif(n, -2, 2)), ncol = 1)

true_pi <- true_pi_fun(X)
m <- sample(80:120, n, replace = TRUE)

y <- rbinom(
  n = n,
  size = m,
  prob = true_pi
)

fit <- fit_BKP(
  X = X,
  y = y,
  m = m,
  Xbounds = Xbounds
)

summary(fit)
plot(fit, engine = "ggplot")

Xnew <- matrix(
  seq(-2, 2, length.out = 100),
  ncol = 1
)

pred <- predict(
  fit,
  Xnew = Xnew
)

head(pred)

Categorical and multinomial responses

Use fit_DKP() when the response is stored as an $n \times q$ matrix Y of nonnegative class counts or frequencies:

dkp_fit <- fit_DKP(
  X = X,
  Y = Y,
  Xbounds = Xbounds
)

summary(dkp_fit)
plot(dkp_fit, engine = "ggplot")

dkp_pred <- predict(
  dkp_fit,
  Xnew = Xnew
)

DKP replaces the pointwise Beta posterior with a pointwise Dirichlet posterior and provides class-specific posterior summaries and classifications.

Scaling to larger datasets

For larger binomial datasets, replace fit_BKP() with fit_TwinBKP():

twin_fit <- fit_TwinBKP(
  X = X,
  y = y,
  m = m,
  Xbounds = Xbounds
)

summary(twin_fit)
plot(twin_fit, engine = "ggplot")

twin_pred <- predict(
  twin_fit,
  Xnew = Xnew
)

TwinBKP combines:

a twinning-selected global subset for broad distributional coverage; and
location-specific nearest neighbours for local refinement.

For categorical or multinomial data, use the analogous fit_TwinDKP() interface:

twindkp_fit <- fit_TwinDKP(
  X = X,
  Y = Y,
  Xbounds = Xbounds
)

Documentation and reproducibility

The statistical foundations, implementation details, and worked examples are available in:

The reproducibility repository contains the manuscript source files, analysis scripts, data-processing code, and materials used to generate the examples and figures in the software paper.