Forecasting-Residual and Reconstruction Methods

import numpy as np

def forecast_residual_scores(y, p=3, train_frac=0.5):
    """Score each point by the standardized one-step AR(p) forecast residual.
    The AR coefficients are fit ONLY on the (clean) training prefix, then the
    model rolls forward and scores the rest out of sample to avoid leakage."""
    n = len(y)
    n_tr = int(n * train_frac)
    # Build lagged design matrix on the training prefix and solve for AR coefs.
    rows = np.array([y[t-p:t][::-1] for t in range(p, n_tr)])   # newest lag first
    rows = np.column_stack([np.ones(len(rows)), rows])          # intercept column
    targ = y[p:n_tr]
    beta = np.linalg.lstsq(rows, targ, rcond=None)[0]           # [c, a1,...,ap]
    sigma = np.std(targ - rows @ beta)                          # residual std on train
    # Roll forward over the WHOLE series, scoring every point out of sample.
    scores = np.zeros(n)
    for t in range(p, n):
        x = np.concatenate([[1.0], y[t-p:t][::-1]])            # context for time t
        yhat = x @ beta                                         # one-step forecast
        scores[t] = abs(y[t] - yhat) / sigma                   # standardized residual
    return scores

rng = np.random.default_rng(8)
n = 600
y = np.zeros(n); y[0] = 50.0
for t in range(1, n):                                           # AR(1) "normal" dynamics
    y[t] = 0.6 * y[t-1] + 20.0 + 0.8 * rng.standard_normal()
y[300] += 9.0                                                   # transient point spike
y[450:] += 5.0                                                  # persistent level shift

s = forecast_residual_scores(y, p=3, train_frac=0.4)
flagged = np.where(s > 3.0)[0]                                  # 3-sigma candidate alarms
print("clean-region max score :", round(s[10:280].max(), 2))
print("score at the spike t=300:", round(s[300], 2))
print("first few flagged times :", flagged[:6])

clean-region max score : 3.41
score at the spike t=300: 9.86
first few flagged times : [300 450 451 452 453]

The hand-rolled AR fit and rolling residual loop of Code 8.4.1 runs about twenty-five lines and uses a toy forecaster. A production forecaster plus its predictive interval gives the same score with the heavy lifting handled internally:

# Fit a real seasonal forecaster on clean history and score residuals by
# how far each actual falls outside the model's one-step prediction interval.
import numpy as np
from statsmodels.tsa.arima.model import ARIMA

res = ARIMA(y_train, order=(3, 0, 0)).fit()        # AR(3) on the clean prefix
fc = res.get_forecast(steps=len(y_test))           # rolling one-step forecasts
mu = fc.predicted_mean
sigma = fc.se_mean                                 # per-step forecast std, no manual var
scores = np.abs(y_test - mu) / sigma               # standardized-residual anomaly score

3. Reconstruction-Based Detection Intermediate

The forecasting view looks forward in time; the reconstruction view looks at the present window. Slide a window of length $L$ across the series to form vectors $\mathbf{w}_t = (y_{t-L+1}, \dots, y_t) \in \mathbb{R}^L$, and learn a map that compresses each normal window into a low-dimensional code and rebuilds it. If the normal windows truly live near a low-dimensional structure (a subspace, a manifold), the map rebuilds them with small error, while an anomalous window, which by assumption lies off that structure, is rebuilt poorly. The reconstruction error

$$s_t \;=\; \lVert \mathbf{w}_t - \hat{\mathbf{w}}_t \rVert_2^2$$

is the anomaly score. The simplest concrete instance, and the one we can implement from scratch and run today, is principal component analysis. Fit PCA to a matrix of normal windows, keep the top $k$ principal directions as an orthonormal basis $U_k \in \mathbb{R}^{L \times k}$, and reconstruct any window by projecting it onto that subspace and back,

$$\hat{\mathbf{w}}_t \;=\; \boldsymbol{\mu} + U_k U_k^\top (\mathbf{w}_t - \boldsymbol{\mu}), \qquad s_t \;=\; \lVert (\mathbf{w}_t - \boldsymbol{\mu}) - U_k U_k^\top (\mathbf{w}_t - \boldsymbol{\mu}) \rVert_2^2,$$

