---
title: "Comparing and Evaluating Changepoint Methods with ggchangepoint"
author: "Youzhi Yu<br><span style='font-size:85%;'>University of Chicago</span>"
bibliography: vignette_reference.bib
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Comparing and Evaluating Changepoint Methods with ggchangepoint}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  message = FALSE,
  warning = FALSE,
  fig.alt = "ggchangepoint plot comparing changepoint detection methods on a time series"
)
library(ggchangepoint)
library(ggplot2)
theme_set(theme_light())

# Optional (Suggests) engines: gate the chunks that need them so the
# vignette builds on any installation.
has_fpop        <- requireNamespace("fpop", quietly = TRUE)
has_wbs         <- requireNamespace("wbs", quietly = TRUE)
has_stepR       <- requireNamespace("stepR", quietly = TRUE)
has_strucchange <- requireNamespace("strucchange", quietly = TRUE)

# The comparison and benchmarking sections run whichever detectors are
# installed. The fallback set comes entirely from the changepoint package (a
# hard dependency), so both branches compare multiple-changepoint detectors
# and the narrative below holds either way.
cmp_methods <- if (has_fpop && has_wbs) {
  c("pelt", "binseg", "fpop", "wbs")
} else {
  c("pelt", "binseg", "segneigh")
}
```

# Abstract

No single changepoint detection method dominates across signal shapes, noise
regimes, and computational budgets: large-scale evaluations find that the
ranking of algorithms is highly dataset-dependent and that default
configurations are often far from optimal [@van2020evaluation;
@truong2020selective]. Sound practice therefore requires running several
detectors, comparing their outputs, and — when ground truth is available —
scoring them with well-defined accuracy metrics. This article presents the
comparison and evaluation toolkit of **ggchangepoint**: visual comparison
across methods (`ggcpt_compare()`), a common tidy representation of competing
segmentations (`ggcpt_compare_table()`), a metrics module implementing
precision/recall/F1 under one-to-one matching, the covering metric, Hausdorff
distance, and the adjusted Rand index (`cpt_metrics()`), multi-annotator
scoring in the style of the Turing Change Point Dataset benchmark
(`cpt_metrics_annotated()`), and visual evaluation against ground truth
(`ggcpt_eval()`). We further discuss three complementary tools for
quantifying the *uncertainty* of a segmentation — engine-native confidence
intervals, bootstrap stability profiles (`cpt_stability()`), and the CROPS
penalty path (`cpt_crops()`) — and close with a small simulation-based
benchmarking workflow built from the package's ground-truth generators.

# Introduction

The changepoint literature offers a wide menu of detectors — penalised
optimal partitioning, binary segmentation and its wild and
narrowest-over-threshold refinements, moving-sum statistics, nonparametric
divergence measures, Bayesian posteriors — each with its own inductive bias
[@aminikhanghahi2017survey; @truong2020selective]. Benchmarks that score many
algorithms on many series find no uniform winner: performance depends on the
kind of change (mean, variance, distribution), the noise (Gaussian,
heavy-tailed, autocorrelated), the number and spacing of changes, and the tuning of
penalties and thresholds [@van2020evaluation]. Two practical consequences
follow. First, an analyst should *compare* several methods on the data at
hand rather than trust one default. Second, when ground truth (or expert
annotation) exists, comparison should be quantitative, using metrics with
agreed conventions.

**ggchangepoint** supports both activities behind one interface. Every
detector returns the same tidy `ggcpt` object, so competing methods are
directly comparable; the comparison module renders them side by side, and the
evaluation module scores them. This article is the package's methods-paper
treatment of that workflow. The companion vignette
`vignette("introduction", package = "ggchangepoint")` documents the full
detection surface; here we take detectors as given and focus on comparison,
evaluation, and uncertainty.

# Problem setup

Let $y_{1:n} = (y_1, \dots, y_n)$ be an ordered sequence. A segmentation with
$m$ changepoints is an ordered set $\tau_{0:m+1}$ with
$0 = \tau_0 < \tau_1 < \cdots < \tau_m < \tau_{m+1} = n$, partitioning the
index set into segments $A_j = \{\tau_{j-1}+1, \dots, \tau_j\}$,
$j = 1, \dots, m+1$. Most offline methods minimise a penalised cost
$$
\sum_{j=1}^{m+1} \mathcal{C}\bigl(y_{(\tau_{j-1}+1):\tau_j}\bigr) + \beta f(m),
$$
where $\mathcal{C}(\cdot)$ is a segment cost and $\beta f(m)$ guards against
over-segmentation [@yao1988estimating; @killick2012pelt].

Throughout the package a changepoint $\tau$ is reported in the **"left"
convention**: $\tau$ is the *last index of the left segment*, the convention
of the **changepoint** package [@killick2012pelt]. Engines that natively
report the first index of the right segment (e.g. the **ecp** family,
@matteson2014nonparametric) are normalised on the way in, so comparisons
across methods are always like for like.

Comparing segmentations is harder than comparing point estimates for three
reasons. First, the *number* of detected changes varies across methods, so a
metric must handle unequal-length sets. Second, a detection a few indices
away from a true change is usually acceptable, so metrics need a tolerance
margin — and a matching rule that prevents one true change from being
"claimed" by several detections. Third, changepoint sets induce
*partitions*, and two very different-looking point sets can induce similar
partitions; good practice therefore reports both point-based metrics
(precision/recall, Hausdorff) and partition-based metrics (covering, Rand)
[@van2020evaluation].

# Visual comparison

We simulate a series with three mean levels (two changepoints) and run
several detectors through the unified dispatcher. `ggcpt_compare()` accepts
the raw series and a vector of method names, runs `cpt_detect()` for each,
and renders the results. The default `layout = "facet"` draws one panel per
method:

```{r compare-facet}
set.seed(2024)
x <- c(rnorm(150, 0), rnorm(150, 3), rnorm(200, 1))
ggcpt_compare(x, methods = cmp_methods)
```

A method that finds *no* changepoints keeps its panel — the series is drawn
with no vertical lines on it — rather than silently disappearing. A method
that ran and found nothing is a result, not a missing value:

```{r compare-nochange}
set.seed(7)
x_null <- rnorm(300)
ggcpt_compare(x_null, methods = c("pelt", "binseg"))
```

The `layout = "overlay"` variant superimposes all methods in a single panel
with colour-coded rules, which makes both kinds of disagreement easier to
see: small differences in an estimated location, and outright differences in
how many changepoints a method reports.

```{r compare-overlay}
ggcpt_compare(x, methods = cmp_methods, layout = "overlay")
```

For a numeric rather than visual comparison, `ggcpt_compare_table()` returns
one tidy tibble with a row per (method, changepoint) pair; a method that
found nothing contributes a single row with `cp = NA`, which is the table's
version of keeping the empty panel:

```{r compare-table}
ggcpt_compare_table(x, methods = cmp_methods)
```

Because every row carries the same columns, this table pipes directly into
**dplyr**/**ggplot2** summaries — counting detections per method, plotting
location agreement, and so on [@wickham2016ggplot2]. When the detectors are
slow, `ggcpt_compare()` honours `future::plan()` (through **future.apply**)
and fits them in parallel.

# Accuracy metrics

When the true changepoints $\mathcal{T} = \{t_1, \dots, t_K\}$ are known,
`cpt_metrics(pred, truth, n, margin)` scores a predicted set
$\mathcal{P} = \{p_1, \dots, p_M\}$ with the following quantities.

**Precision, recall, and F1 under one-to-one matching.** A prediction $p$
*matches* a truth $t$ if $|p - t| \le \texttt{margin}$. Matching is
one-to-one: predictions are scanned in increasing order and each takes the
earliest unmatched truth within the margin, so each truth is claimed by at
most one prediction (for points on a line this greedy rule attains a maximum
matching). With $\mathrm{TP}$ matched pairs,
$$
\mathrm{precision} = \frac{\mathrm{TP}}{|\mathcal{P}|}, \qquad
\mathrm{recall} = \frac{\mathrm{TP}}{|\mathcal{T}|}, \qquad
F_1 = \frac{2\,\mathrm{precision}\cdot\mathrm{recall}}
           {\mathrm{precision} + \mathrm{recall}},
$$
with $F_1 = 0$ by convention when precision and recall are both zero. Because
$\mathrm{TP}$ counts *pairs*, all three quantities lie in $[0, 1]$: three
predictions crowded around one true change score one true positive, not three.

**The covering metric.** Following @van2020evaluation, let $\mathcal{S}$ and
$\mathcal{S}'$ be the partitions induced by the truth and the prediction, and
$J(A, A') = |A \cap A'| / |A \cup A'|$ the Jaccard index of two segments. The
covering of $\mathcal{S}$ by $\mathcal{S}'$ is
$$
\mathrm{cov}(\mathcal{S}, \mathcal{S}') =
  \frac{1}{n} \sum_{A \in \mathcal{S}} |A| \,
  \max_{A' \in \mathcal{S}'} J(A, A'),
$$
a weighted average of the best overlap achieved for each true segment.

**Hausdorff distance.** The worst-case location error, in index units,
$\max\{\max_p \min_t |p - t|,\; \max_t \min_p |p - t|\}$; it is `NA` when
either set is empty, since there is no distance to a nonexistent point.

**Adjusted Rand index.** The chance-corrected agreement of the two induced
segment labellings; 1 for identical partitions, 0 for chance-level
agreement.

**Annotation error and matched location errors.** The absolute difference in
counts $\bigl||\mathcal{P}| - |\mathcal{T}|\bigr|$, and the MAE/RMSE of the
matched location pairs — the latter two `NA` when nothing matched, because an
average over no pairs is not zero error.

```{r metrics-basic}
# perfect detection
cpt_metrics(pred = c(150, 300), truth = c(150, 300), n = 500)

