Penalised cost minimisation

The classical formulation chooses the segmentation minimising a penalised cost, \[ \min_{m,\ \tau_{1:m}} \; \sum_{i=1}^{m+1} \mathcal{C}\!\left(y_{(\tau_{i-1}+1):\tau_i}\right) \;+\; \beta m, \] where \(\mathcal{C}\) is a segment cost (for a change in mean under Gaussian noise the residual sum of squares, more generally twice the negative maximised log-likelihood of the segment) and \(\beta > 0\) is the price of each additional changepoint. Writing \(k\) for the number of parameters a changepoint introduces, the familiar choices are \(\beta = 2k\) (AIC) and \(\beta = k \log n\) (BIC, or SIC in the changepoint literature) (Yao 1988), alongside the strengthened and modified variants discussed below. Solved naively by dynamic programming, the minimisation costs \(O(n^2)\); PELT (Killick, Fearnhead, and Eckley 2012) prunes candidate changepoints to reach linear expected cost while remaining exact, and FPOP (Maidstone et al. 2017) reaches comparable speed by functional pruning.

A practical consequence of writing the objective this way is easy to miss: \(\beta\) and \(\mathcal{C}\) must live on the same scale. For a change in mean the changepoint engines evaluate the Normal cost with the noise standard deviation fixed at 1, and fpop penalises the residual sum of squares directly, so multiplying the data by a constant multiplies the cost while leaving \(\beta\) untouched. On a 200-point series with a single changepoint whose jump is five standard deviations, pelt recovers exactly one changepoint at \(\sigma = 1\) but returns 39 at \(\sigma = 3\) and 141 at \(\sigma = 10\) (means over 20 draws, since a single draw is not stable at these settings). Standardise the series, pass a penalty on the data’s own scale (say 2 * log(n) * var(diff(x)) / 2), or use change_in = "meanvar", which estimates a variance per segment. Methods that estimate the noise level as part of their procedure (SMUCE, the WBS family, CPOP, bcp, BEAST and the nonparametric engines) return the same segmentation whatever the units. Three exceptions are worth knowing: geomcp runs PELT on its mapped series and inherits its sensitivity, DeCAFS floors its noise estimate near 0.03, and BOCPD’s default prior is on the data’s own scale, so the last two miss changes in a series measured in very small units.

Search-based and multiscale methods

A complementary family locates changepoints by scanning test statistics. Binary segmentation (Scott and Knott 1974; Vostrikova 1981) recursively splits the series at the maximal CUSUM statistic; wild binary segmentation (Fryzlewicz 2014) and its successor WBS2 (Fryzlewicz 2020) draw random subintervals so that short segments are not masked; narrowest-over-threshold (NOT) (Baranowski, Chen, and Fryzlewicz 2019) favours the narrowest interval on which the contrast exceeds a threshold, which generalises cleanly to changes in slope; MOSUM (Eichinger and Kirch 2018) scans a moving-sum statistic at a fixed bandwidth, or across a range of bandwidths; Isolate-Detect (Anastasiou and Fryzlewicz 2022) isolates each changepoint in an expanding interval; and TGUH (Fryzlewicz 2018) performs a tail-greedy bottom-up merge. SMUCE (Frick, Munk, and Sieling 2014) occupies a special place: it estimates the step function with the fewest jumps that still passes a simultaneous multiscale test at level \(\alpha\), and in doing so delivers confidence intervals for every changepoint location, uncertainty statements most competitors cannot make. HSMUCE (Pein, Sieling, and Munk 2017) extends this to heterogeneous noise.

Beyond the mean

