Package {focus}

Online Changepoint Detection in Univariate and Multivariate Data Streams

Description

The focus package provides high-performance implementations of online and offline changepoint detection algorithms for univariate and multivariate data streams. The package exposes optimized C++ backends (FOCuS and md-FOCuS style algorithms) with R interfaces for both real-time monitoring (sequential updates) and batch/offline analysis.

Details

focus implements efficient changepoint detection for a range of statistical models and use-cases. The package supports multiple distributional families (Gaussian, Poisson, Bernoulli, Gamma) and a non-parametric NPFOCuS variant. It provides two primary modes:

Offline mode (focus_offline)

Processes all observations in C++ for maximum efficiency. Useful for benchmarking, for computing full statistic trajectories, or for batch processing. By default the offline call stops when a threshold is exceeded; use threshold = Inf to compute statistics for all observations.

Online / sequential mode (detector_create, detector_update, get_statistics)

Create an online detector object and update it one observation at a time from R. This provides a flexible streaming API and access to candidate segments, at the cost of more R/C calls per observation.

Main features:

• Multiple distributions: Gaussian, Poisson, Bernoulli, Gamma, and non-parametric NPFOCuS.

• Univariate and multivariate observations.

• Known or unknown pre-change parameters (generalised likelihood-ratio and Page–CUSUM style tests).

• One-sided and two-sided detection.

• Flexible interface: statistical cost functions are independent from the detector candidate-management strategy.

• High-performance C++ backend designed for integration with R and direct use from C++ projects.

Important functions:

• focus_offline(Y, threshold, type, family, ...): run complete offline detection; returns full statistic trajectories, detected changepoints, candidate segments, thresholds and metadata.

• detector_create(type, ...): construct an online detector object.

• detector_update(det, y): update the detector with a new observation.

• get_statistics(det, family, theta0 = NULL, shape = NULL): compute statistics for the current detector state.

• Inspection utilities: detector_info_n(), detector_info_sn(), detector_cands_len(), detector_candidates().

The detector, the statistics, the offline results and the projection index sets are S3 objects of class "focus_detector", "focus_statistics", "focus_offline" and "focus_projections", respectively, with print methods and, for detectors and offline results, summary methods; offline results also have a plot method. See focus-methods.

Author(s)

Gaetano Romano [aut, cre], g.romano@lancaster.ac.uk\ Kes Ward [aut], k.ward4@lancaster.ac.uk\ Yuntang Fan [aut], y.yuntang@lancaster.ac.uk\ Guillem Rigaill [aut], guillem.rigaill@inrae.fr\ Vincent Runge [aut], vincent.runge@univ-evry.fr\ Paul Fearnhead [aut], p.fearnhead@lancaster.ac.uk\ Idris A. Eckley [aut], i.eckley@lancaster.ac.uk\ Maintainer: Gaetano Romano g.romano@lancaster.ac.uk

References

Pishchagina, L., G. Romano, P. Fearnhead, V. Runge, and G. Rigaill (2025). Online Multivariate Changepoint Detection: Leveraging Links with Computational Geometry. JRSS B. doi:10.1093/jrsssb/qkaf046.

Romano, G., I. A. Eckley, and P. Fearnhead (2024). A Log-Linear Nonparametric Online Changepoint Detection Algorithm Based on Functional Pruning. IEEE Transactions on Signal Processing, 72.

Romano, G., I. A. Eckley, P. Fearnhead, and G. Rigaill (2023). Fast Online Changepoint Detection via Functional Pruning CUSUM Statistics. Journal of Machine Learning Research, 24(81):1–36.

Ward, K., G. Romano, I. Eckley, and P. Fearnhead (2024). A Constant-Per-Iteration Likelihood Ratio Test for Online Changepoint Detection for Exponential Family Models. Statistics and Computing, 34(3).

Note: The package bundles and links to code from Qhull (C. B. Barber et al., The Geometry Center). Qhull functions are used for convex-hull computations in multivariate pruning. See ‘inst/COPYRIGHTS/’ for the Qhull license and attribution details.

See Also

focus_offline, detector_create, detector_update, get_statistics, detector_candidates, generate_projection_indexes

