## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)

## -----------------------------------------------------------------------------
library(gp3bayes)
pupil_advanced_capabilities()
pupil_advanced_compatibility_table()

## -----------------------------------------------------------------------------
sim <- simulate_advanced_pupil_timecourse(
  n_participants = 12,
  trials_per_participant = 4,
  time_points = 31,
  family = "gaussian",
  ar = 0.45,
  heteroskedastic_strength = 0.35,
  missing_fraction = 0.04,
  seed = 2026
)

sim
plot_advanced_pupil_simulation(sim)

## -----------------------------------------------------------------------------
temporal_audit <- audit_pupil_temporal_dependence(sim$data)
temporal_audit
plot_pupil_temporal_dependence(temporal_audit)

## -----------------------------------------------------------------------------
distribution <- specify_pupil_distribution(
  family = "gaussian",
  residual_scale = "condition_time"
)

gp <- create_pupil_gp_spec(
  kernel = "matern32",
  basis = "approximate",
  k = 25
)

spec <- specify_advanced_pupil_timecourse_model(
  prepared = sim$data,
  temporal_structure = "gaussian_process",
  distribution = distribution,
  gp_spec = gp,
  autocorrelation = "none",
  covariates = c("baseline_pupil", "luminance"),
  predictive_target = "future_segment"
)

spec
pupil_advanced_specification_table(spec)
plot_pupil_model_complexity(spec)

## ----eval=FALSE---------------------------------------------------------------
# translated <- translate_advanced_pupil_model_to_brms(spec)
# translated
# 
# fit <- fit_advanced_pupil_model_backend(
#   spec,
#   backend = "cmdstanr",
#   chains = 4,
#   iter = 2000,
#   warmup = 1000,
#   cores = 2,
#   seed = 2026
# )
# 
# diagnose_advanced_pupil_fit(fit)
# trajectory <- predict_advanced_pupil_trajectory(fit)
# plot_advanced_pupil_trajectory(trajectory)
# 
# sigma <- estimate_pupil_residual_scale(fit)
# plot_pupil_residual_scale(sigma)

