Example scripts: manual calibration

script <- system.file("examples", "binary-options.R", package = "vennabers")
file.show(script)
source(script)

# After installation: source(system.file("examples", "binary-options.R", package = "vennabers"))
library(vennabers)

# Precomputed probabilities from a fixed model on its separate calibration set.
# No calibration or test truths were used to train that model or preprocessing.
truth <- c(0, 1, 0, 1)
calibration_probability <- c(.501, .502, .509, .510)
new_probability <- c(.499, .505, .511)
inductive <- venn_abers(truth, calibration_probability, precision = 2)
for (loss in c("log", "brier")) {
  prediction <- predict(inductive, new_probability, loss = loss)
  classes <- predict(inductive, new_probability, loss = loss, type = "class")
  stopifnot(identical(classes, prediction$probability >= .5))
  print(prediction)
}
# Precision affects calibration only; choosing a new precision requires a fit.
# The new values above are deliberately not rounded to two decimal places.

# To select class 0 as the event, reverse truth AND the supplied probabilities.
# Round only after selecting the event; rounded results need not be complements.
other_event <- venn_abers(truth == 0, 1 - calibration_probability, precision = 2)
print(predict(other_event, 1 - new_probability))
# Inductive bounds are VA endpoints, not confidence intervals.

# Cross uses held-out probabilities from each corresponding retained fold model.
# These illustrative folds have unequal sizes; each still receives equal weight.
# See cross-venn-abers.R for complete model fitting and held-out data generation.
cross <- cross_venn_abers(
  truth = c(0, 1, 1, 1), probability = c(.501, .201, .501, .801),
  fold = c("A", "B", "B", "B"), precision = 2
)
# Each column comes from the model that excluded its named calibration fold.
# All columns refer to the same new observations, in the same order.
new_probability <- cbind(B = c(.201, .501, .901), A = c(.199, .501, .801))
test_truth <- c(0, 1, 1)  # Used only for assessment, after fitting.

local({
  plot_file <- tempfile(fileext = ".pdf")
  on.exit(unlink(plot_file), add = TRUE)
  grDevices::pdf(plot_file)
  on.exit(grDevices::dev.off(), add = TRUE, after = FALSE)
  for (loss in c("log", "brier")) {
    details <- predict(cross, new_probability, loss = loss, type = "details")
    reconstructed <- vapply(seq_len(nrow(new_probability)), function(row) {
      f <- details$folds[details$folds$row == row, ]
      if (loss == "log") {
        plogis(mean(log(f$upper)) - mean(log1p(-f$lower)))
      } else mean(f$probability)
    }, 0)
    stopifnot(max(abs(reconstructed - details$prediction$probability)) < 1e-10,
              identical(predict(cross, new_probability, loss = loss, type = "class"),
                        details$prediction$probability >= .5))
    print(details)
    print(c(brier = brier_score(test_truth, reconstructed),
            log = log_loss(test_truth, reconstructed)))
    calibration_plot(test_truth, reconstructed)
  }
})
# Classes use the combined probability with TRUE at .5, never a fold vote.
# Fold details concern new observations; there is no aggregate Cross interval.
# No improvement in test scores or perfect Cross calibration is promised.

# Run from the checkout root after installation:
# Rscript inst/examples/cross-venn-abers.R
library(vennabers)
source(system.file("examples", "assertions.R", package = "vennabers"))

set.seed(20260906)
observations <- data.frame(x1 = rnorm(900L), x2 = rnorm(900L))
observations$truth <- rbinom(900L, 1L, stats::plogis(
  -0.4 + observations$x1 - 0.6 * observations$x2
))  # The event is 1 (equivalently TRUE).

# Set aside test observations before any model fitting or calibration.
test_rows <- sample.int(nrow(observations), 300L)
test <- observations[test_rows, ]
development <- observations[-test_rows, ]
fold_ids <- c("A", "B", "C")
fold <- sample(rep(fold_ids, length.out = nrow(development)))
models <- setNames(vector("list", length(fold_ids)), fold_ids)
out_of_fold_probability <- numeric(nrow(development))