Examples

  ## Univariate Gaussian change in mean
  set.seed(123)
  Y <- c(rnorm(500, mean = 0), rnorm(500, mean = 2))

  ## Online (sequential) example
  det <- detector_create(type = "univariate")
  stat_trace <- numeric(length(Y))
  threshold <- 20
  for (i in seq_along(Y)) {
    detector_update(det, Y[i])
    r <- get_statistics(det, family = "gaussian")
    stat_trace[i] <- r$stat
    if (!is.null(r$stat) && r$stat > threshold) {
      cat("Online detection at", i, "estimate tau =", r$changepoint, "\n")
      break
    }
  }
  summary(det, family = "gaussian")

  ## Offline (batch) example on the same data
  # Run offline detection (C++ loop). Use threshold=Inf to compute full trajectory.
  res <- focus_offline(Y, threshold = 20, type = "univariate", family = "gaussian")
  if (!is.null(res$detection_time)) {
    cat("Detection at time:", res$detection_time, "\n")
  }
  summary(res)
  plot(res, data = Y)

  ## Multivariate detection
  set.seed(123)
  d <- 3
  Y_multi <- rbind(
    matrix(rnorm(300 * d, mean = 0), ncol = d),
    matrix(rnorm(100 * d, mean = 1), ncol = d)
  )
  det_mv <- detector_create(type = "multivariate")
  for (i in seq_len(nrow(Y_multi))) {
    detector_update(det_mv, Y_multi[i, ])
    r <- get_statistics(det_mv, family = "gaussian")
    if (!is.null(r$stat) && r$stat > 30) {
      cat("Multivariate detection at", i, "estimate tau =", r$changepoint, "\n")
      break
    }
  }
  summary(det_mv)

  ## Non-parametric detection (NP-FOCuS)
  # Heavy-tailed data with a change in location
  set.seed(123)
  Y_np <- c(rt(500, df = 2), rt(200, df = 2) + 1)
  # Quantiles of the pre-change distribution
  quants <- qt(c(0.1, 0.25, 0.5, 0.75, 0.9), df = 2)
  det_np <- detector_create(type = "npfocus", quantiles = quants)
  threshold_np <- c(50, 20)  # for the sum and the max statistics
  for (i in seq_along(Y_np)) {
    detector_update(det_np, Y_np[i])
    r <- get_statistics(det_np, family = "npfocus")
    if (!is.null(r$stat) && any(r$stat > threshold_np)) {
      cat("Non-parametric detection at", i, "\n")
      break
    }
  }
  r

  ## Autoregressive detection: AR(3) with a change in mean
  set.seed(123)
  rho <- c(0.5, -0.3, 0.2)  # AR coefficients
  Y_ar <- arima.sim(list(ar = rho), n = 700) + rep(c(0, 1.5), c(500, 200))
  det_ar <- detector_create(type = "arp", rho = rho, mu0_arp = 0)
  for (i in seq_along(Y_ar)) {
    detector_update(det_ar, Y_ar[i])
    r <- get_statistics(det_ar, family = "arp")
    if (!is.null(r$stat) && r$stat > 30) {
      cat("AR(3) detection at", i, "estimate tau =", r$changepoint, "\n")
      break
    }
  }
  det_ar

Get the Candidate Segments

Description

Returns detailed information about all candidate changepoint segments currently tracked by the detector.

Usage

detector_candidates(det_ptr, ...)

## S3 method for class 'focus_detector'
detector_candidates(det_ptr, ...)

Arguments

Details

Each row represents a candidate segment from time tau to the current time. The sufficient statistics in st are used to efficiently compute test statistics without reprocessing the data. detector_candidates() is an S3 generic, with a method for "focus_detector" objects.

Value

A data frame (tibble) with one row per candidate and columns:

Examples

set.seed(1)
det <- detector_create(type = "univariate")
for (y in rnorm(50)) detector_update(det, y)
head(detector_candidates(det))

Get the Number of Candidate Segments

Description

Returns the number of candidate changepoint segments currently tracked by the detector.

Usage

detector_cands_len(det_ptr, ...)

## S3 method for class 'focus_detector'
detector_cands_len(det_ptr, ...)

Arguments

Details

