simulate: Simulate from a Fitted BKP or DKP Model

Description

Generates random draws from the posterior predictive distribution of a fitted BKP or DKP model at specified input locations.

For BKP models, posterior samples are generated from Beta distributions characterizing success probabilities. Optionally, binary class labels can be derived by applying a user-specified classification threshold.

For DKP models, posterior samples are generated from Dirichlet distributions characterizing class probabilities. If training responses are single-label (i.e., one-hot encoded), class labels may additionally be assigned using the maximum a posteriori (MAP) rule.

Usage

# S3 method for BKP
simulate(object, nsim = 1, seed = NULL, Xnew = NULL, threshold = NULL, ...)
# S3 method for DKP
simulate(object, nsim = 1, seed = NULL, Xnew = NULL, ...)

Value

A list with the following components:

samples: For BKP: A numeric matrix of size nrow(Xnew) × nsim, where each column corresponds to one posterior draw of success probabilities.
For DKP: A numeric array of dimension nsim × q × nrow(Xnew), containing simulated class probabilities from Dirichlet posteriors, where q is the number of classes.
mean: For BKP: A numeric vector of posterior mean success probabilities at each Xnew.
For DKP: A numeric matrix of dimension nrow(Xnew) × q, containing posterior mean class probabilities.
class: For BKP: An integer matrix of dimension nrow(Xnew) × nsim, indicating simulated binary class labels (0/1), returned when threshold is specified.
For DKP: An integer matrix of dimension nrow(Xnew) × nsim, where each entry corresponds to a MAP-predicted class label, returned only when training data is single-label.
X: The training input matrix used to fit the BKP/DKP model.
Xnew: The new input locations at which simulations are generated.
threshold: The classification threshold used for generating binary class labels (if provided).

Arguments

object: An object of class "BKP" or "DKP", typically returned by fit_BKP or fit_DKP.
nsim: Number of posterior samples to generate (default = 1).
seed: Optional integer seed for reproducibility.
Xnew: A numeric matrix or vector of new input locations at which simulations are generated.
threshold: Classification threshold for binary decisions (BKP only). When specified, posterior draws exceeding the threshold are classified as 1, and those below as 0. The default is NULL.
...: Additional arguments (currently unused).

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

## -------------------- BKP --------------------
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
model <- fit_BKP(X, y, m, Xbounds=Xbounds)

# Simulate 5 posterior draws of success probabilities
Xnew <- matrix(seq(-2, 2, length.out = 5), ncol = 1)
simulate(model, Xnew = Xnew, nsim = 5)

# Simulate binary classifications (threshold = 0.5)
simulate(model, Xnew = Xnew, nsim = 5, threshold = 0.5)

## -------------------- DKP --------------------
set.seed(123)

# Define true class probability function (3-class)
true_pi_fun <- function(X) {
  p1 <- 1/(1+exp(-3*X))
  p2 <- (1 + exp(-X^2) * cos(10 * (1 - exp(-X)) / (1 + exp(-X)))) / 2
  return(matrix(c(p1/2, p2/2, 1 - (p1+p2)/2), nrow = length(p1)))
}

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

# Generate multinomial responses
Y <- t(sapply(1:n, function(i) rmultinom(1, size = m[i], prob = true_pi[i, ])))

# Fit DKP model
model <- fit_DKP(X, Y, Xbounds = Xbounds)

# Simulate 5 draws from posterior Dirichlet distributions at new point
Xnew <- matrix(seq(-2, 2, length.out = 5), ncol = 1)
simulate(model, Xnew = Xnew, nsim = 5)