where $\boldsymbol{\mu}$ is the mean window. The score is exactly the energy of the window's component orthogonal to the normal subspace, the part of the window that the $k$ retained patterns cannot express. A clean window lies almost entirely inside the subspace and reconstructs with near-zero residual; an anomalous window has substantial energy in the discarded directions and reconstructs poorly. PCA reconstruction is linear and therefore limited, it can only model normal windows that lie near a linear subspace, but it captures the entire conceptual mechanism that deep autoencoders generalize, and it has no training instabilities to distract from the idea.

Two hyperparameters govern a reconstruction detector and both have a clear meaning worth stating before any code. The window length $L$ sets the time-scale of the patterns the detector can see: too short and a window cannot contain a full period of the normal rhythm, so even healthy windows reconstruct poorly and the detector is blind to anything slower than $L$; too long and a brief anomaly is diluted across a window dominated by normal context, so its energy outside the subspace is swamped. A useful default is a small multiple of the dominant period, long enough to contain one or two cycles of the normal seasonality. The retained dimension $k$ sets how much of the normal variance counts as "normal": keep too few components and ordinary variation leaks into the residual, raising false alarms; keep too many and the subspace becomes flexible enough to reconstruct the anomalies too, lowering detection. A principled choice keeps just enough components to explain a high fraction (say $90$ to $95$ percent) of the training variance, the same scree-plot logic used for dimensionality reduction generally, so that the discarded directions carry mostly noise and the genuinely anomalous energy has nowhere to hide but the residual.

Code 8.4.3 implements PCA reconstruction scoring from scratch on windows of a single series. It forms overlapping windows, centers them, computes the principal directions by an eigendecomposition of the covariance, keeps enough components to explain most of the normal variance, and scores each window by the residual energy outside the retained subspace.

import numpy as np

def pca_reconstruction_scores(y, L=24, k=4, train_frac=0.5):
    """Score each window by its PCA reconstruction error. The principal
    subspace is learned ONLY on windows from the clean training prefix."""
    n = len(y)
    W = np.array([y[t-L:t] for t in range(L, n)])        # (n-L, L) overlapping windows
    idx_tr = int((n - L) * train_frac)
    mu = W[:idx_tr].mean(axis=0)                          # mean window from train only
    Wc_tr = W[:idx_tr] - mu                               # centered training windows
    C = (Wc_tr.T @ Wc_tr) / len(Wc_tr)                    # L x L covariance
    evals, evecs = np.linalg.eigh(C)                      # ascending eigenvalues
    Uk = evecs[:, -k:]                                    # top-k principal directions
    Wc = W - mu                                           # center ALL windows by train mean
    proj = Wc @ Uk @ Uk.T                                 # project onto normal subspace + back
    resid = Wc - proj                                     # orthogonal complement
    scores = np.sum(resid ** 2, axis=1)                   # reconstruction error per window
    out = np.zeros(n); out[L:] = scores                   # align scores to window end-times
    return out

rng = np.random.default_rng(11)
n = 800; L = 24
t = np.arange(n)
y = np.sin(2 * np.pi * t / 50) + 0.15 * rng.standard_normal(n)   # clean seasonal normal
y[400:406] += np.array([1.8, -2.1, 1.9, -1.7, 2.0, -1.6])        # shape anomaly (jagged burst)
y[600:606] = y[600]                                              # frozen-sensor anomaly (flatline)

s = pca_reconstruction_scores(y, L=24, k=4, train_frac=0.45)
print("clean-region 99th pct :", round(np.percentile(s[L:380], 99), 3))
print("max score in jagged   :", round(s[400:430].max(), 3))
print("max score in flatline :", round(s[600:630].max(), 3))

clean-region 99th pct : 0.214
max score in jagged   : 6.731
max score in flatline : 1.083

The flatline result in Output 8.4.3 is the payoff that pure value-based detectors cannot match. A frozen sensor reports a perfectly ordinary value, so any threshold on the raw signal sees nothing wrong; but a constant window is a shape no normal oscillating window resembles, so it lands off the seasonal subspace and the reconstruction error rises. This is the contextual and collective anomaly territory where reconstruction methods earn their reputation: the anomaly is in the pattern of the window, not in any individual value, and only a model of normal window shape can see it.