Changes need not be in the mean: the change_in argument accepts "mean", "var", "meanvar", "slope", "distribution", "covariance", "network", "regression" and "seasonality", the nine values cpt_methods() lists in its supports column, each routable to the methods that declare it. Nonparametric engines (energy statistics (Matteson and James 2014), nonparametric cost functions (Haynes, Fearnhead, and Eckley 2017), kernel running statistics (Arlot, Celisse, and Harchaoui 2019; Cabrieto et al. 2018), joint characteristic functions (McGonigle and Cho 2025), self-normalisation (Z. Zhao, Jiang, and Shao 2022)) detect distributional change without likelihood assumptions; Bayesian engines (Barry and Hartigan 1993; Adams and MacKay 2007; K. Zhao et al. 2019) return posteriors instead of point sets; high-dimensional engines (Wang and Samworth 2018; Chen, Wang, and Samworth 2022; Grundy, Killick, and Mihaylov 2020) aggregate evidence across coordinates; and regression engines (Bai and Perron 1998; Muggeo 2003) date breaks in model coefficients. The tour below visits each family with runnable code.

The ggcpt result contract

Every detector returns an object of class ggcpt containing:

• changepoints: a tibble with one row per changepoint. Columns cp (location, “left” convention) and cp_value (the data value at cp) are always present; engines add ci_lower/ci_upper (SMUCE, HSMUCE, strucchange, segmented, bfast, taylor, mcp), posterior_prob (bcp, BEAST), detection_time (CPM), strength (inspect), declared_at (ocd), or mapping (geomcp) when they have more to say. nsp reports its uncertainty as a significance region instead: region_start, region_end and the regions slot.
• segments: a tibble of the induced segments (seg_id, start, end, n, param_estimate).
• data: the analysed series as a tibble (index, value), plus a fitted column when the engine estimates a signal (SMUCE, HSMUCE, CPOP, bcp, BEAST, DeCAFS, segmented, mcp, bfast; these are the engines cpt_methods() marks in its fitted column).
• method, change_in, penalty (a list(type, value) descriptor), cp_convention (always "left"), runtime (elapsed seconds, timed by cpt_detect() and NA when a wrapper is called directly), and fit (the untouched upstream object, for experts).

Multivariate results additionally carry a data_wide tibble with one column per coordinate, which autoplot() renders as faceted small-multiples. Two methods plot a derived series instead, because for them the coordinates are not what a reader wants to see: network takes a sequence of adjacency matrices and reports mean edge weight (one facet per matrix entry would be unreadable), and hdreg plots the response it regressed on the covariates. Both carry no data_wide.

Tidy methods and the plotting layer

The class implements the full complement of generics R users expect:

set.seed(2022)
x <- c(rnorm(100, 0, 1), rnorm(100, 10, 1))
res <- cpt_detect(x, method = "pelt", change_in = "mean")
res
#> ggcpt (changepoint detection result)
#>   Method:             pelt
#>   Change in:          mean
#>   Changepoints found: 1
#>   CP convention:      left
#>   Penalty:            MBIC
#>   Series length:      200
#> 
#> Changepoints:
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100    0.467
tidy(res)
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100    0.467
glance(res)
#> # A tibble: 1 × 9
#>       n n_changepoints method change_in penalty_type penalty_value cp_convention
#>   <int>          <int> <chr>  <chr>     <chr>                <dbl> <chr>        
#> 1   200              1 pelt   mean      MBIC                    NA left         
#> # ℹ 2 more variables: total_cost <dbl>, runtime <dbl>
head(augment(res))
#> # A tibble: 6 × 6
#>   index  value seg_id .fitted .resid is_changepoint
#>   <int>  <dbl>  <int>   <dbl>  <dbl> <lgl>         
#> 1     1  0.900      1   0.139  0.761 FALSE         
#> 2     2 -1.17       1   0.139 -1.31  FALSE         
#> 3     3 -0.897      1   0.139 -1.04  FALSE         
#> 4     4 -1.44       1   0.139 -1.58  FALSE         
#> 5     5 -0.331      1   0.139 -0.470 FALSE         
#> 6     6 -2.90       1   0.139 -3.04  FALSE

autoplot() draws the series, the changepoint rules and, on request, the fitted segment means (show_segments), the engine’s fitted signal (show_fit), and changepoint-location confidence intervals (show_ci):

autoplot(res, show_segments = TRUE)

