## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  collapse = TRUE, comment = "#>"
)
have_cmdstan <- requireNamespace("cmdstanr", quietly = TRUE) &&
  isTRUE(try(cmdstanr::cmdstan_version(), silent = TRUE) != "")

## ----bspline-api, eval=FALSE--------------------------------------------------
# library(gdpar)
# 
# # Polynomial baseline (default)
# W_poly <- W_basis(type = "polynomial", dim = 2L)
# 
# # B-spline with internal knots given explicitly, cubic by default
# W_bs <- W_basis(
#   type           = "bspline",
#   knots          = c(-0.5, 0.0, 0.5),
#   degree         = 3L,
#   boundary_knots = c(-1.5, 1.5)
# )

## ----bspline-example, eval=FALSE----------------------------------------------
# set.seed(2026L)
# n <- 80L
# theta_ref_true <- 0.5
# # A sigmoidal modulating effect that the polynomial baseline cannot
# # represent without raising dim_W substantially.
# x_var <- runif(n, -2, 2)
# W_true <- function(z) 1.2 / (1 + exp(-3 * z)) - 0.6
# y <- theta_ref_true + 0.4 * (x_var - mean(x_var)) +
#   W_true(theta_ref_true) * (x_var - mean(x_var)) * 0.7 +
#   rnorm(n, sd = 0.2)
# d <- data.frame(y = y, x = x_var)
# 
# library(gdpar)
# fit_bs <- gdpar(
#   y ~ a(x) + W(),
#   data   = d,
#   family = gdpar_family("gaussian"),
#   W      = W_basis(type = "bspline",
#                    knots = c(-1, 0, 1),
#                    degree = 3L,
#                    boundary_knots = c(-2.2, 2.2)),
#   chains = 2L, iter_warmup = 400L, iter_sampling = 400L,
#   refresh = 0L
# )
# 
# co <- coef(fit_bs)
# co

## ----het-default-api, eval=FALSE----------------------------------------------
# fit <- gdpar(
#   gdpar_bf(y ~ a(x1), sigma ~ a(x2)),
#   data   = d,
#   family = gdpar_family("gaussian"),
#   chains = 2L, iter_warmup = 400L, iter_sampling = 400L
# )

## ----het-api, eval=FALSE------------------------------------------------------
# fit_het <- gdpar(
#   gdpar_bf(y ~ a(x1), sigma ~ a(x2)),
#   data   = d,
#   family = list(
#     mu    = gdpar_family("gaussian"),
#     sigma = gdpar_family("beta")
#   ),
#   chains = 2L, iter_warmup = 400L, iter_sampling = 400L
# )

## ----het-example, eval=FALSE--------------------------------------------------
# set.seed(818L)
# n <- 100L
# x1 <- rnorm(n); x2 <- rnorm(n)
# mu_true   <- 0.3 + 0.6 * (x1 - mean(x1))
# sigma_eta <- 0.5 + 0.3 * (x2 - mean(x2))
# # Beta-distributed sigma in (0, 1) via inverse-logit of sigma_eta
# sigma_p <- 1 / (1 + exp(-sigma_eta))
# y <- rbeta(n, shape1 = 2 + 5 * sigma_p, shape2 = 5 - 4 * sigma_p)
# d <- data.frame(y = y, x1 = x1, x2 = x2)
# 
# fit_het <- gdpar(
#   gdpar_bf(y ~ a(x1), sigma ~ a(x2)),
#   data   = d,
#   family = list(
#     mu    = gdpar_family("beta"),
#     sigma = gdpar_family("beta")
#   ),
#   chains = 2L, iter_warmup = 400L, iter_sampling = 400L,
#   refresh = 0L
# )
# 
# co <- coef(fit_het)
# str(co, max.level = 2L)

## ----combined-api, eval=FALSE-------------------------------------------------
# fit_combo <- gdpar(
#   gdpar_bf(y ~ a(x1) + W(), sigma ~ a(x2)),
#   data   = d,
#   family = list(
#     mu    = gdpar_family("gaussian"),
#     sigma = gdpar_family("beta")
#   ),
#   W      = W_basis(type = "bspline", knots = c(-1, 0, 1),
#                    degree = 3L, boundary_knots = c(-2.5, 2.5)),
#   chains = 2L, iter_warmup = 400L, iter_sampling = 400L
# )