# near misses within the margin still match one-to-one
cpt_metrics(pred = c(148, 305), truth = c(150, 300), n = 500, margin = 5)

# three predictions around one truth: one TP, precision 1/3
cpt_metrics(pred = c(148, 150, 152), truth = c(150), n = 500, margin = 5)
```

Three edge-case conventions deserve emphasis, because getting them wrong
silently corrupts benchmark averages:

- **Both sets empty.** Predicting "no changepoints" for a series with no
  changepoints is *exactly right*, so precision, recall, F1, covering, and
  the Rand index all equal 1. (Hausdorff distance and the matched-location
  errors remain `NA`: there are no pairs to measure.)
- **Empty prediction, non-empty truth.** An empty changepoint set still
  induces a perfectly well-defined partition — the trivial one with a single
  segment — so covering scores it by segment overlap rather than awarding an
  automatic 0. Precision and recall are 0, and the adjusted Rand index is 0
  because the trivial partition agrees with the truth only at chance level.
- **Out-of-range indices.** Changepoint locations must lie in
  $\{1, \dots, n-1\}$ under the left convention; anything outside is dropped
  with a warning rather than corrupting the partition construction, and the
  metrics are then computed on what remains.

```{r metrics-edges}
# a correct "no change" answer is rewarded
cpt_metrics(pred = integer(0), truth = integer(0), n = 300)