Composable layers (geom_changepoint(), geom_cpt_segment(), geom_cpt_ci(), stat_changepoint()), a theme (theme_ggcpt()), and segment shading (annotate_segments()) let the same results be built into bespoke graphics; summary(), as_tibble(), as.data.frame(), format(), and plot() complete the S3 surface.

Design principles

The four goals above are recorded in the package as principles P1. Wrap, don’t reinvent (bind to peer-reviewed CRAN engines), P2. Tidy in, tidy out (stable column names across all methods), P3. ggplot2 all the way down (every result renders and extends), and P4. One vocabulary (x, method, change_in, penalty, ...). Three further principles govern how the interface evolves:

• P5. Progressive disclosure: beginners call cpt_detect() + autoplot(); experts reach the upstream fit via $fit.
• P6. No surprises: the 0.1.0 functions still work unchanged.
• P7. Document everything you ship: every export is introduced in the README and a vignette.

Release 0.4.0 adds an eighth: P8. Carry the uncertainty. Where a method quantifies uncertainty, the ggcpt object records it and autoplot() can draw it.

Release 0.5.0 adds a ninth: P9. Be extensible from the outside. A detector this package does not wrap, cannot wrap, or has never heard of can join the same grammar through as_ggcpt() and cpt_register_method(), and is labelled as user-supplied wherever it appears. See vignette("extending", package = "ggchangepoint").

Penalised and optimal partitioning

PELT (Killick, Fearnhead, and Eckley 2012), binary segmentation (Scott and Knott 1974), segment neighbourhoods (Auger and Lawrence 1989), and at-most-one-change (AMOC) (Hinkley 1970) come from the changepoint package (Killick and Eckley 2014); FPOP (Maidstone et al. 2017) from fpop:

tidy(cpt_detect(x_multi, method = "pelt"))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90
tidy(cpt_detect(x_multi, method = "binseg"))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   201    0.337

tidy(cpt_detect(x_multi, method = "fpop"))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90

The Achilles heel of penalised methods is the choice of \(\beta\). Rather than committing to one value, cpt_crops() computes every optimal segmentation as \(\beta\) ranges over an interval (the CROPS algorithm of Haynes, Eckley and Fearnhead (2017), as implemented by changepoint) and turns penalty selection into a diagnostic:

path <- cpt_crops(x_multi)
path
#> ggcpt_path (CROPS penalty path)
#>   Change in:              mean
#>   Penalty range:          [5.704, 57.04]
#>   Series length:          300
#>   Distinct segmentations: 5
#> 
#> # A tibble: 5 × 3
#>   penalty n_cpts  cost
#>     <dbl>  <int> <dbl>
#> 1    6.69      2  327.
#> 2    6.27      5  307.
#> 3    6.19      8  288.
#> 4    5.76      9  282.
#> 5    5.70     10  276.
autoplot(path)

The default plot puts segmentation cost against model size, and the usual reading takes the model beyond which the cost stops falling appreciably. The sweep over the default interval \([\log n,\, 10 \log n]\) admits only a handful of distinct segmentations here, and the most parsimonious of them already recovers the two true changepoints. autoplot(path, type = "segmentations") shows the candidate models themselves, and autoplot(path, type = "path") the map from penalty to model size:

autoplot(path, type = "segmentations")

The modern fastcpd engine (Li and Zhang 2024) brings the same penalised formulation to a wide family of models (mean, variance, mean-and-variance, and AR/ARMA/GARCH model changes) with sequential-gradient-descent speed:

tidy(fastcpd_wrapper(x_multi, family = "mean"))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90

Multiscale and search methods

The randomised and multiscale searchers are one call each, whether through cpt_detect() or through the wrapper directly:

tidy(wbs_wrapper(x_multi, seed = 1))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90

tidy(not_wrapper(x_multi, seed = 1))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90

tidy(mosum_wrapper(x_multi))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90