for (id in fold_ids) {
  held_out <- fold == id
  # Exclusion must cover the complete learned pipeline, including any learned
  # preprocessing or tuning. This example needs neither.
  models[[id]] <- stats::glm(truth ~ x1 + x2,
    data = development[!held_out, ], family = stats::binomial()
  )
  # Restore held-out probabilities to their original development positions.
  out_of_fold_probability[held_out] <- predict(models[[id]],
    newdata = development[held_out, ], type = "response"
  )
}

# Retain these models: each predicts the same test rows in the same order.
# A single all-data refit cannot replace them. Names cannot prove model
# provenance, training exclusion, unique observations, or semantic alignment.
test_probability <- vapply(models, function(model) {
  predict(model, newdata = test[c("x1", "x2")], type = "response")
}, numeric(nrow(test)))

# The package owns per-fold calibration and equal-weight log-loss aggregation.
# Model fitting and fold creation remain caller responsibilities.
calibrator <- cross_venn_abers(
  development$truth, out_of_fold_probability, fold
)
prediction <- predict(calibrator, test_probability)
assert_conditions(
  is.data.frame(prediction), nrow(prediction) == nrow(test),
  identical(names(prediction), "probability"),
  is.double(prediction$probability), all(is.finite(prediction$probability)),
  all(prediction$probability >= 0 & prediction$probability <= 1)
)

# Test truth is used only for evaluation. No score improvement is promised.
scores <- c(
  brier_score = brier_score(test$truth, prediction$probability),
  log_loss = log_loss(test$truth, prediction$probability)
)
stopifnot(all(is.finite(scores)))
print(scores)

# A temporary PDF also works without an interactive display. Close and remove
# it even on error. For an interactive plot, call calibration_plot() directly.
local({
  plot_file <- tempfile(fileext = ".pdf")
  on.exit(unlink(plot_file), add = TRUE)
  grDevices::pdf(plot_file)
  on.exit(grDevices::dev.off(), add = TRUE, after = FALSE)
  summary <- calibration_plot(test$truth, prediction$probability)
  stopifnot(sum(summary$count) == nrow(test))
})

# Cross output is a regularized event probability, with no public aggregate
# interval. Individual fold diagnostics are available with type = "details".
# The combined probability inherits no IVAP perfect-calibration
# guarantee and makes no confidence-interval coverage claim.

# After package installation: Rscript inst/examples/inductive-ovo.R
library(vennabers)

# Illustrative caller-supplied outputs from one already fixed multiclass model.
# Calibration and assessment observations are separate from that model's
# training, preprocessing, and tuning; the package does not train the model.
classes <- c("cat", "dog", "bird")
truth <- factor(rep(classes, 2), levels = classes)
p <- rbind(c(.7, .2, .1), c(.2, .6, .2), c(.1, .2, .7),
           c(.5, .3, .2), c(.3, .5, .2), c(.2, .3, .5))
colnames(p) <- classes
fit <- venn_abers(truth, p, precision = 2, method = "ovo")
print(fit)
new_p <- rbind(c(.6, .3, .1), c(.1, .4, .5), c(.1, .8, .1))
colnames(new_p) <- classes
assessment <- factor(c("cat", "bird", "dog"), levels = classes)

