regscan

Regression-based scan statistics for detecting interval anomalies in smoothly varying 1-D signals.

Given a signal and a function family F, the score of an interval I = [a, b] is

S(I) = 1 - (SR_I + SR_O) / SR_A

where SR_A, SR_I and SR_O are the sums of squared residuals from fitting F to the whole signal, to the inside of I, and to the outside. The score approaches 0 when splitting the signal explains nothing and approaches 1 when it explains everything. A scan returns the highest-scoring interval.

Installation

pip install regscan

Requires Python 3.10 or later. The only dependencies are numpy and numba.

The kernel scan is a tight scalar loop, which numba compiles well and numpy vectorizes poorly, so it is JIT-compiled. The first call in a process pays a one-off compilation cost; numba caches the result on disk thereafter.

Quick start

import numpy as np
import regscan

t = np.linspace(0, 1, 300)
x = 4.7 + 0.10 * t + 0.05 * np.random.default_rng(0).normal(0, 1, 300)
x[140:170] -= 0.25                      # the anomaly

res = regscan.scan(x, method="nwkr_gaussian", w=12)
res.score, res.a, res.b                 # 0.248, 140, 169
res.width_frac, res.at_edge()           # 0.10, False

On the same signal with no anomaly planted, nwkr_gaussian scores 0.036 while mean scores 0.185. The constant family is responding to the ramp, not to an anomaly.

Notation

n is fixed by the data. w and r are the parameters worth thinking about, and both may be passed explicitly to regscan.scan.

Methods

regscan.available() lists them at run time.

F_KR is the method this package exists for. A weak family such as F_0 or F_1 cannot represent a curved background, so it lowers the residual by splitting the interval wherever the curvature is worst, which flags smooth structure as an anomaly. A kernel fit follows that structure instead, so the structure enters SR_A, SR_I and SR_O alike and cancels out of the score.

Choosing w and r

w sets the scale of structure the fit can follow. Set it too small and the kernel reproduces the anomaly itself, which cancels it from the score; set it too large and the fit cannot follow the background, which is the failure mode of the weaker families.

r is the longest interval the scan will consider. An anomaly wider than r cannot be returned at all, and cost falls linearly with r, so it is the parameter to reach for when the approximate width of the feature is known.

regscan.scan(x, method="nwkr_gaussian", w=12, r=90)

When omitted, w defaults to max(3, n // 16) and r to 3 * w. These defaults exist so that a scan runs unattended; they are not recommendations.

Sweeping w across a range and checking whether (a, b) holds steady is a cheap way to distinguish a resolved feature from an artifact of the bandwidth. A real interval stays put, whereas one that tracks w is measuring the kernel.

Configuration

from regscan import ScanConfig

cfg = ScanConfig(
    kernel="laplace",     # "gaussian" or "laplace"
    buffer=24,            # exclude this many samples at each end
    min_width=0.01,       # interval length bounds, as a fraction of n
    max_width=0.25,
)
regscan.scan(x, method="nwkr_gaussian", w=16, config=cfg)

ScanConfig is immutable and passed explicitly. Nothing is stored in module globals, so scanning several methods in one process cannot leak settings between them.

buffer excises samples rather than merely forbidding interval placement there, so the fit, SR_A and the outside residuals are all computed on the trimmed signal. Because a non-zero buffer suppresses detections at the ends of the signal, give every family the same value when comparing them.

max_width is worth capping. As an interval approaches n/2 the inside and outside become comparable in size, the statistic stops discriminating, and the maximizer drifts toward whatever split best absorbs slow curvature.

Super-resolution

F_KR can block-mean the signal, scan the shorter version, then search the original samples around the winning blocks to recover exact endpoints.

from regscan import ScanConfig

regscan.scan(x, method="nwkr_gaussian", w=100,
             config=ScanConfig(super_resolution=4))     # or "auto"

At n = 1600 with w = 100:

A factor of 4 runs 39 times faster and returns the same interval. A factor of 8 does not, which is the trade-off: the coarse pass locates each endpoint only to within a block of factor samples, and a wrong block gives a wrong answer.

Refinement therefore searches one block either side of each coarse endpoint. This matters more than it may appear. Block means smooth an anomaly's edges, so the coarse pass selects a neighboring block often enough that confining the search to the winning block alone recovered the exact interval in only 20 of 40 trials at factor 4. Including the neighbors recovered all 40.

It remains an approximation. Verify against super_resolution=1 on a sample of your own data before relying on it.

Passing "auto" derives the factor from the signal length: 1 below 450 samples, then doubling at 900, 1800 and so on. sr_cap bounds it, which matters for narrow features, since an interval must survive decimation to be found. The default is 1, which scans every sample exactly.

Performance

One full scan with w = n // 16, after JIT warm-up:

At n = 800, F_KR is six times faster than F_1 despite fitting a far richer model. That is the practical case for it: the polynomial family pays O(d³) for every candidate interval, whereas the kernel family pays O(r) to extend one.

Each family reaches its stated complexity by carrying state rather than refitting. mean scores an interval in O(1) from prefix sums of y and y². poly_deg1 uses prefix sums of the moments t^p and t^p y, so a fit costs O(d³) to solve regardless of how many samples the interval spans. F_KR grows an interval one sample at a time and updates in O(r): the inside buffer and sse_in, the nin and din arrays holding the inside points' kernel contribution to every index, and sse_out, obtained from those by subtraction from the all-points totals. Since sse_out is adjusted rather than recomputed, it is refreshed exactly on a fixed cadence to keep floating-point error from accumulating.

Super-resolution reduces these times substantially again on longer signals.

Comparing scores

Scores are comparable within a family but not across families. Each family divides by its own SR_A, and a kernel fit has a smaller SR_A than a constant fit before any interval is chosen, so the same interval scores differently under F_0 and F_KR. Compare families by rank, by whether they agree on the interval, or by the contrast between anomalous and clean signals, rather than by absolute score.

Citation

Rakib et al., Efficient Regression Models for Scan Statistics, arXiv:2608.22201 (2026).

@misc{rakib2026efficientregressionmodelsscan,
      title={Efficient Regression Models for Scan Statistics},
      author={Gazi Abdur Rakib and Tristan Ashton and Ryan A. Loomis and Brian S. Mason and Eric J. Murphy and Ci Xue and Jeff M. Phillips},
      year={2026},
      eprint={2608.22201},
      archivePrefix={arXiv},
      primaryClass={stat.ME},
      url={https://arxiv.org/abs/2608.22201},
}

License

BSD 3-Clause. See LICENSE.