Detectors

The online detector (Adams & MacKay 2007), the offline detector (Fearnhead 2006), and the helpers that turn the run-length posterior R into changepoints. All of them can be imported from the top-level package.

Core Bayesian changepoint detection algorithms.

This module implements both online and offline Bayesian changepoint detection algorithms using PyTorch for efficient computation and GPU acceleration.

offline_changepoint_detection(data, prior_function, likelihood_model, truncate=float('-inf'), device=None)

Offline Bayesian changepoint detection using dynamic programming.

Computes the exact posterior distribution over changepoint locations using the algorithm described in Fearnhead (2006).

Parameters:

Name	Type	Description	Default
data	Tensor	Time series data of shape [T] or [T, D] where T is time and D is dimensions.	required
prior_function	callable	Log prior probability mass of a segment length: prior_function(l) returns log P(length = l) for l = 1 .. T (Fearnhead's g). Use const_prior, geometric_prior or negative_binomial_prior with functools.partial. The mass on lengths 1 .. T - 1 must be below 1 (always true for a proper distribution; for const_prior this means p * (T - 1) < 1).	required
likelihood_model	BaseLikelihood	Likelihood model for computing segment probabilities.	required
truncate	float	Deprecated; default -inf (exact sum). A finite value reproduces the truncation rule of versions up to 1.1.0 and of the NumPy original: the sum over segment ends is cut at the first term that falls truncate nats below the running sum. That rule assumed the terms decay monotonically after a peak; they do not (for a segment start t, ends inside the true segment can be far less likely than the true end, then the sequence rises again), and with multivariate likelihoods the cut can discard the dominant term and return changepoint "probabilities" far above 1. Since the segment likelihoods are computed for every end in one vectorized call, the rule also saves no work. Kept only so old results can be reproduced.	float('-inf')
device	str, torch.device, or None	Device to place tensors on; defaults to the likelihood's device (the CPU unless the likelihood was built elsewhere).	None

Returns:

Name	Type	Description
Q	Tensor	Log evidence for data[t:] for each time t. Shape: [T].
P	Tensor	Log likelihood of segment [t, s] with no changepoints. Shape: [T, T].
Pcp	Tensor	Log probability of j-th changepoint at time t. Shape: [T-1, T-1]. Row j sums to the probability that there are more than j changepoints. Once a row's sum falls below exp(-1000), far under the smallest positive float64, the rows after it are left -inf without being computed: their true values are smaller still, and 0 in probability space either way.

Examples:

>>> import torch
>>> from functools import partial
>>> from bayesian_changepoint_detection import (
...     offline_changepoint_detection, const_prior, StudentT
... )
>>>
>>> data = torch.randn(100)
>>> prior_func = partial(const_prior, p=0.01)
>>> likelihood = StudentT()
>>> Q, P, Pcp = offline_changepoint_detection(data, prior_func, likelihood)
>>>
>>> # Get changepoint probabilities
>>> changepoint_probs = torch.exp(Pcp).sum(0)
>>> detected_changepoints = torch.where(changepoint_probs > 0.5)[0]

Notes

The backward recursion for Q and P takes O(T^2) time and memory. The changepoint table Pcp takes O(J T^2) time, where J is the number of rows computed: every row j holds the probability that there are more than j changepoints, which can only decrease with j, and the rows after it drops below exp(-1000) are skipped. J is about the largest plausible number of changepoints plus a margin (for three clear changes, about 190 rows whatever T), and at most T - 1.

Model (Fearnhead 2006, section 2): segment lengths are i.i.d. with mass function g, except the last segment, whose length is only known to be at least what is observed, P(length >= l) = 1 - G(l - 1) with G(l) = sum_{i <= l} g(i). Segments are independent given the changepoints. With P[t, s] the log marginal likelihood of data[t:s+1], the backward recursion (eq. 2 of the paper, 0-indexed) is

Q[t] = sum_{s=t}^{T-2} P[t, s] Q[s+1] g(s+1-t)  +  P[t, T-1] (1 - G(T-1-t))

and the changepoint posteriors are Pcp[0, t] = P[0, t] Q[t+1] g(t+1) / Q[0] for the first changepoint (first segment data[0:t+1] has length t + 1) and, for the j-th, a sum over the previous changepoint s of Pcp[j-1, s] P[s+1, t] Q[t+1] g(t-s) / Q[s+1]. Versions up to 1.0.x evaluated g at length minus one in the first row, paired g with the wrong segment length in the later rows, and included a "length 0" term in G; none of that is visible with const_prior, all of it is with the geometric and negative binomial priors.

References

Fearnhead, P. (2006). Exact and efficient Bayesian inference for multiple changepoint problems. Statistics and Computing, 16(2), 203-213.

online_changepoint_detection(data, hazard_function, likelihood_model, device=None)

Online Bayesian changepoint detection with run length filtering.

Processes data sequentially, maintaining a posterior distribution over run lengths (time since last changepoint) as described in Adams & MacKay (2007).

Parameters:

Name	Type	Description	Default
data	Tensor	Time series data of shape [T] or [T, D] where T is time and D is dimensions.	required
hazard_function	callable	Function that takes run length tensor and returns hazard probabilities. Should accept torch.Tensor of run lengths and return torch.Tensor of same shape.	required
likelihood_model	BaseLikelihood	Online likelihood model that maintains sufficient statistics.	required
device	str, torch.device, or None	Device to place tensors on.	None

Returns:

Name	Type	Description
R	Tensor	Run length posterior. R[r, t] is P(run length = r | x_0..x_{t-1}), i.e. column t is the posterior after t observations; column 0 is the prior (all mass at run length 0). Shape: [T+1, T+1].
map_run_lengths	Tensor	argmax of each column of R: the most likely run length after t observations. Shape: [T+1], dtype long. A changepoint shows up as a drop in this sequence; get_map_changepoints turns the drops into segment start indices, and changepoint_probabilities gives a calibrated probability per position at a chosen detection lag.

Examples:

>>> import torch
>>> from functools import partial
>>> from bayesian_changepoint_detection import (
...     online_changepoint_detection, constant_hazard, StudentT,
...     get_map_changepoints, changepoint_probabilities,
... )
>>>
>>> _ = torch.manual_seed(0)
>>> data = torch.cat([torch.randn(80), torch.randn(80) + 5])
>>> hazard_func = partial(constant_hazard, 100)  # Expected run length = 100
>>> likelihood = StudentT(alpha=0.1, beta=0.01, kappa=1, mu=0, device="cpu")
>>> R, map_run_lengths = online_changepoint_detection(
...     data, hazard_func, likelihood, device="cpu"
... )
>>> get_map_changepoints(R)
tensor([80])
>>> changepoint_probabilities(R, lag=10)[80] > 0.85
tensor(True)

Notes

This algorithm has O(T^2) time complexity but is naturally online and can process streaming data. The run length distribution is normalized at each step for numerical stability.

Why the second return value is a run length and not a probability. Under the Adams & MacKay recursion the posterior probability of run length 0 after each observation, R[0, t], is the hazard evaluated under the previous run-length distribution; with a constant hazard it is identically 1/lam and carries no information about the data. The evidence for a changepoint at position tau accumulates in the following columns as mass at run length k in column tau + k. Versions 1.0.x returned the un-normalized R[0, t] under the name changepoint_probs; that quantity could not detect changepoints. This version restores the pre-1.0 return value (the MAP run length) and adds changepoint_probabilities for the lagged probability.

References

Adams, R. P., & MacKay, D. J. (2007). Bayesian online changepoint detection. arXiv preprint arXiv:0710.3742.

changepoint_probabilities(R, lag=10)

Probability that a new segment started at each position, judged lag observations later.

result[tau] = R[lag, tau + lag]: the posterior probability, after observing x_0 .. x_{tau+lag-1}, that the current run length is exactly lag, which is the event "the segment containing the latest point began at tau". This is the quantity the original notebook plotted as R[Nw, Nw:] and the natural online detector with a fixed decision delay.

Parameters:

Name	Type	Description	Default
R	Tensor	Run length posterior from online_changepoint_detection.	required
lag	int	Detection delay in observations (default 10). lag=0 gives the run-length-0 posterior, which under a constant hazard equals the hazard rate for every tau >= 1 (and 1 at tau = 0, the prior) and is therefore uninformative; use lag >= 1.	10

Returns:

Type	Description
Tensor	Shape [T + 1 - lag]; entry tau refers to data index tau. The last lag positions cannot be judged yet and are not returned.

Examples:

>>> probs = changepoint_probabilities(R, lag=10)
>>> detected = torch.where(probs > 0.5)[0]

get_map_changepoints(R, threshold=None, min_separation=0)

Segment start indices implied by the MAP run-length path.

After t observations the MAP run length r_t = argmax R[:, t] says the current segment began at data index t - r_t. Whenever that implied start moves forward (the MAP run length drops), a changepoint is reported at the new start. This is the classic BOCPD decision rule and is what the pre-1.0 versions of this library exposed as maxes.

Parameters:

Name	Type	Description	Default
R	Tensor	Run length posterior from online_changepoint_detection.	required
threshold	float	Deprecated and ignored. Earlier versions thresholded R[0, :], which is not a changepoint signal (see online_changepoint_detection). For a thresholded probability use changepoint_probabilities.	None
min_separation	int	When the posterior is split between two nearby starts the MAP path can flip between them and both get reported. Starts closer than this many observations to an earlier reported start are dropped (default 0: report every distinct start).	0

Returns:

Type	Description
Tensor	Sorted data indices at which a new segment starts (long). Index 0 is never reported.

Examples:

>>> R, map_run_lengths = online_changepoint_detection(data, hazard_func, likelihood)
>>> get_map_changepoints(R)

compute_run_length_posterior(data, hazard_function, likelihood_model, device=None)

Compute the full run length posterior distribution.

This is a convenience function that returns just the run length posterior from online changepoint detection.

Parameters:

Name	Type	Description	Default
data	Tensor	Time series data.	required
hazard_function	callable	Hazard function for changepoint prior.	required
likelihood_model	BaseLikelihood	Online likelihood model.	required
device	str, torch.device, or None	Device to place tensors on.	None

Returns:

Type	Description
Tensor	Run length posterior distribution R[r, t].

Examples:

>>> posterior = compute_run_length_posterior(data, hazard_func, likelihood)
>>> # Most likely run length at each time
>>> map_run_lengths = torch.argmax(posterior, dim=0)

viterbi_changepoints(data, hazard_function, likelihood_model, device=None)

Most probable run-length path (Viterbi / max-product) under the BOCPD model.

online_changepoint_detection marginalizes over paths and returns the posterior of the run length at each step. This function instead keeps, for every run length, only the single best path leading to it, and returns the jointly most probable sequence of run lengths, i.e. the MAP segmentation of the series under the same model (hazard prior on segment boundaries, conjugate predictive likelihood within a segment).

Parameters:

Name	Type	Description	Default
data	Tensor	Time series of shape [T] or [T, D].	required
hazard_function	callable	Maps a tensor of run lengths to changepoint probabilities.	required
likelihood_model	BaseLikelihood	Fresh online likelihood model (it is consumed by this call).	required
device	str, torch.device, or None	Device to place tensors on; defaults to the likelihood's device.	None

Returns:

Name	Type	Description
run_lengths	Tensor	Shape [T + 1], dtype long. run_lengths[t] is the run length on the best path after t observations, with the same meaning as the row index of R: 0 means the segment ending with data[t - 1] is closed and a new one starts at data[t].
changepoints	Tensor	Data indices at which a new segment starts on the best path (t >= 1 with run_lengths[t] == 0), dtype long; the same convention as get_map_changepoints.

Examples:

>>> import torch
>>> from functools import partial
>>> from bayesian_changepoint_detection import (
...     viterbi_changepoints, constant_hazard, StudentT,
... )
>>> _ = torch.manual_seed(0)
>>> data = torch.cat([torch.randn(80), torch.randn(80) + 5])
>>> run_lengths, changepoints = viterbi_changepoints(
...     data, partial(constant_hazard, 100), StudentT(0.1, 0.01, 1, 0, device="cpu"),
...     device="cpu",
... )
>>> changepoints
tensor([80])

Notes

Same recursion as the forward pass with max in place of sum:

V[r + 1, t + 1] = V[r, t] + log p(x_t | r) + log(1 - H(r))
V[0, t + 1]     = max_r V[r, t] + log p(x_t | r) + log H(r)

in log space, vectorized over r at each step (O(T) per observation, O(T^2) total, like the forward pass). Versions 1.0.x summed over r in the second line, which is neither the forward pass nor Viterbi, and looped over r in Python.