The FOCuS algorithm maintains a set of candidate segments that could potentially contain changepoints. This number grows with time but is controlled by the pruning parameters. detector_cands_len() is an S3 generic, with a method for "focus_detector" objects.

Value

Integer. Number of candidate segments.

Examples

set.seed(1)
det <- detector_create(type = "univariate")
for (y in rnorm(50)) detector_update(det, y)
detector_cands_len(det)

Create a FOCuS Changepoint Detector

Description

Creates an online (sequential) changepoint detector object that provides a step-by-step interface to the FOCuS algorithm. Each call to detector_update() adds new data, and get_statistics() computes the current test statistic and detection result.

Usage

detector_create(
  type,
  dim_indexes = NULL,
  quantiles = NULL,
  pruning_mult = 2L,
  pruning_offset = 1L,
  side = "right",
  anomaly_intensity = NULL,
  rho = NULL,
  mu0_arp = NULL
)

Arguments

Details

The detector maintains sufficient statistics internally and uses pruning to efficiently track candidate changepoints. The pruning_mult and pruning_offset parameters control the pruning strategy.

AutoRegressive Process (ARP). When type = "arp", the rho parameter must be provided as a numeric vector of AR coefficients (lag-1, lag-2, ..., lag-p). The detector then computes statistics optimal for detecting changepoints in AR(p) processes. Use get_statistics(family = "arp") to retrieve the test statistics. The optional mu0_arp parameter specifies the pre-change mean and is tied to the pruning logic: if provided, it enables more efficient pruning by allowing the algorithm to filter candidates based on the known pre-change parameter.

High-dimensional multivariate detectors. For high-dimensional multivariate detection, computing the full convex hull is prohibitive, as the expected number of candidates grows as \log(n)^p, where n is the number of observations and p the number of dimensions. The set of candidate changepoints can then be approximated by computing the hull on lower-dimensional projections: dim_indexes specifies which dimensions to use for each projection. Use generate_projection_indexes() to generate systematic projection sets.

NPFOCuS. For non-parametric detection, create the detector with detector_create(type = "npfocus", quantiles = ...) and compute the statistics with get_statistics(family = "npfocus"). No other family can be used with this detector type.

Value

An object of class "focus_detector": an external pointer to the state of the C++ detector, which is updated in place. It should be passed to the other detector functions, such as detector_update() and get_statistics(). print and summary methods are available, see focus-methods. As external pointers, detectors cannot be saved and restored across R sessions.

Examples

# Univariate detector
det <- detector_create(type = "univariate")
det <- detector_update(det, 0.5)
det <- detector_update(det, 1.2)
det
r <- get_statistics(det, family = "gaussian")
r

## Online (sequential) example
# Generate data with a changepoint
set.seed(123)
Y <- c(rnorm(500, mean = 0), rnorm(500, mean = 1))
det <- detector_create(type = "univariate")
stat_trace <- numeric(length(Y))
threshold <- 20
for (i in seq_along(Y)) {
  detector_update(det, Y[i])
  r <- get_statistics(det, family = "gaussian")
  stat_trace[i] <- r$stat
  if (!is.null(r$stat) && r$stat > threshold) {
    cat("Online detection at", i, "estimate tau =", r$changepoint, "\n")
    plot(stat_trace[1:i], type = "l", ylab = "Test Statistic", xlab = "Time")
    break
  }
}

# Multivariate detector with projections
dim_indexes <- list(c(0,1), c(1,2), c(0,2))  # 0-based indices
det_mv <- detector_create(type = "multivariate", dim_indexes = dim_indexes)
detector_update(det_mv, c(0.5, 1.2, -0.3))

# Nonparametric detector
quants <- qnorm(c(0.25, 0.5, 0.75))
det_np <- detector_create(type = "npfocus", quantiles = quants)

# One-sided univariate detector
det_one_sided <- detector_create(type = "univariate_one_sided", side = "left")

Get the Number of Observations Processed

Description

Returns the total number of observations processed by the detector.

Usage

detector_info_n(det_ptr, ...)

## S3 method for class 'focus_detector'
detector_info_n(det_ptr, ...)

Arguments

Details

detector_info_n() is an S3 generic, with a method for "focus_detector" objects.

