Bayesian simultaneous credible bands for the polynomial model
Overview
This demo illustrates how to use the BSCB package to construct and evaluate Bayesian simultaneous credible bands (BSCB) and Bayesian pointwise credible bands (BPCB) for polynomial regression. These methods are demonstrated:
- BSCB-C: BSCB under the Normal-Gamma conjugate prior
- BSCB-I-J: BSCB under the independent Jeffreys prior
- BPCB-I-J: BPCB under the independent Jeffreys prior
1. Simulate Data
We simulate data from a quadratic regression model \(Y_i = \theta_0 + \theta_1 x_i + \theta_2 x_i^2 + \varepsilon_i\), where \(\varepsilon_i \sim N(0, \sigma^2)\).
Notably, \(L\) should set to \(L=500,000\), whereas here \(L=5000\) just to show a quick example.
# Simulate data from a quadratic model using a D-optimal covariate design
theta_true <- c(-6, -3, 0.25)
alpha <- 0.05
a <- -0.5
b <- 0.5
L <- 5000
p <- 2
n <- 20
e_sd <- 0.2
sim_data <- generate_simulation_data(
p = p,
n = n,
e_sd = e_sd,
theta_true = theta_true,
a = a,
b = b,
replication = 1,
design_index = 2,
center_index = 1
)
#> [1] Call of the function:
#> OptimalDesign::od_KL(Fx = design_matrix, N = n, crit = "D", K = 7,
#> L = 19, t.max = 5)
#> [1] Call of the function:
#> od_REX(Fx = Fx, crit = "D", track = FALSE)
#> [1] Running od_D_KL for cca 5 seconds starting at 2026-07-02 17:51:33.59641.
#> [1] The problem size is n=300000, m=3, N=20
#> [1] Setting K=7, L=19
#> [1] od_D_KL Time: 0.2 Value: 1.562209 Efficiency: 0.590478
#> [1] od_D_KL Time: 1 Value: 2.638816 Efficiency: 0.99741
#> [1] od_D_KL Time: 2 Value: 2.638816 Efficiency: 0.99741
#> [1] od_D_KL Time: 3 Value: 2.638816 Efficiency: 0.99741
#> [1] od_D_KL Time: 4 Value: 2.638816 Efficiency: 0.99741
#> [1] od_D_KL finished after 5.01 seconds at 2026-07-02 17:51:38.60922
#> [1] with 8 restarts and 357 exchanges.
X <- sim_data$X
x <- sim_data$X[,2]
Y <- as.numeric(sim_data$Y.list[[1]])2. Construct the Bands
BSCB-C: BSCB under the Normal-Gamma conjugate prior
The Normal-Gamma conjugate prior supports three hyperparameter
specifications: "empirical" (empirical Bayes),
"unit_info" (unit-information prior), and
"g_prior" (Zellner’s g-prior). The critical constant \(\lambda\) is estimated via \(L\) Monte Carlo draws.
# --- BSCB-C: Bayesian simultaneous credible bands under the Normal-Gamma conjugate prior ---
fit_c <- compute_bscb_conjugate(
X = X,
Y = Y,
alpha = alpha,
a = a,
b = b,
L = L,
theta_true = theta_true,
hyperparameter = "g_prior", # "empirical", "unit_info", or "g_prior"
optimize_type = "P" # "P" = polyroot (recommended)
)
#> Computing lambda via Monte Carlo sampling...
#> The critical constant lambda = 2.647649
cat("Critical constant (BSCB-C):", fit_c$lambda, "\n")
#> Critical constant (BSCB-C): 2.647649
cat("Posterior mean of theta:\n")
#> Posterior mean of theta:
print(round(fit_c$mu_star, 4))
#> [1] -6.0383 -3.0610 0.4922BSCB-I-J: BSCB under the independent Jeffreys prior
BSCB-I-J is equivalent to the FSCB in Liu et al.(2008) for finite sample sizes.
# --- BSCB-I-J: Bayesian simultaneous credible bands under the Independent Jeffreys prior ---
fit_j <- compute_bscb_ind_jeffreys(
X = X,
Y = Y,
alpha = alpha,
a = a,
b = b,
theta_true = theta_true,
L = L
)
#> Computing lambda via Monte Carlo sampling...
#> The critical constant lambda = 2.654988
cat("Critical constant (BSCB-J):", fit_j$lambda, "\n")
#> Critical constant (BSCB-J): 2.654988BSCB-H-C: BSCB under the normal-half-Cauchy prior(0,2) implemented via HMC
If you would like to produce BSCB-H-C, you can use the following codes.
# mod <- instantiate::stan_package_model(
# name = "HMC_model",
# package = "BSCB",
# compile = TRUE
# )# --- BSCB-H-C: BSCB under the normal-half-Cauchy prior(0,2) implemented via HMC ---
# 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
# )
#
# cat("Critical constant (BSCB-H-C):", fit_h$lambda, "\n")BPCB-I-J: BPCB under the independent Jeffreys prior
BPCB-I-J is constructed by connecting confidence intervals at each individual points in the covariate domain. It’s also equivalent to the FPCB.
3. Evaluate Bands over a Grid and Plot
library(ggplot2)
x_seq <- seq(-0.5, 0.5, length.out = 500)
y_true <- as.numeric(cbind(1, x_seq, x_seq^2) %*% theta_true)
df_obs <- data.frame(x = x, Y = as.numeric(Y))
# Collect all band boundaries into a single data frame
df_bands <- data.frame(
x = rep(x_seq, 3),
lower = c(as.numeric(fit_c$lower_bound(x_seq)),
as.numeric(fit_j$lower_bound(x_seq)),
as.numeric(fit_p$lower_bound(x_seq))),
upper = c(as.numeric(fit_c$upper_bound(x_seq)),
as.numeric(fit_j$upper_bound(x_seq)),
as.numeric(fit_p$upper_bound(x_seq))),
method = rep(c("BSCB-C-G",
"BSCB-I-J",
"BPCB-I-J"),
each = length(x_seq))
)
df_true <- data.frame(x = x_seq, y = y_true)
band_colours <- c(
"BSCB-C-G" = "#4DAF4A",
"BSCB-I-J" = "#E41A1C",
"BPCB-I-J" = "#377EB8"
)
band_linetypes <- c(
"BSCB-C-G" = "F1",
"BSCB-I-J" = "dotdash",
"BPCB-I-J" = "solid"
)
ggplot() +
# Shaded credible regions
geom_ribbon(
data = df_bands,
mapping = aes(x = x, ymin = lower, ymax = upper,
fill = method),
alpha = 0.10
) +
# Band boundaries
geom_line(
data = df_bands,
mapping = aes(x = x, y = lower,
colour = method,
linetype = method),
linewidth = 0.8
) +
geom_line(
data = df_bands,
mapping = aes(x = x, y = upper,
colour = method,
linetype = method),
linewidth = 0.8
) +
# True regression curve
geom_line(
data = df_true,
mapping = aes(x = x, y = y),
colour = "navyblue",
linewidth = 0.7,
linetype = "solid"
) +
# Observed data
geom_point(
data = df_obs,
mapping = aes(x = x, y = Y),
colour = "gray50",
size = 1.5
) +
scale_colour_manual(
values = band_colours,
name = "Method",
breaks = c("BSCB-C-G", "BSCB-I-J", "BPCB-I-J")) +
scale_fill_manual(
values = band_colours,
name = "Method",
breaks = c("BSCB-C-G", "BSCB-I-J", "BPCB-I-J")) +
scale_linetype_manual(
values = band_linetypes,
name = "Method",
breaks = c("BSCB-C-G", "BSCB-I-J", "BPCB-I-J")) +
labs(
title = "95% BSCB-C-G, BSCB-I-J and BPCB-I-J for Quadratic Regression",
x = "x",
y = "y"
) +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5))
4. Evaluate Coverage
Empirical Simultaneous Coverage Rate (ESCR)
coverage_ESCR() returns 1 if the true regression
function lies within the band at all points in \([a, b]\), and 0 otherwise.
escr_j <- coverage_ESCR(fit_j, optimize_type = "P", verbose = TRUE)
escr_c <- coverage_ESCR(fit_c, optimize_type = "P", verbose = TRUE)
escr_p <- coverage_ESCR(fit_p, optimize_type = "P", verbose = TRUE)
cat("ESCR (BSCB-I-J):", escr_j, "\n")
#> ESCR (BSCB-I-J): 1
cat("ESCR (BSCB-C-G):", escr_c, "\n")
#> ESCR (BSCB-C-G): 1
cat("ESCR (BPCB-I-J):", escr_p, "\n")
#> ESCR (BPCB-I-J): 1Posterior Simultaneous Coverage Probability (PSCP)
coverage_PSCP() estimates the proportion of posterior
draws for which \(\sup_{x \in [a,b]} T(x) \leq
\lambda\).
pscp_j <- coverage_PSCP(fit_j, draw_num = 10000,
optimize_type = "P", verbose = TRUE)
#> PSCP = 0.9499
pscp_c <- coverage_PSCP(fit_c, draw_num = 10000,
optimize_type = "P", verbose = TRUE)
#> PSCP = 0.954
pscp_p <- coverage_PSCP(fit_p, draw_num = 10000,
optimize_type = "P", verbose = TRUE)
#> PSCP = 0.8606
cat("PSCP (BSCB-I-J):", round(pscp_j, 4), "\n")
#> PSCP (BSCB-I-J): 0.9499
cat("PSCP (BSCB-C-G):", round(pscp_c, 4), "\n")
#> PSCP (BSCB-C-G): 0.954
cat("PSCP (BPCB-I-J):", round(pscp_p, 4), "\n")
#> PSCP (BPCB-I-J): 0.8606Summary Table
summary_tab <- data.frame(
Method = c("BSCB-C-G", "BSCB-I-J"),
Lambda = round(c(fit_c$lambda, fit_j$lambda), 4),
ESCR = c(escr_c, escr_j),
PSCP = round(c(pscp_c, pscp_j), 4)
)
knitr::kable(summary_tab, caption = "Coverage summary for one simulated dataset")| Method | Lambda | ESCR | PSCP |
|---|---|---|---|
| BSCB-C-G | 2.6476 | 1 | 0.9540 |
| BSCB-I-J | 2.6550 | 1 | 0.9499 |