Nothing
R/monitor.R
In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'
Defines functions
control_chart
check_baseline_dependence
autoplot.ggcpt_delay
print.ggcpt_delay
glance.ggcpt_delay
tidy.ggcpt_delay
cpt_delay
cpt_replay
autoplot.ggcpt_monitor
print.ggcpt_monitor
tidy.ggcpt_monitor
alarms.ggcpt_monitor
alarms
monitor_relearn
monitor_width
cpt_update
chart_state
cpt_monitor
Documented in
alarms
alarms.ggcpt_monitor
autoplot.ggcpt_delay
autoplot.ggcpt_monitor
cpt_delay
cpt_monitor
cpt_replay
cpt_update
glance.ggcpt_delay
print.ggcpt_delay
print.ggcpt_monitor
tidy.ggcpt_delay
tidy.ggcpt_monitor
# ---------------------------------------------------------------------------
# Theme G: streaming and online monitoring.
#
# cpm, ocd and BOCPD are *online* methods that 0.4.0 ran in batch-replay mode
# and reported as if they were retrospective. The output the online
# literature cares about -- the alarm sequence, the detection delay, the
# average run length -- was not first-class, and there was no way to feed new
# observations to an existing fit. A monitor is a genuinely different object
# from a segmentation, so it gets its own class rather than another slot.
#
# `edetector` is the one method in this package implemented here rather than
# wrapped. That is a deliberate exception to the wrap-don't-implement rule
# (see ?cpt_monitor), taken because no R package implements e-detectors and
# the construction is short enough to be auditable in one screen.
# ---------------------------------------------------------------------------
#' A stateful sequential changepoint monitor
#'
#' Creates a detector that consumes observations as they arrive and raises
#' alarms, rather than segmenting a series that is already complete. Feed it
#' with \code{\link{cpt_update}()}, read its alarm log with
#' \code{\link{alarms}()}, and score it with \code{\link{cpt_delay}()},
#' because for an online method "did you find the location?" is the wrong
#' question and "how long did you take, and how often do you false-alarm?"
#' is the right one.
#'
#' @param method Which sequential detector:
#' \describe{
#' \item{\code{"edetector"}}{(default) a mixture Shiryaev-Roberts
#' e-detector; see the section below.}
#' \item{\code{"cpm"}}{\pkg{cpm}'s sequential change-point model, tuned
#' by \code{ARL0}.}
#' \item{\code{"ocd"}}{\pkg{ocd}'s high-dimensional online detector.
#' \strong{Multivariate only}: it tracks a projection of the whole
#' vector and needs at least two coordinates, so it refuses a single
#' series rather than falling back to a univariate statistic.}
#' \item{\code{"cusum"}}{a two-sided tabular CUSUM control chart
#' (Page 1954) on the standardised stream, with reference value
#' \code{k} and its decision interval calibrated to \code{arl0}.
#' Implemented in this package; see the section on control charts.}
#' \item{\code{"ewma"}}{an EWMA control chart (Roberts 1959) with
#' smoothing \code{lambda} and limits calibrated to \code{arl0}.
#' Implemented in this package.}
#' }
#' @param baseline A numeric vector (or, for \code{"ocd"}, a matrix with rows
#' as time points) of pre-change training data used to estimate the
#' in-control mean and scale. Required for \code{"edetector"} and
#' \code{"ocd"}; optional for \code{"cpm"}, which has its own start-up
#' period.
#' @param alpha Target false-alarm probability for \code{"edetector"}: the
#' average run length under the null is at least \code{1 / alpha}, by
#' optional stopping on the martingale \eqn{M_t - t} (see Details).
#' Defaults to \code{0.01}.
#' @param arl0 Target in-control average run length for \code{"cpm"},
#' \code{"cusum"} and \code{"ewma"}. Defaults to \code{500}.
#' @param cpm_type Statistic for \code{"cpm"}. Defaults to
#' \code{"Mann-Whitney"}.
#' @param patience Target patience (average run length) for \code{"ocd"}.
#' Defaults to \code{5000}.
#' @param deltas Shift sizes, in baseline standard deviations, mixed over by
#' \code{"edetector"}. Defaults to \code{c(0.5, 1, 2)}, each taken in both
#' directions. The mixture is a uniform average, so adding shifts costs
#' power at the ones already there rather than inflating the false-alarm
#' rate; mix over the range you think the change could fall in, not over
#' everything.
#' @param reset After an alarm, restart the detector (\code{TRUE}, the
#' default) or leave it running? Restarting is what makes a monitor report
#' \emph{repeated} changes rather than latching on the first one.
#' @param relearn How many observations after an alarm are used to re-learn
#' the in-control baseline, during which no further alarm can fire.
#' Defaults to \code{20}. This matters more than it looks: a real change
#' is \emph{persistent}, so a detector that restarts against the stale
#' pre-change baseline alarms again on the very next observation and keeps
#' alarming for the rest of the series: the monitor reports one change as
#' hundreds, and \code{\link{cpt_delay}()} then counts them all as false
#' alarms. Set \code{relearn = 0} to switch the behaviour off and see
#' every threshold crossing.
#' @param thresh Threshold rule for \code{"ocd"}: \code{"MC"} (default)
#' calibrates by Monte Carlo against \code{patience}, or supply a numeric
#' vector of three thresholds. The Monte Carlo calibration is the
#' expensive part of building an \code{"ocd"} monitor (a minute or more
#' at the default \code{patience}), so pass thresholds directly when you
#' already have them, or lower \code{mc_reps} while exploring.
#' @param mc_reps Monte Carlo repetitions for the \code{"ocd"} threshold.
#' @param ... Additional arguments passed to the engine's constructor
#' (\code{"cpm"} and \code{"ocd"}; the three native monitors take none).
#' @param k CUSUM reference value, in baseline standard deviations: half
#' the shift the chart is tuned to detect fastest. Defaults to \code{0.5}
#' (a one-standard-deviation shift).
#' @param lambda EWMA smoothing weight on the newest observation, in (0, 1].
#' Defaults to \code{0.2}; smaller values detect small shifts sooner and
#' large ones later.
#' @param limit Optional control limit, overriding the calibration to
#' \code{arl0}: the CUSUM decision interval \eqn{h}, or the EWMA limit
#' multiplier \eqn{L}, both in baseline standard deviations.
#'
#' @section Control charts:
#' \code{"cusum"} and \code{"ewma"} are the two charts statistical process
#' control runs, and they are sequential changepoint detectors under
#' another name. Each standardises the stream by the baseline's mean and
#' standard deviation (Phase I; the baseline should be in control) and
#' then monitors it (Phase II). The CUSUM keeps \eqn{S_t^+ =
#' \max(0, S_{t-1}^+ + z_t - k)} and \eqn{S_t^- = \max(0, S_{t-1}^- - z_t -
#' k)} and alarms when either reaches \eqn{h}; the EWMA keeps \eqn{z_t =
#' \lambda x_t + (1 - \lambda) z_{t-1}} and alarms when it leaves \eqn{\pm L
#' \sqrt{\lambda / (2 - \lambda)}}. The limit is solved for \code{arl0} by
#' the Markov-chain approximation of the run length (Brook and Evans 1972;
#' Lucas and Saccucci 1990), which reproduces the published tables: a CUSUM
#' with \code{k = 0.5} and \eqn{h = 4} or 5 has in-control run lengths of
#' 168 and 465, and an EWMA with \code{lambda = 0.1} and \eqn{L = 2.814}
#' one of 500. Measured through this function with a 5,000-observation
#' baseline and \code{arl0 = 100}, the first alarm on in-control data came
#' after 99.4 observations on average for the CUSUM and 99.8 for the EWMA
#' (400 runs each, standard errors about 5). The calibration treats the
#' baseline's mean and standard deviation as exact, so a short baseline
#' shortens the in-control run length: a few hundred observations is a
#' reasonable minimum. Rational subgroups, run rules and capability indices are not
#' modelled; for a formal quality-control programme use a dedicated
#' package such as \pkg{qcc}. \code{autoplot(plot_type = "control_chart")}
#' draws the chart.
#'
#' @section The e-detector, and why it is implemented rather than wrapped:
#' Shin, Ramdas and Rinaldo (2023) give a nonparametric sequential framework
#' with non-asymptotic control of the average run length, and it has no R
#' implementation. The construction used here is the mixture Shiryaev-Roberts
#' e-detector for a sub-Gaussian shift. For each candidate shift
#' \eqn{\delta} the increment is the likelihood ratio
#' \eqn{e_t^{(\delta)} = \exp(\delta (X_t - \mu_0)/\sigma^2 -
#' \delta^2/(2\sigma^2))}, which has unit mean under the null; the running
#' statistic is \eqn{R_t^{(\delta)} = (1 + R_{t-1}^{(\delta)})
#' e_t^{(\delta)}}; the shifts are combined by \emph{averaging},
#' \eqn{M_t = K^{-1} \sum_\delta R_t^{(\delta)}}; and an alarm is raised
#' the first time \eqn{M_t \ge 1/\alpha}.
#'
#' The averaging is not a detail. Under the null \eqn{M_t - t} is a
#' mean-zero martingale, so optional stopping at the alarm time \eqn{\tau}
#' gives \eqn{E_\infty[\tau] = E_\infty[M_\tau] \ge 1/\alpha}: a
#' finite-sample lower bound on the in-control average run length, with no
#' asymptotics and no calibration run. A convex combination of e-detectors
#' is an e-detector; a \emph{maximum} of them is not, and taking one
#' silently multiplies the false-alarm rate by roughly the number of shifts
#' mixed over. Apart from it and the two textbook control charts below,
#' every detector in this package wraps published, separately maintained
#' code; the three native monitors are labelled as such wherever they
#' appear.
#'
#' @return A \code{ggcpt_monitor} object: a list carrying the detector's
#' state between calls, so it is meant to be passed to
#' \code{\link{cpt_update}()} and re-assigned rather than read field by
#' field. The parts worth reading are
#' \describe{
#' \item{\code{alarms}}{the tibble \code{\link{alarms}()} and
#' \code{\link[=tidy.ggcpt_monitor]{tidy}()} return: one row per
#' alarm, with \code{time}, \code{statistic} and
#' \code{threshold}. Empty (zero rows) until something fires.}
#' \item{\code{t}}{how many observations the monitor has consumed
#' through \code{\link{cpt_update}()}, \strong{not} counting the
#' baseline; an alarm's \code{time} is on this same clock. A monitor
#' built by \code{\link{cpt_replay}()} also carries \code{offset},
#' the number of leading observations it trained on, so
#' \code{time + offset} is a position in the replayed series
#' (\code{\link{cpt_delay}()} applies it for you).}
#' \item{\code{data}}{every monitored observation so far, baseline
#' excluded (the first coordinate, for a multivariate monitor),
#' which is what \code{autoplot()} draws.}
#' \item{\code{method}, \code{alpha}, \code{n_baseline},
#' \code{reset}, \code{relearn}}{the settings this call fixed,
#' kept so a later \code{cpt_update()} cannot silently disagree
#' with them.}
#' \item{\code{state}}{the detector's internal statistics. Engine
#' internals: their shape differs by \code{method} and is not part
#' of the interface.}
#' }
#' @references
#' \insertRef{shin2023edetectors}{ggchangepoint}
#'
#' \insertRef{page1954continuous}{ggchangepoint}
#'
#' \insertRef{roberts1959ewma}{ggchangepoint}
#'
#' \insertRef{brook1972arl}{ggchangepoint}
#'
#' \insertRef{lucas1990ewma}{ggchangepoint}
#' @seealso \code{\link{cpt_update}()}, \code{\link{alarms}()},
#' \code{\link{cpt_replay}()}, \code{\link{cpt_delay}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_monitor("edetector", baseline = rnorm(100))
#' mon <- cpt_update(mon, rnorm(50)) # still in control
#' mon <- cpt_update(mon, rnorm(50, 3)) # a change arrives
#' alarms(mon)
cpt_monitor <- function(method = c("edetector", "cpm", "ocd", "cusum",
"ewma"),
baseline = NULL, alpha = 0.01, arl0 = 500,
cpm_type = "Mann-Whitney", patience = 5000,
deltas = c(0.5, 1, 2), reset = TRUE, relearn = 20,
thresh = "MC", mc_reps = 100, ..., k = 0.5,
lambda = 0.2, limit = NULL) {
method <- cpt_match_arg(method)
# The three detectors are calibrated in different currencies and each
# ignores the others': setting `arl0` on an e-detector, or `alpha` on
# `cpm`, changes nothing at all. The documentation and the vignette both
# say so, which does not help the user who set one and is now reading an
# unchanged answer. Only an argument the caller supplied is reported --
# the defaults supply all of them.
governs <- list(alpha = "edetector", deltas = "edetector",
arl0 = c("cpm", "cusum", "ewma"), cpm_type = "cpm",
patience = "ocd", thresh = "ocd", mc_reps = "ocd",
k = "cusum", lambda = "ewma", limit = c("cusum", "ewma"))
supplied <- c(alpha = !missing(alpha), deltas = !missing(deltas),
arl0 = !missing(arl0), cpm_type = !missing(cpm_type),
patience = !missing(patience), thresh = !missing(thresh),
mc_reps = !missing(mc_reps), k = !missing(k),
lambda = !missing(lambda), limit = !is.null(limit))
tunes <- vapply(governs, function(g) method %in% g, logical(1))
ignored <- names(governs)[supplied[names(governs)] & !tunes]
if (length(ignored) > 0) {
cpt_warn("`", paste(ignored, collapse = "`, `"), "` ",
if (length(ignored) > 1) "do" else "does", " not affect `method ",
"= \"", method, "\"`, which is tuned by `",
paste(names(governs)[tunes], collapse = "`, `"),
"`. See ?cpt_monitor.", class = "argument_ignored")
}
# `...` reaches cpm::processStream() or ocd::ChangepointDetector(), and
# this function renames three of their arguments and pins two more.
# Passing an engine's own name gave R's raw "formal argument \"cpmType\"
# matched by multiple actual arguments" -- `cpmType`, `ARL0`, `dim`,
# `MC_reps` and `beta`, measured across the two engines. `thresh`,
# `patience` and `startup` cannot reach `...` at all, because they are
# formals here and R matches them there first. `cpt_replay()` forwards
# `...` to this function, so it is covered by the same guard.
reject_renamed_args(list(...), method)
if (identical(method, "ocd")) {
reject_managed_args(list(...), "ocd", c(
beta = paste("the monitor calibrates against a unit-scale change;",
"use `ocd_wrapper()` directly to set it")))
}
validate_scalar(alpha, "alpha", min = 0, max = 1,
min_open = TRUE, max_open = TRUE)
validate_scalar(relearn, "relearn", min = 0)
validate_flag(reset, "reset")
# Three of ten arguments were checked. The rest reached their engine (or,
# for `deltas`, this function's own arithmetic) unvalidated: `deltas = 0`
# makes every `inc` exactly 1, so `R` grows LINEARLY whatever the data
# does and the monitor alarms at t = 1/alpha on pure noise.
if (method %in% c("cusum", "ewma") && length(list(...)) > 0) {
extra <- names(list(...))
cpt_warn("`method = \"", method, "\"` is native to this package and ",
"takes no engine arguments, so ",
if (length(extra) && all(nzchar(extra))) {
paste0("`", paste(extra, collapse = "`, `"), "` ")
} else {
paste0(length(list(...)), " unnamed argument(s) ")
}, "would be ignored. It is tuned by `arl0`, `",
if (method == "cusum") "k" else "lambda", "` and `limit`.",
class = "argument_ignored")
}
if (identical(method, "edetector")) {
if (!is.numeric(deltas) || length(deltas) == 0L ||
any(!is.finite(deltas)) || all(deltas == 0)) {
cpt_abort("`deltas` must be one or more finite non-zero numbers: the ",
"shift ", "sizes, in baseline standard deviations, the ",
"e-detector mixes ", "over. A zero shift contributes a ",
"likelihood ratio of exactly 1 ", "at every observation, so ",
"the statistic grows on nothing.", class = "bad_argument")
}
if (any(deltas == 0)) {
cpt_warn("`deltas` contains 0, which contributes a likelihood ratio of ",
"exactly 1 at every observation and only dilutes the mixture. ",
"Dropping it.", class = "dropped_input")
deltas <- deltas[deltas != 0]
}
# `...` reaches the engine in the cpm and ocd branches and has nowhere
# to go in this one, so it used to vanish -- the same silent-discard
# ggcptplot_internal() warns about.
extra <- names(list(...))
if (length(list(...)) > 0) {
cpt_warn("`method = \"edetector\"` is native to this package and takes ",
"no engine arguments, so ",
if (length(extra) && all(nzchar(extra))) {
paste0("`", paste(extra, collapse = "`, `"), "` ")
} else {
paste0(length(list(...)), " unnamed argument(s) ")
},
"would be ignored. It is tuned by `alpha` and `deltas`.",
class = "argument_ignored")
}
}
if (identical(method, "ocd")) {
validate_scalar(patience, "patience", min = 1)
validate_scalar(mc_reps, "mc_reps", min = 1)
}
# Normalise the series ONCE, before the switch. Each branch used to coerce
# `baseline` its own way -- as.numeric() for edetector and cpm, a matrix
# for ocd -- so a factor reached two of the three detectors as its level
# codes. Guarding here means every method enters through the same door,
# which is the only way this stays true when a fourth branch is added.
if (!is.null(baseline)) {
if (is.matrix(baseline) || is.data.frame(baseline)) {
# `edetector` and `cpm` monitor ONE series, and length() on a
# data.frame counts its columns -- so a 60-row, 2-column baseline was
# refused for having "at least 5 pre-change observations", naming a
# count of 2 for 60 observations. Refuse the shape instead, since a
# multi-column baseline for a univariate detector is outside the
# documented contract either way.
if (!identical(method, "ocd")) {
reject_multicolumn(baseline, "baseline",
paste0("`method = \"", method, "\"` monitors one ",
"series; use `method = \"ocd\"` for a ",
"multivariate stream."))
baseline <- as_mv_matrix(baseline, arg = "baseline")[, 1]
}
} else {
baseline <- coerce_series_values(baseline, arg = "baseline")
}
}
state <- switch(method,
edetector = {
if (is.null(baseline) || NROW(baseline) < 5) {
cpt_abort("`edetector` estimates the in-control mean and scale from ",
"`baseline`, so it needs at least 5 pre-change observations.",
class = "short_series")
}
b <- as.numeric(baseline)
# Check finiteness first: sd() of a vector holding an NA is itself NA,
# so the variability test below would blame a flat baseline for what is
# really a missing value.
if (anyNA(b) || any(!is.finite(b))) {
stop_nonfinite(b, "baseline", stream_na_hint(b))
}
sd0 <- stats::sd(b)
if (!is.finite(sd0) || sd0 <= 0) {
cpt_abort("`baseline` has zero variability, so the e-detector's scale ",
"is ", "undefined.", class = "input_error")
}
deltas <- sort(unique(c(-abs(deltas), abs(deltas)))) * sd0
list(mu0 = mean(b), sd0 = sd0, deltas = deltas,
R = rep(0, length(deltas)),
threshold = 1 / alpha)
},
cpm = {
need_pkg("cpm")
# The same three guards cpm_wrapper() applies, from the same helpers:
# this branch used to call makeChangePointModel() directly, so a
# withheld `cpm_type`, a missing FET `lambda` or an off-grid `arl0`
# each reached the user as an error from inside cpm or base R that
# named no argument at all (`no applicable method for '@' applied to
# an object of class "NULL"`, because cpm *prints* its complaint and
# returns NULL). See cpm_check_type() for what each one looks like.
validate_scalar(arl0, "arl0")
cpm_type <- cpm_check_type(cpm_type, list(...))
engine_dots <- cpm_engine_dots(cpm_type, list(...))
cpm_out <- utils::capture.output(
m <- do.call(cpm::makeChangePointModel,
c(list(cpmType = cpm_type, ARL0 = arl0), engine_dots))
)
cpm_check_printed_error(cpm_out, arl0, engine_dots)
if (length(cpm_out)) cat(cpm_out, sep = "\n")
# A change cpm detects INSIDE the baseline left the model tripped:
# changeDetected() was already TRUE when the loop ended, so the first
# monitored observation inherited the alarm. Measured before this
# was fixed: an alarm at stream index 1 in 2 of 6 clean iid runs and
# 6 of 6 AR(1) runs, and never without a baseline. The model is now
# restarted after each baseline detection, as cpt_update() does after
# an alarm, and the caller is told the baseline was not in control.
n_tripped <- 0L
if (!is.null(baseline)) {
for (v in as.numeric(baseline)) {
m <- cpm::processObservation(m, v)
if (isTRUE(cpm::changeDetected(m))) {
n_tripped <- n_tripped + 1L
m <- cpm::cpmReset(m)
}
}
}
if (n_tripped > 0L) {
cpt_warn("cpm detected ", n_tripped, " change(s) inside `baseline`, ",
"so the baseline is not in control. The monitor was ",
"restarted after each, and only the stretch after the last ",
"one informs it; check the baseline, or build the monitor ",
"without one.", class = "assumption",
data = list(n_changes = n_tripped))
}
list(model = m, cpm_type = cpm_type, arl0 = arl0)
},
cusum = ,
ewma = chart_state(method, baseline, arl0 = arl0, k = k,
lambda = lambda, limit = limit),
ocd = {
need_pkg("ocd")
if (is.null(baseline)) {
cpt_abort("`ocd` needs `baseline` to estimate the pre-change mean and ",
"standard deviation.", class = "bad_argument")
}
B <- if (is.matrix(baseline) || is.data.frame(baseline)) {
as_mv_matrix(baseline, arg = "baseline")
} else {
matrix(baseline, ncol = 1)
}
p <- ncol(B)
if (p < 2) {
cpt_abort("Method `ocd` is high-dimensional and needs at least two ",
"coordinates, but `baseline` has ", p, ". Use method = ",
"\"edetector\" or \"cpm\" for one series.",
class = "wrong_dimension")
}
det <- ocd::ChangepointDetector(dim = p, method = "ocd",
thresh = thresh, beta = 1,
patience = patience,
MC_reps = mc_reps, ...)
det <- ocd::setBaselineMean(det, colMeans(B))
det <- ocd::setBaselineSD(det, apply(B, 2, stats::sd))
det <- ocd::setStatus(det, "monitoring")
list(detector = det, p = p, patience = patience,
mean = colMeans(B),
sd = { sdv <- apply(B, 2, stats::sd)
sdv[!is.finite(sdv) | sdv == 0] <- 1
sdv })
}
)
if (!is.null(baseline)) check_baseline_dependence(baseline, method)
structure(
list(method = method, state = state, reset = reset,
relearn = as.integer(relearn),
n_baseline = if (is.null(baseline)) 0L else NROW(baseline),
t = 0L,
relearn_left = 0L,
relearn_buffer = NULL,
data = numeric(0),
trace = tibble::tibble(time = integer(), statistic = numeric(),
lower = numeric()),
alarms = tibble::tibble(time = integer(), statistic = numeric(),
threshold = numeric()),
alpha = alpha),
class = "ggcpt_monitor"
)
}
# Internal: the state of a CUSUM or EWMA chart, standardised by the
# baseline and calibrated to arl0 unless a limit is given.
#' @noRd
chart_state <- function(method, baseline, arl0, k, lambda, limit) {
if (is.null(baseline) || NROW(baseline) < 5) {
cpt_abort("`", method, "` standardises the stream by the in-control ",
"mean and standard deviation of `baseline`, so it needs at ",
"least 5 pre-change observations.", class = "short_series")
}
b <- as.numeric(baseline)
if (anyNA(b) || any(!is.finite(b))) {
stop_nonfinite(b, "baseline", stream_na_hint(b))
}
sd0 <- stats::sd(b)
if (!is.finite(sd0) || sd0 <= 0) {
cpt_abort("`baseline` has zero variability, so the chart's scale is ",
"undefined.", class = "input_error")
}
check_chart_arl0(arl0)
if (!is.null(limit)) validate_scalar(limit, "limit", min = 0,
min_open = TRUE)
if (method == "cusum") {
validate_scalar(k, "k", min = 0)
h <- limit %||% cusum_h_for(arl0, k)
return(list(mu0 = mean(b), sd0 = sd0, k = k, h = h, sp = 0, sn = 0,
threshold = h, arl0 = arl0,
calibrated = is.null(limit)))
}
validate_scalar(lambda, "lambda", min = 0, max = 1, min_open = TRUE)
L <- limit %||% ewma_L_for(arl0, lambda)
list(mu0 = mean(b), sd0 = sd0, lambda = lambda, L = L, z = 0,
threshold = L * sqrt(lambda / (2 - lambda)), arl0 = arl0,
calibrated = is.null(limit))
}
#' Feed observations to a monitor
#'
#' Pushes new data through a \code{\link{cpt_monitor}()}, updating its
#' internal state and appending to its alarm log. The monitor is returned, so
#' the idiom is \code{mon <- cpt_update(mon, new_obs)}.
#'
#' @param monitor A \code{ggcpt_monitor} object.
#' @param new_obs New observations: a numeric vector, or a matrix with rows
#' as time points for a multivariate monitor.
#' @return The updated \code{ggcpt_monitor}: the same structure as the
#' input (see \code{\link{cpt_monitor}()}), with \code{t} advanced by
#' the number of new observations (rows, for a matrix) and any new rows
#' appended to \code{alarms}.
#' Re-assign it: the monitor is a value, not a mutable object, so
#' \code{cpt_update(mon, y)} without assignment discards the update.
#' @seealso \code{\link{cpt_monitor}()}, \code{\link{alarms}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_monitor("edetector", baseline = rnorm(100))
#' mon <- cpt_update(mon, c(rnorm(50), rnorm(50, 3)))
#' alarms(mon)
cpt_update <- function(monitor, new_obs) {
if (!inherits(monitor, "ggcpt_monitor")) {
cpt_abort("`monitor` must be a ggcpt_monitor from cpt_monitor().",
class = "bad_argument")
}
X <- if (is.matrix(new_obs) || is.data.frame(new_obs)) {
as_mv_matrix(new_obs, arg = "new_obs")
} else {
# a streaming caller passes observations one at a time, so this is the
# most likely place for a factor level code to enter unnoticed
matrix(coerce_series_values(new_obs, arg = "new_obs"), ncol = 1)
}
if (nrow(X) == 0) return(monitor)
if (anyNA(X) || any(!is.finite(X))) {
stop_nonfinite(X, "new_obs", stream_na_hint(X))
}
# A monitor is stateful and dimensioned by its baseline. Feeding it a
# different width is a mistake, not a coercion: `ocd` would consume the
# wrong coordinates and the univariate detectors would silently read
# column 1 and ignore the rest.
expected <- monitor_width(monitor)
if (!is.na(expected) && ncol(X) != expected) {
cpt_abort("This monitor was built on ", expected, " coordinate(s) but ",
"`new_obs` has ", ncol(X), ". A monitor cannot change width ",
"once it is running.", class = "wrong_dimension")
}
st <- monitor$state
new_alarms <- list()
# Grown once at the end, not re-allocated per observation. cpt_replay()
# hands the whole series to a single cpt_update() call, so
# `monitor$data <- c(monitor$data, xt[1])` inside the loop made replaying
# 10,000 observations do on the order of 50 million element copies.
seen <- numeric(nrow(X))
trace_stat <- rep(NA_real_, nrow(X))
trace_low <- rep(NA_real_, nrow(X))
for (i in seq_len(nrow(X))) {
monitor$t <- monitor$t + 1L
xt <- X[i, ]
# Coordinate one only, for every method: `$data` is the series
# autoplot() draws, and a multivariate `ocd` monitor is drawn as its
# first coordinate (see @return).
seen[i] <- xt[1]
# Re-learning window after an alarm: collect the new regime and stay
# silent, then adopt it as the in-control baseline. Without this a
# persistent change makes the detector alarm on every subsequent
# observation (see the `relearn` argument of cpt_monitor()).
if (monitor$relearn_left > 0L) {
monitor$relearn_buffer <- rbind(monitor$relearn_buffer, xt)
monitor$relearn_left <- monitor$relearn_left - 1L
if (monitor$relearn_left == 0L) {
st <- monitor_relearn(monitor$method, st, monitor$relearn_buffer)
monitor$relearn_buffer <- NULL
}
next
}
hit <- FALSE
stat <- NA_real_
thr <- NA_real_
if (monitor$method == "edetector") {
inc <- exp(st$deltas * (xt[1] - st$mu0) / st$sd0^2 -
st$deltas^2 / (2 * st$sd0^2))
st$R <- (1 + st$R) * inc
# An e-detector's mixture grows multiplicatively, so under `reset =
# FALSE` with `relearn = 0` -- the "keep accumulating" configuration
# the arguments offer -- it passes .Machine$double.xmax a few hundred
# observations after a real change and the next product is `Inf`. The
# alarm rule below reads a non-finite statistic as "no alarm", so the
# monitor went permanently silent exactly when the evidence was
# overwhelming: measured 82 alarms and then 1418 observations of
# silence on a shifted stream. Saturating instead keeps the statistic
# ordered against the threshold, which is all the rule needs, and --
# unlike `Inf`, which is absorbing -- it still decays when the stream
# returns to its baseline (measured: 8.99e307 down to 1.8e149 over
# 2000 in-control observations, where `Inf` would have stayed `Inf`).
st$R <- pmin(st$R, .Machine$double.xmax / 2)
# AVERAGE, not max: a convex combination of e-detectors is an
# e-detector and keeps E[M_t] = t, which is what the average-run-length
# bound rests on. A maximum over K shifts crosses the threshold roughly
# K times as often under the null -- measured here at ~1.8x the bound
# with six shifts before this was corrected.
stat <- mean(st$R)
thr <- st$threshold
hit <- is.finite(stat) && stat >= thr
if (hit && monitor$reset) st$R <- rep(0, length(st$deltas))
} else if (monitor$method == "cusum") {
z <- (xt[1] - st$mu0) / st$sd0
st$sp <- max(0, st$sp + z - st$k)
st$sn <- max(0, st$sn - z - st$k)
trace_low[i] <- -st$sn
stat <- max(st$sp, st$sn)
thr <- st$h
hit <- stat >= thr
trace_stat[i] <- st$sp
if (hit && monitor$reset) {
st$sp <- 0
st$sn <- 0
}
} else if (monitor$method == "ewma") {
z <- (xt[1] - st$mu0) / st$sd0
st$z <- st$lambda * z + (1 - st$lambda) * st$z
stat <- st$z
thr <- st$threshold
hit <- abs(stat) >= thr
trace_stat[i] <- st$z
if (hit && monitor$reset) st$z <- 0
} else if (monitor$method == "cpm") {
st$model <- cpm::processObservation(st$model, xt[1])
stat <- suppressWarnings(as.numeric(
utils::tail(cpm::getStatistics(st$model), 1)
))
if (length(stat) == 0) stat <- NA_real_
hit <- isTRUE(cpm::changeDetected(st$model))
if (hit && monitor$reset) {
st$model <- cpm::cpmReset(st$model)
}
} else {
# ocd monitors a standardised stream, so the baseline enters here
# rather than through the detector object.
z <- (as.numeric(xt) - st$mean) / st$sd
utils::capture.output(st$detector <- ocd::getData(st$detector, z))
# `normalisedStatistics()` already divides each statistic by its own
# threshold, so the comparison point is 1. Reporting the raw
# thresholds alongside the normalised statistic would print an alarm
# at 1.10 against a threshold of 16.4, which reads as a bug.
stat <- suppressWarnings(max(as.numeric(
unlist(ocd::normalisedStatistics(st$detector))
), na.rm = TRUE))
thr <- 1
hit <- !identical(ocd::status(st$detector), "monitoring")
if (hit && monitor$reset) {
st$detector <- ocd::reset(st$detector)
st$detector <- ocd::setStatus(st$detector, "monitoring")
}
}
if (!monitor$method %in% c("cusum", "ewma")) trace_stat[i] <- stat
if (isTRUE(hit)) {
new_alarms[[length(new_alarms) + 1L]] <-
tibble::tibble(time = monitor$t, statistic = stat, threshold = thr)
if (monitor$relearn > 0L) {
monitor$relearn_left <- monitor$relearn
monitor$relearn_buffer <- NULL
}
}
}
monitor$state <- st
monitor$data <- c(monitor$data, seen)
# The statistic at every observation (NA while re-learning), so a chart
# can be drawn; a monitor built before 0.7.0 has no trace to extend.
steps <- tibble::tibble(time = monitor$t - rev(seq_len(nrow(X))) + 1L,
statistic = trace_stat, lower = trace_low)
monitor$trace <- rbind(monitor$trace %||% steps[0, ], steps)
if (length(new_alarms) > 0) {
monitor$alarms <- rbind(monitor$alarms, do.call(rbind, new_alarms))
}
monitor
}
# Internal: how many coordinates this monitor was built for. `edetector`
# and `cpm` read a single series; `ocd` carries the dimension it was
# constructed with.
#' @noRd
monitor_width <- function(monitor) {
switch(monitor$method,
ocd = monitor$state$p,
edetector = 1L,
cpm = 1L,
cusum = 1L,
ewma = 1L,
NA_integer_
)
}
# Internal: adopt a post-alarm window as the new in-control baseline. Only
# the detectors that hold an explicit baseline need this; cpm re-derives its
# own reference from the stream after a reset.
#' @noRd
monitor_relearn <- function(method, st, buffer) {
if (is.null(buffer) || nrow(buffer) == 0) return(st)
if (method == "edetector") {
v <- as.numeric(buffer[, 1])
sd0 <- stats::sd(v)
if (is.finite(sd0) && sd0 > 0) {
st$deltas <- st$deltas / st$sd0 * sd0
st$sd0 <- sd0
}
st$mu0 <- mean(v)
st$R <- rep(0, length(st$deltas))
} else if (method %in% c("cusum", "ewma")) {
v <- as.numeric(buffer[, 1])
sd0 <- stats::sd(v)
if (is.finite(sd0) && sd0 > 0) st$sd0 <- sd0
st$mu0 <- mean(v)
if (method == "cusum") {
st$sp <- 0
st$sn <- 0
} else {
st$z <- 0
}
} else if (method == "ocd") {
st$mean <- colMeans(buffer)
sdv <- apply(buffer, 2, stats::sd)
sdv[!is.finite(sdv) | sdv == 0] <- 1
st$sd <- sdv
}
st
}
#' The alarm log of a monitor
#'
#' Every observation at which a sequential monitor crossed its threshold, in
#' the order they fired. A monitor that never fired returns a zero-row tibble
#' rather than \code{NULL}, so the result is always safe to \code{rbind()} or
#' plot.
#'
#' @param x A \code{ggcpt_monitor} object.
#' @param ... Ignored.
#' @return A tibble with one row per alarm: \code{time} (the observation at
#' which it fired, counted from the first one fed to the monitor),
#' \code{statistic} and \code{threshold}.
#' @seealso \code{\link{cpt_monitor}()}, \code{\link{cpt_delay}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_monitor("edetector", baseline = rnorm(100))
#' mon <- cpt_update(mon, c(rnorm(50), rnorm(50, 3)))
#' alarms(mon)
alarms <- function(x, ...) {
UseMethod("alarms")
}
#' @rdname alarms
#' @export
alarms.ggcpt_monitor <- function(x, ...) x$alarms
#' @rdname cpt_monitor
#' @export
tidy.ggcpt_monitor <- function(x, ...) {
alarms(x)
}
#' @rdname cpt_monitor
#' @param x A \code{ggcpt_monitor} object (for \code{print()}).
#' @export
print.ggcpt_monitor <- function(x, ...) {
cat("ggcpt_monitor (", x$method,
if (x$method %in% c("edetector", "cusum", "ewma")) {
", native implementation"
} else "",
")\n", sep = "")
st <- x$state
if (identical(x$method, "cusum")) {
cat(" Chart: k = ", format(st$k), ", h = ", format(signif(st$h, 4)),
if (isTRUE(st$calibrated)) paste0(" (in-control ARL ",
format(st$arl0), ")"),
"\n", sep = "")
} else if (identical(x$method, "ewma")) {
cat(" Chart: lambda = ", format(st$lambda), ", L = ",
format(signif(st$L, 4)),
if (isTRUE(st$calibrated)) paste0(" (in-control ARL ",
format(st$arl0), ")"),
"\n", sep = "")
}
cat(" Baseline observations: ", x$n_baseline, "\n", sep = "")
cat(" Monitored observations: ", x$t, "\n", sep = "")
cat(" Alarms: ", nrow(x$alarms), "\n", sep = "")
cat(" Restart after alarm: ", x$reset,
if (x$relearn > 0L) paste0(" (re-learning ", x$relearn,
" observations)") else "",
"\n", sep = "")
if (nrow(x$alarms) > 0) {
cat("\n")
print(x$alarms, n = 10)
}
invisible(x)
}
#' @rdname cpt_monitor
#' @param object A \code{ggcpt_monitor} object (for \code{autoplot()}).
#' @param plot_type \code{"timeline"} (the monitored series with the alarms
#' marked), \code{"control_chart"} (the detection statistic at every
#' observation against its limits, with the alarms flagged: for a CUSUM
#' the upper and lower sums against \eqn{\pm h}, for an EWMA the
#' smoothed statistic between its limits), \code{"statistic"} (the
#' statistic at each alarm only, against its threshold) or
#' \code{"runlength"} (the gaps between alarms, which estimate the run
#' length). For a multivariate monitor the timeline draws
#' the first coordinate: the alarms are shared, so the rules are right
#' whichever coordinate is shown, but the line is one of several.
#' @export
autoplot.ggcpt_monitor <- function(object,
plot_type = c("timeline", "statistic",
"runlength",
"control_chart"),
...) {
plot_type <- cpt_match_arg(plot_type)
al <- object$alarms
d <- tibble::tibble(time = seq_along(object$data), value = object$data)
if (plot_type == "control_chart") return(control_chart(object))
if (plot_type == "timeline") {
p <- ggplot2::ggplot(d, ggplot2::aes(time, value)) +
ggplot2::geom_line(colour = "grey40") +
ggplot2::labs(x = "Observation", y = "Value",
title = paste0("Monitoring timeline (", object$method,
")"),
subtitle = paste0(nrow(al), " alarm(s) in ", object$t,
" observations"))
if (nrow(al) > 0) {
p <- p + ggplot2::geom_vline(xintercept = al$time, colour = "#D55E00",
linetype = "dashed", linewidth = 0.5)
}
return(p)
}
if (plot_type == "statistic") {
if (nrow(al) == 0) {
# cpt_replay() used to be offered here as the way to a "full
# statistic trace". It returns the same object, which records the
# statistic at alarms only, so the advice led nowhere.
cpt_abort("No alarm has fired, so there is no statistic to draw: a ",
"monitor records its statistic only when an alarm fires. ",
"`plot_type = \"timeline\"` shows the monitored series.",
class = "capability_absent")
}
return(
ggplot2::ggplot(al, ggplot2::aes(time, statistic)) +
ggplot2::geom_point(size = 2, colour = "#0072B2") +
ggplot2::geom_hline(ggplot2::aes(yintercept = threshold),
linetype = "dotted", colour = "grey30",
na.rm = TRUE) +
ggplot2::labs(x = "Observation", y = "Detection statistic",
title = "Statistic at each alarm")
)
}
# runlength
if (nrow(al) < 2) {
cpt_abort("At least two alarms are needed to show run lengths; this ",
"monitor ", "has ", nrow(al), ".", class = "capability_absent")
}
rl <- tibble::tibble(run_length = diff(c(0L, al$time)))
ggplot2::ggplot(rl, ggplot2::aes(run_length)) +
ggplot2::geom_histogram(bins = max(5, min(30, nrow(rl))),
fill = "#0072B2", colour = "white") +
ggplot2::geom_vline(xintercept = mean(rl$run_length), colour = "#D55E00",
linetype = "dashed") +
ggplot2::labs(x = "Run length (observations between alarms)",
y = "Count",
title = "Run-length distribution",
subtitle = paste0("Mean run length ",
format(mean(rl$run_length), digits = 4)))
}
#' Replay a series through a sequential detector
#'
#' Runs a whole series through a monitor in one call, so the alarm timeline
#' of an online method can be studied retrospectively without hand-rolling
#' the update loop. This is the honest version of what 0.4.0's \code{cpm},
#' \code{ocd} and \code{bocpd} wrappers do implicitly: it reports
#' \emph{alarms}, which arrive after the change, rather than pretending they
#' are estimated changepoint locations.
#'
#' @param x A numeric vector, or a matrix with rows as time points.
#' @param method Sequential detector; see \code{\link{cpt_monitor}()}.
#' @param baseline Number of leading observations used as pre-change training
#' data (an integer), or an explicit baseline vector. Defaults to
#' \code{min(100, floor(n / 4))}.
#' @param ... Additional arguments passed to \code{\link{cpt_monitor}()}.
#' @return A \code{ggcpt_monitor} that has already consumed the series:
#' the same object \code{\link{cpt_monitor}()} returns, described
#' there, with \code{alarms} filled in for whatever fired.
#' @seealso \code{\link{cpt_delay}()}, \code{\link{cpt_monitor}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_replay(c(rnorm(200), rnorm(200, 3)), method = "edetector")
#' alarms(mon)
cpt_replay <- function(x, method = c("edetector", "cpm", "ocd", "cusum",
"ewma"),
baseline = NULL, ...) {
method <- cpt_match_arg(method)
# coerce_series_values()/as_mv_matrix(), not a bare as.numeric(): a factor
# coerces to its LEVEL CODES, which for labels like "10", "2", "30" is the
# alphabetical order 3, 1, 2 rather than the numbers -- so a replay ran on
# a series the user never supplied. cpt_detect() has refused that since
# 0.4.0, and as_mv_matrix() already refuses a factor *column*; only the
# vector path was open.
X <- if (is.matrix(x) || is.data.frame(x)) as_mv_matrix(x) else {
matrix(coerce_series_values(x), ncol = 1)
}
n <- nrow(X)
if (n < 10) {
cpt_abort("A replay needs at least 10 observations; got ", n, ".",
class = "short_series")
}
# Without this, a non-finite value was caught downstream by whichever of
# cpt_monitor() or cpt_update() happened to receive the slice containing
# it -- so `cpt_replay(x)` reported a problem with `baseline` for an NA at
# position 20 and with `new_obs` for one at 95, named an argument the
# caller never passed, and counted "1 of 45 values" against the baseline
# slice rather than the series. Checked here, the message describes `x`.
if (anyNA(X) || any(!is.finite(X))) {
stop_nonfinite(X, hint = stream_na_hint(X))
}
if (is.null(baseline)) baseline <- min(100L, floor(n / 4))
# `baseline` is documented two ways -- a count of leading observations, or
# an explicit series -- so it needs the same type guard as `x` above before
# either reading is applied. The explicit-series branch used a bare
# as.numeric(), and it hands cpt_monitor() an already-numeric vector, so
# cpt_monitor()'s own guard could never see the factor.
if (!is.matrix(baseline) && !is.data.frame(baseline)) {
baseline <- coerce_series_values(baseline, arg = "baseline")
}
if (length(baseline) == 1L) {
# A lone number is the count. It used to fall through to the
# explicit-series branch whenever it was fractional or out of range,
# which turned it into a ONE-point baseline: cpt_replay(x, baseline =
# 500) on a 180-point series then failed with "needs at least 5
# pre-change observations" and never mentioned the 500.
if (is.na(baseline) || baseline %% 1 != 0 || baseline < 1 ||
baseline >= n) {
cpt_abort("`baseline` as a single number is the count of leading ",
"observations to train on, so it must be a whole number in ",
"1..", n - 1L, "; the series has ", n, " observations and ",
"`baseline` is ", format(baseline), ". Pass an explicit ",
"baseline series if you meant a value rather than a count.",
class = "bad_argument")
}
n_base <- as.integer(baseline)
base_data <- X[seq_len(n_base), , drop = FALSE]
rest <- X[seq.int(n_base + 1L, n), , drop = FALSE]
} else {
base_data <- if (is.matrix(baseline) || is.data.frame(baseline)) {
as_mv_matrix(baseline, arg = "baseline")
} else {
matrix(baseline, ncol = ncol(X))
}
n_base <- 0L
rest <- X
}
mon <- cpt_monitor(method,
baseline = if (ncol(base_data) == 1) {
as.numeric(base_data)
} else {
base_data
},
...)
mon <- cpt_update(mon, rest)
mon$offset <- n_base
mon
}
#' Detection delay and false-alarm rate
#'
#' Scores an online result the way the sequential literature does: how long
#' after each true change did the first alarm arrive, and how many alarms
#' were raised with no change behind them. \code{\link{cpt_metrics}()} is the
#' wrong tool for an online detector: it asks whether the \emph{location}
#' was recovered, which a sequential procedure never claims. It will not
#' stop you: \code{cpt_metrics()} takes bare integer vectors and never
#' sees which detector produced them, so it cannot know. (An earlier
#' version of this sentence said it warns. It does not, and given that
#' signature it could not.)
#'
#' @param object A \code{ggcpt_monitor}, or a tibble of alarms with a
#' \code{time} column.
#' @param truth Integer vector of true changepoint positions, on the same
#' clock as the alarms. When \code{object} is a monitor built by
#' \code{\link{cpt_replay}()} the baseline offset is applied automatically.
#' @param max_delay Alarms further than this after a change are treated as
#' false alarms rather than late detections. Defaults to \code{Inf}.
#' @return A \code{ggcpt_delay} object: a list with \code{per_change} (one
#' row per true change: \code{truth}, \code{alarm}, \code{delay},
#' \code{detected}), \code{false_alarms}, and the summary statistics
#' \code{mean_delay}, \code{median_delay}, \code{n_false_alarms} and
#' \code{arl}: monitored observations per false alarm. The baseline a
#' \code{\link{cpt_replay}()} monitor trained on is not counted, because
#' no alarm can fire there; \code{n_obs} is the length of the stream in
#' series positions, baseline included.
#' @seealso \code{\link{cpt_monitor}()}, \code{\link{cpt_replay}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_replay(c(rnorm(200), rnorm(200, 3)), method = "edetector")
#' cpt_delay(mon, truth = 200)
cpt_delay <- function(object, truth, max_delay = Inf) {
# Compared straight against alarm times, so `NA` or a negative value
# silently counted every alarm as false, a vector recycled, and a string
# compared as text.
if (!is.numeric(max_delay) || length(max_delay) != 1L || is.na(max_delay) ||
max_delay < 0) {
cpt_abort("`max_delay` must be a single non-negative number (`Inf` for no ",
"limit); got ", paste(format(max_delay), collapse = ", "), ".",
class = "bad_argument")
}
if (missing(truth)) {
cpt_abort("`truth` is required: detection delay is measured from the true ",
"changepoint(s), so there is nothing to measure without them. Pass ",
"the location(s) as an integer vector.", class = "bad_argument")
}
al <- if (inherits(object, "ggcpt_monitor")) {
a <- object$alarms
off <- object$offset %||% 0L
a$time <- a$time + off
a
} else if (is.data.frame(object) && "time" %in% names(object)) {
tibble::as_tibble(object)
} else {
cpt_abort("`object` must be a ggcpt_monitor or a tibble with a `time` ",
"column.", class = "bad_argument")
}
truth <- as_cp_locations(truth, "truth", sort = TRUE)
if (length(truth) == 0 || anyNA(truth) || any(truth < 1)) {
cpt_abort("`truth` must be one or more positive changepoint locations.",
class = "bad_argument")
}
n_obs <- if (inherits(object, "ggcpt_monitor")) {
object$t + (object$offset %||% 0L)
} else {
max(c(al$time, truth), na.rm = TRUE)
}
n_monitored <- if (inherits(object, "ggcpt_monitor")) object$t else n_obs
# A truth past the end of the stream can never be detected, so scoring it
# would report a miss that says more about the argument than the detector.
if (inherits(object, "ggcpt_monitor") && any(truth > n_obs)) {
cpt_abort("`truth` runs past the end of the stream: ",
paste(utils::head(truth[truth > n_obs], 5), collapse = ", "),
" vs ", n_obs, " observations seen.", class = "bad_argument")
}
used <- rep(FALSE, nrow(al))
rows <- lapply(seq_along(truth), function(i) {
k <- truth[i]
next_change <- if (i < length(truth)) truth[i + 1] else Inf
cand <- which(!used & al$time >= k & al$time < next_change &
al$time - k <= max_delay)
if (length(cand) == 0) {
return(tibble::tibble(truth = k, alarm = NA_integer_,
delay = NA_real_, detected = FALSE))
}
j <- cand[1]
used[j] <<- TRUE
tibble::tibble(truth = k, alarm = as.integer(al$time[j]),
delay = as.numeric(al$time[j] - k), detected = TRUE)
})
per_change <- do.call(rbind, rows)
false_alarms <- al[!used, , drop = FALSE]
structure(
list(
per_change = per_change,
false_alarms = false_alarms,
mean_delay = mean(per_change$delay, na.rm = TRUE),
median_delay = stats::median(per_change$delay, na.rm = TRUE),
n_false_alarms = nrow(false_alarms),
n_detected = sum(per_change$detected),
n_changes = length(truth),
# Monitored observations, not the whole stream: a replay's training
# baseline cannot raise an alarm, so counting it inflated the ARL by
# n_baseline / n_false_alarms (300 monitored observations with three
# false alarms read as 133, not 100), in the direction that flatters
# the calibration the ARL exists to check.
arl = if (nrow(false_alarms) > 0) {
n_monitored / nrow(false_alarms)
} else {
Inf
},
n_obs = n_obs
),
class = "ggcpt_delay"
)
}
#' @rdname cpt_delay
#' @export
tidy.ggcpt_delay <- function(x, ...) {
x$per_change
}
#' @rdname cpt_delay
#' @export
glance.ggcpt_delay <- function(x, ...) {
tibble::tibble(n_changes = x$n_changes, n_detected = x$n_detected,
mean_delay = x$mean_delay, median_delay = x$median_delay,
n_false_alarms = x$n_false_alarms, arl = x$arl,
n_obs = x$n_obs)
}
#' @rdname cpt_delay
#' @param x A \code{ggcpt_delay} object.
#' @param ... Ignored.
#' @export
print.ggcpt_delay <- function(x, ...) {
cat("ggcpt_delay\n")
cat_field("True changes", x$n_changes)
cat_field("Detected", x$n_detected)
cat_field("Mean delay",
if (is.nan(x$mean_delay)) "-" else format(x$mean_delay, digits = 4))
cat_field("Median delay",
if (is.na(x$median_delay)) "-" else
format(x$median_delay, digits = 4))
cat_field("False alarms", x$n_false_alarms)
cat_field("Average run length",
if (is.infinite(x$arl)) "no false alarms" else
format(x$arl, digits = 4))
cat("\n")
print(x$per_change, n = 10)
invisible(x)
}
#' @rdname cpt_delay
#' @param object A \code{ggcpt_delay} object (for \code{autoplot()}).
#' @export
autoplot.ggcpt_delay <- function(object, ...) {
d <- object$per_change
if (all(!d$detected)) {
cpt_abort("No change was detected, so there are no delays to draw.",
class = "capability_absent")
}
ggplot2::ggplot(d, ggplot2::aes(factor(truth), delay)) +
ggplot2::geom_col(ggplot2::aes(fill = detected), width = 0.6,
na.rm = TRUE) +
ggplot2::scale_fill_manual(values = c(`TRUE` = "#0072B2",
`FALSE` = "grey70"),
name = "Detected") +
ggplot2::labs(x = "True changepoint", y = "Detection delay (observations)",
title = "Detection delay",
subtitle = paste0(object$n_false_alarms,
" false alarm(s); average run length ",
if (is.infinite(object$arl)) {
"infinite"
} else {
format(object$arl, digits = 4)
}))
}
# Internal: the one moment the streaming API can test its independence
# assumption before a wrong answer exists. A monitor is handed an in-control
# baseline before any decision; measured (ยง206), under AR(1) noise with
# rho = 0.7 the e-detector raised 3.9 times and cpm 10 times their declared
# false-alarm rates. Warned once, at construction, rather than after the
# pager has been muted.
#' @noRd
check_baseline_dependence <- function(baseline, method) {
B <- if (is.matrix(baseline) || is.data.frame(baseline)) {
as.matrix(baseline)
} else {
matrix(as.numeric(baseline), ncol = 1)
}
if (nrow(B) < 20L) return(invisible(NULL))
deps <- lapply(seq_len(ncol(B)), function(j) {
residual_dependence(B[, j], integer(0))
})
p <- vapply(deps, function(d) if (is.null(d)) NA_real_ else d$p_value,
numeric(1))
acf1 <- vapply(deps, function(d) if (is.null(d)) NA_real_ else d$acf1,
numeric(1))
worst <- which.min(p)
if (!length(worst) || is.na(p[worst]) || p[worst] >= 0.05) {
return(invisible(NULL))
}
cpt_warn("`baseline` is autocorrelated (Ljung-Box p = ",
format(signif(p[worst], 2)), ", lag-1 autocorrelation ",
format(round(acf1[worst], 2)),
if (ncol(B) > 1L) paste0(", coordinate ", worst) else "",
"). The `", method, "` monitor's false-alarm rate assumes ",
"independent observations; measured under AR(1) noise with ",
"rho = 0.7, the monitors raised four to ten times their declared ",
"rate. Expect alarms sooner than `arl0` or `alpha` says, or model ",
"the dependence before monitoring (for example, monitor the ",
"residuals of an AR fit).", class = "dependence",
data = list(p_value = p[worst], acf1 = acf1[worst]))
invisible(NULL)
}
# Internal: the control chart of a monitor: its statistic at every
# observation against its limits, alarms flagged.
#' @noRd
control_chart <- function(object) {
tr <- object$trace
if (is.null(tr) || !nrow(tr)) {
cpt_abort("This monitor has consumed no observations yet (or was built ",
"before 0.7.0, which recorded no trace), so there is no chart ",
"to draw. Feed it with cpt_update().",
class = "capability_absent")
}
st <- object$state
al <- object$alarms
lim <- switch(object$method,
cusum = c(-st$h, st$h),
ewma = c(-st$threshold, st$threshold),
edetector = c(NA_real_, st$threshold),
ocd = c(NA_real_, 1),
c(NA_real_, NA_real_)
)
ylab <- switch(object$method,
cusum = "CUSUM (upper, lower)",
ewma = "EWMA (standardised)",
edetector = "E-detector",
"Statistic"
)
p <- ggplot2::ggplot(tr, ggplot2::aes(x = time, y = statistic)) +
ggplot2::geom_line(colour = "#0072B2", na.rm = TRUE)
if (any(!is.na(tr$lower))) {
p <- p + ggplot2::geom_line(ggplot2::aes(y = lower), colour = "#0072B2",
na.rm = TRUE)
}
lims <- lim[!is.na(lim)]
if (length(lims)) {
p <- p + ggplot2::geom_hline(yintercept = lims, linetype = "dashed",
colour = "#D55E00")
}
if (nrow(al)) {
hit <- tr[tr$time %in% al$time, , drop = FALSE]
yv <- if (object$method == "cusum") {
ifelse(abs(hit$lower) > hit$statistic, hit$lower, hit$statistic)
} else hit$statistic
p <- p + ggplot2::geom_point(data = tibble::tibble(time = hit$time,
statistic = yv),
colour = "#D55E00", size = 2.2,
na.rm = TRUE)
}
p + ggplot2::labs(x = "Observation", y = ylab,
title = paste0("Control chart (", object$method, ")"),
subtitle = paste0(nrow(al), " alarm(s) in ", object$t,
" observations; dashed: control ",
"limits")) +
theme_ggcpt()
}
Try the ggchangepoint package in your browser
Any scripts or data that you put into this service are public.
ggchangepoint documentation built on Nov. 5, 2026, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions control_chart check_baseline_dependence autoplot.ggcpt_delay print.ggcpt_delay glance.ggcpt_delay tidy.ggcpt_delay cpt_delay cpt_replay autoplot.ggcpt_monitor print.ggcpt_monitor tidy.ggcpt_monitor alarms.ggcpt_monitor alarms monitor_relearn monitor_width cpt_update chart_state cpt_monitor
Documented in alarms alarms.ggcpt_monitor autoplot.ggcpt_delay autoplot.ggcpt_monitor cpt_delay cpt_monitor cpt_replay cpt_update glance.ggcpt_delay print.ggcpt_delay print.ggcpt_monitor tidy.ggcpt_delay tidy.ggcpt_monitor
# ---------------------------------------------------------------------------
# Theme G: streaming and online monitoring.
#
# cpm, ocd and BOCPD are *online* methods that 0.4.0 ran in batch-replay mode
# and reported as if they were retrospective. The output the online
# literature cares about -- the alarm sequence, the detection delay, the
# average run length -- was not first-class, and there was no way to feed new
# observations to an existing fit. A monitor is a genuinely different object
# from a segmentation, so it gets its own class rather than another slot.
#
# `edetector` is the one method in this package implemented here rather than
# wrapped. That is a deliberate exception to the wrap-don't-implement rule
# (see ?cpt_monitor), taken because no R package implements e-detectors and
# the construction is short enough to be auditable in one screen.
# ---------------------------------------------------------------------------
#' A stateful sequential changepoint monitor
#'
#' Creates a detector that consumes observations as they arrive and raises
#' alarms, rather than segmenting a series that is already complete. Feed it
#' with \code{\link{cpt_update}()}, read its alarm log with
#' \code{\link{alarms}()}, and score it with \code{\link{cpt_delay}()},
#' because for an online method "did you find the location?" is the wrong
#' question and "how long did you take, and how often do you false-alarm?"
#' is the right one.
#'
#' @param method Which sequential detector:
#' \describe{
#' \item{\code{"edetector"}}{(default) a mixture Shiryaev-Roberts
#' e-detector; see the section below.}
#' \item{\code{"cpm"}}{\pkg{cpm}'s sequential change-point model, tuned
#' by \code{ARL0}.}
#' \item{\code{"ocd"}}{\pkg{ocd}'s high-dimensional online detector.
#' \strong{Multivariate only}: it tracks a projection of the whole
#' vector and needs at least two coordinates, so it refuses a single
#' series rather than falling back to a univariate statistic.}
#' \item{\code{"cusum"}}{a two-sided tabular CUSUM control chart
#' (Page 1954) on the standardised stream, with reference value
#' \code{k} and its decision interval calibrated to \code{arl0}.
#' Implemented in this package; see the section on control charts.}
#' \item{\code{"ewma"}}{an EWMA control chart (Roberts 1959) with
#' smoothing \code{lambda} and limits calibrated to \code{arl0}.
#' Implemented in this package.}
#' }
#' @param baseline A numeric vector (or, for \code{"ocd"}, a matrix with rows
#' as time points) of pre-change training data used to estimate the
#' in-control mean and scale. Required for \code{"edetector"} and
#' \code{"ocd"}; optional for \code{"cpm"}, which has its own start-up
#' period.
#' @param alpha Target false-alarm probability for \code{"edetector"}: the
#' average run length under the null is at least \code{1 / alpha}, by
#' optional stopping on the martingale \eqn{M_t - t} (see Details).
#' Defaults to \code{0.01}.
#' @param arl0 Target in-control average run length for \code{"cpm"},
#' \code{"cusum"} and \code{"ewma"}. Defaults to \code{500}.
#' @param cpm_type Statistic for \code{"cpm"}. Defaults to
#' \code{"Mann-Whitney"}.
#' @param patience Target patience (average run length) for \code{"ocd"}.
#' Defaults to \code{5000}.
#' @param deltas Shift sizes, in baseline standard deviations, mixed over by
#' \code{"edetector"}. Defaults to \code{c(0.5, 1, 2)}, each taken in both
#' directions. The mixture is a uniform average, so adding shifts costs
#' power at the ones already there rather than inflating the false-alarm
#' rate; mix over the range you think the change could fall in, not over
#' everything.
#' @param reset After an alarm, restart the detector (\code{TRUE}, the
#' default) or leave it running? Restarting is what makes a monitor report
#' \emph{repeated} changes rather than latching on the first one.
#' @param relearn How many observations after an alarm are used to re-learn
#' the in-control baseline, during which no further alarm can fire.
#' Defaults to \code{20}. This matters more than it looks: a real change
#' is \emph{persistent}, so a detector that restarts against the stale
#' pre-change baseline alarms again on the very next observation and keeps
#' alarming for the rest of the series: the monitor reports one change as
#' hundreds, and \code{\link{cpt_delay}()} then counts them all as false
#' alarms. Set \code{relearn = 0} to switch the behaviour off and see
#' every threshold crossing.
#' @param thresh Threshold rule for \code{"ocd"}: \code{"MC"} (default)
#' calibrates by Monte Carlo against \code{patience}, or supply a numeric
#' vector of three thresholds. The Monte Carlo calibration is the
#' expensive part of building an \code{"ocd"} monitor (a minute or more
#' at the default \code{patience}), so pass thresholds directly when you
#' already have them, or lower \code{mc_reps} while exploring.
#' @param mc_reps Monte Carlo repetitions for the \code{"ocd"} threshold.
#' @param ... Additional arguments passed to the engine's constructor
#' (\code{"cpm"} and \code{"ocd"}; the three native monitors take none).
#' @param k CUSUM reference value, in baseline standard deviations: half
#' the shift the chart is tuned to detect fastest. Defaults to \code{0.5}
#' (a one-standard-deviation shift).
#' @param lambda EWMA smoothing weight on the newest observation, in (0, 1].
#' Defaults to \code{0.2}; smaller values detect small shifts sooner and
#' large ones later.
#' @param limit Optional control limit, overriding the calibration to
#' \code{arl0}: the CUSUM decision interval \eqn{h}, or the EWMA limit
#' multiplier \eqn{L}, both in baseline standard deviations.
#'
#' @section Control charts:
#' \code{"cusum"} and \code{"ewma"} are the two charts statistical process
#' control runs, and they are sequential changepoint detectors under
#' another name. Each standardises the stream by the baseline's mean and
#' standard deviation (Phase I; the baseline should be in control) and
#' then monitors it (Phase II). The CUSUM keeps \eqn{S_t^+ =
#' \max(0, S_{t-1}^+ + z_t - k)} and \eqn{S_t^- = \max(0, S_{t-1}^- - z_t -
#' k)} and alarms when either reaches \eqn{h}; the EWMA keeps \eqn{z_t =
#' \lambda x_t + (1 - \lambda) z_{t-1}} and alarms when it leaves \eqn{\pm L
#' \sqrt{\lambda / (2 - \lambda)}}. The limit is solved for \code{arl0} by
#' the Markov-chain approximation of the run length (Brook and Evans 1972;
#' Lucas and Saccucci 1990), which reproduces the published tables: a CUSUM
#' with \code{k = 0.5} and \eqn{h = 4} or 5 has in-control run lengths of
#' 168 and 465, and an EWMA with \code{lambda = 0.1} and \eqn{L = 2.814}
#' one of 500. Measured through this function with a 5,000-observation
#' baseline and \code{arl0 = 100}, the first alarm on in-control data came
#' after 99.4 observations on average for the CUSUM and 99.8 for the EWMA
#' (400 runs each, standard errors about 5). The calibration treats the
#' baseline's mean and standard deviation as exact, so a short baseline
#' shortens the in-control run length: a few hundred observations is a
#' reasonable minimum. Rational subgroups, run rules and capability indices are not
#' modelled; for a formal quality-control programme use a dedicated
#' package such as \pkg{qcc}. \code{autoplot(plot_type = "control_chart")}
#' draws the chart.
#'
#' @section The e-detector, and why it is implemented rather than wrapped:
#' Shin, Ramdas and Rinaldo (2023) give a nonparametric sequential framework
#' with non-asymptotic control of the average run length, and it has no R
#' implementation. The construction used here is the mixture Shiryaev-Roberts
#' e-detector for a sub-Gaussian shift. For each candidate shift
#' \eqn{\delta} the increment is the likelihood ratio
#' \eqn{e_t^{(\delta)} = \exp(\delta (X_t - \mu_0)/\sigma^2 -
#' \delta^2/(2\sigma^2))}, which has unit mean under the null; the running
#' statistic is \eqn{R_t^{(\delta)} = (1 + R_{t-1}^{(\delta)})
#' e_t^{(\delta)}}; the shifts are combined by \emph{averaging},
#' \eqn{M_t = K^{-1} \sum_\delta R_t^{(\delta)}}; and an alarm is raised
#' the first time \eqn{M_t \ge 1/\alpha}.
#'
#' The averaging is not a detail. Under the null \eqn{M_t - t} is a
#' mean-zero martingale, so optional stopping at the alarm time \eqn{\tau}
#' gives \eqn{E_\infty[\tau] = E_\infty[M_\tau] \ge 1/\alpha}: a
#' finite-sample lower bound on the in-control average run length, with no
#' asymptotics and no calibration run. A convex combination of e-detectors
#' is an e-detector; a \emph{maximum} of them is not, and taking one
#' silently multiplies the false-alarm rate by roughly the number of shifts
#' mixed over. Apart from it and the two textbook control charts below,
#' every detector in this package wraps published, separately maintained
#' code; the three native monitors are labelled as such wherever they
#' appear.
#'
#' @return A \code{ggcpt_monitor} object: a list carrying the detector's
#' state between calls, so it is meant to be passed to
#' \code{\link{cpt_update}()} and re-assigned rather than read field by
#' field. The parts worth reading are
#' \describe{
#' \item{\code{alarms}}{the tibble \code{\link{alarms}()} and
#' \code{\link[=tidy.ggcpt_monitor]{tidy}()} return: one row per
#' alarm, with \code{time}, \code{statistic} and
#' \code{threshold}. Empty (zero rows) until something fires.}
#' \item{\code{t}}{how many observations the monitor has consumed
#' through \code{\link{cpt_update}()}, \strong{not} counting the
#' baseline; an alarm's \code{time} is on this same clock. A monitor
#' built by \code{\link{cpt_replay}()} also carries \code{offset},
#' the number of leading observations it trained on, so
#' \code{time + offset} is a position in the replayed series
#' (\code{\link{cpt_delay}()} applies it for you).}
#' \item{\code{data}}{every monitored observation so far, baseline
#' excluded (the first coordinate, for a multivariate monitor),
#' which is what \code{autoplot()} draws.}
#' \item{\code{method}, \code{alpha}, \code{n_baseline},
#' \code{reset}, \code{relearn}}{the settings this call fixed,
#' kept so a later \code{cpt_update()} cannot silently disagree
#' with them.}
#' \item{\code{state}}{the detector's internal statistics. Engine
#' internals: their shape differs by \code{method} and is not part
#' of the interface.}
#' }
#' @references
#' \insertRef{shin2023edetectors}{ggchangepoint}
#'
#' \insertRef{page1954continuous}{ggchangepoint}
#'
#' \insertRef{roberts1959ewma}{ggchangepoint}
#'
#' \insertRef{brook1972arl}{ggchangepoint}
#'
#' \insertRef{lucas1990ewma}{ggchangepoint}
#' @seealso \code{\link{cpt_update}()}, \code{\link{alarms}()},
#' \code{\link{cpt_replay}()}, \code{\link{cpt_delay}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_monitor("edetector", baseline = rnorm(100))
#' mon <- cpt_update(mon, rnorm(50)) # still in control
#' mon <- cpt_update(mon, rnorm(50, 3)) # a change arrives
#' alarms(mon)
cpt_monitor <- function(method = c("edetector", "cpm", "ocd", "cusum",
"ewma"),
baseline = NULL, alpha = 0.01, arl0 = 500,
cpm_type = "Mann-Whitney", patience = 5000,
deltas = c(0.5, 1, 2), reset = TRUE, relearn = 20,
thresh = "MC", mc_reps = 100, ..., k = 0.5,
lambda = 0.2, limit = NULL) {
method <- cpt_match_arg(method)
# The three detectors are calibrated in different currencies and each
# ignores the others': setting `arl0` on an e-detector, or `alpha` on
# `cpm`, changes nothing at all. The documentation and the vignette both
# say so, which does not help the user who set one and is now reading an
# unchanged answer. Only an argument the caller supplied is reported --
# the defaults supply all of them.
governs <- list(alpha = "edetector", deltas = "edetector",
arl0 = c("cpm", "cusum", "ewma"), cpm_type = "cpm",
patience = "ocd", thresh = "ocd", mc_reps = "ocd",
k = "cusum", lambda = "ewma", limit = c("cusum", "ewma"))
supplied <- c(alpha = !missing(alpha), deltas = !missing(deltas),
arl0 = !missing(arl0), cpm_type = !missing(cpm_type),
patience = !missing(patience), thresh = !missing(thresh),
mc_reps = !missing(mc_reps), k = !missing(k),
lambda = !missing(lambda), limit = !is.null(limit))
tunes <- vapply(governs, function(g) method %in% g, logical(1))
ignored <- names(governs)[supplied[names(governs)] & !tunes]
if (length(ignored) > 0) {
cpt_warn("`", paste(ignored, collapse = "`, `"), "` ",
if (length(ignored) > 1) "do" else "does", " not affect `method ",
"= \"", method, "\"`, which is tuned by `",
paste(names(governs)[tunes], collapse = "`, `"),
"`. See ?cpt_monitor.", class = "argument_ignored")
}
# `...` reaches cpm::processStream() or ocd::ChangepointDetector(), and
# this function renames three of their arguments and pins two more.
# Passing an engine's own name gave R's raw "formal argument \"cpmType\"
# matched by multiple actual arguments" -- `cpmType`, `ARL0`, `dim`,
# `MC_reps` and `beta`, measured across the two engines. `thresh`,
# `patience` and `startup` cannot reach `...` at all, because they are
# formals here and R matches them there first. `cpt_replay()` forwards
# `...` to this function, so it is covered by the same guard.
reject_renamed_args(list(...), method)
if (identical(method, "ocd")) {
reject_managed_args(list(...), "ocd", c(
beta = paste("the monitor calibrates against a unit-scale change;",
"use `ocd_wrapper()` directly to set it")))
}
validate_scalar(alpha, "alpha", min = 0, max = 1,
min_open = TRUE, max_open = TRUE)
validate_scalar(relearn, "relearn", min = 0)
validate_flag(reset, "reset")
# Three of ten arguments were checked. The rest reached their engine (or,
# for `deltas`, this function's own arithmetic) unvalidated: `deltas = 0`
# makes every `inc` exactly 1, so `R` grows LINEARLY whatever the data
# does and the monitor alarms at t = 1/alpha on pure noise.
if (method %in% c("cusum", "ewma") && length(list(...)) > 0) {
extra <- names(list(...))
cpt_warn("`method = \"", method, "\"` is native to this package and ",
"takes no engine arguments, so ",
if (length(extra) && all(nzchar(extra))) {
paste0("`", paste(extra, collapse = "`, `"), "` ")
} else {
paste0(length(list(...)), " unnamed argument(s) ")
}, "would be ignored. It is tuned by `arl0`, `",
if (method == "cusum") "k" else "lambda", "` and `limit`.",
class = "argument_ignored")
}
if (identical(method, "edetector")) {
if (!is.numeric(deltas) || length(deltas) == 0L ||
any(!is.finite(deltas)) || all(deltas == 0)) {
cpt_abort("`deltas` must be one or more finite non-zero numbers: the ",
"shift ", "sizes, in baseline standard deviations, the ",
"e-detector mixes ", "over. A zero shift contributes a ",
"likelihood ratio of exactly 1 ", "at every observation, so ",
"the statistic grows on nothing.", class = "bad_argument")
}
if (any(deltas == 0)) {
cpt_warn("`deltas` contains 0, which contributes a likelihood ratio of ",
"exactly 1 at every observation and only dilutes the mixture. ",
"Dropping it.", class = "dropped_input")
deltas <- deltas[deltas != 0]
}
# `...` reaches the engine in the cpm and ocd branches and has nowhere
# to go in this one, so it used to vanish -- the same silent-discard
# ggcptplot_internal() warns about.
extra <- names(list(...))
if (length(list(...)) > 0) {
cpt_warn("`method = \"edetector\"` is native to this package and takes ",
"no engine arguments, so ",
if (length(extra) && all(nzchar(extra))) {
paste0("`", paste(extra, collapse = "`, `"), "` ")
} else {
paste0(length(list(...)), " unnamed argument(s) ")
},
"would be ignored. It is tuned by `alpha` and `deltas`.",
class = "argument_ignored")
}
}
if (identical(method, "ocd")) {
validate_scalar(patience, "patience", min = 1)
validate_scalar(mc_reps, "mc_reps", min = 1)
}
# Normalise the series ONCE, before the switch. Each branch used to coerce
# `baseline` its own way -- as.numeric() for edetector and cpm, a matrix
# for ocd -- so a factor reached two of the three detectors as its level
# codes. Guarding here means every method enters through the same door,
# which is the only way this stays true when a fourth branch is added.
if (!is.null(baseline)) {
if (is.matrix(baseline) || is.data.frame(baseline)) {
# `edetector` and `cpm` monitor ONE series, and length() on a
# data.frame counts its columns -- so a 60-row, 2-column baseline was
# refused for having "at least 5 pre-change observations", naming a
# count of 2 for 60 observations. Refuse the shape instead, since a
# multi-column baseline for a univariate detector is outside the
# documented contract either way.
if (!identical(method, "ocd")) {
reject_multicolumn(baseline, "baseline",
paste0("`method = \"", method, "\"` monitors one ",
"series; use `method = \"ocd\"` for a ",
"multivariate stream."))
baseline <- as_mv_matrix(baseline, arg = "baseline")[, 1]
}
} else {
baseline <- coerce_series_values(baseline, arg = "baseline")
}
}
state <- switch(method,
edetector = {
if (is.null(baseline) || NROW(baseline) < 5) {
cpt_abort("`edetector` estimates the in-control mean and scale from ",
"`baseline`, so it needs at least 5 pre-change observations.",
class = "short_series")
}
b <- as.numeric(baseline)
# Check finiteness first: sd() of a vector holding an NA is itself NA,
# so the variability test below would blame a flat baseline for what is
# really a missing value.
if (anyNA(b) || any(!is.finite(b))) {
stop_nonfinite(b, "baseline", stream_na_hint(b))
}
sd0 <- stats::sd(b)
if (!is.finite(sd0) || sd0 <= 0) {
cpt_abort("`baseline` has zero variability, so the e-detector's scale ",
"is ", "undefined.", class = "input_error")
}
deltas <- sort(unique(c(-abs(deltas), abs(deltas)))) * sd0
list(mu0 = mean(b), sd0 = sd0, deltas = deltas,
R = rep(0, length(deltas)),
threshold = 1 / alpha)
},
cpm = {
need_pkg("cpm")
# The same three guards cpm_wrapper() applies, from the same helpers:
# this branch used to call makeChangePointModel() directly, so a
# withheld `cpm_type`, a missing FET `lambda` or an off-grid `arl0`
# each reached the user as an error from inside cpm or base R that
# named no argument at all (`no applicable method for '@' applied to
# an object of class "NULL"`, because cpm *prints* its complaint and
# returns NULL). See cpm_check_type() for what each one looks like.
validate_scalar(arl0, "arl0")
cpm_type <- cpm_check_type(cpm_type, list(...))
engine_dots <- cpm_engine_dots(cpm_type, list(...))
cpm_out <- utils::capture.output(
m <- do.call(cpm::makeChangePointModel,
c(list(cpmType = cpm_type, ARL0 = arl0), engine_dots))
)
cpm_check_printed_error(cpm_out, arl0, engine_dots)
if (length(cpm_out)) cat(cpm_out, sep = "\n")
# A change cpm detects INSIDE the baseline left the model tripped:
# changeDetected() was already TRUE when the loop ended, so the first
# monitored observation inherited the alarm. Measured before this
# was fixed: an alarm at stream index 1 in 2 of 6 clean iid runs and
# 6 of 6 AR(1) runs, and never without a baseline. The model is now
# restarted after each baseline detection, as cpt_update() does after
# an alarm, and the caller is told the baseline was not in control.
n_tripped <- 0L
if (!is.null(baseline)) {
for (v in as.numeric(baseline)) {
m <- cpm::processObservation(m, v)
if (isTRUE(cpm::changeDetected(m))) {
n_tripped <- n_tripped + 1L
m <- cpm::cpmReset(m)
}
}
}
if (n_tripped > 0L) {
cpt_warn("cpm detected ", n_tripped, " change(s) inside `baseline`, ",
"so the baseline is not in control. The monitor was ",
"restarted after each, and only the stretch after the last ",
"one informs it; check the baseline, or build the monitor ",
"without one.", class = "assumption",
data = list(n_changes = n_tripped))
}
list(model = m, cpm_type = cpm_type, arl0 = arl0)
},
cusum = ,
ewma = chart_state(method, baseline, arl0 = arl0, k = k,
lambda = lambda, limit = limit),
ocd = {
need_pkg("ocd")
if (is.null(baseline)) {
cpt_abort("`ocd` needs `baseline` to estimate the pre-change mean and ",
"standard deviation.", class = "bad_argument")
}
B <- if (is.matrix(baseline) || is.data.frame(baseline)) {
as_mv_matrix(baseline, arg = "baseline")
} else {
matrix(baseline, ncol = 1)
}
p <- ncol(B)
if (p < 2) {
cpt_abort("Method `ocd` is high-dimensional and needs at least two ",
"coordinates, but `baseline` has ", p, ". Use method = ",
"\"edetector\" or \"cpm\" for one series.",
class = "wrong_dimension")
}
det <- ocd::ChangepointDetector(dim = p, method = "ocd",
thresh = thresh, beta = 1,
patience = patience,
MC_reps = mc_reps, ...)
det <- ocd::setBaselineMean(det, colMeans(B))
det <- ocd::setBaselineSD(det, apply(B, 2, stats::sd))
det <- ocd::setStatus(det, "monitoring")
list(detector = det, p = p, patience = patience,
mean = colMeans(B),
sd = { sdv <- apply(B, 2, stats::sd)
sdv[!is.finite(sdv) | sdv == 0] <- 1
sdv })
}
)
if (!is.null(baseline)) check_baseline_dependence(baseline, method)
structure(
list(method = method, state = state, reset = reset,
relearn = as.integer(relearn),
n_baseline = if (is.null(baseline)) 0L else NROW(baseline),
t = 0L,
relearn_left = 0L,
relearn_buffer = NULL,
data = numeric(0),
trace = tibble::tibble(time = integer(), statistic = numeric(),
lower = numeric()),
alarms = tibble::tibble(time = integer(), statistic = numeric(),
threshold = numeric()),
alpha = alpha),
class = "ggcpt_monitor"
)
}
# Internal: the state of a CUSUM or EWMA chart, standardised by the
# baseline and calibrated to arl0 unless a limit is given.
#' @noRd
chart_state <- function(method, baseline, arl0, k, lambda, limit) {
if (is.null(baseline) || NROW(baseline) < 5) {
cpt_abort("`", method, "` standardises the stream by the in-control ",
"mean and standard deviation of `baseline`, so it needs at ",
"least 5 pre-change observations.", class = "short_series")
}
b <- as.numeric(baseline)
if (anyNA(b) || any(!is.finite(b))) {
stop_nonfinite(b, "baseline", stream_na_hint(b))
}
sd0 <- stats::sd(b)
if (!is.finite(sd0) || sd0 <= 0) {
cpt_abort("`baseline` has zero variability, so the chart's scale is ",
"undefined.", class = "input_error")
}
check_chart_arl0(arl0)
if (!is.null(limit)) validate_scalar(limit, "limit", min = 0,
min_open = TRUE)
if (method == "cusum") {
validate_scalar(k, "k", min = 0)
h <- limit %||% cusum_h_for(arl0, k)
return(list(mu0 = mean(b), sd0 = sd0, k = k, h = h, sp = 0, sn = 0,
threshold = h, arl0 = arl0,
calibrated = is.null(limit)))
}
validate_scalar(lambda, "lambda", min = 0, max = 1, min_open = TRUE)
L <- limit %||% ewma_L_for(arl0, lambda)
list(mu0 = mean(b), sd0 = sd0, lambda = lambda, L = L, z = 0,
threshold = L * sqrt(lambda / (2 - lambda)), arl0 = arl0,
calibrated = is.null(limit))
}
#' Feed observations to a monitor
#'
#' Pushes new data through a \code{\link{cpt_monitor}()}, updating its
#' internal state and appending to its alarm log. The monitor is returned, so
#' the idiom is \code{mon <- cpt_update(mon, new_obs)}.
#'
#' @param monitor A \code{ggcpt_monitor} object.
#' @param new_obs New observations: a numeric vector, or a matrix with rows
#' as time points for a multivariate monitor.
#' @return The updated \code{ggcpt_monitor}: the same structure as the
#' input (see \code{\link{cpt_monitor}()}), with \code{t} advanced by
#' the number of new observations (rows, for a matrix) and any new rows
#' appended to \code{alarms}.
#' Re-assign it: the monitor is a value, not a mutable object, so
#' \code{cpt_update(mon, y)} without assignment discards the update.
#' @seealso \code{\link{cpt_monitor}()}, \code{\link{alarms}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_monitor("edetector", baseline = rnorm(100))
#' mon <- cpt_update(mon, c(rnorm(50), rnorm(50, 3)))
#' alarms(mon)
cpt_update <- function(monitor, new_obs) {
if (!inherits(monitor, "ggcpt_monitor")) {
cpt_abort("`monitor` must be a ggcpt_monitor from cpt_monitor().",
class = "bad_argument")
}
X <- if (is.matrix(new_obs) || is.data.frame(new_obs)) {
as_mv_matrix(new_obs, arg = "new_obs")
} else {
# a streaming caller passes observations one at a time, so this is the
# most likely place for a factor level code to enter unnoticed
matrix(coerce_series_values(new_obs, arg = "new_obs"), ncol = 1)
}
if (nrow(X) == 0) return(monitor)
if (anyNA(X) || any(!is.finite(X))) {
stop_nonfinite(X, "new_obs", stream_na_hint(X))
}
# A monitor is stateful and dimensioned by its baseline. Feeding it a
# different width is a mistake, not a coercion: `ocd` would consume the
# wrong coordinates and the univariate detectors would silently read
# column 1 and ignore the rest.
expected <- monitor_width(monitor)
if (!is.na(expected) && ncol(X) != expected) {
cpt_abort("This monitor was built on ", expected, " coordinate(s) but ",
"`new_obs` has ", ncol(X), ". A monitor cannot change width ",
"once it is running.", class = "wrong_dimension")
}
st <- monitor$state
new_alarms <- list()
# Grown once at the end, not re-allocated per observation. cpt_replay()
# hands the whole series to a single cpt_update() call, so
# `monitor$data <- c(monitor$data, xt[1])` inside the loop made replaying
# 10,000 observations do on the order of 50 million element copies.
seen <- numeric(nrow(X))
trace_stat <- rep(NA_real_, nrow(X))
trace_low <- rep(NA_real_, nrow(X))
for (i in seq_len(nrow(X))) {
monitor$t <- monitor$t + 1L
xt <- X[i, ]
# Coordinate one only, for every method: `$data` is the series
# autoplot() draws, and a multivariate `ocd` monitor is drawn as its
# first coordinate (see @return).
seen[i] <- xt[1]
# Re-learning window after an alarm: collect the new regime and stay
# silent, then adopt it as the in-control baseline. Without this a
# persistent change makes the detector alarm on every subsequent
# observation (see the `relearn` argument of cpt_monitor()).
if (monitor$relearn_left > 0L) {
monitor$relearn_buffer <- rbind(monitor$relearn_buffer, xt)
monitor$relearn_left <- monitor$relearn_left - 1L
if (monitor$relearn_left == 0L) {
st <- monitor_relearn(monitor$method, st, monitor$relearn_buffer)
monitor$relearn_buffer <- NULL
}
next
}
hit <- FALSE
stat <- NA_real_
thr <- NA_real_
if (monitor$method == "edetector") {
inc <- exp(st$deltas * (xt[1] - st$mu0) / st$sd0^2 -
st$deltas^2 / (2 * st$sd0^2))
st$R <- (1 + st$R) * inc
# An e-detector's mixture grows multiplicatively, so under `reset =
# FALSE` with `relearn = 0` -- the "keep accumulating" configuration
# the arguments offer -- it passes .Machine$double.xmax a few hundred
# observations after a real change and the next product is `Inf`. The
# alarm rule below reads a non-finite statistic as "no alarm", so the
# monitor went permanently silent exactly when the evidence was
# overwhelming: measured 82 alarms and then 1418 observations of
# silence on a shifted stream. Saturating instead keeps the statistic
# ordered against the threshold, which is all the rule needs, and --
# unlike `Inf`, which is absorbing -- it still decays when the stream
# returns to its baseline (measured: 8.99e307 down to 1.8e149 over
# 2000 in-control observations, where `Inf` would have stayed `Inf`).
st$R <- pmin(st$R, .Machine$double.xmax / 2)
# AVERAGE, not max: a convex combination of e-detectors is an
# e-detector and keeps E[M_t] = t, which is what the average-run-length
# bound rests on. A maximum over K shifts crosses the threshold roughly
# K times as often under the null -- measured here at ~1.8x the bound
# with six shifts before this was corrected.
stat <- mean(st$R)
thr <- st$threshold
hit <- is.finite(stat) && stat >= thr
if (hit && monitor$reset) st$R <- rep(0, length(st$deltas))
} else if (monitor$method == "cusum") {
z <- (xt[1] - st$mu0) / st$sd0
st$sp <- max(0, st$sp + z - st$k)
st$sn <- max(0, st$sn - z - st$k)
trace_low[i] <- -st$sn
stat <- max(st$sp, st$sn)
thr <- st$h
hit <- stat >= thr
trace_stat[i] <- st$sp
if (hit && monitor$reset) {
st$sp <- 0
st$sn <- 0
}
} else if (monitor$method == "ewma") {
z <- (xt[1] - st$mu0) / st$sd0
st$z <- st$lambda * z + (1 - st$lambda) * st$z
stat <- st$z
thr <- st$threshold
hit <- abs(stat) >= thr
trace_stat[i] <- st$z
if (hit && monitor$reset) st$z <- 0
} else if (monitor$method == "cpm") {
st$model <- cpm::processObservation(st$model, xt[1])
stat <- suppressWarnings(as.numeric(
utils::tail(cpm::getStatistics(st$model), 1)
))
if (length(stat) == 0) stat <- NA_real_
hit <- isTRUE(cpm::changeDetected(st$model))
if (hit && monitor$reset) {
st$model <- cpm::cpmReset(st$model)
}
} else {
# ocd monitors a standardised stream, so the baseline enters here
# rather than through the detector object.
z <- (as.numeric(xt) - st$mean) / st$sd
utils::capture.output(st$detector <- ocd::getData(st$detector, z))
# `normalisedStatistics()` already divides each statistic by its own
# threshold, so the comparison point is 1. Reporting the raw
# thresholds alongside the normalised statistic would print an alarm
# at 1.10 against a threshold of 16.4, which reads as a bug.
stat <- suppressWarnings(max(as.numeric(
unlist(ocd::normalisedStatistics(st$detector))
), na.rm = TRUE))
thr <- 1
hit <- !identical(ocd::status(st$detector), "monitoring")
if (hit && monitor$reset) {
st$detector <- ocd::reset(st$detector)
st$detector <- ocd::setStatus(st$detector, "monitoring")
}
}
if (!monitor$method %in% c("cusum", "ewma")) trace_stat[i] <- stat
if (isTRUE(hit)) {
new_alarms[[length(new_alarms) + 1L]] <-
tibble::tibble(time = monitor$t, statistic = stat, threshold = thr)
if (monitor$relearn > 0L) {
monitor$relearn_left <- monitor$relearn
monitor$relearn_buffer <- NULL
}
}
}
monitor$state <- st
monitor$data <- c(monitor$data, seen)
# The statistic at every observation (NA while re-learning), so a chart
# can be drawn; a monitor built before 0.7.0 has no trace to extend.
steps <- tibble::tibble(time = monitor$t - rev(seq_len(nrow(X))) + 1L,
statistic = trace_stat, lower = trace_low)
monitor$trace <- rbind(monitor$trace %||% steps[0, ], steps)
if (length(new_alarms) > 0) {
monitor$alarms <- rbind(monitor$alarms, do.call(rbind, new_alarms))
}
monitor
}
# Internal: how many coordinates this monitor was built for. `edetector`
# and `cpm` read a single series; `ocd` carries the dimension it was
# constructed with.
#' @noRd
monitor_width <- function(monitor) {
switch(monitor$method,
ocd = monitor$state$p,
edetector = 1L,
cpm = 1L,
cusum = 1L,
ewma = 1L,
NA_integer_
)
}
# Internal: adopt a post-alarm window as the new in-control baseline. Only
# the detectors that hold an explicit baseline need this; cpm re-derives its
# own reference from the stream after a reset.
#' @noRd
monitor_relearn <- function(method, st, buffer) {
if (is.null(buffer) || nrow(buffer) == 0) return(st)
if (method == "edetector") {
v <- as.numeric(buffer[, 1])
sd0 <- stats::sd(v)
if (is.finite(sd0) && sd0 > 0) {
st$deltas <- st$deltas / st$sd0 * sd0
st$sd0 <- sd0
}
st$mu0 <- mean(v)
st$R <- rep(0, length(st$deltas))
} else if (method %in% c("cusum", "ewma")) {
v <- as.numeric(buffer[, 1])
sd0 <- stats::sd(v)
if (is.finite(sd0) && sd0 > 0) st$sd0 <- sd0
st$mu0 <- mean(v)
if (method == "cusum") {
st$sp <- 0
st$sn <- 0
} else {
st$z <- 0
}
} else if (method == "ocd") {
st$mean <- colMeans(buffer)
sdv <- apply(buffer, 2, stats::sd)
sdv[!is.finite(sdv) | sdv == 0] <- 1
st$sd <- sdv
}
st
}
#' The alarm log of a monitor
#'
#' Every observation at which a sequential monitor crossed its threshold, in
#' the order they fired. A monitor that never fired returns a zero-row tibble
#' rather than \code{NULL}, so the result is always safe to \code{rbind()} or
#' plot.
#'
#' @param x A \code{ggcpt_monitor} object.
#' @param ... Ignored.
#' @return A tibble with one row per alarm: \code{time} (the observation at
#' which it fired, counted from the first one fed to the monitor),
#' \code{statistic} and \code{threshold}.
#' @seealso \code{\link{cpt_monitor}()}, \code{\link{cpt_delay}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_monitor("edetector", baseline = rnorm(100))
#' mon <- cpt_update(mon, c(rnorm(50), rnorm(50, 3)))
#' alarms(mon)
alarms <- function(x, ...) {
UseMethod("alarms")
}
#' @rdname alarms
#' @export
alarms.ggcpt_monitor <- function(x, ...) x$alarms
#' @rdname cpt_monitor
#' @export
tidy.ggcpt_monitor <- function(x, ...) {
alarms(x)
}
#' @rdname cpt_monitor
#' @param x A \code{ggcpt_monitor} object (for \code{print()}).
#' @export
print.ggcpt_monitor <- function(x, ...) {
cat("ggcpt_monitor (", x$method,
if (x$method %in% c("edetector", "cusum", "ewma")) {
", native implementation"
} else "",
")\n", sep = "")
st <- x$state
if (identical(x$method, "cusum")) {
cat(" Chart: k = ", format(st$k), ", h = ", format(signif(st$h, 4)),
if (isTRUE(st$calibrated)) paste0(" (in-control ARL ",
format(st$arl0), ")"),
"\n", sep = "")
} else if (identical(x$method, "ewma")) {
cat(" Chart: lambda = ", format(st$lambda), ", L = ",
format(signif(st$L, 4)),
if (isTRUE(st$calibrated)) paste0(" (in-control ARL ",
format(st$arl0), ")"),
"\n", sep = "")
}
cat(" Baseline observations: ", x$n_baseline, "\n", sep = "")
cat(" Monitored observations: ", x$t, "\n", sep = "")
cat(" Alarms: ", nrow(x$alarms), "\n", sep = "")
cat(" Restart after alarm: ", x$reset,
if (x$relearn > 0L) paste0(" (re-learning ", x$relearn,
" observations)") else "",
"\n", sep = "")
if (nrow(x$alarms) > 0) {
cat("\n")
print(x$alarms, n = 10)
}
invisible(x)
}
#' @rdname cpt_monitor
#' @param object A \code{ggcpt_monitor} object (for \code{autoplot()}).
#' @param plot_type \code{"timeline"} (the monitored series with the alarms
#' marked), \code{"control_chart"} (the detection statistic at every
#' observation against its limits, with the alarms flagged: for a CUSUM
#' the upper and lower sums against \eqn{\pm h}, for an EWMA the
#' smoothed statistic between its limits), \code{"statistic"} (the
#' statistic at each alarm only, against its threshold) or
#' \code{"runlength"} (the gaps between alarms, which estimate the run
#' length). For a multivariate monitor the timeline draws
#' the first coordinate: the alarms are shared, so the rules are right
#' whichever coordinate is shown, but the line is one of several.
#' @export
autoplot.ggcpt_monitor <- function(object,
plot_type = c("timeline", "statistic",
"runlength",
"control_chart"),
...) {
plot_type <- cpt_match_arg(plot_type)
al <- object$alarms
d <- tibble::tibble(time = seq_along(object$data), value = object$data)
if (plot_type == "control_chart") return(control_chart(object))
if (plot_type == "timeline") {
p <- ggplot2::ggplot(d, ggplot2::aes(time, value)) +
ggplot2::geom_line(colour = "grey40") +
ggplot2::labs(x = "Observation", y = "Value",
title = paste0("Monitoring timeline (", object$method,
")"),
subtitle = paste0(nrow(al), " alarm(s) in ", object$t,
" observations"))
if (nrow(al) > 0) {
p <- p + ggplot2::geom_vline(xintercept = al$time, colour = "#D55E00",
linetype = "dashed", linewidth = 0.5)
}
return(p)
}
if (plot_type == "statistic") {
if (nrow(al) == 0) {
# cpt_replay() used to be offered here as the way to a "full
# statistic trace". It returns the same object, which records the
# statistic at alarms only, so the advice led nowhere.
cpt_abort("No alarm has fired, so there is no statistic to draw: a ",
"monitor records its statistic only when an alarm fires. ",
"`plot_type = \"timeline\"` shows the monitored series.",
class = "capability_absent")
}
return(
ggplot2::ggplot(al, ggplot2::aes(time, statistic)) +
ggplot2::geom_point(size = 2, colour = "#0072B2") +
ggplot2::geom_hline(ggplot2::aes(yintercept = threshold),
linetype = "dotted", colour = "grey30",
na.rm = TRUE) +
ggplot2::labs(x = "Observation", y = "Detection statistic",
title = "Statistic at each alarm")
)
}
# runlength
if (nrow(al) < 2) {
cpt_abort("At least two alarms are needed to show run lengths; this ",
"monitor ", "has ", nrow(al), ".", class = "capability_absent")
}
rl <- tibble::tibble(run_length = diff(c(0L, al$time)))
ggplot2::ggplot(rl, ggplot2::aes(run_length)) +
ggplot2::geom_histogram(bins = max(5, min(30, nrow(rl))),
fill = "#0072B2", colour = "white") +
ggplot2::geom_vline(xintercept = mean(rl$run_length), colour = "#D55E00",
linetype = "dashed") +
ggplot2::labs(x = "Run length (observations between alarms)",
y = "Count",
title = "Run-length distribution",
subtitle = paste0("Mean run length ",
format(mean(rl$run_length), digits = 4)))
}
#' Replay a series through a sequential detector
#'
#' Runs a whole series through a monitor in one call, so the alarm timeline
#' of an online method can be studied retrospectively without hand-rolling
#' the update loop. This is the honest version of what 0.4.0's \code{cpm},
#' \code{ocd} and \code{bocpd} wrappers do implicitly: it reports
#' \emph{alarms}, which arrive after the change, rather than pretending they
#' are estimated changepoint locations.
#'
#' @param x A numeric vector, or a matrix with rows as time points.
#' @param method Sequential detector; see \code{\link{cpt_monitor}()}.
#' @param baseline Number of leading observations used as pre-change training
#' data (an integer), or an explicit baseline vector. Defaults to
#' \code{min(100, floor(n / 4))}.
#' @param ... Additional arguments passed to \code{\link{cpt_monitor}()}.
#' @return A \code{ggcpt_monitor} that has already consumed the series:
#' the same object \code{\link{cpt_monitor}()} returns, described
#' there, with \code{alarms} filled in for whatever fired.
#' @seealso \code{\link{cpt_delay}()}, \code{\link{cpt_monitor}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_replay(c(rnorm(200), rnorm(200, 3)), method = "edetector")
#' alarms(mon)
cpt_replay <- function(x, method = c("edetector", "cpm", "ocd", "cusum",
"ewma"),
baseline = NULL, ...) {
method <- cpt_match_arg(method)
# coerce_series_values()/as_mv_matrix(), not a bare as.numeric(): a factor
# coerces to its LEVEL CODES, which for labels like "10", "2", "30" is the
# alphabetical order 3, 1, 2 rather than the numbers -- so a replay ran on
# a series the user never supplied. cpt_detect() has refused that since
# 0.4.0, and as_mv_matrix() already refuses a factor *column*; only the
# vector path was open.
X <- if (is.matrix(x) || is.data.frame(x)) as_mv_matrix(x) else {
matrix(coerce_series_values(x), ncol = 1)
}
n <- nrow(X)
if (n < 10) {
cpt_abort("A replay needs at least 10 observations; got ", n, ".",
class = "short_series")
}
# Without this, a non-finite value was caught downstream by whichever of
# cpt_monitor() or cpt_update() happened to receive the slice containing
# it -- so `cpt_replay(x)` reported a problem with `baseline` for an NA at
# position 20 and with `new_obs` for one at 95, named an argument the
# caller never passed, and counted "1 of 45 values" against the baseline
# slice rather than the series. Checked here, the message describes `x`.
if (anyNA(X) || any(!is.finite(X))) {
stop_nonfinite(X, hint = stream_na_hint(X))
}
if (is.null(baseline)) baseline <- min(100L, floor(n / 4))
# `baseline` is documented two ways -- a count of leading observations, or
# an explicit series -- so it needs the same type guard as `x` above before
# either reading is applied. The explicit-series branch used a bare
# as.numeric(), and it hands cpt_monitor() an already-numeric vector, so
# cpt_monitor()'s own guard could never see the factor.
if (!is.matrix(baseline) && !is.data.frame(baseline)) {
baseline <- coerce_series_values(baseline, arg = "baseline")
}
if (length(baseline) == 1L) {
# A lone number is the count. It used to fall through to the
# explicit-series branch whenever it was fractional or out of range,
# which turned it into a ONE-point baseline: cpt_replay(x, baseline =
# 500) on a 180-point series then failed with "needs at least 5
# pre-change observations" and never mentioned the 500.
if (is.na(baseline) || baseline %% 1 != 0 || baseline < 1 ||
baseline >= n) {
cpt_abort("`baseline` as a single number is the count of leading ",
"observations to train on, so it must be a whole number in ",
"1..", n - 1L, "; the series has ", n, " observations and ",
"`baseline` is ", format(baseline), ". Pass an explicit ",
"baseline series if you meant a value rather than a count.",
class = "bad_argument")
}
n_base <- as.integer(baseline)
base_data <- X[seq_len(n_base), , drop = FALSE]
rest <- X[seq.int(n_base + 1L, n), , drop = FALSE]
} else {
base_data <- if (is.matrix(baseline) || is.data.frame(baseline)) {
as_mv_matrix(baseline, arg = "baseline")
} else {
matrix(baseline, ncol = ncol(X))
}
n_base <- 0L
rest <- X
}
mon <- cpt_monitor(method,
baseline = if (ncol(base_data) == 1) {
as.numeric(base_data)
} else {
base_data
},
...)
mon <- cpt_update(mon, rest)
mon$offset <- n_base
mon
}
#' Detection delay and false-alarm rate
#'
#' Scores an online result the way the sequential literature does: how long
#' after each true change did the first alarm arrive, and how many alarms
#' were raised with no change behind them. \code{\link{cpt_metrics}()} is the
#' wrong tool for an online detector: it asks whether the \emph{location}
#' was recovered, which a sequential procedure never claims. It will not
#' stop you: \code{cpt_metrics()} takes bare integer vectors and never
#' sees which detector produced them, so it cannot know. (An earlier
#' version of this sentence said it warns. It does not, and given that
#' signature it could not.)
#'
#' @param object A \code{ggcpt_monitor}, or a tibble of alarms with a
#' \code{time} column.
#' @param truth Integer vector of true changepoint positions, on the same
#' clock as the alarms. When \code{object} is a monitor built by
#' \code{\link{cpt_replay}()} the baseline offset is applied automatically.
#' @param max_delay Alarms further than this after a change are treated as
#' false alarms rather than late detections. Defaults to \code{Inf}.
#' @return A \code{ggcpt_delay} object: a list with \code{per_change} (one
#' row per true change: \code{truth}, \code{alarm}, \code{delay},
#' \code{detected}), \code{false_alarms}, and the summary statistics
#' \code{mean_delay}, \code{median_delay}, \code{n_false_alarms} and
#' \code{arl}: monitored observations per false alarm. The baseline a
#' \code{\link{cpt_replay}()} monitor trained on is not counted, because
#' no alarm can fire there; \code{n_obs} is the length of the stream in
#' series positions, baseline included.
#' @seealso \code{\link{cpt_monitor}()}, \code{\link{cpt_replay}()}.
#' @export
#' @examples
#' set.seed(2026)
#' mon <- cpt_replay(c(rnorm(200), rnorm(200, 3)), method = "edetector")
#' cpt_delay(mon, truth = 200)
cpt_delay <- function(object, truth, max_delay = Inf) {
# Compared straight against alarm times, so `NA` or a negative value
# silently counted every alarm as false, a vector recycled, and a string
# compared as text.
if (!is.numeric(max_delay) || length(max_delay) != 1L || is.na(max_delay) ||
max_delay < 0) {
cpt_abort("`max_delay` must be a single non-negative number (`Inf` for no ",
"limit); got ", paste(format(max_delay), collapse = ", "), ".",
class = "bad_argument")
}
if (missing(truth)) {
cpt_abort("`truth` is required: detection delay is measured from the true ",
"changepoint(s), so there is nothing to measure without them. Pass ",
"the location(s) as an integer vector.", class = "bad_argument")
}
al <- if (inherits(object, "ggcpt_monitor")) {
a <- object$alarms
off <- object$offset %||% 0L
a$time <- a$time + off
a
} else if (is.data.frame(object) && "time" %in% names(object)) {
tibble::as_tibble(object)
} else {
cpt_abort("`object` must be a ggcpt_monitor or a tibble with a `time` ",
"column.", class = "bad_argument")
}
truth <- as_cp_locations(truth, "truth", sort = TRUE)
if (length(truth) == 0 || anyNA(truth) || any(truth < 1)) {
cpt_abort("`truth` must be one or more positive changepoint locations.",
class = "bad_argument")
}
n_obs <- if (inherits(object, "ggcpt_monitor")) {
object$t + (object$offset %||% 0L)
} else {
max(c(al$time, truth), na.rm = TRUE)
}
n_monitored <- if (inherits(object, "ggcpt_monitor")) object$t else n_obs
# A truth past the end of the stream can never be detected, so scoring it
# would report a miss that says more about the argument than the detector.
if (inherits(object, "ggcpt_monitor") && any(truth > n_obs)) {
cpt_abort("`truth` runs past the end of the stream: ",
paste(utils::head(truth[truth > n_obs], 5), collapse = ", "),
" vs ", n_obs, " observations seen.", class = "bad_argument")
}
used <- rep(FALSE, nrow(al))
rows <- lapply(seq_along(truth), function(i) {
k <- truth[i]
next_change <- if (i < length(truth)) truth[i + 1] else Inf
cand <- which(!used & al$time >= k & al$time < next_change &
al$time - k <= max_delay)
if (length(cand) == 0) {
return(tibble::tibble(truth = k, alarm = NA_integer_,
delay = NA_real_, detected = FALSE))
}
j <- cand[1]
used[j] <<- TRUE
tibble::tibble(truth = k, alarm = as.integer(al$time[j]),
delay = as.numeric(al$time[j] - k), detected = TRUE)
})
per_change <- do.call(rbind, rows)
false_alarms <- al[!used, , drop = FALSE]
structure(
list(
per_change = per_change,
false_alarms = false_alarms,
mean_delay = mean(per_change$delay, na.rm = TRUE),
median_delay = stats::median(per_change$delay, na.rm = TRUE),
n_false_alarms = nrow(false_alarms),
n_detected = sum(per_change$detected),
n_changes = length(truth),
# Monitored observations, not the whole stream: a replay's training
# baseline cannot raise an alarm, so counting it inflated the ARL by
# n_baseline / n_false_alarms (300 monitored observations with three
# false alarms read as 133, not 100), in the direction that flatters
# the calibration the ARL exists to check.
arl = if (nrow(false_alarms) > 0) {
n_monitored / nrow(false_alarms)
} else {
Inf
},
n_obs = n_obs
),
class = "ggcpt_delay"
)
}
#' @rdname cpt_delay
#' @export
tidy.ggcpt_delay <- function(x, ...) {
x$per_change
}
#' @rdname cpt_delay
#' @export
glance.ggcpt_delay <- function(x, ...) {
tibble::tibble(n_changes = x$n_changes, n_detected = x$n_detected,
mean_delay = x$mean_delay, median_delay = x$median_delay,
n_false_alarms = x$n_false_alarms, arl = x$arl,
n_obs = x$n_obs)
}
#' @rdname cpt_delay
#' @param x A \code{ggcpt_delay} object.
#' @param ... Ignored.
#' @export
print.ggcpt_delay <- function(x, ...) {
cat("ggcpt_delay\n")
cat_field("True changes", x$n_changes)
cat_field("Detected", x$n_detected)
cat_field("Mean delay",
if (is.nan(x$mean_delay)) "-" else format(x$mean_delay, digits = 4))
cat_field("Median delay",
if (is.na(x$median_delay)) "-" else
format(x$median_delay, digits = 4))
cat_field("False alarms", x$n_false_alarms)
cat_field("Average run length",
if (is.infinite(x$arl)) "no false alarms" else
format(x$arl, digits = 4))
cat("\n")
print(x$per_change, n = 10)
invisible(x)
}
#' @rdname cpt_delay
#' @param object A \code{ggcpt_delay} object (for \code{autoplot()}).
#' @export
autoplot.ggcpt_delay <- function(object, ...) {
d <- object$per_change
if (all(!d$detected)) {
cpt_abort("No change was detected, so there are no delays to draw.",
class = "capability_absent")
}
ggplot2::ggplot(d, ggplot2::aes(factor(truth), delay)) +
ggplot2::geom_col(ggplot2::aes(fill = detected), width = 0.6,
na.rm = TRUE) +
ggplot2::scale_fill_manual(values = c(`TRUE` = "#0072B2",
`FALSE` = "grey70"),
name = "Detected") +
ggplot2::labs(x = "True changepoint", y = "Detection delay (observations)",
title = "Detection delay",
subtitle = paste0(object$n_false_alarms,
" false alarm(s); average run length ",
if (is.infinite(object$arl)) {
"infinite"
} else {
format(object$arl, digits = 4)
}))
}
# Internal: the one moment the streaming API can test its independence
# assumption before a wrong answer exists. A monitor is handed an in-control
# baseline before any decision; measured (ยง206), under AR(1) noise with
# rho = 0.7 the e-detector raised 3.9 times and cpm 10 times their declared
# false-alarm rates. Warned once, at construction, rather than after the
# pager has been muted.
#' @noRd
check_baseline_dependence <- function(baseline, method) {
B <- if (is.matrix(baseline) || is.data.frame(baseline)) {
as.matrix(baseline)
} else {
matrix(as.numeric(baseline), ncol = 1)
}
if (nrow(B) < 20L) return(invisible(NULL))
deps <- lapply(seq_len(ncol(B)), function(j) {
residual_dependence(B[, j], integer(0))
})
p <- vapply(deps, function(d) if (is.null(d)) NA_real_ else d$p_value,
numeric(1))
acf1 <- vapply(deps, function(d) if (is.null(d)) NA_real_ else d$acf1,
numeric(1))
worst <- which.min(p)
if (!length(worst) || is.na(p[worst]) || p[worst] >= 0.05) {
return(invisible(NULL))
}
cpt_warn("`baseline` is autocorrelated (Ljung-Box p = ",
format(signif(p[worst], 2)), ", lag-1 autocorrelation ",
format(round(acf1[worst], 2)),
if (ncol(B) > 1L) paste0(", coordinate ", worst) else "",
"). The `", method, "` monitor's false-alarm rate assumes ",
"independent observations; measured under AR(1) noise with ",
"rho = 0.7, the monitors raised four to ten times their declared ",
"rate. Expect alarms sooner than `arl0` or `alpha` says, or model ",
"the dependence before monitoring (for example, monitor the ",
"residuals of an AR fit).", class = "dependence",
data = list(p_value = p[worst], acf1 = acf1[worst]))
invisible(NULL)
}
# Internal: the control chart of a monitor: its statistic at every
# observation against its limits, alarms flagged.
#' @noRd
control_chart <- function(object) {
tr <- object$trace
if (is.null(tr) || !nrow(tr)) {
cpt_abort("This monitor has consumed no observations yet (or was built ",
"before 0.7.0, which recorded no trace), so there is no chart ",
"to draw. Feed it with cpt_update().",
class = "capability_absent")
}
st <- object$state
al <- object$alarms
lim <- switch(object$method,
cusum = c(-st$h, st$h),
ewma = c(-st$threshold, st$threshold),
edetector = c(NA_real_, st$threshold),
ocd = c(NA_real_, 1),
c(NA_real_, NA_real_)
)
ylab <- switch(object$method,
cusum = "CUSUM (upper, lower)",
ewma = "EWMA (standardised)",
edetector = "E-detector",
"Statistic"
)
p <- ggplot2::ggplot(tr, ggplot2::aes(x = time, y = statistic)) +
ggplot2::geom_line(colour = "#0072B2", na.rm = TRUE)
if (any(!is.na(tr$lower))) {
p <- p + ggplot2::geom_line(ggplot2::aes(y = lower), colour = "#0072B2",
na.rm = TRUE)
}
lims <- lim[!is.na(lim)]
if (length(lims)) {
p <- p + ggplot2::geom_hline(yintercept = lims, linetype = "dashed",
colour = "#D55E00")
}
if (nrow(al)) {
hit <- tr[tr$time %in% al$time, , drop = FALSE]
yv <- if (object$method == "cusum") {
ifelse(abs(hit$lower) > hit$statistic, hit$lower, hit$statistic)
} else hit$statistic
p <- p + ggplot2::geom_point(data = tibble::tibble(time = hit$time,
statistic = yv),
colour = "#D55E00", size = 2.2,
na.rm = TRUE)
}
p + ggplot2::labs(x = "Observation", y = ylab,
title = paste0("Control chart (", object$method, ")"),
subtitle = paste0(nrow(al), " alarm(s) in ", object$t,
" observations; dashed: control ",
"limits")) +
theme_ggcpt()
}
Try the ggchangepoint package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.