local({
  plot_file <- tempfile(fileext = ".pdf")
  grDevices::pdf(plot_file)
  on.exit({ grDevices::dev.off(); unlink(plot_file) })
  for (loss in c("log", "brier")) {
    details <- predict(fit, new_p[, rev(classes)], loss = loss, type = "details")
    probability <- predict(fit, new_p, loss = loss)
    predicted_class <- predict(fit, new_p, loss = loss, type = "class")
    reconstructed <- probability
    for (row in seq_len(nrow(new_p))) {
      pairs <- details$pairs[details$pairs$row == row, ]
      for (class in classes) {
        q <- c(pairs$probability[pairs$event == class],
               1 - pairs$probability[pairs$non_event == class])
        reconstructed[row, class] <- 1 / (sum(1 / q) - (length(classes) - 2))
      }
    }
    reconstructed <- reconstructed / rowSums(reconstructed)
    stopifnot(max(abs(reconstructed - probability)) <= 1e-10,
              identical(probability, details$prediction),
              identical(predicted_class,
                        factor(classes[max.col(probability, ties.method = "first")], levels = classes)))
    print(details)
    print(predicted_class)
    print(c(brier = brier_score(assessment, probability), log = log_loss(assessment, probability)))
    groups <- calibration_plot(assessment, probability)
    stopifnot(all(tapply(groups$count, groups$class, sum) == length(assessment)))
    print(groups)
  }
})
# Pair endpoints have only their binary interpretation. Assessment scores use
# the final matrix; Brier sums class errors without dividing by class count.

# After package installation: Rscript inst/examples/inductive-ovr.R
library(vennabers)

# Illustrative probabilities from an already fixed multiclass model.
# Calibration and assessment are separate from its training and tuning.
classes <- c("cat", "dog", "bird")
truth <- factor(rep(classes, 2), levels = classes)
p <- rbind(c(.7, .2, .1), c(.2, .6, .2), c(.1, .2, .7),
           c(.5, .3, .2), c(.3, .5, .2), c(.2, .3, .5))
colnames(p) <- classes
fit <- venn_abers(truth, p, method = "ovr", precision = 2)
print(fit)
new_p <- rbind(c(.6, .3, .1), c(.1, .4, .5), c(.1, .8, .1))
colnames(new_p) <- classes
assessment <- factor(c("cat", "bird", "dog"), levels = classes)

local({
  plot_file <- tempfile(fileext = ".pdf")
  grDevices::pdf(plot_file)
  on.exit({ grDevices::dev.off(); unlink(plot_file) })
  for (loss in c("log", "brier")) {
    details <- predict(fit, new_p[, rev(classes)], loss = loss, type = "details")
    probability <- predict(fit, new_p, loss = loss)
    predicted_class <- predict(fit, new_p, loss = loss, type = "class")
    points <- vapply(classes, function(class) {
      parts <- details$classes[details$classes$event == class, ]
      parts$probability[match(seq_len(nrow(new_p)), parts$row)]
    }, numeric(nrow(new_p)))
    stopifnot(max(abs(points / rowSums(points) - probability)) <= 1e-10,
              identical(probability, details$prediction),
              identical(predicted_class,
                        factor(classes[max.col(probability, ties.method = "first")], levels = classes)))
    print(details)
    print(predicted_class)
    print(c(brier = brier_score(assessment, probability), log = log_loss(assessment, probability)))
    groups <- calibration_plot(assessment, probability)
    stopifnot(all(tapply(groups$count, groups$class, sum) == length(assessment)))
    print(groups)
  }
})
# Endpoints concern each binary class-against-rest problem. Scores and plots
# use the final normalized matrix. Manual Brier is an R extension over Python.

# After package installation: Rscript inst/examples/cross-multiclass.R
library(vennabers)

# Precomputed out-of-fold probabilities from two retained multiclass models.
# Model A excluded rows 1:2; model B excluded rows 3:6, including preprocessing
# and tuning. Retain both models yourself. No model is trained by this package.
classes <- c("cat", "dog", "bird")
truth <- factor(c("cat", "dog", "cat", "bird", "dog", "bird"), levels = classes)
fold <- c("A", "A", "B", "B", "B", "B")
p <- rbind(c(.7, .2, .1), c(.2, .6, .2), c(.5, .3, .2),
           c(.1, .2, .7), c(.3, .5, .2), c(.2, .3, .5))
