Bayesian Simultaneous Credible Bands for Polynomial Regression

Provides functions to construct two-sided Bayesian simultaneous credible bands (BSCBs) for the regression curve in univariate polynomial regression over a finite covariate interval. Six methods are implemented, including Normal-Gamma conjugate priors (with empirical Bayes, unit-information, and g-prior hyperparameter specifications), non-conjugate priors fitted via Hamiltonian Monte Carlo (HMC) using 'cmdstanr', and a non-informative independent Jeffreys prior approach. Also includes functions for computing the empirical simultaneous coverage rate (ESCR) and posterior simultaneous coverage probability (PSCP), enabling performance comparison across methods. The methodology is described in: Yang et al. (2026). "Bayesian simultaneous credible bands for polynomial regression" .

BSCB

Overview

BSCB provides methods for constructing Bayesian simultaneous credible bands (BSCB) and Bayesian pointwise credible bands (BPCB) for polynomial regression models. The package implements the following approaches based on different prior specifications:

• BSCB-C: Normal-Gamma conjugate prior (empirical Bayes, unit-information, or Zellner’s g-prior)
• BSCB-H: Non-conjugate prior implemented via HMC
• BSCB-I-J: Independent Jeffreys prior (objective Bayesian inference)
• BPCB-I-J: Bayesian pointwise credible band under the independent Jeffreys prior

The methodology is based on the following paper:

• Yang, F., Han, Y., Liu, W., & Hall, I. (2026). Bayesian simultaneous credible bands for polynomial regression. arXiv preprint arXiv:2606.28015.

A full demo is available here.

Installation

You can install BSCB on CRAN using:

install.packages("BSCB")

library(BSCB)

Alternatively, you can install the package through github:

# install.packages("devtools")
devtools::install_github("fannyyang73/BSCB")

Requirements

• R (≥ 4.0.0)
• R packages: mvtnorm, MASS, OptimalDesign, instantiate, posterior

Quick Start

library(BSCB)

# Simulate data from a quadratic model
set.seed(123)
n <- 50
x <- seq(-5, 5, length.out = n)
X <- cbind(1, x, x^2)
theta_true <- c(-6, -3, 0.25)
Y <- X %*% theta_true + rnorm(n, sd = 0.2)

# --- BSCB-C: Bayesian simultaneous credible bands under the Normal-Gamma conjugate prior ---
fit_c <- compute_bscb_conjugate(
  X              = X,
  Y              = Y,
  alpha          = 0.05,
  a              = -5,
  b              =  5,
  L              = 500000,
  theta_true     = theta_true,
  hyperparameter = "g_prior",   # "empirical", "unit_info", or "g_prior"
  optimize_type  = "P"          # "P" = polyroot (recommended)
)

# --- BSCB-H: Bayesian simultaneous credible bands under a non-conjugate prior implemented via HMC

mod <- instantiate::stan_package_model(
  name    = "HMC_model",
  package = "BSCB",
  compile = TRUE
)

fit_h <- compute_bscb_hmc(
  X     = X,
  Y     = Y,
  V     = diag(n),
  alpha = alpha,
  a     = a,
  b     = b,
  theta_true = theta_true,
  prior_type = "normal_half_cauchy",
  L     = L,
  draw_num = 10000
)


# --- BSCB-I-J: Bayesian simultaneous credible bands under the Independent Jeffreys prior ---
fit_j <- compute_bscb_ind_jeffreys(
  X     = X,
  Y     = Y,
  alpha = 0.05,
  a     = -5,
  b     =  5,
  L     = 500000,
  theta_true     = theta_true
)

# --- BPCB-I-J: Bayesian pointwise credible bands under the Independent Jeffreys prior ---
fit_p <- compute_bpcb_ind_jeffreys(
  X     = X,
  Y     = Y,
  alpha = 0.05,
  a     = -5,
  b     =  5,
  theta_true     = theta_true
)

# Evaluate bands over a grid and plot
x_seq <- seq(-5, 5, length.out = 500)

plot(x_seq, fit_c$lower_bound(x_seq), type = "l",
     col = "red", lty = 2, lwd = 2,
     ylim = range(c(fit_c$lower_bound(x_seq),
                    fit_c$upper_bound(x_seq), Y)),
     xlab = "x", ylab = "y",
     main = "95% Bayesian Simultaneous Credible Band")
lines(x_seq, fit_c$upper_bound(x_seq), col = "red",  lty = 2, lwd = 2)
lines(x_seq, cbind(1, x_seq, x_seq^2) %*% theta_true,
      col = "blue", lwd = 2)
points(x, Y, pch = 16, col = "gray")
legend("topright",
       legend = c("True curve", "Data", "95% BSCB-C"),
       col    = c("blue", "gray", "red"),
       lty    = c(1, NA, 2),
       pch    = c(NA, 16, NA))
       
# Evaluate coverage
coverage_ESCR(fit_c, optimize_type = "P", verbose = TRUE)
coverage_PSCP(fit_c, draw_num = 10000, optimize_type = "P", verbose = TRUE)

Main Functions

Key Arguments