R/quantile_NBKP.R

In NBKP: Beta Kernel Process Modeling with Negative Binomial Response

#' @rdname quantile.NBKP
#' @title Posterior Quantiles from a Fitted NBKP Model
#'
#' @description Compute posterior quantiles from a fitted \code{NBKP} model.
#' This returns posterior quantiles of the latent mean count from the Gamma posterior distribution.
#'
#' @param x An object of class \code{NBKP}, typically the result of a call to \code{\link{fit_NBKP}}.
#' @param probs Numeric vector of probabilities specifying which posterior quantiles to return.
#' Defaults to \code{c(0.025, 0.5, 0.975)}.
#' @param ... Additional arguments (currently unused).
#'
#' @return A numeric vector (if \code{length(probs) = 1}) or numeric matrix (if \code{length(probs) > 1})
#' of posterior quantiles. Rows correspond to observations, and columns correspond to the requested probabilities.
#'
#' @details For a \code{NBKP} model, posterior quantiles are computed from the Gamma Kernel Process posterior distribution.
#'
#' @seealso \code{\link{fit_NBKP}} for model fitting.
#'
#' @keywords NBKP
#'
#' @references Zhao J, Qing K, Xu J (2025). \emph{BKP: An R Package for Beta
#' Kernel Process Modeling}. arXiv. https://doi.org/10.48550/arxiv.2508.10447.
#'
#' @examples
#' \donttest{
#' set.seed(123)
#'
#' # Define true mean function
#' true_mu_fun <- function(x) {
#'   exp(sin(x) + 0.5)
#' }
#'
#' n <- 30
#' Xbounds <- matrix(c(-2, 2), nrow = 1)
#' X <- tgp::lhs(n = n, rect = Xbounds)
#' true_mu <- true_mu_fun(X)
#' y <- rnbinom(n, size = 1, mu = true_mu)
#'
#' # Fit NBKP model
#' model <- fit_NBKP(X, y, Xbounds = Xbounds)
#'
#' # Extract posterior quantiles
#' quantile(model)
#' }
#'
#' @export
#' @method quantile NBKP
quantile.NBKP <- function(x, probs = c(0.025, 0.5, 0.975), ...) {
  # arguments checking
  if (!is.numeric(probs) || any(probs < 0 | probs > 1)) {
    stop("'probs' must be a numeric vector with all values in [0,1].")
  }

  # Extract posterior gamma parameters
  alpha_n <- x$alpha_n
  beta_n  <- x$beta_n

  if (length(probs) > 1) {
    # Posterior quantiles matrix: rows = observations, cols = probs
    post_q <- t(mapply(function(a, b) stats::qgamma(probs, a, b), alpha_n, beta_n))
    colnames(post_q) <- paste0(probs * 100, "%")
  } else {
    # Single probability: return a vector
    post_q <- mapply(function(a, b) stats::qgamma(probs, a, b), alpha_n, beta_n)
  }
  return(post_q)
}