# empty prediction against one true change at the midpoint: the trivial
# one-segment partition still overlaps half the series, so covering is 0.5
cpt_metrics(pred = integer(0), truth = c(150), n = 300)
```

```{r metrics-outofrange, warning = TRUE}
# index 700 cannot be a changepoint of a length-500 series
cpt_metrics(pred = c(100, 700), truth = c(100, 300), n = 500)
```

# Multi-annotator evaluation

For real data, "ground truth" is often a set of human annotations that
disagree with one another. The Turing Change Point Dataset benchmark
[@van2020evaluation] therefore scores a prediction against *each* annotator
and averages. `cpt_metrics_annotated()` implements this convention: it takes
a list of annotation vectors and returns the averaged precision, recall, F1,
and covering.

```{r annotated}
annotations <- list(
  ann1 = c(150, 300),
  ann2 = c(152, 301),
  ann3 = c(149)          # a third annotator missed the second change
)
cpt_metrics_annotated(pred = c(150, 300), annotations, n = 500, margin = 5)
```

Averaging *per annotator* rather than pooling the annotations matters. The
prediction matches both of the first two annotators completely; against the
third it finds the one change that annotator marked but also reports a second
one, scoring precision $1/2$ and recall $1$. The averages are therefore
precision $5/6$ and recall $1$: the disagreement is charged to precision,
which is the honest place for it, since the extra detection may well be real
and merely unannotated.

# Visual evaluation

`ggcpt_eval()` overlays predictions and ground truth on the series, shades
the tolerance window around each true change, and colours each prediction as
a true positive or false positive, with missed truths drawn as dashed "Miss"
rules. It uses the *same one-to-one matching* as `cpt_metrics()`, so the
picture and the numbers always agree:

```{r eval-plot, fig.alt = "Series with predictions coloured as true positives and false positives, shaded tolerance windows around each true changepoint, and missed truths as dashed rules"}
truth <- c(150, 300)
pred <- c(151, 240)      # one hit, one false alarm, one miss
ggcpt_eval(pred, truth, data_vec = x, margin = 5)
cpt_metrics(pred, truth, n = length(x), margin = 5)
```

One true positive out of two predictions and two truths gives precision,
recall, and F1 all equal to $0.5$ — the single blue rule, the single orange
rule, and the single dashed rule in the plot, counted.

# Uncertainty beyond point sets

A segmentation is a point estimate. Three complementary tools quantify how
much confidence it deserves.

## Engine-native confidence intervals

Some engines deliver genuine confidence statements for changepoint
*locations*. SMUCE [@frick2014smuce] controls, at level $\alpha$, the
probability of overestimating the number of changepoints, and returns a
confidence interval for every location; the Bai–Perron dynamic program
[@bai2003computation; @zeileis2002strucchange] returns break-date intervals
for regression breaks. Four wrapped engines report such intervals — `smuce`,
`hsmuce`, `strucchange`, and `segmented` — as `ci_lower`/`ci_upper` columns on
the `ggcpt` changepoints tibble, which `autoplot(show_ci = TRUE)` draws as
whiskers below the series. `show_fit = TRUE` adds the engine's own fitted
signal, available here and from `decafs`, `cpop`, `segmented`, `bcp`, and
`beast`:

```{r smuce-ci, eval = has_stepR, fig.alt = "Series with the SMUCE step fit overlaid and confidence-interval whiskers for each estimated changepoint location"}
res_smuce <- smuce_wrapper(x)
tidy(res_smuce)
autoplot(res_smuce, show_ci = TRUE, show_fit = TRUE)
```

```{r strucchange-ci, eval = has_strucchange}
res_bp <- strucchange_wrapper(x)
tidy(res_bp)
```

Such intervals are worth reading alongside each other: two engines can agree
on where a change is and still, under their different noise models, disagree
about how tightly the location is pinned down.

## Bootstrap stability

Most engines ship no intervals at all. `cpt_stability()` provides a cheap,
model-agnostic substitute: it fits the detector once, resamples residuals
*within* the fitted segments (preserving the estimated regime structure),
re-runs the detector on each replicate, and reports the proportion of
replicates that re-detect a change within `margin` of each index. Locations
that survive resampling are trustworthy; locations that appear in only a
fraction of replicates are fragile.

```{r stability, fig.alt = "Bootstrap detection-frequency profile along the series index, with the original changepoints marked as dashed rules"}
st <- cpt_stability(x, method = "pelt", B = 30, seed = 1)
st
autoplot(st)
```

Both changepoints here are re-detected in every replicate, so the point
estimate is as stable as this diagnostic can report. The profile is also
informative where the print method is silent: a broad, low plateau marks a
region the detector keeps splitting somewhere without agreeing where.

## Penalty-path sensitivity: CROPS

Penalised methods commit to one penalty $\beta$, and the segmentation can
change qualitatively as $\beta$ moves. CROPS [@haynes2017computationally]
computes *every* optimal segmentation as the penalty ranges over an interval,
at roughly one PELT run per distinct solution — penalty selection as a
diagnostic rather than a guess. `cpt_crops()` returns the full path:

```{r crops}
path <- cpt_crops(x)
path
```

Over the default range $[\log n, 10 \log n]$ this series admits only two
distinct segmentations, and the two-changepoint solution holds over all but
the very bottom of it. That insensitivity is itself the diagnostic: on a
signal this cleanly separated, the answer does not depend on the penalty.
Lowering `pen_min` opens up the rest of the path, where the elbow plot shows
the cost reduction per additional changepoint flattening out and the
segmentation facets show the candidate models being chosen among:

```{r crops-elbow, fig.alt = "CROPS cost elbow: segmentation cost plotted against the number of changepoints for each solution on the penalty path"}
path_wide <- cpt_crops(x, pen_min = 4)
autoplot(path_wide)                          # cost elbow
```

```{r crops-segs, fig.alt = "The candidate CROPS segmentations, faceted by number of changepoints, each panel showing the series with that solution's changepoints"}
autoplot(path_wide, type = "segmentations")
```

Both `cpt_stability()` and `cpt_crops()` operate on a single numeric series
and reject multi-column input; for a panel of series, loop over the columns
(or use `cpt_batch()` for the detection step).

# A benchmarking workflow

The package's simulation module closes the loop: generators with *known*
truth feed the metrics module, so methods can be scored over replications.
`cpt_simulate()` draws series with specified changepoints under Gaussian,
Student-$t$, AR(1), or random-walk noise, and the canonical test signals of
the literature ship as ready-made generators: `signal_blocks()`
[the Donoho–Johnstone blocks signal, @donoho1994ideal], `signal_fms()`,
`signal_mix()`, `signal_teeth()`, and `signal_stairs()`. Every generator
attaches its true changepoints as a `true_changepoints` attribute, which is
exactly the `truth` argument `cpt_metrics()` expects.

```{r signals, fig.height = 4, fig.alt = "Donoho-Johnstone blocks test signal with its eleven true changepoints marked by vertical rules"}
blocks <- signal_blocks(n = 500, seed = 3)
ggplot(blocks, aes(index, value)) +
  geom_line(colour = "grey40") +
  geom_vline(xintercept = attr(blocks, "true_changepoints"),
             colour = "blue", linewidth = 0.3) +
  labs(title = "signal_blocks(): Donoho-Johnstone blocks with true changepoints")
