| Type: | Package |
| Title: | Object-Oriented Interface for Offline Change-Point Detection |
| Version: | 2.0.0 |
| Description: | A collection of efficient implementations of popular offline change-point detection algorithms, featuring a consistent, object-oriented interface for practical use. |
| Encoding: | UTF-8 |
| URL: | https://edelweiss611428.github.io/rupturesRcpp/, https://github.com/edelweiss611428/rupturesRcpp |
| BugReports: | https://github.com/edelweiss611428/rupturesRcpp/issues |
| RoxygenNote: | 7.3.3 |
| License: | MIT + file LICENSE |
| LinkingTo: | Rcpp, RcppArmadillo |
| Imports: | Rcpp, R6, ggplot2, patchwork, methods |
| Suggests: | testthat (≥ 3.0.0), reticulate, binsegRcpp, covr |
| Config/testthat/edition: | 3 |
| Collate: | 'costFuncR6.R' 'DynpR6.R' 'PeltR6.R' 'WindowR6.R' 'binSegR6.R' 'costFactoryR6.R' 'rupturesRcpp-package.R' 'zzz.R' |
| NeedsCompilation: | yes |
| Packaged: | 2026-10-11 14:17:04 UTC; edelweiss |
| Author: | Minh Long Nguyen [aut, cre], Toby Hocking [aut], Charles Truong [aut], Huy Nhat Minh Nguyen [ctb] |
| Maintainer: | Minh Long Nguyen <edelweiss611428@gmail.com> |
| Repository: | CRAN |
| Date/Publication: | 2026-10-11 15:00:02 UTC |
Exact Dynamic Programming (Dynp)
Description
An R6 class implementing exact dynamic programming for offline change-point detection.
Details
Dynp finds the segmentation that globally minimises the total cost for a specified number of
change-points – unlike PELT (exact for a penalty, not a change-point count) and binSeg
(greedy at each split, not guaranteed globally optimal for a fixed count), Dynp is exact for
whatever nBkps is requested. The price is complexity: building the full table costs
O(\text{nBkpsMax} \cdot M^2) where M is the number of (minSize, jump)-admissible
grid points (M \approx n/\text{jump} in the worst case), versus PELT's near-linear pruning
or binSeg's O(n \log n)-ish greedy search. $fit() computes the exact minimal cost for
every change-point count from 0 to nBkpsMax in one pass (see $costPath()). Dynp is
therefore useful when an exact fixed-number-of-change-points solution is required,
particularly when using custom cost functions for which PELT pruning cannot be guaranteed.
Dynp requires a R6 object of class costFunc, exactly as PELT/binSeg/Window do – see
costFunc for the supported cost functions ("L1", "L2", "SIGMA", "VAR", "LinearL2",
"LinearSIGMA", "LinearL1", "Custom").
Some examples are provided below. See the package website for detailed usage!
Methods
$new()Initialises a
Dynpobject.$describe()Describes the
Dynpobject.$fit()Constructs a
Dynpmodule inC++.$eval()Evaluates the cost of a segment.
$predict()Returns the exact optimal breakpoints, for a given
penornBkps.$segments()Returns the cost and parameter estimates of each segment from the latest
$predict().$costPath()Returns the exact minimal cost for every change-point count up to
nBkpsMax.$getHistory()Same data as
$costPath(), as ak/costdata.frame(same shape asbinSeg/Window's$getHistory(), minusadded_bkp– see its own docs for why).$plotElbow()Plots the elbow curve (Total Cost vs. Number of Change-Points).
$plot()Plots change-point segmentation in
ggplotstyle.$clone()Clones the
R6object.
Active bindings
minSizeInteger. Minimum allowed segment length. Can be accessed or modified via
$minSize. ModifyingminSizewill automatically trigger$fit().jumpInteger. Search grid step size. Can be accessed or modified via
$jump. Modifyingjumpwill automatically trigger$fit().nBkpsMaxInteger or
NULL. Upper bound on the number of change-points the exact DP table is built for;$predict(nBkps = ...)treats this as a cap, not a requirement – see$predict(). IfNULL(default),$fit()resolves it tomin(20, floor(n / minSize) - 1), i.e. the maximum feasible value givenminSize, capped at 20. Set explicitly to go beyond 20 (up to the maximum feasible value) or to lower it; the table costsO(\text{nBkpsMax} \cdot M^2)work andO(\text{nBkpsMax} \cdot M)memory. Can be accessed or modified via$nBkpsMax; modifying it will automatically trigger$fit().costFuncR6object of classcostFunc. Can be accessed or modified via$costFunc. ModifyingcostFuncwill automatically trigger$fit().tsMatNumeric matrix. Input time series matrix of size
n \times p. Can be accessed or modified via$tsMat. ModifyingtsMatwill automatically trigger$fit().covariatesNumeric matrix. Input time series matrix having a similar number of observations as
tsMat. Can be accessed or modified via$covariates. Modifyingcovariateswill automatically trigger$fit().
Methods
Public methods
Method new()
Initialises a Dynp object.
Usage
Dynp$new(minSize, jump, nBkpsMax, costFunc)
Arguments
minSizeInteger. Minimum allowed segment length. Default:
1L.jumpInteger. Search grid step size: only positions in {k, 2k, ...} are considered. Default:
1L.nBkpsMaxInteger or
NULL. Upper bound on the number of change-points to build the exact DP table for. Default:NULL(resolved tomin(20, floor(n / minSize) - 1)at$fit()time).costFuncA
R6object of classcostFunc. Should be created viacostFunc$new()to avoid error. Default:costFunc$new("L2").
Returns
Invisibly returns NULL.
Method describe()
Describes a Dynp object.
Usage
Dynp$describe(printConfig = FALSE)
Arguments
printConfigLogical. Whether to print object configurations. Default:
FALSE.
Returns
Invisibly returns a list storing at least the following fields:
minSizeMinimum allowed segment length.
jumpSearch grid step size.
nBkpsMaxThe user-set
nBkpsMax(possiblyNULL).resolvedNBkpsMaxThe
nBkpsMaxactually used by the last$fit()(NULLif not fitted).costFuncThe
costFuncobject.fittedWhether or not
$fit()has been run.tsMatTime series matrix.
covariatesCovariate matrix (if exists).
nNumber of observations.
pNumber of features.
Method fit()
Constructs a C++ module for Dynp and builds the exact DP table.
Usage
Dynp$fit(tsMat = NULL, covariates = NULL)
Arguments
tsMatNumeric matrix. A time series matrix of size
n \times pwhose rows are observations ordered in time. IftsMat = NULL, the method will use the previously assignedtsMat(e.g., set via the active binding$tsMator from a prior$fit(tsMat)). Default:NULL.covariatesNumeric matrix. A time series matrix having a similar number of observations as
tsMat. Required for models involving both dependent and independent variables. Ifcovariates = NULLand no prior covariates were set (i.e.,$covariatesis stillNULL), the model is force-fitted with only an intercept. Default:NULL.
Details
This method constructs a C++ Dynp module and sets private$.fitted to TRUE, enabling the
use of $predict(), $costPath() and $eval(). If $nBkpsMax is NULL, it is resolved here to
min(20, floor(n / minSize) - 1) (with a message); if set higher than floor(n / minSize) - 1,
it is capped to that (with a warning).
Returns
Invisibly returns NULL.
Method eval()
Evaluate the cost of the segment (a,b]
Usage
Dynp$eval(a, b)
Arguments
aInteger. Start index of the segment (exclusive). Must satisfy
start < end.bInteger. End index of the segment (inclusive).
Returns
The segment cost. See costFunc for the cost formulas.
Method predict()
Returns the exact optimal breakpoints, either under a linear penalty or for a specified number of change-points.
Usage
Dynp$predict(pen = 0, nBkps = NULL)
Arguments
penNumeric. Penalty per change-point; the change-point count is chosen by minimising
\text{cost} + \text{pen} \cdot koverk = 0, \dots, \text{nBkpsMax}(same convention asPELT/binSeg/Window). Ignored ifnBkpsis supplied. Default:0.nBkpsInteger. If supplied, takes precedence over
pen: returns the exact optimal segmentation for exactlynBkpschange-points. Treated as an upper bound, not a strict requirement, in the same spirit asbinSeg/Window: ifnBkpsexceedsnBkpsMax, it is capped tonBkpsMaxand a message reports the shortfall. UnlikebinSeg/Window, the result is still exact for whatever count is actually used – capping only ever happens because the DP table wasn't built that far (a configuration choice via$nBkpsMax), not becauseDynpran out of candidates the way a greedy search can. If no valid segmentation exists at all for the (possibly capped) count – only possible whenjumpis large relative tominSize– this still errors, since there is no well-defined smaller-but-still-what-you-asked-for answer to fall back to. Default:NULL.
Details
With nBkps = k, this returns the segmentation into exactly k change-points that globally
minimises the total cost – exact, unlike binSeg$predict() for the same k. With pen,
this instead minimises the penalised cost over every change-point count the DP table covers,
matching how PELT/binSeg/Window already select a count from a penalty.
Both the DP table (via $fit()) and the traceback here only depend on minSize/jump
admissible positions, so $predict() itself is cheap (O(nBkps)) regardless of which mode
is used – the cost was already paid in $fit().
Temporary segment end-points are saved to private$.tmpEndPoints after $predict(), enabling
users to call $plot() without specifying endpoints manually.
Returns
An integer vector of regime end-points. By design, the last element is the number of observations.
Method segments()
Returns the cost and parameter estimates of each segment from the latest $predict().
Usage
Dynp$segments()
Details
Let 0 = c_0 < c_1 < \dots < c_{k+1} = n be the end-points from the latest $predict(). Segment i is
(c_{i-1}, c_i], Cost is c_{(c_{i-1}, c_i]} and Params is its minimiser. The costs sum to
$costPath()[k + 1], the exact minimal total cost for k change-points. With nBkps, k is the
requested count (capped at the resolved nBkpsMax). With pen, k is the count, from 0 up to the resolved
nBkpsMax, that minimises the total cost plus pen times k.
Temporary end-points are cleared by $fit(), so $predict() must be run again after modifying the object via
its active bindings.
Returns
A list with one element per segment (Start, End]. Each element is a list with:
StartStart index of the segment (exclusive, 0-based), same convention as
$eval().EndEnd index of the segment (inclusive).
CostThe segment cost, as returned by
$eval(Start, End).ParamsA named list of the segment parameter estimates, e.g.
list(mean = ...)for"L2". SeecostFactoryfor the fields returned by each cost function.
Method costPath()
Returns the exact minimal cost for every change-point count the DP table covers.
Usage
Dynp$costPath()
Details
This is the data behind the "elbow method" for choosing the number of change-points:
plot it (see $getHistory()/$plotElbow()) and look for where the marginal decrease in cost
flattens out. It comes for free out of $fit() – no extra computation is triggered here.
Returns
A numeric vector of length resolvedNBkpsMax + 1: element k+1 is the exact minimal
total cost of segmenting the series into exactly k change-points, for k = 0, ..., resolvedNBkpsMax.
Inf at position k+1 means no valid segmentation with exactly k change-points exists given
minSize/jump (only possible when jump is large relative to minSize).
Method getHistory()
Retrieves the exact minimal cost for every change-point count the DP table covers,
as a data.frame (same shape as binSeg/Window's $getHistory(), minus added_bkp).
Usage
Dynp$getHistory()
Details
Unlike binSeg/Window, there is no added_bkp column here: binSeg and Window
each build a single nested sequence of change-points, where every k's answer is the
previous k-1's answer plus one more point, so "the breakpoint added at this step" is
well-defined. Dynp's per-k solutions are each independently exact and need not be nested at
all – the optimal segmentation for k change-points can differ completely from the one for
k-1. Use $predict(nBkps = k) to get the full breakpoint set for a given k.
Returns
A data.frame with two columns:
kThe number of change-points, from
0to the resolvednBkpsMax.costThe exact minimal total cost of segmenting the series into exactly
kchange-points (see$costPath()).
Method plotElbow()
Plots the elbow curve (Total Cost vs. Number of Change-Points).
Usage
Dynp$plotElbow(maxK = NULL)
Arguments
maxKInteger. The maximum number of change-points to display on the plot. If
NULL, displays the full history. Default:NULL.
Returns
A ggplot object.
Method plot()
Plots change-point segmentation
Usage
Dynp$plot(
d = 1L,
endPts,
dimNames,
main,
xlab,
tsWidth = 0.25,
tsCol = "#5B9BD5",
bgCol = c("#A3C4F3", "#FBB1BD"),
bgAlpha = 0.5,
ncol = 1L
)Arguments
dInteger vector. Dimensions to plot. Default:
1L.endPtsInteger vector. End points. Default: latest temporary changepoints obtained via
$predict().dimNamesCharacter vector. Feature names matching length of
d. Defaults to"X1", "X2", ....mainCharacter. Main title. Defaults to
"Dynp: d = ...".xlabCharacter. X-axis label. Default:
"Time".tsWidthNumeric. Line width for time series and segments. Default:
0.25.tsColCharacter. Time series color. Default:
"#5B9BD5".bgColCharacter vector. Segment colors, recycled to length of
endPts. Default:c("#A3C4F3", "#FBB1BD").bgAlphaNumeric. Background transparency. Default:
0.5.ncolInteger. Number of columns in facet layout. Default:
1L.
Details
Plots change-point segmentation results. Based on ggplot2. Multiple plots can easily be
horizontally and vertically stacked using patchwork's operators / and |, respectively.
Returns
An object of classes gg and ggplot.
Method clone()
The objects of this class are cloneable with this method.
Usage
Dynp$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.
Author(s)
Minh Long Nguyen edelweiss611428@gmail.com
Toby Dylan Hocking toby.hocking@r-project.org
Charles Truong ctruong@ens-paris-saclay.fr
Huy Nhat Minh Nguyen sleepysnorlax0115@gmail.com
References
Truong, C., Oudre, L., & Vayatis, N. (2020). Selective review of offline change point detection methods. Signal Processing, 167, 107299.
Examples
## L2 example
set.seed(1121)
signals = as.matrix(c(rnorm(100,0,1),
rnorm(100,5,1)))
# Default L2 cost function; nBkpsMax left NULL -> resolved to min(20, maximum feasible value)
DynpObj = Dynp$new(minSize = 1L, jump = 1L)
DynpObj$fit(signals)
DynpObj$predict(nBkps = 1)
DynpObj$plot()
# Same result, reached via a penalty instead of an explicit count
DynpObj$predict(pen = 100)
# The exact cost of every change-point count from 0 to nBkpsMax -- elbow-method data
DynpObj$getHistory()
DynpObj$plotElbow()
Pruned Exact Linear Time (PELT)
Description
An R6 class implementing the PELT algorithm for offline change-point detection.
Details
PELT (Pruned Exact Linear Time) is an efficient algorithm for change point detection that prunes the search space to achieve optimal segmentation in linear time under certain conditions.
PELT requires a R6 object of class costFunc, which can be created via costFunc$new(). Currently, the following cost functions are supported:
-
"L1"and"L2"for (independent) piecewise Gaussian process with constant variance -
"SIGMA": for (independent) piecewise Gaussian process with varying variance -
"VAR": for piecewise Gaussian vector-regressive process with constant noise variance -
"LinearL2": for piecewise linear regression process with constant noise variance -
"LinearSIGMA": for piecewise linear regression process with varying noise covariance -
"LinearL1": for piecewise linear regression process under L1 (least absolute deviations) loss
See $eval() method for more details on computation of cost.
Some examples are provided below. See the package website for detailed usage!
Methods
$new()Initialises a
PELTobject.$describe()Describes the
PELTobject.$fit()Constructs a
PELTmodule inC++.$eval()Evaluates the cost of a segment.
$predict()Performs
PELTgiven a linear penalty value.$segments()Returns the cost and parameter estimates of each segment from the latest
$predict().$plot()Plots change-point segmentation in
ggplotstyle.$clone()Clones the
R6object.
Active bindings
minSizeInteger. Minimum allowed segment length. Can be accessed or modified via
$minSize. ModifyingminSizewill automatically trigger$fit().jumpInteger. Search grid step size. Can be accessed or modified via
$jump. Modifyingjumpwill automatically trigger$fit().costFuncR6object of classcostFunc. Search grid step size. Can be accessed or modified via$costFunc. ModifyingcostFuncwill automatically trigger$fit().tsMatNumeric matrix. Input time series matrix of size
n \times p. Can be accessed or modified via$tsMat. ModifyingtsMatwill automatically trigger$fit().covariatesNumeric matrix. Input time series matrix having a similar number of observations as
tsMat. Can be accessed or modified via$covariates. Modifyingcovariateswill automatically trigger$fit().
Methods
Public methods
Method new()
Initialises a PELT object.
Usage
PELT$new(minSize, jump, costFunc)
Arguments
minSizeInteger. Minimum allowed segment length. Default:
1L.jumpInteger. Search grid step size: only positions in {k, 2k, ...} are considered. Default:
1L.costFuncA
R6object of classcostFunc. Should be created viacostFunc$new()to avoid error. Default:costFunc$new("L2").
Returns
Invisibly returns NULL.
Method describe()
Describes a PELT object.
Usage
PELT$describe(printConfig = FALSE)
Arguments
printConfigLogical. Whether to print object configurations. Default:
FALSE.
Returns
Invisibly returns a list storing at least the following fields:
minSizeMinimum allowed segment length.
jumpSearch grid step size.
costFuncThe
costFuncobject.fittedWhether or not
$fit()has been run.tsMatTime series matrix.
covariatesCovariate matrix (if exists).
nNumber of observations.
pNumber of features.
Method fit()
Constructs a C++ module for PELT.
Usage
PELT$fit(tsMat = NULL, covariates = NULL)
Arguments
tsMatNumeric matrix. A time series matrix of size
n \times pwhose rows are observations ordered in time. IftsMat = NULL, the method will use the previously assignedtsMat(e.g., set via the active binding$tsMator from a prior$fit(tsMat)). Default:NULL.covariatesNumeric matrix. A time series matrix having a similar number of observations as
tsMat. Required for models involving both dependent and independent variables. Ifcovariates = NULLand no prior covariates were set (i.e.,$covariatesis stillNULL), the model is force-fitted with only an intercept. Default:NULL.
Details
This method constructs a C++ PELT module and sets private$.fitted to TRUE, enabling the use of $predict() and $eval().
Returns
Invisibly returns NULL.
Method eval()
Evaluate the cost of the segment (a,b]
Usage
PELT$eval(a, b)
Arguments
aInteger. Start index of the segment (exclusive). Must satisfy
start < end.bInteger. End index of the segment (inclusive).
Details
The segment cost is evaluated as follows:
-
L1 cost function:
c_{L_1}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \tilde{y}_{(a+1)...b} \|_1where
\tilde{y}_{(a+1)...b}is the coordinate-wise median of the segment. Ifa \ge b - 1, return 0. -
L2 cost function:
c_{L_2}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \bar{y}_{(a+1)...b} \|_2^2where
\bar{y}_{(a+1)...b}is the empirical mean of the segment. Ifa \ge b - 1, return 0. -
SIGMA cost function:
c_{\Sigma}(y_{(a+1)...b}) := (b - a)\log \det \hat{\Sigma}_{(a+1)...b}where
\hat{\Sigma}_{(a+1)...b}is the empirical covariance matrix of the segment without Bessel's correction. Here, ifaddSmallDiag = TRUE, a small biasepsilonis added to the diagonal of estimated covariance matrices to improve numerical stability.
By default,addSmallDiag = TRUEandepsilon = 1e-6. In caseaddSmallDiag = TRUE, if the covariance matrix is numerically singular (its log-determinant cannot be computed) or its log-determinant is smaller than the lower boundp*log(epsilon), return(b - a)*p*log(epsilon), otherwise, output an error message. -
VAR(r) cost function:
c_{\mathrm{VAR}}(y_{(a+1)...b}) := \sum_{t = \max(a, r)+1}^{b} \left\| y_t - \hat c - \sum_{j=1}^r \hat A_j y_{t-j} \right\|_2^2where
\hat cand\hat A_jare the OLS estimates of the intercept and VAR coefficients on the segment. The lagged valuesy_{t-j}may come from beforea+1, so only the firstrobservations of the whole series lack a full set of lags. If the system is singular, an approximate (minimum-norm) least-squares solve is used. Ifb-a < p*r+1(i.e., not enough observations), ora > n-r(wherenis the time series length), return 0.
"LinearL2" for piecewise linear regression process with constant noise variance
c_{\text{LinearL2}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_2^2
where \hat{\beta} are OLS estimates on segment (a+1):b. If segment is shorter than the minimum number of
points needed for OLS, return 0.
-
"LinearSIGMA" for piecewise linear regression process with varying noise covariance
c_{\text{LinearSIGMA}}(y_{(a+1):b}) := (b-a)\log \det \hat\Sigma_{(a+1):b}where
\hat\Sigma_{(a+1):b}is the empirical covariance matrix of OLS residualsy - X\hat{\beta}on segment(a+1):b, divided byb-awith no degrees-of-freedom correction for the fitted coefficients, and otherwise estimated the same way as in the SIGMA cost function. -
"LinearL1" for piecewise linear regression process under L1 (least absolute deviations) loss
c_{\text{LinearL1}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_1where
\hat{\beta}is fit column-by-column via IRLS on segment(a+1):b, iterated until convergence (tol) ormaxIter. Unlike the other regression costs, this has noO(1)-per-segment closed form.
Returns
The segment cost.
Method predict()
Performs PELT given a linear penalty value.
Usage
PELT$predict(pen = 0)
Arguments
penNumeric. Penalty per change-point. Default:
0.
Details
The PELT algorithm detects multiple change-points by finding the set of break-points that globally minimises
a penalised cost function. PELT uses dynamic programming combined with a pruning rule to reduce the number of candidate change-points, achieving efficient computation.
Let [c_1, \dots, c_k, c_{k+1}] denote the set of segment end-points with 0 = c_0 < c_1 < c_2 < \dots < c_k < c_{k+1} = n,
where k is the number of detected change-points and n is the total number of data points.
Let c_{(c_{i-1}, c_i]} be the cost of segment (c_{i-1}, c_i].
The total penalised cost is
\text{TotalCost} = \sum_{i=1}^{k+1} c_{(c_{i-1}, c_i]} + \lambda \cdot k,
where \lambda is a linear penalty applied per change-point. PELT finds the set of endpoints that minimises this cost exactly.
The pruning step eliminates candidate change-points that cannot lead to an optimal solution,
allowing PELT to run in linear time with respect to the number of data points.
Temporary segment end-points are saved to private$.tmpEndPoints after $predict(), enabling users to call $plot() without
specifying endpoints manually.
Returns
An integer vector of regime end-points. By design, the last element is the number of observations.
Method segments()
Returns the cost and parameter estimates of each segment from the latest $predict().
Usage
PELT$segments()
Details
Let 0 = c_0 < c_1 < \dots < c_{k+1} = n be the end-points from the latest $predict(pen). Segment i is
(c_{i-1}, c_i], Cost is c_{(c_{i-1}, c_i]} and Params is its minimiser. The costs therefore sum to the
optimal total penalised cost minus \lambda \cdot k.
Temporary end-points are cleared by $fit(), so $predict() must be run again after modifying the object via
its active bindings.
Returns
A list with one element per segment (Start, End]. Each element is a list with:
StartStart index of the segment (exclusive, 0-based), same convention as
$eval().EndEnd index of the segment (inclusive).
CostThe segment cost, as returned by
$eval(Start, End).ParamsA named list of the segment parameter estimates, e.g.
list(mean = ...)for"L2". SeecostFactoryfor the fields returned by each cost function.
Method plot()
Plots change-point segmentation
Usage
PELT$plot(
d = 1L,
endPts,
dimNames,
main,
xlab,
tsWidth = 0.25,
tsCol = "#5B9BD5",
bgCol = c("#A3C4F3", "#FBB1BD"),
bgAlpha = 0.5,
ncol = 1L
)Arguments
dInteger vector. Dimensions to plot. Default:
1L.endPtsInteger vector. End points. Default: latest temporary changepoints obtained via
$predict().dimNamesCharacter vector. Feature names matching length of
d. Defaults to"X1", "X2", ....mainCharacter. Main title. Defaults to
"PELT: d = ...".xlabCharacter. X-axis label. Default:
"Time".tsWidthNumeric. Line width for time series and segments. Default:
0.25.tsColCharacter. Time series color. Default:
"#5B9BD5".bgColCharacter vector. Segment colors, recycled to length of
endPts. Default:c("#A3C4F3", "#FBB1BD").bgAlphaNumeric. Background transparency. Default:
0.5.ncolInteger. Number of columns in facet layout. Default:
1L.
Details
Plots change-point segmentation results. Based on ggplot2. Multiple plots can easily be
horizontally and vertically stacked using patchwork's operators / and |, respectively.
Returns
An object of classes gg and ggplot.
Method clone()
The objects of this class are cloneable with this method.
Usage
PELT$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.
Author(s)
Minh Long Nguyen edelweiss611428@gmail.com
Toby Dylan Hocking toby.hocking@r-project.org
Charles Truong ctruong@ens-paris-saclay.fr
Huy Nhat Minh Nguyen sleepysnorlax0115@gmail.com
References
Truong, C., Oudre, L., & Vayatis, N. (2020). Selective review of offline change point detection methods. Signal Processing, 167, 107299.
Killick, R., Fearnhead, P., & Eckley, I. A. (2012). Optimal detection of change points with a linear computational cost. Journal of the American Statistical Association, 107(500), 1590-1598.
Examples
## L2 example
set.seed(1121)
signals = as.matrix(c(rnorm(100,0,1),
rnorm(100,5,1)))
# Default L2 cost function
PELTObj = PELT$new(minSize = 1L, jump = 1L)
PELTObj$fit(signals)
PELTObj$predict(pen = 100)
PELTObj$plot()
## SIGMA example
set.seed(111)
signals = as.matrix(c(rnorm(100,-5,1),
rnorm(100,-5,10),
rnorm(100,-5,1)))
# L2 cost function
PELTObj = PELT$new(minSize = 1L, jump = 1L)
PELTObj$fit(signals)
# We choose pen = 50.
PELTObj$predict(pen = 50)
PELTObj$plot()
# The standard L2 cost function is not suitable.
# Use the SIGMA cost function.
PELTObj$costFunc = costFunc$new(costFunc = "SIGMA")
PELTObj$predict(pen = 50)
PELTObj$plot()
Slicing Window (Window)
Description
An R6 class implementing slicing window for offline change-point detection.
Details
Slicing window is a scalable, linear-time change-point detection algorithm that selects breakpoints based on local gains computed over sliding windows.
Currently supports the following cost functions:
-
"L1"and"L2"for (independent) piecewise Gaussian process with constant variance -
"SIGMA": for (independent) piecewise Gaussian process with varying variance -
"VAR": for piecewise Gaussian vector-regressive process with constant noise variance -
"LinearL2": for piecewise linear regression process with constant noise variance -
"LinearSIGMA": for piecewise linear regression process with varying noise covariance -
"LinearL1": for piecewise linear regression process under L1 (least absolute deviations) loss
Window requires a R6 object of class costFunc, which can be created via costFunc$new(). Currently, the following cost functions are supported:
-
"L1"and"L2"for (independent) piecewise Gaussian process with constant variance -
"SIGMA": for (independent) piecewise Gaussian process with varying variance -
"VAR": for piecewise Gaussian vector-regressive process with constant noise variance -
"LinearL2": for piecewise linear regression process with constant noise variance -
"LinearSIGMA": for piecewise linear regression process with varying noise covariance -
"LinearL1": for piecewise linear regression process under L1 (least absolute deviations) loss
See $eval() method for more details on computation of cost.
Some examples are provided below. See the package website for detailed usage!
Methods
$new()Initialises a
Windowobject.$describe()Describes the
Windowobject.$fit()Constructs a
Windowmodule inC++.$eval()Evaluates the cost of a segment.
$predict()Performs
Windowgiven a linear penalty value.$segments()Returns the cost and parameter estimates of each segment from the latest
$predict().$getHistory()Retrieves the full cost history and sequentially added breakpoints.
$plotElbow()Plots the elbow curve (Total Cost vs. Number of Change-Points).
$plot()Plots change-point segmentation in
ggplotstyle.$clone()Clones the
R6object.
Active bindings
minSizeInteger. Minimum allowed segment length. Can be accessed or modified via
$minSize. ModifyingminSizewill automatically trigger$fit().radiusInteger. Window radius. Can be accessed or modified via
$radius. Modifyingradiuswill automatically trigger$fit().jumpInteger. Search grid step size. Can be accessed or modified via
$jump. Modifyingjumpwill automatically trigger$fit().costFuncR6object of classcostFunc. Search grid step size. Can be accessed or modified via$costFunc. ModifyingcostFuncwill automatically trigger$fit().tsMatNumeric matrix. Input time series matrix of size
n \times p. Can be accessed or modified via$tsMat. ModifyingtsMatwill automatically trigger$fit().covariatesNumeric matrix. Input time series matrix having a similar number of observations as
tsMat. Can be accessed or modified via$covariates. Modifyingcovariateswill automatically trigger$fit().
Methods
Public methods
Method new()
Initialises a Window object.
Usage
Window$new(minSize, jump, radius, costFunc)
Arguments
minSizeInteger. Minimum allowed segment length. Default:
1L.jumpInteger. Search grid step size: only positions in {k, 2k, ...} are considered. Default:
1L.radiusInteger. Radius of each sliding window. Default:
1L.costFuncA
R6object of classcostFunc. Should be created viacostFunc$new()to avoid error. Default:costFunc$new("L2").
Returns
Invisibly returns NULL.
Method describe()
Describes a Window object.
Usage
Window$describe(printConfig = FALSE)
Arguments
printConfigLogical. Whether to print object configurations. Default:
FALSE.
Returns
Invisibly returns a list storing at least the following fields:
minSizeMinimum allowed segment length.
jumpSearch grid step size.
radiusRadius of each sliding window.
costFuncThe
costFuncobject.fittedWhether or not
$fit()has been run.tsMatTime series matrix.
covariatesCovariate matrix (if exists).
nNumber of observations.
pNumber of features.
Method fit()
Constructs a C++ module for Window.
Usage
Window$fit(tsMat = NULL, covariates = NULL)
Arguments
tsMatNumeric matrix. A time series matrix of size
n \times pwhose rows are observations ordered in time. IftsMat = NULL, the method will use the previously assignedtsMat(e.g., set via the active binding$tsMator from a prior$fit(tsMat)). Default:NULL.covariatesNumeric matrix. A time series matrix having a similar number of observations as
tsMat. Required for models involving both dependent and independent variables. Ifcovariates = NULLand no prior covariates were set (i.e.,$covariatesis stillNULL), the model is force-fitted with only an intercept. Default:NULL..
Details
This method constructs a C++ Window module and sets private$.fitted to TRUE,
enabling the use of $predict() and $eval(). Some precomputations are performed to allow
$predict() to run in linear time with respect to the number of local change-points
(see $predict() for more details).
Returns
Invisibly returns NULL.
Method eval()
Evaluate the cost of the segment (a,b]
Usage
Window$eval(a, b)
Arguments
aInteger. Start index of the segment (exclusive). Must satisfy
start < end.bInteger. End index of the segment (inclusive).
Details
The segment cost is evaluated as follows:
-
L1 cost function:
c_{L_1}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \tilde{y}_{(a+1)...b} \|_1where
\tilde{y}_{(a+1)...b}is the coordinate-wise median of the segment. Ifa \ge b - 1, return 0. -
L2 cost function:
c_{L_2}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \bar{y}_{(a+1)...b} \|_2^2where
\bar{y}_{(a+1)...b}is the empirical mean of the segment. Ifa \ge b - 1, return 0. -
SIGMA cost function:
c_{\Sigma}(y_{(a+1)...b}) := (b - a)\log \det \hat{\Sigma}_{(a+1)...b}where
\hat{\Sigma}_{(a+1)...b}is the empirical covariance matrix of the segment without Bessel's correction. Here, ifaddSmallDiag = TRUE, a small biasepsilonis added to the diagonal of estimated covariance matrices to improve numerical stability.
By default,addSmallDiag = TRUEandepsilon = 1e-6. In caseaddSmallDiag = TRUE, if the covariance matrix is numerically singular (its log-determinant cannot be computed) or its log-determinant is smaller than the lower boundp*log(epsilon), return(b - a)*p*log(epsilon), otherwise, output an error message. -
VAR(r) cost function:
c_{\mathrm{VAR}}(y_{(a+1)...b}) := \sum_{t = \max(a, r)+1}^{b} \left\| y_t - \hat c - \sum_{j=1}^r \hat A_j y_{t-j} \right\|_2^2where
\hat cand\hat A_jare the OLS estimates of the intercept and VAR coefficients on the segment. The lagged valuesy_{t-j}may come from beforea+1, so only the firstrobservations of the whole series lack a full set of lags. If the system is singular, an approximate (minimum-norm) least-squares solve is used. Ifb-a < p*r+1(i.e., not enough observations), ora > n-r(wherenis the time series length), return 0. -
"LinearL2" for piecewise linear regression process with constant noise variance
c_{\text{LinearL2}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_2^2where
\hat{\beta}are OLS estimates on segment(a+1):b. If segment is shorter than the minimum number of points needed for OLS, return 0. -
"LinearSIGMA" for piecewise linear regression process with varying noise covariance
c_{\text{LinearSIGMA}}(y_{(a+1):b}) := (b-a)\log \det \hat\Sigma_{(a+1):b}where
\hat\Sigma_{(a+1):b}is the empirical covariance matrix of OLS residualsy - X\hat{\beta}on segment(a+1):b, divided byb-awith no degrees-of-freedom correction for the fitted coefficients, and otherwise estimated the same way as in the SIGMA cost function. -
"LinearL1" for piecewise linear regression process under L1 (least absolute deviations) loss
c_{\text{LinearL1}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_1where
\hat{\beta}is fit column-by-column via IRLS on segment(a+1):b, iterated until convergence (tol) ormaxIter. Unlike the other regression costs, this has noO(1)-per-segment closed form.
Returns
The segment cost.
Method predict()
Performs Window given a linear penalty value, or a target number of change-points.
Usage
Window$predict(pen = 0, nBkps = NULL)
Arguments
penNumeric. Penalty per change-point. Ignored if
nBkpsis supplied. Default:0.nBkpsInteger. If supplied, takes precedence over
pen: returns thenBkpshighest-gain local maxima found by$fit()(see$getHistory()), i.e.Window's own best answer for that many change-points – not necessarily the globally optimal one for that count (seeDynpfor that guarantee). Treated as an upper bound, not a strict requirement:Windowonly ever finds finitely many local maxima, so if fewer thannBkpsexist, all of them are returned and a message reports the shortfall. Default:NULL.
Details
The algorithm scans the data with a fixed-size window to detect candidate local change-points (lcps)
if the gains of its k_\text{thresh} neighbors to the left and right are all smaller than its gain,
where k_\text{thresh} is defined as
k_\text{thresh} = \max \left( \left\lfloor \frac{\max(2\text{radius}, 2 \cdot \text{minSize})}{2 \cdot \text{jump}} \right\rfloor, 1 \right)
After candidate local change-points and computing the local gains, the algorithm selects the "optimal" set of break-points given the linear penalty threshold.
Let G_i denote the local gain for candidate change-point i, for i = 1, \dots, n_\text{lcps}. The local gains are ordered such that
G_1 \ge G_2 \ge \dots \ge G_{n_\text{lcps}}. Note that it is possible that no local change-points are detected,
for example if the window size is too large.
The total cost for the selected k change-points is then calculated as
\text{TotalCost} = - \sum_{i=1}^{k} G_i + \lambda \cdot k,
where \lambda is a linear penalty applied per change-point. We then optimise over
k = 1, \dots, n_\text{lcps} to minimise the penalised cost function. k = 0 is not a candidate, so at least one change-point is returned whenever a local change-point exists, however large the penalty.
This approach allows detecting multiple change-points in a time series while controlling model complexity through the linear penalty threshold.
In our implementation, scanning the data to detect candidate local change-points and computing their corresponding local gains is already performed in $fit().
Therefore, $predict() runs in linear time with respect to the number of local change-points.
Temporary segment end-points are saved to private$.tmpEndPoints after $predict(), enabling users to call $plot() without
specifying endpoints manually.
Returns
An integer vector of regime end-points. By design, the last element is the number of observations.
Method segments()
Returns the cost and parameter estimates of each segment from the latest $predict().
Usage
Window$segments()
Details
Let 0 = c_0 < c_1 < \dots < c_{k+1} = n be the end-points from the latest $predict(), called with either
pen or nBkps. Segment i is (c_{i-1}, c_i], Cost is c_{(c_{i-1}, c_i]} and Params is its
minimiser. The costs sum to the total cost of this segmentation. For k \ge 1 this generally differs from
$getHistory()$cost[k + 1], which subtracts local gains computed on windows of 2 * radius points instead of
re-evaluating the segments.
Temporary end-points are cleared by $fit(), so $predict() must be run again after modifying the object via
its active bindings.
Returns
A list with one element per segment (Start, End]. Each element is a list with:
StartStart index of the segment (exclusive, 0-based), same convention as
$eval().EndEnd index of the segment (inclusive).
CostThe segment cost, as returned by
$eval(Start, End).ParamsA named list of the segment parameter estimates, e.g.
list(mean = ...)for"L2". SeecostFactoryfor the fields returned by each cost function.
Method getHistory()
Retrieves the full cost history and sequentially added breakpoints.
Usage
Window$getHistory()
Returns
A data.frame with three columns:
kThe number of change-points.
costThe total unpenalised cost of the segmentation.
added_bkpThe breakpoint added at this step to achieve the cost.
Method plotElbow()
Plots the elbow curve (Total Cost vs. Number of Change-Points).
Usage
Window$plotElbow(maxK = NULL)
Arguments
maxKInteger. The maximum number of change-points to display on the plot. If
NULL, displays the full history. Default:NULL.
Returns
A ggplot object.
Method plot()
Plots change-point segmentation
Usage
Window$plot(
d = 1L,
endPts,
dimNames,
main,
xlab,
tsWidth = 0.25,
tsCol = "#5B9BD5",
bgCol = c("#A3C4F3", "#FBB1BD"),
bgAlpha = 0.5,
ncol = 1L
)Arguments
dInteger vector. Dimensions to plot. Default:
1L.endPtsInteger vector. End points. Default: latest temporary changepoints obtained via
$predict().dimNamesCharacter vector. Feature names matching length of
d. Defaults to"X1", "X2", ....mainCharacter. Main title. Defaults to
"Window: d = ...".xlabCharacter. X-axis label. Default:
"Time".tsWidthNumeric. Line width for time series and segments. Default:
0.25.tsColCharacter. Time series color. Default:
"#5B9BD5".bgColCharacter vector. Segment colors, recycled to length of
endPts. Default:c("#A3C4F3", "#FBB1BD").bgAlphaNumeric. Background transparency. Default:
0.5.ncolInteger. Number of columns in facet layout. Default:
1L.
Details
Plots change-point segmentation results. Based on ggplot2. Multiple plots can easily be
horizontally and vertically stacked using patchwork's operators / and |, respectively.
Returns
An object of classes gg and ggplot.
Method clone()
The objects of this class are cloneable with this method.
Usage
Window$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.
Author(s)
Minh Long Nguyen edelweiss611428@gmail.com
Toby Dylan Hocking toby.hocking@r-project.org
Charles Truong ctruong@ens-paris-saclay.fr
Huy Nhat Minh Nguyen sleepysnorlax0115@gmail.com
References
Truong, C., Oudre, L., & Vayatis, N. (2020). Selective review of offline change point detection methods. Signal Processing, 167, 107299.
Examples
## L2 example
set.seed(1121)
signals = as.matrix(c(rnorm(100,0,1),
rnorm(100,5,1)))
# Default L2 cost function
WindowObj = Window$new(minSize = 1L, jump = 1L)
WindowObj$fit(signals)
WindowObj$predict(pen = 100)
WindowObj$plot()
## SIGMA example
set.seed(111)
signals = as.matrix(c(rnorm(100,-5,1),
rnorm(100,-5,10),
rnorm(100,-5,1)))
# L2 cost function
WindowObj = Window$new(minSize = 1L, jump = 1L)
WindowObj$fit(signals)
# We choose pen = 50.
WindowObj$predict(pen = 50)
WindowObj$plot()
# The standard L2 cost function is not suitable.
# Use the SIGMA cost function.
WindowObj$costFunc = costFunc$new(costFunc = "SIGMA")
WindowObj$predict(pen = 50)
WindowObj$plot()
Binary Segmentation (binSeg)
Description
An R6 class implementing binary segmentation for offline change-point detection.
Details
Binary segmentation is a classic algorithm for change-point detection that recursively splits the data at locations that minimise the cost function.
binSeg requires a R6 object of class costFunc, which can be created via costFunc$new(). Currently, the following cost functions are supported:
-
"L1"and"L2"for (independent) piecewise Gaussian process with constant variance -
"SIGMA": for (independent) piecewise Gaussian process with varying variance -
"VAR": for piecewise Gaussian vector-regressive process with constant noise variance -
"LinearL2": for piecewise linear regression process with constant noise variance -
"LinearSIGMA": for piecewise linear regression process with varying noise covariance -
"LinearL1": for piecewise linear regression process under L1 (least absolute deviations) loss
See $eval() method for more details on computation of cost.
Some examples are provided below. See the package website for detailed usage!
Methods
$new()Initialises a
binSegobject.$describe()Describes the
binSegobject.$fit()Constructs a
binSegmodule inC++.$eval()Evaluates the cost of a segment.
$predict()Performs
binSeggiven a linear penalty value.$segments()Returns the cost and parameter estimates of each segment from the latest
$predict().$plot()Plots change-point segmentation in
ggplotstyle.$clone()Clones the
R6object.
Active bindings
minSizeInteger. Minimum allowed segment length. Can be accessed or modified via
$minSize. ModifyingminSizewill automatically trigger$fit().jumpInteger. Search grid step size. Can be accessed or modified via
$jump. Modifyingjumpwill automatically trigger$fit().costFuncR6object of classcostFunc. Search grid step size. Can be accessed or modified via$costFunc. ModifyingcostFuncwill automatically trigger$fit().tsMatNumeric matrix. Input time series matrix of size
n \times p. Can be accessed or modified via$tsMat. ModifyingtsMatwill automatically trigger$fit().covariatesNumeric matrix. Input time series matrix having a similar number of observations as
tsMat. Can be accessed or modified via$covariates. Modifyingcovariateswill automatically trigger$fit().
Methods
Public methods
Method new()
Initialises a binSeg object.
Usage
binSeg$new(minSize, jump, costFunc)
Arguments
minSizeInteger. Minimum allowed segment length. Default:
1L.jumpInteger. Search grid step size: only positions in {k, 2k, ...} are considered. Default:
1L.costFuncA
R6object of classcostFunc. Should be created viacostFunc$new()to avoid error. Default:costFunc$new("L2").
Returns
Invisibly returns NULL.
Method describe()
Describes a binSeg object.
Usage
binSeg$describe(printConfig = FALSE)
Arguments
printConfigLogical. Whether to print object configurations. Default:
FALSE.
Returns
Invisibly returns a list storing at least the following fields:
minSizeMinimum allowed segment length.
jumpSearch grid step size.
costFuncThe
costFuncobject.fittedWhether or not
$fit()has been run.tsMatTime series matrix.
covariatesCovariate matrix (if exists).
nNumber of observations.
pNumber of features.
Method fit()
Constructs a C++ module for binary segmentation.
Usage
binSeg$fit(tsMat = NULL, covariates = NULL)
Arguments
tsMatNumeric matrix. A time series matrix of size
n \times pwhose rows are observations ordered in time. IftsMat = NULL, the method will use the previously assignedtsMat(e.g., set via the active binding$tsMator from a prior$fit(tsMat)). Default:NULL.covariatesNumeric matrix. A time series matrix having a similar number of observations as
tsMat. Required for models involving both dependent and independent variables. Ifcovariates = NULLand no prior covariates were set (i.e.,$covariatesis stillNULL), the model is force-fitted with only an intercept. Default:NULL.
Details
This method constructs a C++ binSeg module and sets private$.fitted to TRUE,
enabling the use of $predict() and $eval(). Some precomputations are performed to allow
$predict() to run in linear time with respect to the number of data points
(see $predict() for more details).
Returns
Invisibly returns NULL.
Method eval()
Evaluate the cost of the segment (a,b]
Usage
binSeg$eval(a, b)
Arguments
aInteger. Start index of the segment (exclusive). Must satisfy
start < end.bInteger. End index of the segment (inclusive).
Details
The segment cost is evaluated as follows:
-
L1 cost function:
c_{L_1}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \tilde{y}_{(a+1)...b} \|_1where
\tilde{y}_{(a+1)...b}is the coordinate-wise median of the segment. Ifa \ge b - 1, return 0. -
L2 cost function:
c_{L_2}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \bar{y}_{(a+1)...b} \|_2^2where
\bar{y}_{(a+1)...b}is the empirical mean of the segment. Ifa \ge b - 1, return 0. -
SIGMA cost function:
c_{\Sigma}(y_{(a+1)...b}) := (b - a)\log \det \hat{\Sigma}_{(a+1)...b}where
\hat{\Sigma}_{(a+1)...b}is the empirical covariance matrix of the segment without Bessel's correction. Here, ifaddSmallDiag = TRUE, a small biasepsilonis added to the diagonal of estimated covariance matrices to improve numerical stability.
By default,addSmallDiag = TRUEandepsilon = 1e-6. In caseaddSmallDiag = TRUE, if the covariance matrix is numerically singular (its log-determinant cannot be computed) or its log-determinant is smaller than the lower boundp*log(epsilon), return(b - a)*p*log(epsilon), otherwise, output an error message. -
VAR(r) cost function:
c_{\mathrm{VAR}}(y_{(a+1)...b}) := \sum_{t = \max(a, r)+1}^{b} \left\| y_t - \hat c - \sum_{j=1}^r \hat A_j y_{t-j} \right\|_2^2where
\hat cand\hat A_jare the OLS estimates of the intercept and VAR coefficients on the segment. The lagged valuesy_{t-j}may come from beforea+1, so only the firstrobservations of the whole series lack a full set of lags. If the system is singular, an approximate (minimum-norm) least-squares solve is used. Ifb-a < p*r+1(i.e., not enough observations), ora > n-r(wherenis the time series length), return 0. -
"LinearL2" for piecewise linear regression process with constant noise variance
c_{\text{LinearL2}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_2^2where
\hat{\beta}are OLS estimates on segment(a+1):b. If segment is shorter than the minimum number of points needed for OLS, return 0. -
"LinearSIGMA" for piecewise linear regression process with varying noise covariance
c_{\text{LinearSIGMA}}(y_{(a+1):b}) := (b-a)\log \det \hat\Sigma_{(a+1):b}where
\hat\Sigma_{(a+1):b}is the empirical covariance matrix of OLS residualsy - X\hat{\beta}on segment(a+1):b, divided byb-awith no degrees-of-freedom correction for the fitted coefficients, and otherwise estimated the same way as in the SIGMA cost function. -
"LinearL1" for piecewise linear regression process under L1 (least absolute deviations) loss
c_{\text{LinearL1}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_1where
\hat{\beta}is fit column-by-column via IRLS on segment(a+1):b, iterated until convergence (tol) ormaxIter. Unlike the other regression costs, this has noO(1)-per-segment closed form.
Returns
The segment cost.
Method predict()
Performs binSeg given a linear penalty value, or a target number of change-points.
Usage
binSeg$predict(pen = 0, nBkps = NULL)
Arguments
penNumeric. Penalty per change-point. Ignored if
nBkpsis supplied. Default:0.nBkpsInteger. If supplied, takes precedence over
pen: returns the firstnBkpschange-pointsbinSeg's greedy splitting added (see$getHistory()), i.e. its own best answer for that many change-points – not necessarily the globally optimal one for that count (seeDynpfor that guarantee). Treated as an upper bound, not a strict requirement: if fewer thannBkpschange-points were found (e.g.minSizeleaves no further valid split), all of them are returned and a message reports the shortfall. Default:NULL.
Details
The algorithm recursively partitions a time series to detect multiple change-points. At each step, the algorithm identifies the segment that, if split, would result in the greatest reduction in total cost. This process continues until no further splits are possible (e.g., each segment is of minimal length or each breakpoint corresponds to a single data point).
Then, the algorithm selects the "optimal" set of break-points given the linear penalty threshold. Let [c_1, \dots, c_k, c_{k+1}] denote the set of segment end-points with 0 = c_0 < c_1 < c_2 < \dots < c_k < c_{k+1} = n,
where k is the number of detected change-points and n is the total number of data points.
and k is the number of change-points. Let c_{(c_{i-1}, c_i]} be the cost of segment (c_{i-1}, c_i].
The total penalised cost is then
\text{TotalCost} = \sum_{i=1}^{k+1} c_{(c_{i-1}, c_i]} + \lambda \cdot k,
where \lambda is a linear penalty applied per change-point. We then optimise over
k to minimise the penalised cost function.
This approach allows detecting multiple change-points in a time series while controlling model complexity through the linear penalty threshold.
In our implementation, the recursive step is carried out during $fit().
Therefore, $predict() runs in linear time with respect to the number of data points.
Temporary segment end-points are saved to private$.tmpEndPoints after $predict(), enabling users to call $plot() without
specifying endpoints manually.
Returns
An integer vector of regime end-points. By design, the last element is the number of observations.
Method segments()
Returns the cost and parameter estimates of each segment from the latest $predict().
Usage
binSeg$segments()
Details
Let 0 = c_0 < c_1 < \dots < c_{k+1} = n be the end-points from the latest $predict(), called with either
pen or nBkps. Segment i is (c_{i-1}, c_i], Cost is c_{(c_{i-1}, c_i]} and Params is its
minimiser. Both modes return the first k change-points of the greedy splitting path, so the costs sum to
$getHistory()$cost[k + 1]. This is the cost of the greedy segmentation, not necessarily the minimal cost for
k change-points (see Dynp).
Temporary end-points are cleared by $fit(), so $predict() must be run again after modifying the object via
its active bindings.
Returns
A list with one element per segment (Start, End]. Each element is a list with:
StartStart index of the segment (exclusive, 0-based), same convention as
$eval().EndEnd index of the segment (inclusive).
CostThe segment cost, as returned by
$eval(Start, End).ParamsA named list of the segment parameter estimates, e.g.
list(mean = ...)for"L2". SeecostFactoryfor the fields returned by each cost function.
Method getHistory()
Retrieves the full cost history and sequentially added breakpoints.
Usage
binSeg$getHistory()
Returns
A data.frame with three columns:
kThe number of change-points.
costThe total unpenalised cost of the segmentation.
added_bkpThe breakpoint added at this step to achieve the cost.
Method plotElbow()
Plots the elbow curve (Total Cost vs. Number of Change-Points).
Usage
binSeg$plotElbow(maxK = NULL)
Arguments
maxKInteger. The maximum number of change-points to display on the plot. If
NULL, displays the full history. Default:NULL.
Returns
A ggplot object.
Method plot()
Plots change-point segmentation
Usage
binSeg$plot(
d = 1L,
endPts,
dimNames,
main,
xlab,
tsWidth = 0.25,
tsCol = "#5B9BD5",
bgCol = c("#A3C4F3", "#FBB1BD"),
bgAlpha = 0.5,
ncol = 1L
)Arguments
dInteger vector. Dimensions to plot. Default:
1L.endPtsInteger vector. End points. Default: latest temporary changepoints obtained via
$predict().dimNamesCharacter vector. Feature names matching length of
d. Defaults to"X1", "X2", ....mainCharacter. Main title. Defaults to
"binSeg: d = ...".xlabCharacter. X-axis label. Default:
"Time".tsWidthNumeric. Line width for time series and segments. Default:
0.25.tsColCharacter. Time series color. Default:
"#5B9BD5".bgColCharacter vector. Segment colors, recycled to length of
endPts. Default:c("#A3C4F3", "#FBB1BD").bgAlphaNumeric. Background transparency. Default:
0.5.ncolInteger. Number of columns in facet layout. Default:
1L.
Details
Plots change-point segmentation results. Based on ggplot2. Multiple plots can easily be
horizontally and vertically stacked using patchwork's operators / and |, respectively.
Returns
An object of classes gg and ggplot.
Method clone()
The objects of this class are cloneable with this method.
Usage
binSeg$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.
Author(s)
Minh Long Nguyen edelweiss611428@gmail.com
Toby Dylan Hocking toby.hocking@r-project.org
Charles Truong ctruong@ens-paris-saclay.fr
Huy Nhat Minh Nguyen sleepysnorlax0115@gmail.com
References
Truong, C., Oudre, L., & Vayatis, N. (2020). Selective review of offline change point detection methods. Signal Processing, 167, 107299.
Hocking, T. D. (2024). Finite Sample Complexity Analysis of Binary Segmentation. arXiv preprint arXiv:2410.08654.
Examples
## L2 example
set.seed(1121)
signals = as.matrix(c(rnorm(100,0,1),
rnorm(100,5,1)))
# Default L2 cost function
binSegObj = binSeg$new(minSize = 1L, jump = 1L)
binSegObj$fit(signals)
binSegObj$predict(pen = 100)
binSegObj$plot()
## SIGMA example
set.seed(111)
signals = as.matrix(c(rnorm(100,-5,1),
rnorm(100,-5,10),
rnorm(100,-5,1)))
# L2 cost function
binSegObj = binSeg$new(minSize = 1L, jump = 1L)
binSegObj$fit(signals)
# We choose pen = 50.
binSegObj$predict(pen = 50)
binSegObj$plot()
# The standard L2 cost function is not suitable.
# Use the SIGMA cost function.
binSegObj$costFunc = costFunc$new(costFunc = "SIGMA")
binSegObj$predict(pen = 50)
binSegObj$plot()
costFactory class
Description
An R6 class for fast segment-cost evaluation and parameter estimation, without a detection algorithm.
Details
costFactory builds the C++ cost module selected by $costFunc at $fit(), and keeps it, so every query reuses
its precomputations, as PELT, binSeg and Window do. $eval() and $get_params() call the module directly:
segments are (a, b] with 0-based a, and all checks are done in C++.
$new() only stores the costFunc object; $fit() validates and attaches the data. $costFunc is an active
binding, so it can be inspected or replaced after construction – if data has already been supplied, replacing it
automatically triggers $fit() again.
$get_params() returns:
-
"L1":median. -
"L2":mean. -
"SIGMA":meanandcov(plusepsilonon the diagonal ifaddSmallDiag = TRUE). -
"VAR","LinearL2"and"LinearL1":coef, intercept first. -
"LinearSIGMA":coefandcov(plusepsilonon the diagonal ifaddSmallDiag = TRUE). -
"Custom":params, the output ofparamFun; an empty list ifparamFunisNULL.
Both covs are biased maximum-likelihood estimates: divided by the segment length b - a, with no Bessel
(n - 1) or degrees-of-freedom correction.
Methods
$new()Initialises a
costFactoryobject.$fit()Constructs the
C++cost module.$eval()Evaluates the cost of a segment.
$get_params()Returns the parameter estimates of a segment.
$segments()Returns the cost and parameter estimates of each segment of a given segmentation.
$clone()Clones the
R6object.
Active bindings
costFuncR6object of classcostFunc. Can be accessed or modified via$costFunc. ModifyingcostFuncwill automatically trigger$fit()if atsMathas already been fitted.
Methods
Public methods
Method new()
Initialises a costFactory object. Does not build the C++ cost module; call $fit() for that.
Usage
costFactory$new(costFunc)
Arguments
costFuncA
R6object of classcostFunc. Should be created viacostFunc$new()to avoid error. Default:costFunc$new("L2").
Returns
Invisibly returns NULL.
Method fit()
Validates the supplied data and constructs the C++ cost module selected by $costFunc.
Usage
costFactory$fit(tsMat = NULL, covariates = NULL)
Arguments
tsMatNumeric matrix. Time series of size
n \times p. IfNULL, the method will use the previously assignedtsMat(i.e., from a prior$fit(tsMat)). Default:NULL.covariatesNumeric matrix with
nrows, used by"LinearL2","LinearSIGMA"and"LinearL1". IfNULLand no priorcovariateswere set, the model is force-fitted with only an intercept. Default:NULL.
Details
This method constructs the C++ cost module and sets private$.fitted to TRUE, enabling the use
of $eval() and $get_params().
Returns
Invisibly returns NULL.
Method eval()
Evaluates the cost of the segment (a, b] with the module's eval().
Usage
costFactory$eval(a, b)
Arguments
aInteger. Start index (exclusive, 0-based).
bInteger. End index (inclusive).
Returns
The segment cost.
Method get_params()
Returns the parameter estimates of the segment (a, b] with the module's get_params().
Usage
costFactory$get_params(a, b)
Arguments
aInteger. Start index (exclusive, 0-based).
bInteger. End index (inclusive).
Returns
A named list; see Details.
Method segments()
Returns the cost and parameter estimates of each segment of a given segmentation.
Usage
costFactory$segments(endPts)
Arguments
endPtsInteger vector. Segment end-points, e.g. from
$predict()ofPELT,binSeg,WindoworDynp. Sorted internally; must be unique, at least1, and end atn.
Returns
A list with one element per segment (Start, End], in the same format as the $segments() of the
segmentation classes. Each element is a list with Start (exclusive, 0-based), End (inclusive), Cost
(as returned by $eval(Start, End)) and Params (as returned by $get_params(Start, End)).
Method clone()
The objects of this class are cloneable with this method.
Usage
costFactory$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.
Author(s)
Minh Long Nguyen edelweiss611428@gmail.com
Huy Nhat Minh Nguyen sleepysnorlax0115@gmail.com
Examples
set.seed(1)
tsMat = cbind(c(rnorm(100, 0), rnorm(100, 5, 5)))
cf = costFactory$new(costFunc$new("L2"))
cf$fit(tsMat)
cf$eval(0, 100)
cf$get_params(0, 100)
cf$segments(c(100, 200))
# `costFunc` is an active binding: swapping it re-fits automatically.
cf$costFunc = costFunc$new("SIGMA")
cf$eval(0, 100)
costFunc class
Description
An R6 class specifying a cost function
Details
Creates an instance of costFunc R6 class, used in initialisation of change-point detection modules. Currently
supports the following cost functions:
-
L1 cost function:
c_{L_1}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \tilde{y}_{(a+1)...b} \|_1where
\tilde{y}_{(a+1)...b}is the coordinate-wise median of the segment. Ifa \ge b - 1, return 0. -
L2 cost function:
c_{L_2}(y_{(a+1)...b}) := \sum_{t = a+1}^{b} \| y_t - \bar{y}_{(a+1)...b} \|_2^2where
\bar{y}_{(a+1)...b}is the empirical mean of the segment. Ifa \ge b - 1, return 0. -
SIGMA cost function:
c_{\Sigma}(y_{(a+1)...b}) := (b - a)\log \det \hat{\Sigma}_{(a+1)...b}where
\hat{\Sigma}_{(a+1)...b}is the empirical covariance matrix of the segment without Bessel's correction. Here, ifaddSmallDiag = TRUE, a small biasepsilonis added to the diagonal of estimated covariance matrices to improve numerical stability.
By default,addSmallDiag = TRUEandepsilon = 1e-6. In caseaddSmallDiag = TRUE, if the covariance matrix is numerically singular (its log-determinant cannot be computed) or its log-determinant is smaller than the lower boundp*log(epsilon), return(b - a)*p*log(epsilon), otherwise, output an error message. -
VAR(r) cost function:
c_{\mathrm{VAR}}(y_{(a+1)...b}) := \sum_{t = \max(a, r)+1}^{b} \left\| y_t - \hat c - \sum_{j=1}^r \hat A_j y_{t-j} \right\|_2^2where
\hat cand\hat A_jare the OLS estimates of the intercept and VAR coefficients on the segment. The lagged valuesy_{t-j}may come from beforea+1, so only the firstrobservations of the whole series lack a full set of lags. If the system is singular, an approximate (minimum-norm) least-squares solve is used. Ifb-a < p*r+1(i.e., not enough observations), ora > n-r(wherenis the time series length), return 0. -
"LinearL2" for piecewise linear regression process with constant noise variance
c_{\text{LinearL2}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_2^2where
\hat{\beta}are OLS estimates on segment(a+1):b. If segment is shorter than the minimum number of points needed for OLS, return 0. -
"LinearSIGMA" for piecewise linear regression process with varying noise covariance
c_{\text{LinearSIGMA}}(y_{(a+1):b}) := (b-a)\log \det \hat\Sigma_{(a+1):b}where
\hat\Sigma_{(a+1):b}is the empirical covariance matrix of OLS residualsy - X\hat{\beta}on segment(a+1):b, divided byb-awith no degrees-of-freedom correction for the fitted coefficients, and otherwise estimated the same way as in the SIGMA cost function (including theaddSmallDiag/epsilonstabilisation and lower-bound fallback). -
"LinearL1" for piecewise linear regression process under L1 (least absolute deviations) loss
c_{\text{LinearL1}}(y_{(a+1):b}) := \sum_{t=a+1}^b \| y_t - X_t \hat{\beta} \|_1where
\hat{\beta}is fit column-by-column via Iteratively Reweighted Least Squares (IRLS), iterated until the change in the fit's cost is withintolormaxIteriterations are reached. Unlike the other regression cost functions, this has noO(1)-per-segment closed form, since IRLS must be re-run on each queried segment. -
"Custom" for a user-defined cost function supplied from R
c_{\text{Custom}}(y_{(a+1):b}) := \text{evalFun}(y_{(a+1):b}, a, b)where
evalFunis a user-supplied function called on the raw segment matrix, plus the segment's own(a,b]bounds – the latter letevalFunalign the segment against any externally-captured, position-indexed data (e.g. a weight vector or exogenous seriesevalFuncloses over) without the package needing to know that data exists. Because each call crosses back into R, this is substantially slower per call than the built-in costs above; see$evalFunand$paramFun.
If active binding $costFunc is modified (via assignment operator), the default parameters will be used.
Methods
$new()Initialises a
costFuncobject.$pass()Describes the
costFuncobject.$clone()Clones the
costFuncobject.
Active bindings
costFuncCharacter. Cost function. Can be accessed or modified via
$costFunc. IfcostFuncis modified and required parameters are missing, the default parameters are used.pVARInteger. Vector autoregressive order. Can be accessed or modified via
$pVAR.addSmallDiagLogical. Whether to add a bias value to the diagonal of estimated covariance matrices to stabilise matrix operations. Can be accessed or modified via
$addSmallDiag.epsilonDouble. A bias value added to the diagonal of estimated covariance matrices to stabilise matrix operations. Can be accessed or modified via
$epsilon.interceptLogical. Whether to include the intercept in regression problems. Can be accessed or modified via
$intercept.tolDouble. IRLS convergence tolerance: iteration stops once the change in the fit's cost falls below
tol. Can be accessed or modified via$tol.maxIterInteger. Maximum number of IRLS iterations. Can be accessed or modified via
$maxIter.evalFunFunction. Required for
costFunc = "Custom". A user-defined cost function, called asevalFun(segment, a, b), wheresegmentis the numeric matrix of rows(a+1):bfor the queried segment(a,b](0-indexed, same convention as$eval(a, b)).aandbletevalFunalignsegmentagainst externally-captured, position-indexed data it closes over (e.g.externalSeries[(a+1):b]), which the package itself never needs to see. Must return a single numeric value. Can be accessed or modified via$evalFun.paramFunFunction or
NULL. Optional forcostFunc = "Custom". A user-defined function called asparamFun(segment, a, b)(same convention asevalFun), used by$get_params()to report segment-level estimates. IfNULL(default),$get_params()returns an empty list for"Custom". Can be accessed or modified via$paramFun.
Methods
Public methods
Method new()
Initialises a costFunc object.
Usage
costFunc$new(costFunc, ...)
Arguments
costFuncCharacter. Cost function. Supported values include
"L2","VAR", and"SIGMA". Default:L2....Optional named parameters required by specific cost functions.
If any required parameters are missing or null, default values will be used.For
"L1"and"L2", there is no extra parameter.For
"SIGMA", supported parameters are:addSmallDiagLogical. If
TRUE, add a small value to the diagonal of estimated covariance matrices to stabilise matrix operations. Default:TRUE.epsilonDouble. If
addSmallDiag = TRUE, a small positive value added to the diagonal of estimated covariance matrices to stabilise matrix operations. Default:1e-6.
For
"VAR",pVARis required:pVARInteger. Vector autoregressive order. Must be a positive integer. Default:
1L.
For
"LinearL2",interceptis required:interceptLogical. Whether to include the intercept in regression problems. Default:
TRUE.
For
"LinearSIGMA", supported parameters are:interceptLogical. Whether to include the intercept in regression problems. Default:
TRUE.addSmallDiagLogical. If
TRUE, add a small value to the diagonal of estimated residual covariance matrices to stabilise matrix operations. Default:TRUE.epsilonDouble. If
addSmallDiag = TRUE, a small positive value added to the diagonal of estimated residual covariance matrices to stabilise matrix operations. Default:1e-6.
For
"LinearL1", supported parameters are:interceptLogical. Whether to include the intercept in regression problems. Default:
TRUE.tolDouble. IRLS convergence tolerance: iteration stops once the change in the fit's cost falls below
tol. Default:1e-6.maxIterInteger. Maximum number of IRLS iterations. Default:
1000L.
For
"Custom", supported parameters are:evalFunFunction. Required. See
$evalFunfor details.paramFunFunction or
NULL. Optional. See$paramFunfor details.
Method pass()
Returns a list of configuration parameters to initialise detection modules.
Usage
costFunc$pass()
Method clone()
The objects of this class are cloneable with this method.
Usage
costFunc$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.
Author(s)
Minh Long Nguyen edelweiss611428@gmail.com
Examples
## L2 costFunc (default)
costFuncObj = costFunc$new()
costFuncObj$pass()
## SIGMA costFunc
costFuncObj = costFunc$new(costFunc = "SIGMA")
costFuncObj$pass()
# Modify active bindings
costFuncObj$epsilon = 10^-5
costFuncObj$pass()
costFuncObj$costFunc = "VAR"
costFuncObj$pass()