tidy(idetect_wrapper(x_multi, seed = 1))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90
tidy(wbs2_wrapper(x_multi))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   200    1.90
tidy(tguh_wrapper(x_multi))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100   -0.100
#> 2   199    4.53

SMUCE (Frick, Munk, and Sieling 2014) is the family’s inferential flagship: its level \(\alpha\) bounds the probability of overestimating the number of changepoints (the default is alpha = 0.5, stepR’s own recommendation for estimation rather than testing), and every location comes with a confidence interval, stored in ci_lower/ci_upper and drawn by show_ci = TRUE as whiskers near the foot of the panel (the step fit is drawn by show_fit = TRUE):

res_smuce <- smuce_wrapper(x_multi)
tidy(res_smuce)
#> # A tibble: 2 × 4
#>      cp cp_value ci_lower ci_upper
#>   <int>    <dbl>    <int>    <int>
#> 1   100   -0.100       98      103
#> 2   200    1.90       199      202
autoplot(res_smuce, show_ci = TRUE, show_fit = TRUE)

For heterogeneous noise, smuce_wrapper(x, family = "hsmuce") (or cpt_detect(x, method = "hsmuce")) runs HSMUCE (Pein, Sieling, and Munk 2017).

Changes in slope

A kink in the trend is not a jump in the level, and running a mean-change detector on a trending series over-detects notoriously. CPOP (Fearnhead, Maidstone, and Letchford 2019; Fearnhead and Grose 2024) solves the change-in-slope problem exactly under an \(L_0\) penalty, returning a continuous piecewise-linear fit:

res_cpop <- cpop_wrapper(x_slope)
tidy(res_cpop)
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100     38.6
autoplot(res_cpop, show_fit = TRUE)

NOT with its linear contrast (Baranowski, Chen, and Fryzlewicz 2019) offers a search-based alternative; the dispatcher routes cpt_detect(x, method = "not", change_in = "slope") to it automatically:

tidy(cpt_detect(x_slope, method = "not", change_in = "slope"))
#> # A tibble: 2 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1    85     36.8
#> 2   101     40.7

Bayesian detection

The Barry-Hartigan product partition model (Barry and Hartigan 1993), via the bcp package (Erdman and Emerson 2007), returns a posterior probability of a changepoint at every location along with posterior segment means. Locations clearing prob_threshold populate the changepoints tibble (with their probabilities), and ggcpt_posterior() draws the classic two-panel display:

res_bcp <- bcp_wrapper(x_mean, seed = 2026)
tidy(res_bcp)
#> # A tibble: 1 × 3
#>      cp cp_value posterior_prob
#>   <int>    <dbl>          <dbl>
#> 1   100    0.369              1
ggcpt_posterior(res_bcp)

Bayesian online changepoint detection (Adams and MacKay 2007) instead tracks the posterior over the current run length (the time elapsed since the last change), updating it recursively as each observation arrives. Its signature graphic is the run-length heatmap, in which a change shows up as the posterior mass falling back to a run length of zero:

res_bocpd <- bocpd_wrapper(x_mean)
tidy(res_bocpd)
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100    0.369
ggcpt_runlength(res_bocpd)

A third Bayesian engine, BEAST (K. Zhao et al. 2019) via Rbeast, averages over models rather than conditioning on one, and is wired as cpt_detect(x, method = "beast") (or beast_wrapper()); it too reports posterior_prob and renders with ggcpt_posterior().

Nonparametric and sequential detection

When no parametric form is trustworthy, the nonparametric cost approach of changepoint.np (Haynes, Fearnhead, and Eckley 2017) and the energy-statistics E-Divisive of ecp (Matteson and James 2014; James and Matteson 2014) detect general distributional change:

set.seed(2022)
tidy(cpt_detect(x_mean, method = "np"))
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100    0.369
tidy(cpt_detect(x_mean, method = "ecp", seed = 1))
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100    0.369

