📈 Bayesian Changepoint Detection

Find the points where a time series changes regime, with calibrated posterior probabilities instead of a threshold. Online (Adams & MacKay 2007) and offline (Fearnhead 2006) Bayesian changepoint detection on PyTorch tensors, with conjugate Normal-Gamma and Normal-Wishart likelihoods for univariate and multivariate series.

✨ Features

• 🔭 Online detection: the run-length posterior after every observation (Adams & MacKay 2007), for streams and for measuring how quickly a change would have been noticed
• 🔍 Offline detection: the exact posterior probability of a changepoint at every position given the whole series (Fearnhead 2006)
• 🎯 Calibrated outputs: probabilities you can threshold, MAP segment starts, and the single most probable segmentation (viterbi_changepoints)
• 📐 Conjugate likelihoods: Student-t predictive for univariate data (unknown mean and variance), multivariate-t for vector data (unknown mean and covariance), independent-features and covariance-only variants, Gamma-Poisson for counts, and Normal with known variance for mean changes at a known noise level
• 🧮 Verified mathematics: closed forms checked against scipy and against exhaustive enumeration of segmentations; every pinned number in the test suite says where it comes from
• ⚡ Vectorized recursions: both detectors are O(T²) with the inner work on tensors, not Python loops; 1 000 points offline in under 4 s on a laptop CPU
• 🖥️ Runs where your tensors are: CPU, CUDA or Apple MPS through one device argument, with measured guidance on when an accelerator is not worth it
• 🪶 One dependency: torch; NumPy, SciPy and Matplotlib are only needed for the tests and examples

🚀 Quick Start

Installation

pip install bayesian-changepoint

Using uv:

uv add bayesian-changepoint

The import name is bayesian_changepoint_detection, whatever the distribution is called:

import bayesian_changepoint_detection

To run the examples and notebooks, add the plot extra (pip install "bayesian-changepoint[plot]", or ".[plot]" from a clone).

Package names

bayesian-changepoint is the distribution name from 1.1.0 on. The same project was published before as bayescd (0.4, April 2022) and, earlier, as bayesian-changepoint-detection (0.2.dev1). Both are frozen at those releases and neither gets updates; if you have one installed, replace it:

pip uninstall bayescd bayesian-changepoint-detection
pip install bayesian-changepoint

The code is the same project and lives in the same repository; only the name on PyPI changed. 1.1.0 is a rewrite on PyTorch relative to 0.4 and changes the online API relative to 1.0.x — the CHANGELOG lists every breaking change.

System Requirements

• Python 3.9 or higher
• PyTorch 2.0 or higher (installed automatically). For a CUDA build of PyTorch, install it first following https://pytorch.org/get-started/locally/; the CPU build is enough for everything in this README.

📖 Usage

Online detection

The online detector processes the series one point at a time and keeps the posterior over the run length, the number of observations since the last change. Two helpers turn that posterior into changepoints.

from functools import partial

import torch

from bayesian_changepoint_detection import (
    StudentT,
    changepoint_probabilities,
    constant_hazard,
    get_map_changepoints,
    online_changepoint_detection,
)

torch.manual_seed(42)
data = torch.cat([
    torch.randn(50) + 0,  # first segment: mean 0
    torch.randn(50) + 3,  # second segment: mean 3
    torch.randn(50) + 0,  # third segment: mean 0
])

hazard = partial(constant_hazard, 250)  # prior: one change every ~250 points
likelihood = StudentT(alpha=0.1, beta=0.01, kappa=1, mu=0)  # unknown mean and variance

R, map_run_lengths = online_changepoint_detection(data, hazard, likelihood)

# Index of the first point of each new segment on the MAP run-length path.
# min_separation merges starts closer than that many points when the
# posterior hesitates between neighbors.
print(get_map_changepoints(R, min_separation=10))  # tensor([ 50, 100])

# Or a probability per position, judged `lag` observations later.
probs = changepoint_probabilities(R, lag=10)  # probs[t] refers to data index t
print(torch.where(probs[1:] > 0.5)[0] + 1)  # tensor([ 50, 100])

R[r, t] is P(run length = r | first t observations). Why not simply threshold R[0, :]? Under a constant hazard the posterior probability of run length 0 is the hazard rate at every step, whatever the data say; the evidence for a change at t shows up in the following columns, as mass at run length k in column t + k. changepoint_probabilities reads exactly that. viterbi_changepoints(data, hazard, likelihood) returns the single most probable run-length path instead, i.e. the MAP segmentation.

Streaming

