Time-series anomaly detection

Deep dive

How it works

Parts:

• Baseline — expected value for this phase of the cycle: per-phase median over the last D cycles, STL seasonal plus trend, or a forecast.
• Scale — rolling median and MAD of recent residuals (robust statistics).
• Score and threshold — z = (r - median) / (1.4826 * MAD), flag when |z| > z*.
• Incident builder — merges flags separated by at most g quiet points; carries start, end and peak score.

Operations (streaming: only past data is used):

1. Look up the baseline for the current phase (same time of day; same weekday if weekly).
2. Residual = value minus baseline; push it into the residual window.
3. Window median and MAD (recompute or update incrementally), current point excluded.
4. Score z; flag if |z| > z*.
5. Extend the open incident if the last flag is within g points, else open one; close it after g quiet points and deliver it once (peak score, duration, direction).

x --> [baseline by phase, median of D past cycles] --> r = x - base
r --> [rolling median / MAD of last w residuals] --> z --> |z|>z* --> flags --> [merge gaps<=g] --> incidents

Anomaly types (point, contextual, collective: Chandola et al. 2009). A spike is a point anomaly. A seasonal/contextual one is normal in value but wrong for its phase; only the baseline catches it. A dip or level shift is collective: a run of same-sign residuals. A variance change gives many moderate residuals, so also alert on the rolling MAD. A shift that never returns is a change point; the baseline absorbs it after about D/2 cycles.

Thresholds. "3 sigma" is folklore: per-point alarm probability times points per day is your pager load. Pick the per-point rate a (N series at T points per day give N*T*a false flags per day), then z* = Phi^-1(1 - a/2); 5 gives a = 5.7e-7. Real residuals are heavy-tailed and autocorrelated: measure the realised rate on clean history and calibrate z* to it.

Contaminated history. Past anomalies sit in the baseline data. Per-phase medians, MAD and STL's robust weights tolerate a minority; mean-based fits do not. Exclude known incidents and holidays when possible.

Many series. At thousands of series even a small a pages every hour. Use Benjamini-Hochberg across series, k-consecutive or m-of-n flag rules, and grouping of co-occurring incidents into one page.

S-H-ESD (Hochenbaum et al. 2017): seasonal decomposition, then generalized ESD (Rosner 1983) with median and MAD instead of mean and standard deviation; usable when up to 50 % of the data is anomalous.

Code

Std-only. Seeded synthetic metric (daily season, noise) with an injected spike, level shift and dip; the asserts require all three as incidents and at most one false incident.

use std::time::Instant;

const P: usize = 144; // samples per day (10-minute data)
const DAYS: usize = 7; // trailing days in the per-phase baseline
const W: usize = 288; // residuals in the rolling scale window
const Z_THR: f64 = 5.0; // two-sided Gaussian tail 5.7e-7 per point
const MAX_GAP: usize = 5; // unflagged points tolerated inside one incident

struct Rng(u64);
impl Rng {
    fn next(&mut self) -> f64 { // xorshift64*, uniform in (0,1)
        self.0 ^= self.0 >> 12; self.0 ^= self.0 << 25; self.0 ^= self.0 >> 27;
        ((self.0.wrapping_mul(0x2545F4914F6CDD1D) >> 11) as f64 + 0.5) / (1u64 << 53) as f64
    }
    fn normal(&mut self) -> f64 { // Box-Muller
        (-2.0 * self.next().ln()).sqrt() * (std::f64::consts::TAU * self.next()).cos()
    }
}

fn median(v: &mut [f64]) -> f64 { // selection, O(n); NaNs must be removed upstream
    let (n, mid) = (v.len(), v.len() / 2);
    let (lo, m, _) = v.select_nth_unstable_by(mid, |a, b| a.total_cmp(b));
    if n % 2 == 1 { return *m; }
    0.5 * (lo.iter().cloned().fold(f64::MIN, f64::max) + *m)
}