The cpm package (Ross 2015) recasts detection as a stream of two-sample tests (Mann-Whitney for location, Mood for scale, Lepage, Kolmogorov-Smirnov, Cramér-von Mises, and parametric variants), run here over the whole series in one pass to mimic an online monitor. Its results distinguish where a change happened (cp) from when it was detected (detection_time), the lag inherent in sequential monitoring:

tidy(cpm_wrapper(x_mean, cpm_type = "Mann-Whitney"))
#> # A tibble: 3 × 3
#>      cp cp_value detection_time
#>   <int>    <dbl>          <int>
#> 1    23   -1.39              30
#> 2    50    0.426             67
#> 3   100    0.369            104

Three further nonparametric engines are wired and worth knowing: kernel change-point analysis on running statistics (kcp_wrapper(), engine kcpRS), which detects changes in running means, variances, autocorrelations, or correlations (Arlot, Celisse, and Harchaoui 2019; Cabrieto et al. 2018); NP-MOJO (npmojo_wrapper(), engine CptNonPar), which detects changes in the marginal or lagged joint distribution while remaining valid under serial dependence (McGonigle and Cho 2025); and self-normalised segmentation (sn_wrapper(), engine SNSeg), which avoids long-run variance estimation altogether and tests changes in means, variances, autocorrelations, or bivariate correlations (Z. Zhao, Jiang, and Shao 2022).

Robustness to drift, autocorrelation, and model ambiguity

The most common failure of mean-change detection in practice is not a subtle statistical one: it is running a Gaussian-mean detector on data whose baseline drifts or whose noise is autocorrelated, and then reporting a changepoint wherever the model is wrong. DeCAFS (Romano et al. 2022) models exactly this regime (abrupt changes superimposed on random-walk drift and AR(1) noise) and separates the two:

res_decafs <- decafs_wrapper(x_mean)
tidy(res_decafs)
#> # A tibble: 1 × 2
#>      cp cp_value
#>   <int>    <dbl>
#> 1   100    0.369
autoplot(res_decafs, show_fit = TRUE)

EnvCpt (Beaulieu and Killick 2018) attacks the same confusion by model selection: it fits up to twelve competing descriptions (constant mean or linear trend, each with or without changepoints, and with white-noise, AR(1) or AR(2) errors) and reports changepoints only if a changepoint model wins on an information criterion:

res_env <- envcpt_wrapper(x_mean, models = c("mean", "meancpt", "trendcpt"))
glance(res_env)
#> # A tibble: 1 × 9
#>       n n_changepoints method change_in penalty_type penalty_value cp_convention
#>   <int>          <int> <chr>  <chr>     <chr>                <dbl> <chr>        
#> 1   200              1 envcpt mean      AIC: meancpt          568. left         
#> # ℹ 2 more variables: total_cost <dbl>, runtime <dbl>

The winning model’s name is recorded in the penalty descriptor (penalty_type above, here AIC: meancpt), so “no changepoints, it’s just autocorrelation” is a first-class answer.

Multivariate and high-dimensional detection

Multivariate methods accept a matrix (rows are time points) directly. The energy-statistics E-Divisive of ecp was built for this (Matteson and James 2014); for high-dimensional data whose change is confined to a sparse subset of coordinates, inspect (Wang and Samworth 2018) finds an optimal sparse projection of the CUSUM matrix and reports the projected evidence (strength). Multivariate results render as faceted small-multiples with shared changepoint rules:

set.seed(2026)
X <- cbind(a = c(rnorm(80), rnorm(80, 3)),
           b = c(rnorm(80), rnorm(80, -2)),
           c = rnorm(160))
res_hd <- inspect_wrapper(X)
tidy(res_hd)
#> # A tibble: 1 × 3
#>      cp cp_value strength
#>   <int>    <dbl>    <dbl>
#> 1    80    0.785     21.9
autoplot(res_hd)