online_changepoint_detection needs the whole series and returns the (T+1)² matrix R. For a stream of unknown length, feed observations one at a time to OnlineChangepointDetector; it keeps only the current run-length posterior. With max_run_length the memory and the time per observation stay bounded: run lengths above the bound are dropped and the posterior renormalized, which is exact until the bound is reached and a close approximation afterwards when segments are shorter than the bound.

from bayesian_changepoint_detection import OnlineChangepointDetector

detector = OnlineChangepointDetector(
    hazard, StudentT(alpha=0.1, beta=0.01, kappa=1, mu=0), max_run_length=500
)
starts = []
for x in data:  # any iterable: a socket, a file, a generator
    detector.update(x)
    # P(a segment started 10 observations ago), as changepoint_probabilities
    if detector.t > 10 and detector.changepoint_probability(lag=10) > 0.5:
        starts.append(detector.t - 10)
print(starts)  # [50, 100]

Without max_run_length each posterior equals the corresponding column of R from online_changepoint_detection (up to float32 rounding).

Direction and size of each change

segment_statistics summarizes the segments between changepoints: mean, standard deviation, and for every changepoint the change in mean and a Welch z-score (issue #42).

from bayesian_changepoint_detection import segment_statistics

stats = segment_statistics(data, get_map_changepoints(R, min_separation=10))
print(stats.means.tolist())         # about [0.096, 3.173, -0.165]
print(stats.mean_changes.tolist())  # about [3.078, -3.339]: up, then down

The z-score ignores that the changepoints were found in the same data, so it overstates significance; treat it as a rough guide. Offline positions are the last index of the old segment: pass positions + 1.

Offline detection

The offline detector sees the whole series and returns, for every position, the posterior probability that a segment ends there. It is usually sharper than the online detector; use it for retrospective analysis.

from bayesian_changepoint_detection import const_prior, offline_changepoint_detection
from bayesian_changepoint_detection.offline_likelihoods import StudentT as OfflineStudentT

prior = partial(const_prior, p=1 / (len(data) + 1))  # flat prior on segment length
Q, P, changepoint_log_probs = offline_changepoint_detection(data, prior, OfflineStudentT())

changepoint_probs = torch.exp(changepoint_log_probs).sum(0)  # P(a segment ends at t)
print(torch.where(changepoint_probs > 0.5)[0])  # tensor([49, 99])

The two detectors use different index conventions: online reports the first point of the new segment (50), offline the last point of the old one (49). See the FAQ.

Multivariate data

Pass a [T, d] tensor, one row per observation, and a multivariate likelihood; everything else is the same. All three detectors check their input first: data must be [T] or [T, d], non-empty, real and finite, and match the likelihood's dims. A transposed [d, T] tensor is rejected with a hint rather than read as d observations.

from bayesian_changepoint_detection import MultivariateT

dims = 3
mv_data = torch.cat([
    torch.randn(50, dims) + torch.tensor([0.0, 0.0, 0.0]),
    torch.randn(50, dims) + torch.tensor([2.0, -1.0, 1.0]),
    torch.randn(50, dims) + torch.tensor([0.0, 0.0, 0.0]),
])

R, _ = online_changepoint_detection(mv_data, hazard, MultivariateT(dims=dims))
print(get_map_changepoints(R, min_separation=10))  # tensor([ 48, 100])

The first start lands two points early on this draw: the lag-10 posterior puts 0.53 on 48, 0.13 on 49 and 0.26 on 50, and the MAP path takes the mode. Read changepoint_probabilities when the exact position matters.

Devices

Every likelihood and both detectors take a device argument. The default is the CPU; name a device on the likelihood ("cuda", "mps", or "auto" for the first available) to opt into an accelerator, and the detectors follow it. On a laptop the CPU is the faster choice for the online detector (measured: 6–30x faster than MPS), and the offline detector always runs on the CPU under MPS because it needs float64. How the argument is resolved, what has been measured, how to time your own workload and how much memory the tables need: docs/devices.md.

API at a glance

All public functions have NumPy-style docstrings with the formulas and the paper they come from; the API reference on the documentation site is generated from them.

🏗️ Architecture

bayesian_changepoint_detection/
├── __init__.py             # Public API and __version__ (from package metadata)
├── bayesian_models.py      # The two detectors, viterbi_changepoints, and the R helpers
├── streaming.py            # OnlineChangepointDetector: the online recursion one observation at a time
├── segments.py             # segment_statistics: mean, spread and direction of each change
├── online_likelihoods.py   # Online StudentT and MultivariateT: per-run-length predictive densities
├── offline_likelihoods.py  # Offline StudentT, MultivariateT, IndependentFeatures, FullCovariance: segment marginals
├── priors.py               # const_prior, geometric_prior, negative_binomial_prior (segment-length priors)
├── hazard_functions.py     # constant_hazard, negative_binomial_hazard
├── device.py               # get_device, get_device_info, to_tensor, ensure_tensor
└── generate_data.py        # Synthetic series with known changepoints, for tests and examples

Supporting directories: tests/ (the suite, see below), examples/ (scripts and two notebooks, all run in CI), docs/ (pages whose code blocks are executed by the tests), benchmarks/ (timings across released versions, see Performance).

🧪 Development

Setup Development Environment

# Clone repository
git clone https://github.com/hildensia/bayesian_changepoint_detection.git
cd bayesian_changepoint_detection

# Install with development dependencies
uv venv && source .venv/bin/activate
uv pip install -e ".[dev]"
# ...or, without uv:  python -m venv .venv && source .venv/bin/activate && pip install -e ".[dev]"

# Install pre-commit hooks (ruff lint + format on staged files)
pre-commit install

Running Tests

# Run all tests (about 15 s on a CPU)
pytest

# Only the tests that check the mathematics against independent references
pytest -m math

# Only the tests that pin current behavior (contracts, edge cases, devices, goldens)
pytest -m behavior

# With coverage
pytest --cov=bayesian_changepoint_detection --cov-report=term-missing

# One file
pytest tests/test_online_detection.py -v

Every test carries exactly one of the markers math and behavior; collection fails otherwise. Tests pass device="cpu" explicitly even though it is the default, so that a test never lands on an accelerator by accident (the suite is much slower there). The Python blocks in this README and in docs/ are executed as part of the suite.

Code Quality

# Lint with ruff
ruff check .

# Format code
ruff format .

# Type checking (configured, advisory: not enforced in CI)
mypy bayesian_changepoint_detection

ruff check and ruff format --check are enforced in CI, together with the test suite on Python 3.9–3.13, the example scripts, and a build job that installs the wheel into a clean environment. See CONTRIBUTING.md for the workflow and the review process.

Building

# Build sdist and wheel
python -m build

# Check the metadata PyPI will see
twine check --strict dist/*

📊 Example Output

examples/simple_example.py runs both detectors on a 150-point series with changes at 50 and 100 and saves a figure:

============================================================
Bayesian Changepoint Detection - Simple Example
============================================================
Generated data with 150 points
True changepoints at: [50, 100]
Running online changepoint detection...
✓ Online detection completed
  Segment starts on the MAP path: [50, 100]
  Max lag-10 changepoint probability (t > 0): 0.9117
Running offline changepoint detection...
✓ Offline detection completed
  Max changepoint probability: 0.9322

Detected changepoints:
  Online method: [50, 100]...
  Offline method: [49, 99]...
Creating visualization...
✓ Visualization saved as 'changepoint_detection_results.png'

============================================================
✅ Example completed successfully!
============================================================

Other scripts in examples/: basic_usage.py (400 points, four segments, both detectors), multivariate_example.py, gpu_acceleration.py (device selection and CPU/GPU comparison), benchmark_offline.py (offline timing at several lengths), and the notebooks Example_Code.ipynb and Multivariate_Example.ipynb. The scripts, and the notebooks' code cells (examples/run_notebooks.py), run in CI on every push.

⚡ Performance

Measured with benchmarks/performance.py, which runs the same series and parameters against this version, the previous release (1.1.0), the first PyTorch release (1.0.0) and the original NumPy implementation (0.4), each in a fresh process, and scores every run against the true changepoints (F1, margin 5). Apple M1, CPU, 4 threads, PyTorch 2.14, Python 3.12; median of up to five runs. Raw results, with more sizes: benchmarks/results/2026-09-23-apple-m1-cpu.json.

1.2.0 finds every change in each of these series (F1 1.00) except the streaming one (F1 0.94, 199 changes). In short: the offline detector is 31–51x faster than the NumPy original, 206–360x faster than 1.0.0 (which also misses changes) and faster than 1.1.0 by a factor that grows with length (1.5x at 250 points, 6x at 1 000, 12x at 2 000); the online detector runs at the speed of the NumPy original (both are a Python loop over time), and its multivariate model is correct only since 1.1.0.

Detection quality on real data, measured with benchmarks/tcpd.py on the Turing Change Point Dataset (van den Burg and Williams, 2020; 30 annotated real series, TCPDBench's F1 and covering metrics, higher is better):

The online model reproduces the published BOCPD (identical F1 on 27 of the 31 series both score, which add the 2-D run_log to these 30; also with tuned settings, 0.887 against 0.890); the offline detector beats it without tuning. For a finished series, use the offline detector, or read the online posterior with viterbi_changepoints. The online readout get_map_changepoints reports changes as data arrive, without hindsight, and scores 0.571 F1 on the same series.

Complexity: the online recursion is O(T²) in time and memory (the run-length posterior R is (T+1)² float32); OnlineChangepointDetector with max_run_length=K is O(K) per observation. The offline recursion is O(T²) (vectorized per start point), and so, in practice, is the table of changepoint locations Pcp: it is O(J T²) for J rows, and rows stop once the probability of that many changepoints drops below exp(-1000) (about 190 rows for a series with three clear changes, whatever its length; up to 19x faster than 1.1.0 at 4 000 points); memory is about 16 T² bytes (see docs/devices.md).

Accelerators: see the FAQ; MPS is slower than the CPU on all of these, CUDA is unmeasured (issue #43). Only measured numbers appear in this README.

❓ FAQ

Which detector should I use, online or offline?

online_changepoint_detection (Adams & MacKay 2007) processes the series one point at a time and, after each point, gives the posterior over how long the current segment has lasted. Use it for streams, or when you want to know how quickly a change would have been noticed. offline_changepoint_detection (Fearnhead 2006) sees the whole series and returns the posterior probability of a changepoint at each position, using data on both sides of it. Use it for retrospective analysis; it is usually sharper. Both cost O(T²).

The two detectors report the same change at indices one apart. Why?

Different conventions, both documented in the docstrings:

• Online (get_map_changepoints, changepoint_probabilities, viterbi_changepoints): the index of the first point of the new segment. A series whose first 80 points come from one regime reports 80.
• Offline (Pcp[j, t], and torch.exp(Pcp).sum(0)[t]): the probability that a segment ends at t, i.e. the last point of the old regime. The same series reports 79.

So offline index + 1 == online index.

Does the scale of my data matter? (issue #34)

Yes. The priors are on the mean and variance of the data, so their hyperparameters have units, and rescaling the data without rescaling them changes the model. For the univariate Normal-Gamma model (online StudentT with alpha, beta, kappa, mu; offline StudentT with alpha0, beta0, kappa0, mu0):

Multiplying the data by c is equivalent to using mu * c and beta * c² with kappa and alpha unchanged. With beta / alpha far from the actual within-segment variance, or mu far from the data, the first points of every segment look surprising and the detector over- or under-reacts.

Practical choices: standardize the data (subtract a typical level, divide by a typical within-segment standard deviation, ideally estimated on a calibration window rather than on the whole series), or set mu to the expected level and beta = alpha * expected_variance. The values in the examples (alpha=0.1, beta=0.01, kappa=1, mu=0) encode "around zero, variance about 0.1, but I am not sure": with df = 2 * alpha = 0.2 the predictive is extremely heavy-tailed, which is why they still work on roughly unit-scale data.

The multivariate classes work the same way but parametrize the prior on the covariance differently. Online MultivariateT takes scale, the Wishart scale W on the precision: to encode a prior covariance C pass scale = inv(C) / dof (default I / dof, unit prior covariance). Offline MultivariateT takes Psi0, the inverse-Wishart scale on the covariance side (Psi0 = inv(W)): the same prior covariance C is Psi0 = dof0 * C, and the default dof0 * I is the same unit prior covariance as online. mu/mu0 are in data units in both.

Data far from zero (timestamps, prices, counters) is handled without loss of precision: the offline likelihoods compute their statistics on data centered on its mean, and the online ones keep their state relative to the first observation, so an offset of 1e8 gives the same result as the same series around 0. Pass such data as float64 (a float64 tensor or NumPy array): a float32 value near 1e8 is only resolved to steps of 8, before the detector ever sees it.

How do I make the detector more or less sensitive? (issue #31)

In order of importance:

1. The hazard, i.e. the expected segment length. constant_hazard(lam) puts prior probability 1 / lam on a change at every step. Larger lam means fewer detections, more confidence needed, slightly longer delay; smaller lam means more, earlier, and more false alarms. This is the main knob and it is about the data, not the model: set it near the segment length you expect. If segments shorter than some minimum are implausible, negative_binomial_hazard(k, p) (mean length k / p) puts little probability on changes soon after the last one.
2. How much you trust the prior versus the first points of a new segment. kappa (for the mean) and alpha (for the variance) act as pseudo-counts. Small values let a few points establish a new regime quickly; larger values make the detector wait for more evidence. beta and mu should describe the data (previous question) rather than be used as sensitivity knobs.
3. How you read the output. changepoint_probabilities(R, lag) trades delay for confidence: a larger lag gives a more decisive probability, lag observations later. get_map_changepoints(R, min_separation=k) drops starts closer than k points to an earlier one, for when the posterior hesitates between neighboring points.

Offline, the equivalent of the hazard is the segment-length prior: const_prior(p=1/(T+1)) is the flat default; geometric_prior(p=1/L) encodes an expected segment length L; negative_binomial_prior allows a peaked length distribution, and negative_binomial_hazard with the same k and p is its online counterpart. Leave truncate at its default: the sum is exact and the legacy truncation can drop the dominant term.

My data are not normally distributed. Can I still use this? (issue #36)

Every likelihood here assumes that within a segment the observations are independent draws from one distribution. The Gaussian ones detect changes in the mean and/or the (co)variance; Poisson detects changes in the rate of counts:

When the data are not Gaussian the detector still runs, and the question is what the misspecification does to it:

• Heavy tails or outliers: single extreme points look like the start of a new segment. The Student-t predictive already tolerates some of this; a larger lam or kappa helps, and so does a transform (log for positive, right-skewed quantities such as latencies or prices).
• Counts: use online_likelihoods.Poisson / offline_likelihoods.Poisson (non-negative integers only). Counts that vary more than a Poisson allows (overdispersion) will show extra changepoints; a square-root or Anscombe transform with StudentT is the alternative.
• Bounded data: a transform (logit for proportions) usually gets you close enough.
• Autocorrelation or slow drift: the model has no notion of dynamics within a segment, so a drift is reported as a sequence of small changes. Differencing, or modeling residuals from a trend, is the usual fix.
• Changes in something other than mean or variance (e.g. in autocorrelation) are not detected.

In short: use it when "piecewise stationary with Gaussian-ish noise" is a reasonable description after a transform, and check on a segment you trust that the residuals look plausible.

Why is it slow on my laptop with a GPU?

Since 1.2.0 the default device is the CPU. Earlier versions selected CUDA or Apple MPS automatically when present, but the online recursion is a sequential loop over small tensors, and each step on an accelerator pays a launch cost. Measured on an Apple M-series laptop (PyTorch 2.14), CPU against MPS:

The offline detector needs float64 and always runs on the CPU when MPS is selected. On 1.1.0 or earlier, pass device="cpu" to both the likelihood and the detector. Opt into an accelerator only after measuring on your hardware; CUDA has not been benchmarked (issue #43).

🤝 Contributing

Contributions are welcome. Please see the Contributing Guidelines for the development setup, the conventions (including the rule that a test pinning a number says where the number comes from) and the review process.

1. Fork the repository
2. Create a feature branch (git checkout -b feat/amazing-feature)
3. Commit your changes (git commit -m 'feat: add amazing feature')
4. Push to the branch (git push origin feat/amazing-feature)
5. Open a Pull Request

Project documentation

📄 License

This project is licensed under the MIT License - see the LICENSE file for details.

🔗 Related Resources

• Ryan P. Adams and David J. C. MacKay (2007). Bayesian Online Changepoint Detection. arXiv:0710.3742. https://arxiv.org/abs/0710.3742 — the online algorithm.
• Paul Fearnhead (2006). Exact and Efficient Bayesian Inference for Multiple Changepoint Problems. Statistics and Computing 16(2), 203–213. https://doi.org/10.1007/s11222-006-8450-8 — the offline algorithm.
• Xiang Xuan and Kevin Murphy (2007). Modeling Changing Dependency Structure in Multivariate Time Series. ICML 2007, 1055–1062. https://doi.org/10.1145/1273496.1273629 — the multivariate likelihoods.
• Kevin P. Murphy (2007). Conjugate Bayesian analysis of the Gaussian distribution. Technical note. https://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf — the Normal-Gamma and Normal-Wishart closed forms used in the likelihoods.

🙏 Acknowledgements

• Johannes Kulick wrote the original NumPy implementation (2014–2022), published as bayesian-changepoint-detection and then bayescd, and owns this repository.
• Esteban Carisimo did the PyTorch rewrite, the vectorized recursions, the verified likelihoods and the current maintenance.

Citation

If you use this library in your research, please cite it (GitHub's "Cite this repository" button reads CITATION.cff):

@software{bayesian_changepoint_detection,
  title   = {Bayesian Changepoint Detection: A PyTorch Implementation},
  author  = {Kulick, Johannes and Carisimo, Esteban},
  url     = {https://github.com/hildensia/bayesian_changepoint_detection},
  year    = {2026},
  version = {1.2.0}
}

The algorithms are due to Adams & MacKay (2007) and Fearnhead (2006); please cite those papers as well.