/// Incidents (start, end inclusive, peak |z|) from a stream; uses only past data at each step.
fn detect(x: &[f64]) -> Vec<(usize, usize, f64)> {
    let (mut res, mut tmp) = (vec![0.0; x.len()], Vec::with_capacity(W));
    let mut flags: Vec<(usize, f64)> = Vec::new();
    for t in DAYS * P..x.len() {
        // 1. expected value: median of the same phase over the last DAYS days
        let mut same: [f64; DAYS] = std::array::from_fn(|d| x[t - (d + 1) * P]);
        res[t] = x[t] - median(&mut same);
        // 2. scale from the previous W residuals, current point excluded
        let from = t.saturating_sub(W).max(DAYS * P);
        if t - from < W / 2 { continue; } // warm-up
        tmp.clear(); tmp.extend_from_slice(&res[from..t]);
        let c = median(&mut tmp);
        tmp.iter_mut().for_each(|v| *v = (*v - c).abs());
        let s = 1.4826 * median(&mut tmp).max(1e-9); // MAD -> sigma; floor avoids /0
        // 3. robust z and threshold
        let z = (res[t] - c) / s;
        if z.abs() > Z_THR { flags.push((t, z.abs())); }
    }
    // 4. merge flags closer than MAX_GAP into incidents
    let mut out: Vec<(usize, usize, f64)> = Vec::new();
    for (t, z) in flags {
        match out.last_mut() {
            Some(l) if t - l.1 <= MAX_GAP + 1 => { l.1 = t; l.2 = l.2.max(z); }
            _ => out.push((t, t, z)),
        }
    }
    out
}

fn synth(n: usize, seed: u64) -> Vec<f64> {
    let mut r = Rng(seed);
    (0..n).map(|t| 100.0 + 40.0 * (std::f64::consts::TAU * t as f64 / P as f64).sin() + 3.0 * r.normal()).collect()
}

fn main() {
    let mut x = synth(30 * P, 42);
    let truth = [(12 * P + 50, 12 * P + 50), (18 * P + 20, 18 * P + 49), (25 * P + 100, 25 * P + 107)];
    x[truth[0].0] += 30.0; // spike (10 sigma)
    x[truth[1].0..=truth[1].1].iter_mut().for_each(|v| *v += 24.0); // level shift, 5 h
    x[truth[2].0..=truth[2].1].iter_mut().for_each(|v| *v -= 30.0); // dip, 80 min
    let inc = detect(&x);
    for i in &inc { println!("incident {}..{} peak |z| {:.1}", i.0, i.1, i.2); }
    for t in truth { assert!(inc.iter().any(|i| i.0 <= t.1 && t.0 <= i.1), "missed {:?}", t); }
    let false_inc = inc.iter().filter(|i| !truth.iter().any(|t| i.0 <= t.1 && t.0 <= i.1)).count();
    assert!(false_inc <= 1, "false incidents: {false_inc}");
    println!("found {}/3, false incidents {false_inc}", 3);

    let big = synth(1_000_000, 7); // throughput, clean data
    let t0 = Instant::now();
    let n = detect(&big).len();
    let dt = t0.elapsed().as_secs_f64();
    println!("1M points: {:.2} s, {:.0} points/s, {n} incidents", dt, 1e6 / dt);
}

Output: three incidents (peaks |z| 8.8, 7.9, 9.5) at the injected places, no false incident. A clean 1 M-point series gives 2 flagged incidents, against 0.57 expected for ideal Gaussian residuals with a known scale.

Complexity and constants

• Per point: a D-element selection, then two selections over w residuals: O(w); memory 8 bytes per history value plus w residuals.
• Measured on this Mac (Apple M4 Pro), release build, one thread, D = 7, w = 288, 1 M points: 0.95-0.99 s, about 1.0 M points/s; the scale step dominates.

Efficiency