colnames(p) <- classes
# The same three untouched assessment observations, in the same row order,
# predicted by each corresponding retained model. List order is immaterial.
new_p <- list(B = rbind(c(.6, .3, .1), c(.1, .4, .5), c(.1, .8, .1)),
              A = rbind(c(.5, .4, .1), c(.2, .2, .6), c(.2, .7, .1)))
new_p <- lapply(new_p, function(x) { colnames(x) <- classes; x })
assessment <- factor(c("cat", "bird", "dog"), levels = classes)

local({
  plot_file <- tempfile(fileext = ".pdf")
  grDevices::pdf(plot_file)
  on.exit({ grDevices::dev.off(); unlink(plot_file) })
  for (method in c("ovo", "ovr")) {
    fit <- cross_venn_abers(truth, p, fold, method = method, precision = 2)
    print(fit)
    for (loss in c("log", "brier")) {
      details <- predict(fit, new_p, loss = loss, type = "details")
      probability <- predict(fit, new_p, loss = loss)
      predicted_class <- predict(fit, new_p, loss = loss, type = "class")
      # Reconstruct from the fold endpoints/points, then combine subproblems.
      parts <- details[[2L]]
      problems <- unique(parts[c(if (method == "ovo") "non_event", "event")])
      merged <- vapply(seq_len(nrow(problems)), function(i) {
        selected <- parts$event == problems$event[i]
        if (method == "ovo") selected <- selected & parts$non_event == problems$non_event[i]
        vapply(seq_len(nrow(probability)), function(row) {
          f <- parts[selected & parts$row == row, ]
          if (loss == "log") {
            plogis(mean(log(f$upper)) - mean(log1p(-f$lower)))
          } else mean(f$probability)
        }, 0)
      }, numeric(nrow(probability)))
      recovered <- if (method == "ovr") merged else vapply(classes, function(class) {
        q <- cbind(1 - merged[, problems$non_event == class, drop = FALSE],
                   merged[, problems$event == class, drop = FALSE])
        1 / (1 + rowSums((1 - q) / q))
      }, numeric(nrow(probability)))
      recovered <- recovered / rowSums(recovered)
      stopifnot(max(abs(recovered - probability)) <= 1e-10,
                identical(probability, details$prediction),
                identical(predicted_class,
                          factor(classes[max.col(probability, ties.method = "first")],
                                 levels = classes)))
      print(details)
      print(predicted_class)
      print(c(brier = brier_score(assessment, probability),
              log = log_loss(assessment, probability)))
      groups <- calibration_plot(assessment, probability, bins = 1)
      stopifnot(max(abs(groups$mean_probability - colMeans(probability))) <= 1e-10,
                all(groups$count == length(assessment)), all(groups$event_rate == 1 / 3))
      print(groups)
    }
  }
})
# One shared fold plan serves all pairs/classes. Fold A lacks bird but every
# OvO pair remains nonempty. Empty pairs fail with the pair and fold identified.
# Endpoints describe binary subproblems, with no aggregate confidence interval.
# Manual Cross multiclass is an R extension; Python has no matching public
# manual whole-procedure API. Numerical tests reconstruct binary components.

# After package installation: Rscript inst/examples/one-hot-classes.R
library(vennabers)

# One-hot matrices represent class decisions the predictor has already made.
# venn_abers_one_hot() never chooses a class from probabilities itself.

# Binary: the logical class decision is TRUE (the event) at probability >= 0.5,
# so an exact 0.5 tie sets the TRUE column.
binary <- venn_abers(c(FALSE, TRUE, FALSE, TRUE), c(.2, .2, .8, .8))
new_scores <- c(.2, .5, .8)
probability <- predict(binary, new_scores)$probability
binary_class <- predict(binary, new_scores, type = "class")
binary_one_hot <- venn_abers_one_hot(binary_class)
print(cbind(probability, binary_one_hot))
stopifnot(probability[[2L]] == .5,
          identical(colnames(binary_one_hot), c("FALSE", "TRUE")),
          identical(binary_one_hot[, "TRUE"], as.double(binary_class)))

