## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----setup--------------------------------------------------------------------
library(bayesiansurpriser)
library(sf)

## ----uniform------------------------------------------------------------------
model_uniform <- bs_model_uniform()
print(model_uniform)

## ----baserate-----------------------------------------------------------------
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

# Create base rate model from birth counts
model_baserate <- bs_model_baserate(nc$BIR74)
print(model_baserate)

## ----gaussian-----------------------------------------------------------------
# Fit from data
model_gaussian_fit <- bs_model_gaussian()

# Fixed parameters
model_gaussian_fixed <- bs_model_gaussian(mu = 100, sigma = 20, fit_from_data = FALSE)
print(model_gaussian_fixed)

## ----sampled------------------------------------------------------------------
# Sample 20% of data for KDE
model_sampled <- bs_model_sampled(sample_frac = 0.2, kernel = "gaussian")
print(model_sampled)

## ----funnel-------------------------------------------------------------------
# Create funnel model with sample sizes
model_funnel <- bs_model_funnel(nc$BIR74, type = "count")
print(model_funnel)

## ----model-space--------------------------------------------------------------
# Default model space (uniform + baserate + funnel)
space_default <- default_model_space(nc$BIR74)
print(space_default)

## ----custom-space-------------------------------------------------------------
# Custom model space with explicit prior
space_custom <- model_space(
  bs_model_uniform(),
  bs_model_baserate(nc$BIR74),
  bs_model_gaussian(),
  bs_model_funnel(nc$BIR74),
  prior = c(0.1, 0.4, 0.2, 0.3),
  names = c("Uniform", "Population", "Normal", "Funnel")
)
print(space_custom)

## ----comparison---------------------------------------------------------------
# Only uniform model - highlights any non-uniform distribution
result_uniform <- surprise(nc, observed = SID74,
                          models = model_space(bs_model_uniform()))

# Uniform + baserate - highlights deviation from population proportion
result_pop <- surprise(nc, observed = SID74, expected = BIR74,
                      models = model_space(
                        bs_model_uniform(),
                        bs_model_baserate(nc$BIR74)
                      ))

# Full default - accounts for sampling variance too
result_full <- surprise(nc, observed = SID74, expected = BIR74)

# Compare max surprise values
cat("Uniform only, max surprise:", max(result_uniform$surprise), "\n")
cat("With base rate, max surprise:", max(result_pop$surprise), "\n")
cat("Full model space, max surprise:", max(result_full$surprise), "\n")

## ----posteriors---------------------------------------------------------------
space <- default_model_space(nc$BIR74)
updated_space <- bayesian_update(space, nc$SID74)

cat("Prior:", updated_space$prior, "\n")
cat("Posterior:", updated_space$posterior, "\n")