predict: Posterior Prediction for BKP or DKP Models

Description

Generates posterior predictive summaries from a fitted Beta Kernel Process (BKP) or Dirichlet Kernel Process (DKP) model at new input locations. Supports prediction of posterior mean, variance, credible intervals, and classification labels (where applicable).

Usage

# S3 method for BKP
predict(object, Xnew = NULL, CI_level = 0.95, threshold = 0.5, ...)
# S3 method for DKP
predict(object, Xnew = NULL, CI_level = 0.95, ...)

Value

A list containing posterior predictive summaries:

X

The original training input locations.

Xnew

The new input locations for prediction (same as Xnew if provided).

alpha_n, beta_n

Posterior shape parameters for each location:

BKP: Vectors of posterior shape parameters (alpha_n, beta_n) for each input location.
DKP: Posterior shape parameter matrix alpha_n (rows = input locations, columns = classes).

mean

Posterior mean prediction:

BKP: Posterior mean success probability at each location.
DKP: Matrix of posterior mean class probabilities (rows = inputs, columns = classes).

variance

Posterior predictive variance:

BKP: Variance of success probability.
DKP: Matrix of predictive variances for each class.

lower

Lower bound of the posterior credible interval:

BKP: Lower bound (e.g., 2.5th percentile for 95% CI).
DKP: Matrix of lower bounds for each class.

upper

Upper bound of the posterior credible interval:

BKP: Upper bound (e.g., 97.5th percentile for 95% CI).
DKP: Matrix of upper bounds for each class.

class

Predicted label:

BKP: Binary class (0 or 1) based on posterior mean and threshold, only if m = 1.
DKP: Predicted class label with highest posterior mean probability.

CI_level

The specified credible interval level.

Arguments

object: An object of class "BKP" or "DKP", typically returned by fit_BKP or fit_DKP.
Xnew: A numeric vector or matrix of new input locations at which to generate predictions. If NULL, predictions are returned for the training data.
CI_level: Numeric between 0 and 1 specifying the credible level for posterior intervals (default 0.95 for 95% credible interval).
threshold: Numeric between 0 and 1 specifying the classification threshold for binary predictions based on posterior mean (used only for BKP; default is 0.5).
...: Additional arguments passed to generic predict methods (currently not used; included for S3 method consistency).

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 Examples ======================= #
# ============================================================== #

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

# Prediction on training data
predict(model1)

# Prediction on new data
Xnew = matrix(seq(-2, 2, length = 10), ncol=1) #new data points
predict(model1, Xnew = Xnew)


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

# Prediction on training data
predict(model2)

# Prediction on new data
x1 <- seq(Xbounds[1,1], Xbounds[1,2], length.out = 10)
x2 <- seq(Xbounds[2,1], Xbounds[2,2], length.out = 10)
Xnew <- expand.grid(x1 = x1, x2 = x2)
predict(model2, Xnew = Xnew)

# ============================================================== #
# ========================= DKP Examples ======================= #
# ============================================================== #

#-------------------------- 1D Example ---------------------------
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
model1 <- fit_DKP(X, Y, Xbounds = Xbounds)

# Prediction on training data
predict(model1)

# Prediction on new data
Xnew = matrix(seq(-2, 2, length = 10), ncol=1) #new data points
predict(model1, Xnew)


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

# Define latent function and transform to 3-class probabilities
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)- m)/s
  p1 <- pnorm(f) # Transform to probability
  p2 <- sin(pi * X[,1]) * sin(pi * X[,2])
  return(matrix(c(p1/2, p2/2, 1 - (p1+p2)/2), nrow = length(p1)))
}

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(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
model2 <- fit_DKP(X, Y, Xbounds = Xbounds)

# Prediction on training data
predict(model2)

# Prediction on new data
x1 <- seq(Xbounds[1,1], Xbounds[1,2], length.out = 10)
x2 <- seq(Xbounds[2,1], Xbounds[2,2], length.out = 10)
Xnew <- expand.grid(x1 = x1, x2 = x2)
predict(model2, Xnew)