# Multiclass: columns follow the stored factor levels, even for classes that
# no row predicts, and an exact tie keeps the first stored level.
classes <- c("cat", "dog", "bird")
truth <- factor(rep(classes, 2), levels = classes)
p <- rbind(c(.7, .2, .1), c(.2, .6, .2), c(.1, .2, .7),
           c(.5, .3, .2), c(.3, .5, .2), c(.2, .3, .5))
colnames(p) <- classes
fit <- venn_abers(truth, p, method = "ovo")
new_p <- rbind(c(.6, .3, .1), c(.1, .4, .5), c(1, 1, 1) / 3)
colnames(new_p) <- classes
multiclass_class <- predict(fit, new_p, type = "class")
multiclass_one_hot <- venn_abers_one_hot(multiclass_class)
print(multiclass_one_hot)
# The third row is an exact tie, so its one sits in the first level "cat".
stopifnot(length(unique(predict(fit, new_p)[3L, ])) == 1L,
          identical(multiclass_one_hot[3L, ], c(cat = 1, dog = 0, bird = 0)))
stopifnot(identical(colnames(multiclass_one_hot), classes),
          all(rowSums(multiclass_one_hot) == 1),
          identical(multiclass_class, factor(classes[max.col(multiclass_one_hot)], levels = classes)))

# Missing decisions give missing rows and keep the row order.
print(venn_abers_one_hot(factor(c("bird", NA, "cat"), levels = classes)))

# Workflow, Grid, parsnip and frozen predictions pass their `.pred_class` table:
#   venn_abers_one_hot(predict(fitted_workflow, new_data, type = "class"))
# Schema-encoded tables keep the encoded class names; decode first for originals:
#   venn_abers_one_hot(venn_abers_decode(prediction, schema)$class)

# After installation: source(system.file("examples", "inductive-regression.R",
#                                       package = "vennabers"))
library(vennabers)
source(system.file("examples", "assertions.R", package = "vennabers"))

# Independent IID observations; roles are fixed before examining their outcomes.
set.seed(20260907)
x <- runif(240, -3, 3)
data <- data.frame(id = seq_along(x), x = x, truth = 2 + x + .4 * x^2 + rnorm(240))
training <- data[1:120, ]
calibration <- data[121:180, ]
assessment <- data[181:240, ]
stopifnot(!anyDuplicated(c(training$id, calibration$id, assessment$id)))

# Only training data determine this model; no preprocessing or tuning uses the
# calibration or assessment observations. The caller retains the fixed model.
model <- lm(truth ~ x, data = training)
calibration_prediction <- predict(model, calibration)
fit <- venn_abers_regression(calibration$truth, calibration_prediction, epsilon = .1)
uncalibrated <- predict(model, assessment)
result <- predict(fit, uncalibrated)
print(list(requested = fit$requested, m = fit$m, epsilon = fit$epsilon))

# Assessment truths enter only here. These are point-prediction scores in target
# units; the VA bounds concern regularized conditional expectations, not coverage
# of individual outcomes. Neither improvement nor a general selector guarantee
# for this particular point merger is promised.
scores <- vapply(list(uncalibrated = uncalibrated, calibrated = result$prediction),
                 function(prediction) {
                   c(RMSE = sqrt(mean((assessment$truth - prediction)^2)),
                     MAE = mean(abs(assessment$truth - prediction)))
                 }, c(RMSE = 0, MAE = 0))
assert_conditions(all(is.finite(scores)), nrow(result) == nrow(assessment),
          all(result$lower <= result$prediction), all(result$prediction <= result$upper),
          abs(scores["RMSE", "calibrated"]^2 -
              sum((assessment$truth - result$prediction)^2) / nrow(assessment)) < 1e-12,
          abs(scores["MAE", "calibrated"] -
              sum(abs(assessment$truth - result$prediction)) / nrow(assessment)) < 1e-12)