Value

Integer. Number of observations processed (current time index).

Examples

det <- detector_create(type = "univariate")
for (y in c(0.5, 1.2, -0.3)) detector_update(det, y)
detector_info_n(det)

Get the Cumulative Sum Statistic

Description

Returns the current cumulative sum statistic maintained by the detector.

Usage

detector_info_sn(det_ptr, ...)

## S3 method for class 'focus_detector'
detector_info_sn(det_ptr, ...)

Arguments

Details

detector_info_sn() is an S3 generic, with a method for "focus_detector" objects.

Value

Numeric vector. Cumulative sum statistic. For univariate detectors, a scalar (length-1 vector). For multivariate detectors, a vector of length equal to the number of dimensions.

Examples

det <- detector_create(type = "multivariate")
detector_update(det, c(0.5, 1.2))
detector_update(det, c(-0.3, 0.4))
detector_info_sn(det)

Update a Detector with New Observations

Description

Adds a new observation to the detector's internal state, updates the sufficient statistics and prunes the set of candidate changepoints.

Usage

detector_update(det_ptr, y, lambda = 1)

Arguments

Value

The same "focus_detector" object that was passed in. The detector is updated in place (no copy is made), so the return value is provided only for convenience, for example to chain calls with the native pipe operator, as in det |> detector_update(y) |> get_statistics(family = "gaussian").

Examples

# Univariate example
det <- detector_create(type = "univariate")
detector_update(det, 0.5)
detector_update(det, 1.2)

# Multivariate example
det_mv <- detector_create(type = "multivariate")
detector_update(det_mv, c(0.5, 1.2, -0.3))

## Online (sequential) example
# Generate data with a changepoint
set.seed(123)
Y <- c(rnorm(500, mean = 0), rnorm(500, mean = 1))
det <- detector_create(type = "univariate")
stat_trace <- numeric(length(Y))
threshold <- 20
for (i in seq_along(Y)) {
  detector_update(det, Y[i])
  r <- get_statistics(det, family = "gaussian")
  stat_trace[i] <- r$stat
  if (!is.null(r$stat) && r$stat > threshold) {
    cat("Online detection at", i, "estimate tau =", r$changepoint, "\n")
    plot(stat_trace[1:i], type = "l", ylab = "Test Statistic", xlab = "Time")
    break
  }
}

Methods for Detector, Statistics, Offline Result and Projection Objects

Description

Methods for the S3 classes returned by the main functions of the package: "focus_detector" objects, created by detector_create(); "focus_statistics" objects, returned by get_statistics(); "focus_offline" objects, returned by focus_offline(); and "focus_projections" objects, returned by generate_projection_indexes().

Usage

## S3 method for class 'focus_detector'
print(x, ...)

## S3 method for class 'focus_statistics'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'focus_offline'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'focus_offline'
summary(object, ...)

## S3 method for class 'summary.focus_offline'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'focus_offline'
plot(
  x,
  data = NULL,
  type = "l",
  lty = 1,
  col = NULL,
  xlab = "Time",
  ylab = "Statistic",
  main = NULL,
  ...
)

## S3 method for class 'focus_detector'
summary(object, family = NULL, theta0 = NULL, shape = NULL, ...)

## S3 method for class 'summary.focus_detector'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'focus_projections'
print(x, n = 10L, ...)

## S3 method for class 'focus_projections'
x[i]

## S3 method for class 'focus_projections'
as.matrix(x, ...)

## S3 method for class 'focus_offline'
autoplot(object, data = NULL, ...)

Arguments

Details

The print methods give a compact description of the object: for a detector, its type, the number of observations processed and the number of candidate changepoints currently stored; for the statistics, the current time, changepoint estimate and test statistic(s); for an offline result, the detector type and family, the threshold(s) and the detection, if any; for projection index sets, the (0-based) indices of each projection, one projection per row.

The summary method for "focus_detector" objects additionally reports the cumulative sums (see detector_info_sn()) and the locations of the candidate changepoints, at the time it is called, and, if a family is given, the current statistics.

The summary method for "focus_offline" objects additionally reports, for each test statistic, its maximum over time, the time at which the maximum is attained, the changepoint estimate at that time and the threshold.