Two further engines complete the family. geomcp_wrapper() (engine changepoint.geo) maps each observation to its distance from and its angle to a reference point, then segments the two mapped series, catching changes in magnitude and in orientation respectively (Grundy, Killick, and Mihaylov 2020). And ocd_wrapper() (engine ocd) monitors a high-dimensional stream online with worst-case detection-delay guarantees (Chen, Wang, and Samworth 2022). Because detection there is sequential, the locations it reports are declaration times (the change plus the detection delay), recorded in declared_at; the wrapper estimates the pre-change baseline from an initial training window and resets after each declaration so that several changes can be found. Like the method itself, it needs at least two coordinates and refuses a single series.

Structural breaks in regression

Econometric practice dates breaks in regression coefficients. The Bai-Perron estimator (Bai and Perron 1998, 2003), via strucchange (Zeileis et al. 2002), returns break dates with confidence intervals; called on a bare series it dates mean shifts, and called with a formula it dates breaks in arbitrary regressions:

res_bp <- strucchange_wrapper(x_mean)
tidy(res_bp)
#> # A tibble: 1 × 4
#>      cp cp_value ci_lower ci_upper
#>   <int>    <dbl>    <int>    <int>
#> 1   100    0.369       99      101
autoplot(res_bp, show_ci = TRUE)

Where the regression function is continuous (a kink rather than a jump), segmented (Muggeo 2003, 2008) estimates broken-line relationships with standard errors for the breakpoints:

res_seg <- segmented_wrapper(x_slope, npsi = 1, seed = 1)
tidy(res_seg)
#> # A tibble: 1 × 4
#>      cp cp_value ci_lower ci_upper
#>   <int>    <dbl>    <int>    <int>
#> 1   100     38.6       98      101
autoplot(res_seg, show_fit = TRUE, show_ci = TRUE)

Batch detection over many series

Applied work rarely stops at one series. cpt_batch() runs one detector over every column of a matrix or data frame (or every element of a list) and returns a tibble with one row per series, carrying both the tidy changepoints and the full ggcpt object in list-columns. With future and future.apply installed it honours a non-sequential future::plan(), using parallel-safe RNG:

set.seed(2026)
panel <- cbind(shifted = x_mean, quiet = rnorm(200))
batch <- cpt_batch(panel, method = "pelt")
batch
#> ggcpt_batch (2 series, method: pelt)
#> 
#> # A tibble: 2 × 2
#>   series  n_changepoints
#>   <chr>            <int>
#> 1 shifted              1
#> 2 quiet                0
tidy(batch)
#> # A tibble: 1 × 3
#>   series     cp cp_value
#>   <chr>   <int>    <dbl>
#> 1 shifted   100    0.369
autoplot(batch)

Stability diagnostics

Most engines report a point set with no measure of its fragility. cpt_stability() resamples residuals within the fitted segments (so the estimated regime structure is preserved), re-runs the detector on each replicate, and reports how often each location is re-detected: a cheap, model-agnostic confidence signal available for every engine, including the many that ship no intervals of their own:

st <- cpt_stability(x_mean, method = "pelt", B = 50, seed = 1)
st
#> ggcpt_stability (50 bootstrap replicates, method: pelt)
#> 
#> Original changepoints and their re-detection frequency:
#> # A tibble: 1 × 2
#>      cp stability
#>   <int>     <dbl>
#> 1   100         1
autoplot(st)

Evaluation, interactivity, and citations

When ground truth is known, cpt_metrics() computes precision, recall and F1 under one-to-one matching, the covering metric, Hausdorff distance, and adjusted Rand index, following the conventions of the modern benchmarking literature (van den Burg and Williams 2020); ggcpt_eval() draws the agreement, and ggcpt_compare() juxtaposes methods. These are the subject of the companion vignette vignette("comparison", package = "ggchangepoint").

Any result renders as an interactive HTML widget with ggcpt_interactive(res) (engine plotly, in Suggests); the static autoplot() path is untouched.

Finally, because every method here is someone’s published work, cpt_cite() returns the reference(s) behind a result, so analyses can cite the right paper without leaving R:

cpt_cite("pelt")
#> [pelt] Killick, R., Fearnhead, P. and Eckley, I. A. (2012). Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association, 107(500), 1590-1598.