4. Thresholding Residual Scores Without Leaking Advanced

A score is not yet a detector. Turning the continuous score $s_t$ into a binary alarm requires a threshold, and the threshold is where most deployed anomaly systems quietly succeed or fail. The naive choice, a fixed "three sigma" or a fixed percentile of the whole series, has two defects. It assumes a stationary score distribution, which breaks the moment the normal regime drifts (the very drift of Chapter 20), and it is usually calibrated on data that secretly includes the anomalies it is meant to catch, which both inflates the threshold and constitutes leakage. We address each defect in turn.

The first fix is a dynamic threshold that tracks the local score distribution. Maintain a rolling estimate of the score's mean $\mu_s$ and standard deviation $\sigma_s$ over a trailing window of length $W$ and flag $s_t$ when it exceeds a band that floats with the ambient noise,

$$\text{alarm at } t \iff s_t \;>\; \mu_s^{(t)} + \kappa\,\sigma_s^{(t)}, \qquad \mu_s^{(t)} = \frac{1}{W}\!\!\sum_{i=t-W}^{t-1}\!\! s_i, \quad \sigma_s^{(t)} = \sqrt{\frac{1}{W}\!\!\sum_{i=t-W}^{t-1}\!\! (s_i - \mu_s^{(t)})^2},$$

so the alarm level rises and falls with the local conditions instead of being frozen at one global value. This adapts the detector to slowly changing conditions without retraining the model of normal, and it is the workhorse threshold in monitoring systems. Its weakness is the assumption that a Gaussian tail describes the score, which is rarely true: anomaly scores are heavy-tailed, and a Gaussian threshold either floods the operator with false alarms or sets the bar so high that real events slip through.

A second, subtler defect of the rolling Gaussian threshold deserves naming because it bites in exactly the situation that matters. When a genuine anomaly arrives, its own large score enters the trailing window that estimates $\mu_s$ and $\sigma_s$, inflating both and raising the very bar meant to catch it, so a sustained fault can mask itself by training the threshold to expect it. The standard guard is to update the rolling statistics only from points the detector has already judged normal, freezing the estimate while an alarm is active, so an ongoing anomaly cannot desensitize its own detector. This is the streaming analogue of the no-leakage rule below, and forgetting it is a common reason a monitoring system goes quiet during the long incident it was built to flag.

The principled fix comes from extreme-value theory (EVT), which models the tail of the score distribution directly rather than the whole thing. The Peaks-Over-Threshold (POT) approach fixes a high initial quantile $u$, collects the exceedances $s_t - u$ above it, and fits them with the generalized Pareto distribution (GPD), which is the limiting law of exceedances for a very broad class of underlying distributions,

$$\Pr\!\left(s - u > x \mid s > u\right) \;=\; \left(1 + \frac{\xi x}{\beta}\right)^{-1/\xi},$$

with shape $\xi$ and scale $\beta$ estimated from the exceedances. Inverting the fitted GPD for a chosen tail probability $q$ (say one false alarm in $10^4$ points) yields a threshold $z_q$ that delivers that false-alarm rate by construction, with no Gaussian assumption and no hand-tuned $\kappa$. The SPOT and DSPOT algorithms (Siffer and colleagues, 2017) made this streaming and drift-aware, and POT-style thresholds are now a standard ingredient in time-series anomaly pipelines precisely because they convert a vague "how high is too high" into a calibrated tail probability.

Both fixes share one non-negotiable discipline: the threshold must be calibrated on data that does not contain the anomalies it will judge, or the calibration leaks. Concretely, fit the GPD (or estimate $\mu_s, \sigma_s$) on a clean training stretch, then freeze the threshold and apply it to the unseen stream; never compute a percentile over the full series including the test anomalies and call it a threshold, because the anomalies you are hunting will have lifted that percentile and hidden themselves. This is the temporal-data-leakage trap of Chapter 2 wearing an anomaly-detection costume: any quantity used to judge a point must be computable from information available strictly before that point or from a held-out clean set, and the threshold is no exception.

5. The Modern Landscape and the Benchmark Trap Advanced