The plot method for "focus_offline" objects draws the trace of the test statistic(s) over time, with the finite threshold(s) as dashed horizontal lines and, if a detection occurred, the estimated changepoint as a dotted vertical line. If data are given, they are drawn above the trace, with one panel per dimension. If ggplot2 is installed, autoplot() returns the same plot as a "ggplot" object.

Projection index sets can be subset with [ (e.g., with head()), keeping their class, and as.matrix returns them as a matrix with one row per projection.

Value

The print and plot methods return x invisibly. The summary method for detectors returns an object of class "summary.focus_detector", a list with the reported elements. The summary method for offline results returns an object of class "summary.focus_offline": a list with the elements type, family, shape, n, detection_time and detected_changepoint of the offline result, the number of final candidates n_candidates, and statistics, a data frame with one row per test statistic. The [ method returns a "focus_projections" object and as.matrix an integer matrix.

Examples

set.seed(123)
Y <- c(rnorm(100, mean = 0), rnorm(100, mean = 2))

# Online interface
det <- detector_create(type = "univariate")
for (y in Y[1:120]) detector_update(det, y)
det
summary(det, family = "gaussian")
get_statistics(det, family = "gaussian")

# Offline interface
res <- focus_offline(Y, threshold = 20, type = "univariate",
                     family = "gaussian")
res
summary(res)
plot(res)
plot(res, data = Y)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::autoplot(res, data = Y)
}

# Projection index sets
proj <- generate_projection_indexes(d = 6, p = 2)
proj
head(proj, 2)
as.matrix(proj)

Run a FOCuS Detector in Offline Batch Mode

Description

Processes all data at once and returns detection results and trajectories. This is the most efficient way to run changepoint detection when all data is available upfront.

Usage

focus_offline(
  Y,
  threshold,
  type = "univariate",
  family = "gaussian",
  theta0 = NULL,
  dim_indexes = NULL,
  quantiles = NULL,
  pruning_mult = 2L,
  pruning_offset = 1L,
  side = "right",
  shape = NULL,
  anomaly_intensity = NULL,
  rho = NULL,
  mu0_arp = NULL
)

Arguments

Details

This function runs the complete detection algorithm in C++ for maximum efficiency. It processes observations sequentially and stops at the first detection (when any statistic exceeds its threshold).

When more than one statistic is computed (e.g., the sum and maximum statistics of "npfocus"), a detection occurs when any statistic exceeds its threshold.

Value

An object of class "focus_offline", with print, summary and plot methods (see focus-methods): a list with components

Examples

# Univariate Gaussian detection
set.seed(123)
Y <- c(rnorm(100, mean = 0), rnorm(100, mean = 2))
result <- focus_offline(Y, threshold = 10, type = "univariate",
                       family = "gaussian")
result
summary(result)

# Plot the trace of the statistic, the threshold and the detection
plot(result)

# Plot the data above the trace of the statistic
plot(result, data = Y)

# Poisson detection
Y_poisson <- c(rpois(100, lambda = 2), rpois(100, lambda = 5))
result_poisson <- focus_offline(Y_poisson, threshold = 10,
                                type = "univariate",
                                family = "poisson")

Generate Projection Index Sets

Description

Generates projection index sets for high-dimensional multivariate detectors using circular combinations.

Usage

generate_projection_indexes(d, p)

Arguments

Details

This function generates systematic projection sets for use with multivariate detectors. The circular combination approach ensures good coverage of the dimensional space while keeping the number of projections manageable.

Value

An object of class "focus_projections": a list of integer vectors, each holding the 0-based column indices of one projection, with the number of dimensions and the projection size stored in the attributes "d" and "p". It can be passed directly as the dim_indexes argument of detector_create() and focus_offline(). A print method shows one projection per row, subsetting with [ (e.g., head()) keeps the class, and as.matrix() returns the indices as a matrix with one row per projection; see focus-methods.

Examples

# Generate 2-dimensional projections from 5 dimensions
proj <- generate_projection_indexes(d = 5, p = 2)
proj
as.matrix(proj)

# Use with multivariate detector
det <- detector_create(type = "multivariate", dim_indexes = proj)
detector_update(det, c(0.5, 1.2, -0.3, 0.8, 0.1))

# Offline comparison of the full detector and approximation
set.seed(42)
d <- 6

# Create data: changepoint at t=1000
Y_multi <- rbind(
    matrix(rnorm(1000 * d, mean = -1, 1), ncol = d),
    matrix(rnorm(500 * d, mean = 1.2), ncol = d)
)

# Full multivariate detection
system.time(
res_multi <- focus_offline(Y_multi, threshold = Inf,
                           type = "multivariate", family = "gaussian")
)

# Low-dimensional projection approximation
dim_indexes <- generate_projection_indexes(6, 2)
system.time(
res_multi_approx <- focus_offline(Y_multi, threshold = Inf,
                                  type = "multivariate", family = "gaussian",
                                  dim_indexes = dim_indexes)
)

# Verify similarity
all.equal(res_multi$stat, res_multi_approx$stat)

Compute the Current Changepoint Statistics

Description

Computes the current changepoint test statistic and detection result based on all observations processed so far.

Usage

get_statistics(det_ptr, family, theta0 = NULL, shape = NULL)

Arguments

Details

The function computes a log-likelihood ratio test statistic comparing the null hypothesis (no change) against the alternative (a change at the optimal location). The statistic is typically compared against a threshold to determine if a changepoint should be declared.

Gamma family. When family = "gamma" a positive shape parameter must be provided; otherwise an error is raised. Passing shape for a non-gamma family raises a warning and the parameter is ignored.

NPFOCuS. For non-parametric detection, the detector must be created with detector_create(type = "npfocus", quantiles = ...). NPFOCuS returns two statistics (sum and maximum over the quantiles) as a vector; in the offline interface, stat is a matrix with two columns.

AutoRegressive Process (ARP). For ARP detection, use family = "arp" with a detector created via detector_create(type = "arp", rho = ...). The AR coefficients and the (optional) pre-change mean mu0_arp are set at detector creation, so theta0 is ignored (with a warning) for this family.

Value

An object of class "focus_statistics", with a print method (see focus-methods): a list with components

The family is stored in the attribute "family".

Examples

## Online (sequential) example
# Generate data with a changepoint
set.seed(123)
Y <- c(rnorm(500, mean = 0), rnorm(500, mean = 1))
det <- detector_create(type = "univariate")
stat_trace <- numeric(length(Y))
threshold <- 20
for (i in seq_along(Y)) {
  detector_update(det, Y[i])
  r <- get_statistics(det, family = "gaussian")
  stat_trace[i] <- r$stat
  if (!is.null(r$stat) && r$stat > threshold) {
    cat("Online detection at", i, "estimate tau =", r$changepoint, "\n")
    plot(stat_trace[1:i], type = "l", ylab = "Test Statistic", xlab = "Time")
    break
  }
}

## Note that multiple models can be tested simultaneously on the same detector
# as the statistics is independent of the detector state.
# For example, testing both Gaussian and Poisson costs
set.seed(2024)
# Generate Poisson count data with a rate change
Y_counts <- c(rpois(500, lambda = 10), rpois(500, lambda = 15))
# Compute full trajectories for comparison
det2 <- detector_create(type = "univariate")
stat_gaussian <- numeric(length(Y_counts))
stat_poisson <- numeric(length(Y_counts))

for (i in seq_along(Y_counts)) {
  detector_update(det2, Y_counts[i])
  stat_gaussian[i] <- get_statistics(det2, family = "gaussian")$stat
  stat_poisson[i] <- get_statistics(det2, family = "poisson", theta0 = 10)$stat
}

# Plot comparison
oldpar <- par(mfrow = c(1, 2))
plot(stat_gaussian, type = "l", main = "Gaussian Statistic on Poisson Data",
     xlab = "Time", ylab = "Statistic", lwd = 2, col = "blue")
abline(v = 500, col = "green", lty = 3, lwd = 2)

plot(stat_poisson, type = "l", main = "Poisson Statistic on Poisson Data",
     xlab = "Time", ylab = "Statistic", lwd = 2, col = "red")
abline(v = 500, col = "green", lty = 3, lwd = 2)
par(oldpar)