## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  collapse = TRUE, comment = "#>",
  eval = FALSE
)

## ----packages-----------------------------------------------------------------
# library(gdpar)
# # These two are Suggests; gdpar_eb() and gdpar_compare_eb_fb() require
# # them at runtime.
# library(cmdstanr)
# library(posterior)

## ----data---------------------------------------------------------------------
# set.seed(20260526L)
# n  <- 150L
# df <- data.frame(x = stats::rnorm(n))
# df$y_scalar <- 1.0 + 0.4 * df$x + stats::rnorm(n, sd = 0.3)
# # Multivariate (p = 2) outcome for Path A.
# df$y_p2 <- cbind(
#   1.0 + 0.4 * df$x + stats::rnorm(n, sd = 0.3),
#   -0.5 + 0.2 * df$x + stats::rnorm(n, sd = 0.4)
# )
# # Same dataset is fine for Path B (K > 1, p = 1) and Path C (K > 1,
# # p > 1) by reusing y_scalar / y_p2 with a multi-slot family below.

## ----eb-base------------------------------------------------------------------
# fit_eb <- gdpar_eb(
#   formula        = y_scalar ~ x,
#   family         = gdpar_family("gaussian"),
#   amm            = amm_spec(a = ~ x),
#   data           = df,
#   iter_warmup    = 500L,
#   iter_sampling  = 500L,
#   chains         = 2L,
#   refresh        = 0L,
#   seed           = 1L,
#   laplace_control = list(multi_start_M = 5L)
# )
# print(fit_eb)

## ----eb-base-summary----------------------------------------------------------
# summary(fit_eb)

## ----eb-path-A----------------------------------------------------------------
# fit_eb_A <- gdpar_eb(
#   formula        = y_p2 ~ x,
#   family         = gdpar_family_multi("gaussian", p = 2L),
#   amm            = amm_spec(p = 2L, dims = dimwise(a = ~ x)),
#   data           = df,
#   iter_warmup    = 500L,
#   iter_sampling  = 500L,
#   chains         = 2L,
#   refresh        = 0L,
#   seed           = 2L
# )
# print(fit_eb_A)

## ----eb-path-B----------------------------------------------------------------
# fs <- gdpar_bf(y_scalar ~ a(x), sigma ~ a(x))
# fit_eb_B <- gdpar_eb(
#   formula        = fs,
#   family         = gdpar_family("gaussian"),
#   data           = df,
#   iter_warmup    = 500L,
#   iter_sampling  = 500L,
#   chains         = 2L,
#   refresh        = 0L,
#   seed           = 3L,
#   skip_id_check  = TRUE
# )
# print(fit_eb_B)

## ----eb-path-C----------------------------------------------------------------
# y_p2_int <- matrix(
#   rnbinom(n * 2L, size = 5, mu = exp(0.5 + 0.2 * df$x)),
#   n, 2L
# )
# df$y_p2_int <- y_p2_int
# fit_eb_C <- gdpar_eb(
#   formula        = y_p2_int ~ x,
#   family         = gdpar_family("neg_binomial_2"),
#   amm            = list(
#     mu  = amm_spec(p = 2L, dims = dimwise(a = ~ x)),
#     phi = amm_spec(p = 2L, dims = dimwise(a = ~ x))
#   ),
#   data           = df,
#   iter_warmup    = 500L,
#   iter_sampling  = 500L,
#   chains         = 2L,
#   refresh        = 0L,
#   seed           = 4L,
#   skip_id_check  = TRUE
# )
# print(fit_eb_C)

## ----compare-eb-fb------------------------------------------------------------
# fit_fb <- gdpar(
#   formula        = y_scalar ~ x,
#   family         = gdpar_family("gaussian"),
#   amm            = amm_spec(a = ~ x),
#   data           = df,
#   iter_warmup    = 500L,
#   iter_sampling  = 500L,
#   chains         = 2L,
#   refresh        = 0L,
#   seed           = 1L
# )
# cmp <- gdpar_compare_eb_fb(fit_eb, fit_fb, level = 0.95,
#                             tv_bins = 30L)
# print(cmp)
# summary(cmp)