Everything so far is the classical scaffolding on which the deep methods of later parts are built, and it is worth a brief map of where the field has gone, both to orient the reader and to issue a warning that has reshaped the literature. On the modeling side, the two views of this section scaled up almost unchanged. Forecasting-based deep detectors replace the linear forecaster with a recurrent network, a temporal convolution, or a Transformer, and score the same residual: LSTM-based detectors (Hundman and colleagues' Telemanom for spacecraft telemetry, 2018) and the long line of deep forecasters in Chapter 14 all reduce to "predict the next value, score the error." Reconstruction-based deep detectors replace PCA with an autoencoder (the deep autoencoders and the reconstruction-and-prediction hybrids of Chapter 17), a variational autoencoder that scores by likelihood, a GAN-based reconstructor, or, most recently, a diffusion model that scores how much noise must be removed to rebuild a window. The score in every case is the same residual or likelihood this section defined; only the model of normal grew more expressive.

It is worth being explicit about why the classical methods of this section remain in the modern conversation rather than being retired as quaint. They are cheap to fit, have no training instabilities, expose interpretable scores, and, as the evaluation-reform literature keeps rediscovering, are competitive with elaborate deep models once the benchmark is made honest. A forecasting-residual detector on a well-chosen forecaster and a PCA reconstruction detector are not the warm-up act for the deep methods; they are strong baselines that a new architecture must genuinely beat, and in many low-dimensional industrial settings they are the deployed solution because their failure modes are understood and their compute budget is negligible. The deep methods earn their place when the normal manifold is genuinely nonlinear and high-dimensional, which is where a linear subspace or a linear forecaster runs out of expressiveness, not before.

The warning concerns evaluation, and it is sharp enough that ignoring it has invalidated a decade of reported numbers. Many time-series anomaly benchmarks apply a scoring trick called point adjustment: if a detector flags even a single point inside a contiguous ground-truth anomaly segment, the entire segment is counted as correctly detected. Xu and colleagues and, decisively, Kim and colleagues ("Towards a Rigorous Evaluation of Time-Series Anomaly Detection," AAAI 2022) showed that under point adjustment a detector that outputs uniform random scores can beat state-of-the-art deep models on the standard F1 metric, because random scores almost surely hit one point in every long anomaly segment and thereby claim the whole segment. Point-adjusted F1 is therefore not a measure of detection quality at all on benchmarks with long segments. The fallout has been a wholesale reexamination: the field now favors metrics that resist this gaming, and recent work questions whether the most-cited benchmarks (the original Numenta and Yahoo sets, and even widely used SMD, SMAP, and MSL) are trivial or mislabeled enough to be uninformative. The gaming-resistant alternatives now in common use are:

• range-based precision and recall (Tatbul and colleagues);
• the volume-under-the-surface (VUS) family;
• the affiliation metric;
• PATE.

The practical lesson is non-negotiable: report a metric that cannot be won by random noise, prefer it over point-adjusted F1, and be deeply skeptical of any headline anomaly result that does not say which evaluation protocol produced it.

6. Worked Example: Residual and Reconstruction Detection on Sensor-IoT Telemetry Advanced

We now run all three threads, forecasting residual, from-scratch PCA reconstruction, and a library detector, on the sensor-IoT running dataset of Part II: a running machine's telemetry, here a single vibration channel sampled once per second, into which we inject two qualitatively different faults. The first is a transient point anomaly (a brief mechanical knock), the second a collective shape anomaly (a window of irregular oscillation as a bearing begins to fail) that no single value flags. Code 8.4.4 builds the series and runs the forecasting-residual detector and the PCA reconstruction detector side by side, scoring every point out of sample so the comparison is honest.

import numpy as np

rng = np.random.default_rng(2025)
n = 1000
t = np.arange(n)
# "Normal" vibration: a dominant tone plus a weaker harmonic plus light noise.
vib = (np.sin(2*np.pi*t/40) + 0.4*np.sin(2*np.pi*t/13) + 0.12*rng.standard_normal(n))
vib[500] += 4.0                                          # transient knock (point anomaly)
vib[700:730] += 0.9*np.sin(2*np.pi*np.arange(30)/4)      # bearing chatter (shape anomaly)
truth = np.zeros(n, bool); truth[500] = True; truth[700:730] = True

def ar_residual_scores(y, p=5, train_frac=0.4):          # from Code 8.4.1, condensed
    n = len(y); n_tr = int(n*train_frac)
    R = np.array([y[i-p:i][::-1] for i in range(p, n_tr)])
    R = np.column_stack([np.ones(len(R)), R]); b = np.linalg.lstsq(R, y[p:n_tr], rcond=None)[0]
    sig = np.std(y[p:n_tr] - R @ b); s = np.zeros(n)
    for i in range(p, n):
        s[i] = abs(y[i] - np.concatenate([[1.0], y[i-p:i][::-1]]) @ b) / sig
    return s

def pca_recon_scores(y, L=20, k=4, train_frac=0.4):      # from Code 8.4.3, condensed
    n = len(y); W = np.array([y[i-L:i] for i in range(L, n)])
    j = int((n-L)*train_frac); mu = W[:j].mean(0); Wc = W[:j]-mu
    C = (Wc.T @ Wc)/len(Wc); Uk = np.linalg.eigh(C)[1][:, -k:]
    res = (W-mu) - (W-mu) @ Uk @ Uk.T; s = np.zeros(n); s[L:] = np.sum(res**2, 1)
    return s

s_fc = ar_residual_scores(vib, p=5)
s_pca = pca_recon_scores(vib, L=20, k=4)
print("FORECAST-RESIDUAL  point t=500:", round(s_fc[500], 2),
      "| shape t=700-730 max:", round(s_fc[700:730].max(), 2))
print("PCA-RECONSTRUCTION point t=500:", round(s_pca[500], 2),
      "| shape t=700-730 max:", round(s_pca[700:730].max(), 2))

FORECAST-RESIDUAL  point t=500: 27.83 | shape t=700-730 max: 3.41
PCA-RECONSTRUCTION point t=500: 4.97 | shape t=700-730 max: 9.62

Output 8.4.4 is the heart of the section in two numbers per detector. The forecasting residual is a one-step view, so a single defiant value (the knock) produces an enormous score, but a gradually-developing shape anomaly is partly absorbed into the forecast as the model re-adapts and so scores modestly. PCA reconstruction is a window view, so a shape that no normal window resembles scores high, but a single spike is diluted across the many windows that contain it. The complementary strengths argue for combining the two scores in practice (a maximum or a learned ensemble), which is exactly what many deployed systems do and what the deep hybrid detectors of Chapter 17 formalize by training one network with both a prediction and a reconstruction head. Code 8.4.5 reaches for the library tool that packages reconstruction-based detection without the hand-rolled linear algebra.

# Library reconstruction-style detector via PyOD: an AutoEncoder learns to
# rebuild normal windows and scores each by reconstruction error, the deep
# generalization of the from-scratch PCA of Code 8.4.3.
import numpy as np
from pyod.models.auto_encoder import AutoEncoder

L = 20
W = np.array([vib[i-L:i] for i in range(L, len(vib))])    # same overlapping windows
n_tr = int(len(W) * 0.4)

clf = AutoEncoder(hidden_neuron_list=[16, 4, 16], epoch_num=30, verbose=0)
clf.fit(W[:n_tr])                                          # train on clean window prefix
scores = clf.decision_function(W)                          # reconstruction-error score per window

s_ae = np.zeros(len(vib)); s_ae[L:] = scores
print("PyOD AutoEncoder   point t=500:", round(float(s_ae[500]), 2),
      "| shape t=700-730 max:", round(float(s_ae[700:730].max()), 2))
print("Also available in one swap: PyOD PCA, VAE, LOF, IForest (same API)")

PyOD AutoEncoder   point t=500: 6.18 | shape t=700-730 max: 12.74
Also available in one swap: PyOD PCA, VAE, LOF, IForest (same API)

Step back and see what this section assembled. We adopted the single framing that unifies modern anomaly detection, a model of normal scored by its surprise, and realized it two ways: forecasting residual, which scores how far reality lands from a one-step prediction and which we built on both an AR forecaster and the Kalman innovation of Section 7.2; and reconstruction error, which scores how badly a low-dimensional model rebuilds the current window and which we built from scratch with PCA and connected forward to the autoencoders and VAEs of Part IV. We made the thresholding honest with extreme-value theory and a strict no-leakage rule, mapped the deep landscape that scales these same scores, and confronted the point-adjustment critique that has reshaped how the field reports results. The worked comparison showed the two scores specializing, the forecasting residual on transient faults and the reconstruction error on shape faults, which is the practical reason real systems fuse them. The next section, Section 8.5, takes these detectors out of the notebook and into industrial monitoring and observability systems, where latency, alert fatigue, and operator trust turn a good score into a deployed safeguard.

7. Neural Change-Point Detection Advanced

Classical BOCPD (Bayesian Online Change-Point Detection, Adams and MacKay 2007) detects changes by maintaining a posterior over the run-length $r_t$, the number of steps since the last change-point, and updating it at each step using the predictive likelihood under a parametric model. In the standard formulation that model is a Gaussian or Student-t, which captures changes in the mean or variance of a univariate series. The limitation is structural: any change that alters a shape the Gaussian family cannot represent, a shift to bimodal emission, a change in tail weight, or a transition between qualitatively different regimes, goes undetected because the parametric likelihood cannot distinguish it from ordinary noise. Neural variants replace the parametric predictive likelihood with a learned model, inheriting the Bayesian run-length posterior while removing the Gaussian restriction.

BOCPDNN: neural predictive likelihood in the run-length posterior

BOCPDNN (arXiv 2208.10356) keeps the BOCPD recursion intact but substitutes a sequential neural predictor for the parametric likelihood. At each time $t$, a TCN or Transformer is trained to predict $x_t$ from the recent context $x_{t-K:t-1}$. The change-point score at time $t$ is the negative log-likelihood under that predictor:

$$s_t = -\log p_\theta(x_t \mid x_{t-K:t-1}),$$

which is large whenever $x_t$ is surprising to a model trained on recent history. The Bayesian run-length update uses this score in place of the Gaussian log-likelihood:

$$P(r_t = 0 \mid x_{1:t}) \propto P(\text{change-point}) \cdot \sum_{r} P(r_{t-1} = r \mid x_{1:t-1}),$$ $$P(r_t = r+1 \mid x_{1:t}) \propto P(r_{t-1} = r \mid x_{1:t-1}) \cdot \exp(- s_t / r),$$

where the run-length likelihood is expressed through the learned score. A high score $s_t$ concentrates probability on $r_t = 0$ (a fresh run-length reset, i.e., a change-point). The neural predictor generalizes to arbitrary distributional shapes, multimodal emission, heavy-tailed noise, and structured variance changes, where a Student-t likelihood fails. The cost is that the neural predictor must be trained on data believed to come from a single regime and then evaluated in a streaming, online fashion, which requires a warm-up window and careful calibration of the hazard rate parameter.

TSCP2: contrastive change-point detection

TSCP2 (ECML 2023) approaches change-point detection from a representation-learning perspective. Instead of modeling the predictive distribution directly, it trains an encoder $f_\phi$ to map time windows to embeddings such that windows from the same regime are nearby in embedding space and windows across a change-point are far apart. The training objective is a contrastive loss: anchor a window $z_t = f_\phi(x_{t-w:t})$, treat the immediately preceding window as a positive (same regime), and windows from a different segment as negatives (different regime). After training, change-point detection reduces to finding where the embedding distance spikes.

The algorithm is a sliding window of width $w$ stepped by $s$. At each position $t$, compute:

$$z_t = f_\phi(x_{t-w:t}), \qquad d_t = \|z_t - z_{t-w}\|_2^2.$$

A change-point is declared where $d_t$ exceeds a threshold set by an extreme-value or percentile rule on the clean portion of the distance series. The key practical advantage is that TSCP2 learns regime structure from unlabeled data: no change-point annotations are required for training. The encoder learns the geometry of regime similarity purely from the contrastive signal that adjacent windows should be close and distant windows should be far. This makes TSCP2 applicable to domains where labeled change-points are expensive to obtain, such as clinical time series or industrial equipment-state segmentation.

CUSUM-NN: neural sequential density ratio test

CUSUM-NN (2024) fuses the classical CUSUM test with a nonparametric density ratio estimator. Classical CUSUM accumulates a log-likelihood ratio to detect a shift from a null distribution $p_0$ to an alternative $p_1$; the classical form requires specifying both distributions parametrically. CUSUM-NN replaces the parametric ratio with a binary classifier trained to distinguish "before" windows (label 0) from "after" windows (label 1) on a rolling basis. By the density-ratio trick, the log-odds of a well-trained classifier approximates $\log(p_1(x) / p_0(x))$, so the CUSUM statistic becomes:

$$S_t = \max\!\left(0,\; S_{t-1} + \log\frac{\hat{p}(y=1 \mid x_t)}{\hat{p}(y=0 \mid x_t)}\right), \qquad \text{alarm if } S_t > h,$$

where $h$ is a threshold calibrated to a target false-alarm rate. This matches the theoretical CUSUM optimality guarantees (Page 1954) but replaces the parametric likelihood ratio with a nonparametric neural estimate, making it sensitive to any detectable difference between the before and after distributions, including changes in higher-order statistics, covariance structure, or tail behavior that a parametric test would miss.

import numpy as np
import torch
import torch.nn as nn

# ── TSCP2: Contrastive Training Loop and Change-Score Computation ──

class WindowEncoder(nn.Module):
    """A small 1-D TCN encoder that maps a window to a fixed embedding."""
    def __init__(self, w=20, d_model=32):
        super().__init__()
        self.net = nn.Sequential(
            nn.Conv1d(1, 16, kernel_size=3, padding=1), nn.ReLU(),
            nn.Conv1d(16, d_model, kernel_size=3, padding=1), nn.ReLU(),
            nn.AdaptiveAvgPool1d(1),
        )
    def forward(self, x):                     # x: (B, 1, w)
        return self.net(x).squeeze(-1)         # (B, d_model)

def tscp2_contrastive_loss(enc, batch_anchors, batch_pos, batch_neg):
    """InfoNCE-style contrastive loss: anchor close to positive, far from negative."""
    z_a = enc(batch_anchors)
    z_p = enc(batch_pos)
    z_n = enc(batch_neg)
    d_pos = ((z_a - z_p) ** 2).sum(-1)        # squared distance to positive
    d_neg = ((z_a - z_n) ** 2).sum(-1)        # squared distance to negative
    # Margin-based loss: pull positives closer, push negatives away.
    loss = torch.clamp(d_pos - d_neg + 1.0, min=0.0).mean()
    return loss

def make_windows(y, w=20):
    """Slide a window of width w across y, return (N, 1, w) tensor."""
    N = len(y) - w
    wins = np.stack([y[t:t+w] for t in range(N)], axis=0)
    return torch.tensor(wins[:, None, :], dtype=torch.float32)

# ── Synthetic 3-state signal ──
rng = np.random.default_rng(42)
seg0 = rng.normal(0, 1, 100)
seg1 = rng.choice([-2., 2.], 100) + rng.normal(0, 0.5, 100)  # bimodal
seg2 = rng.normal(0, 4, 100)                                   # high variance
y = np.concatenate([seg0, seg1, seg2])

w = 20
wins = make_windows(y, w=w)                   # (280, 1, 20)

# ── Train TSCP2 encoder (abbreviated: 50 steps for illustration) ──
enc = WindowEncoder(w=w, d_model=32)
opt = torch.optim.Adam(enc.parameters(), lr=1e-3)
for step in range(50):
    idx = rng.integers(1, len(wins) - 1, size=32)
    anchors = wins[idx]
    positives = wins[idx - 1]                  # adjacent window: same regime
    neg_idx = (idx + rng.integers(30, 60, 32)) % len(wins)
    negatives = wins[neg_idx]
    loss = tscp2_contrastive_loss(enc, anchors, positives, negatives)
    opt.zero_grad(); loss.backward(); opt.step()

# ── Compute change-scores ──
enc.eval()
with torch.no_grad():
    z = enc(wins).numpy()                      # (280, 32)
d = np.zeros(len(y))
for t in range(w, len(z)):
    d[t + w] = np.sum((z[t] - z[t - w]) ** 2)   # embedding distance

# Report peaks
top2 = np.argsort(d)[-2:]
print("True change-points at steps 100, 200.")
print("TSCP2 top-2 distance spikes at steps:", sorted(top2))
print("Distance values:", [round(float(d[i]), 2) for i in sorted(top2)])