## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(probcal)
library(dplyr)

## ----data---------------------------------------------------------------------
set.seed(1001)
iris_binary <- iris |>
  filter(Species != "setosa") |>
  mutate(y = as.integer(Species == "virginica")) |>
  group_by(y) |>
  mutate(
    split = sample(rep(
      c("train", "calibration", "test"),
      times = c(25, 12, 13)
    ))
  ) |>
  ungroup()

iris_binary |>
  count(split, y)

## ----classifier---------------------------------------------------------------
train <- iris_binary |>
  filter(split == "train")

calibration <- iris_binary |>
  filter(split == "calibration")

test <- iris_binary |>
  filter(split == "test")

classifier <- glm(
  y ~ Sepal.Length + Sepal.Width,
  data = train,
  family = binomial()
)

calibration <- calibration |>
  mutate(raw_p = predict(classifier, calibration, type = "response"))

test <- test |>
  mutate(raw_p = predict(classifier, test, type = "response"))

## ----calibrators--------------------------------------------------------------
beta_fit <- cal_beta(calibration$raw_p, calibration$y)
platt_fit <- cal_platt(calibration$raw_p, calibration$y)

test <- test |>
  mutate(
    beta = predict(beta_fit, raw_p),
    platt = predict(platt_fit, raw_p)
  )

## ----metrics------------------------------------------------------------------
metric_table <- bind_rows(
  test |>
    summarise(method = "raw", ece = ece(raw_p, y, bins = 5),
              mce = mce(raw_p, y, bins = 5), ace = ace(raw_p, y, bins = 5)),
  test |>
    summarise(method = "beta", ece = ece(beta, y, bins = 5),
              mce = mce(beta, y, bins = 5), ace = ace(beta, y, bins = 5)),
  test |>
    summarise(method = "platt", ece = ece(platt, y, bins = 5),
              mce = mce(platt, y, bins = 5), ace = ace(platt, y, bins = 5))
) |>
  mutate(across(where(is.numeric), function(x) round(x, 3)))

metric_table

## ----diagram, fig.width = 6, fig.height = 5, fig.alt = "Reliability diagram for beta-calibrated iris probabilities, with binned points compared to the diagonal line."----
reliability_diagram(test$beta, test$y, bins = 5)