• Cost: about 1 microsecond per point (measured, Apple M4 Pro, one core, w = 288), compute-bound in the two selections; 1 M one-minute series need about 1 core-second per minute of data (arithmetic).
• Levers:
  • Refresh median/MAD every k points — k = 12 gave about 9 M points/s, 9x faster, same 3/3 incidents (measured here) — slowly drifting scale.
  • Incremental running median (two heaps or order-statistic tree) — O(w) to O(log w) per update (standard) — large windows.
  • Cheap pre-filter before the full score — fewer full evaluations (no figure) — fleets of mostly quiet series.
• Naive → standard → state of the art:
  • Global mean + 3 sigma, no seasonality · O(1) per point, floods or misses
  • Per-phase median + rolling MAD, merged incidents · about 1 microsecond per point (this page)
  • Forecast or decomposition model, calibrated threshold, FDR, incident grouping · model cost per point, fewer pages

Implementation notes

• Missing points are not zeros: skip them and alert on missing data separately.
• Weekly cycles: take D weeks of the same weekday; give holidays their own phase or exclude them.
• Integer counts: MAD can be 0 (floor the scale as the code does); take sqrt or log of skewed counts and latencies first.
• No alerts during warm-up (D cycles plus w/2 points); alert one-sided where only one tail matters.

Pitfalls

• Scale computed with the current point in the window, or from the standard deviation: the anomaly hides itself.
• Paging on every flagged point of a long incident instead of merging.
• A baseline that adapts so fast it absorbs a level shift within a few cycles.
• Ignoring that thousands of series multiply the false-alarm rate.
• Evaluating with point adjustment (below).

Labels and evaluation. Score events: an incident overlapping a labelled event is a hit, other incidents are false alarms; report detection delay. Tatbul et al. (2018) define range-based precision and recall. Point adjustment (a labelled segment counts as detected once one of its points is flagged) can make a random score look state-of-the-art (Kim et al. 2022). Wu and Keogh (2021) report four flaws in common benchmarks: triviality (one line of code solves many), unrealistic anomaly density, mislabelled ground truth, run-to-failure bias (anomalies near the end); they propose the UCR Time Series Anomaly Archive. Beat this page's recipe before claiming progress.

Variants and alternatives

In production

Cloud monitors model an expected range per metric and alarm on departures (Amazon CloudWatch anomaly detection); Twitter built S-H-ESD for service metrics. The same chain runs on payment volumes, factory sensors, utility meters and patient vitals.

Implementations

Sources

• Jordan Hochenbaum, Owen S. Vallis, Arun Kejariwal, Automatic Anomaly Detection in the Cloud Via Statistical Learning, arXiv 2017. paper · code
• Renjie Wu, Eamonn J. Keogh, Current Time Series Anomaly Detection Benchmarks are Flawed and are Creating the Illusion of Progress, arXiv 2020. paper
• Siwon Kim, Kukjin Choi, Hyun-Soo Choi, Byunghan Lee, Sungroh Yoon, Towards a Rigorous Evaluation of Time-series Anomaly Detection, AAAI 2022. paper
• Nesime Tatbul, Tae Jun Lee, Stan Zdonik, Mejbah Alam, Justin Gottschlich, Precision and Recall for Time Series, NeurIPS 2018. paper
• Varun Chandola, Arindam Banerjee, Vipin Kumar, Anomaly detection: A survey, ACM Computing Surveys 41(3), 2009. paper
• Christophe Leys, Christophe Ley, Olivier Klein, Philippe Bernard, Laurent Licata, Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median, J. Experimental Social Psychology 49(4), 2013. paper
• Bernard Rosner, Percentage Points for a Generalized ESD Many-Outlier Procedure, Technometrics 25(2), 1983. paper
• Yoav Benjamini, Yosef Hochberg, Controlling the False Discovery Rate, J. Royal Statistical Society B 57(1), 1995. paper
• Amazon CloudWatch User Guide, Using CloudWatch anomaly detection. docs