```

A minimal benchmark: three methods, two signal-to-noise regimes, ten
replications each, scored by precision, recall, F1, and covering. The two
regimes differ only in the size of the jumps — 2.5 noise standard deviations
against 1.0 — so any difference between the panels is attributable to
difficulty alone. (Everything stays sequential and small here; `cpt_batch()`
runs one method over many series and honours `future::plan()` for parallel
execution in larger studies.)

```{r benchmark}
methods <- cmp_methods[1:3]
n_rep <- 10
regimes <- list(
  "high SNR (jump 2.5 sd)" = c(0, 2.5, 0.5),
  "low SNR (jump 1.0 sd)"  = c(0, 1.0, 0.3)
)

results <- do.call(rbind, lapply(names(regimes), function(g) {
  do.call(rbind, lapply(seq_len(n_rep), function(r) {
    dat <- cpt_simulate(400, changepoints = c(130, 260), change_in = "mean",
                        params = regimes[[g]], sd = 1, seed = 100 + r)
    truth <- attr(dat, "true_changepoints")
    do.call(rbind, lapply(methods, function(m) {
      fit <- cpt_detect(dat$value, method = m)
      score <- cpt_metrics(fit$changepoints$cp, truth, n = 400, margin = 5)
      cbind(tibble::tibble(regime = g, rep = r, method = m),
            score[, c("precision", "recall", "f1", "covering")])
    }))
  }))
}))

