BSCB

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

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)