print(scores)

# After installation: source(system.file("examples", "cross-regression.R",
#                                       package = "vennabers"))
library(vennabers)

# IID observations; fold membership and the separate assessment role are fixed
# before examining outcomes. Unequal folds still receive equal aggregate weight.
set.seed(20260908)
x <- runif(240, -3, 3)
data <- data.frame(id = seq_along(x), x = x, truth = 2 + x + .4 * x^2 + rnorm(240))
development <- data[1:180, ]
assessment <- data[181:240, ]
fold <- sample(rep(c("A", "B", "C"), c(50, 60, 70)))
ids <- unique(fold)
models <- setNames(vector("list", length(ids)), ids)
out_of_fold <- rep(NA_real_, nrow(development))
newdata <- matrix(NA_real_, nrow(assessment), length(ids), dimnames = list(NULL, ids))
for (id in ids) {
  held_out <- fold == id
  training <- development[!held_out, ]
  calibration <- development[held_out, ]
  stopifnot(!any(training$id %in% c(calibration$id, assessment$id)),
            !any(calibration$id %in% assessment$id))
  # No learned preprocessing or tuning outside this training set. The caller
  # retains this exact model for both its calibration and assessment predictions.
  models[[id]] <- lm(truth ~ x, data = training)
  out_of_fold[held_out] <- predict(models[[id]], calibration)
  newdata[, id] <- predict(models[[id]], assessment)
}
fit <- cross_venn_abers_regression(development$truth, out_of_fold, fold, epsilon = .1)
details <- predict(fit, newdata, type = "details")
print(lapply(fit$calibrators, function(x) x[c("observations", "requested", "m", "epsilon")]))
result <- details$prediction
for (column in names(result)) {
  reconstructed <- Reduce(`+`, lapply(details$folds, `[[`, column)) / length(ids)
  stopifnot(all(abs(result[[column]] - reconstructed) <= 1e-10 + 1e-10 * abs(reconstructed)))
}

# Matching names cannot prove model provenance or semantic row identity. Each
# column above uses the same assessment rows in the same order. A model refitted
# on all development data cannot replace these fold models. Predict trains nothing.
# Assessment truths enter only in scoring. Aggregated bounds are descriptive
# fold summaries, with no additional coverage or selector guarantee; improvement
# in point-prediction scores is not promised.
scores <- vapply(list(uncalibrated = rowMeans(newdata), calibrated = result$prediction),
                 function(prediction) {
                   c(RMSE = sqrt(mean((assessment$truth - prediction)^2)),
                     MAE = mean(abs(assessment$truth - prediction)))
                 }, c(RMSE = 0, MAE = 0))
stopifnot(all(is.finite(scores)), nrow(result) == nrow(assessment),
          all(result$lower <= result$prediction), all(result$prediction <= result$upper))
for (column in colnames(scores)) {
  residual <- assessment$truth - if (column == "calibrated") result$prediction else rowMeans(newdata)
  stopifnot(abs(scores["RMSE", column]^2 - sum(residual^2) / length(residual)) < 1e-12,
            abs(scores["MAE", column] - sum(abs(residual)) / length(residual)) < 1e-12)
}
print(scores)

# After installation: source(system.file("examples", "calibration-graphics.R",
#                                       package = "vennabers"))
library(vennabers)
source(system.file("examples", "assertions.R", package = "vennabers"))

# Fit diagnoses draw the stored state of a calibrator and need no new truth
# values. Assessment graphics need held-out pairs of truth and prediction.
set.seed(20261021)
observations <- data.frame(x = rnorm(900L))
observations$truth <- rbinom(900L, 1L, stats::plogis(0.3 + 1.5 * observations$x))
training <- observations[1:300, ]
calibration <- observations[301:600, ]
assessment <- observations[601:900, ]