# average over replications, within regime
res_summary <- aggregate(cbind(precision, recall, f1, covering) ~ method + regime,
                         data = results, FUN = mean)
knitr::kable(res_summary, digits = 3,
             caption = "Mean accuracy over 10 replications per regime (margin = 5).")
```

Two regimes are enough to reproduce the claim this article opened with. In the
high-SNR panel the problem is close to solved and the methods are hard to tell
apart; shrinking the jump to one noise standard deviation lowers all four
metrics for every method and opens a visible gap between them, so an ordering
read off the easy regime need not survive into the hard one. Accuracy is a
property of the (method, signal, noise) triple rather than of the method
alone, which is why this loop is worth running on the series at hand instead
of adopting a ranking from the literature.

The same skeleton scales to the studies of @van2020evaluation: more
generators (heavy-tailed and autocorrelated noise via
`cpt_simulate(noise = "t")` and `noise = "ar1"`), more methods (everything in
`cpt_methods()`), more replications, and parallel execution via
`future::plan(multisession)`.

# Discussion

Comparison and evaluation are not afterthoughts in changepoint analysis; they
are how an analyst earns confidence in a segmentation. The design of
**ggchangepoint**'s toolkit follows three principles. First, *a common
representation makes comparison trivial*: because every engine returns the
same `ggcpt` contract in the same location convention, `ggcpt_compare()` and
the metrics module work for all 31 wired methods without special cases.
Second, *metrics must agree with their pictures*: `ggcpt_eval()` and
`cpt_metrics()` share one matching routine, so a plotted true positive is a
counted true positive. Third, *uncertainty deserves first-class treatment*:
engine-native intervals, bootstrap stability, and penalty paths give three
independent views of how much a reported changepoint should be trusted, and
all three render directly with **ggplot2** [@wickham2016ggplot2].

For the detection surface itself — the dispatcher, the engine wave, the
Bayesian displays, and the multivariate tools — see
`vignette("introduction", package = "ggchangepoint")` and the feature tour in
`vignette("ggchangepoint", package = "ggchangepoint")`.

# References