model <- stats::glm(truth ~ x, data = training, family = stats::binomial())
calibration_probability <- predict(model, calibration, type = "response")
assessment_probability <- predict(model, assessment, type = "response")
calibrator <- venn_abers(calibration$truth, calibration_probability)
probability <- predict(calibrator, assessment_probability)$probability

# A temporary PDF also works without an interactive display. Close and remove
# it even on error. For interactive plots, call the functions directly.
local({
  plot_file <- tempfile(fileext = ".pdf")
  on.exit(unlink(plot_file), add = TRUE)
  grDevices::pdf(plot_file)
  on.exit(grDevices::dev.off(), add = TRUE, after = FALSE)

  # Fit diagnosis: the calibration map and the Venn-Abers width upper - lower
  # at every stored knot. The width is neither a confidence level nor a
  # nominal coverage.
  map <- plot(calibrator, which = c("map", "width"))
  stored <- predict(calibrator, map$score)
  assert_conditions(
    identical(names(map), c("score", "lower", "upper", "probability", "width")),
    identical(map$width, map$upper - map$lower),
    identical(map$probability, stored$probability),
    identical(range(map$score), c(0, 1))
  )
  print(head(map))

  # Cross fits overlay the map of every fold calibrator; the table keeps one
  # block per fold, not the aggregated Cross probability.
  fold <- rep(c("A", "B", "C"), length.out = nrow(calibration))
  cross <- cross_venn_abers(calibration$truth, calibration_probability, fold)
  fold_maps <- plot(cross, loss = "brier")
  stopifnot(identical(unique(fold_maps$fold), c("A", "B", "C")))

  # Assessment graphic: the CORP reliability diagram of the held-out pairs.
  # The default consistency band shows where the estimate of perfectly
  # calibrated probabilities would fall; a seed makes it reproducible without
  # changing the caller's random stream.
  state <- .Random.seed
  corp <- corp_plot(assessment$truth, probability, seed = 1)
  assert_conditions(
    identical(.Random.seed, state),
    identical(corp, corp_plot(assessment$truth, probability, seed = 1)),
    identical(corp$probability, sort(unique(probability))),
    sum(corp$count) == nrow(assessment),
    all(diff(corp$estimate) >= 0),
    all(corp$band_lower >= 0 & corp$band_upper <= 1)
  )
  print(head(corp))

  # The decomposition splits the Brier score of the same pairs into
  # miscalibration, discrimination and uncertainty.
  parts <- corp_decomposition_plot(assessment$truth, probability)
  assert_conditions(
    abs(parts$brier - brier_score(assessment$truth, probability)) < 1e-12,
    abs(parts$brier - (parts$mcb - parts$dsc + parts$unc)) < 1e-12
  )
  print(parts)

  # Regression: the map over the uncalibrated prediction with reference lines
  # at the moderated truths a and b, then held-out pairs against the
  # calibrated points.
  x <- runif(600L, 0, 10)
  regression_data <- data.frame(x = x, truth = x + sin(x) + rnorm(600L, sd = 0.5))
  regression_model <- stats::lm(truth ~ x, data = regression_data[1:200, ])
  uncalibrated <- predict(regression_model, regression_data)
  regression <- venn_abers_regression(regression_data$truth[201:400],
                                      uncalibrated[201:400])
  regression_map <- plot(regression, which = c("map", "width"))
  stopifnot(identical(regression_map$width,
                      regression_map$upper - regression_map$lower))
  point <- predict(regression, uncalibrated[401:600])$prediction
  held_out_truth <- regression_data$truth[401:600]
  held_out_pairs <- regression_plot(held_out_truth, point,
                                    which = c("observed", "residual"))
  stopifnot(identical(held_out_pairs$residual, held_out_truth - point))
})

# CORP bands are pointwise resampling bands for the assessment pairs given.
# They are not Venn-Abers lower and upper probabilities, and neither graphic
# makes a statement about Venn-Abers validity.