Moving Percentile and Median, in depth: rank definitions, sorted windows, two heaps with lazy deletion, Fenwick ranks and when to approximate

Define the window and the percentile

Fix the window first. A count window holds the last k values. A time window holds everything from the last T seconds, so its size varies. Most of this article uses count windows, and the last section adapts the structures to time windows.

Then fix the definition of a percentile, because there are several and they disagree on small windows. Nearest rank returns the value at rank ceil(q k) in sorted order, counting from 1. It always returns an observed value, which suits latency, where a made-up 117.5 ms is not something any request experienced. Linear interpolation (the default in many numeric libraries) places a fractional rank q(k - 1) and blends the two neighbours. It is smoother but returns values nobody observed. For the median of an even-sized window, the usual convention is the mean of the two middle values. Pick one definition, write it down, and test the edges: k = 1, k = 2, q = 1.0 and duplicate values.

With k = 10 and q = 0.9, nearest rank picks the 9th smallest value; interpolation blends the 9th and 10th, which on a tail with one 310 ms outlier differ by tens of milliseconds.

Why a moving median, and why tails need big windows

A small latency trace shows why people want medians. The values in milliseconds are [12, 15, 11, 240, 14, 13, 16, 12, 310, 15, 14], two slow requests in an otherwise steady stream. A window of 5 gives these outputs (measured):

Each outlier lifts the mean to four times the typical value for k steps. The median does not move, because one value cannot shift the middle of five. That is the breakdown point: a median of k values tolerates up to (k - 1) / 2 outliers per window. The same robustness is a weakness for tails. The moving p90 over a window of 10 on this trace is 240 both times, because with k = 10 the p90 is just the second-largest value. A tail percentile needs a window large enough that ceil(q k) leaves several values above it, about 10 / (1 - q) for a stable estimate. That is roughly 100 values for p90 and 1,000 for p99.

Exact method 1: keep the window sorted

The baseline is to sort every window: O(k log k) per step, O(n k log k) in total. The first real improvement keeps the window sorted between steps. Insert the new value by binary search and delete the departing one the same way.

import bisect, math
from fractions import Fraction

def sorted_window_quantile(xs, k, q):
    """Nearest-rank q-quantile of each full window, via a sorted list."""
    win, out = [], []
    r = max(1, math.ceil(Fraction(str(q)) * k))   # exact ceil(q*k)
    for i, x in enumerate(xs):
        bisect.insort(win, x)
        if i >= k:
            win.pop(bisect.bisect_left(win, xs[i - k]))
        if i >= k - 1:
            out.append(win[r - 1])
    return out

The binary searches are O(log k), but insertion and deletion shift list elements, so each step is O(k). In practice that shift is a single memmove of contiguous pointers and very fast, so for windows up to about a thousand this is often the quickest exact method in Python. It also answers any percentile, or several at once, in O(1) after the update. Note the rank is computed with exact rational arithmetic. In floating point, 0.29 * 100 is 28.999999999999996, so truncating or rounding a float product can land on the wrong rank, and p99.9 can silently become p99.

Exact method 2: two heaps with lazy deletion

For the median specifically, two heaps do each step in O(log k). The max-heap lo holds the smaller half, the min-heap hi the larger, and the invariant is that lo has the same number of live elements as hi, or one more. The median is then the top of lo for odd k, and the mean of the two tops for even k.

The difficulty is deletion. Heaps cannot remove an arbitrary element cheaply. The standard trick is lazy deletion: record the departing value in a dead table, adjust the live count of the heap it belongs to, and physically pop it only when it reaches a top. Because the median only ever reads the tops, buried dead entries are harmless.

import heapq
from collections import defaultdict

class SlidingMedian:
    """Exact median of the last k values. Two heaps plus lazy deletion."""
    def __init__(self, k):
        self.k = k
        self.lo, self.hi = [], []          # lo: max-heap via negation, hi: min-heap
        self.n_lo = self.n_hi = 0          # live counts, excluding pending deletions
        self.dead = defaultdict(int)       # value -> pending deletions

    def _prune(self, heap, sign):
        while heap and self.dead[sign * heap[0]]:
            self.dead[sign * heap[0]] -= 1
            heapq.heappop(heap)

    def _rebalance(self):
        if self.n_lo > self.n_hi + 1:
            heapq.heappush(self.hi, -heapq.heappop(self.lo))
            self.n_lo -= 1; self.n_hi += 1
            self._prune(self.lo, -1)
        elif self.n_lo < self.n_hi:
            heapq.heappush(self.lo, -heapq.heappop(self.hi))
            self.n_hi -= 1; self.n_lo += 1
            self._prune(self.hi, 1)

    def add(self, x):
        if not self.lo or x <= -self.lo[0]:
            heapq.heappush(self.lo, -x); self.n_lo += 1
        else:
            heapq.heappush(self.hi, x); self.n_hi += 1
        self._rebalance()

    def remove(self, x):
        self.dead[x] += 1
        if x <= -self.lo[0]:
            self.n_lo -= 1
            if x == -self.lo[0]: self._prune(self.lo, -1)
        else:
            self.n_hi -= 1
            if self.hi and x == self.hi[0]: self._prune(self.hi, 1)
        self._rebalance()

    def median(self):
        if self.k % 2: return float(-self.lo[0])
        return (-self.lo[0] + self.hi[0]) / 2

Two subtleties. First, remove decides which heap a value lives in by comparing with lo's top. That is correct because every value in lo is at most that top and every value in hi is at least it, and duplicates equal to the top can be counted against lo without harm. Second, the lazy entries are not bounded by k. A value buried deep in a heap can stay there long after it left the window, and on adversarial input, such as a slowly rising trend, heaps can grow toward O(n). Production code tracks the physical heap size and rebuilds both heaps from the live window when it exceeds, say, 2k. The rebuild costs O(k) and happens at most once every k steps, so the amortised cost stays O(log k).

The heap approach generalises to any fixed percentile by keeping lo at ceil(q k) live elements instead of half. Each structure then serves one q. If you need p50, p95 and p99 together, the sorted window or the Fenwick tree is simpler. For the heap mechanics themselves, see priority queues via binary heaps.

Exact method 3: a Fenwick tree over value ranks

When values come from a small, known domain, such as latencies in whole milliseconds capped at a few thousand, or quantised sensor readings, a Fenwick tree over value ranks gives every percentile in O(log V), where V is the number of distinct values, independent of k. The tree stores a count per value. Adding or removing a sample is a point update, and the r-th smallest value is found by descending the tree's implicit binary structure.

class Fenwick:
    def __init__(self, n): self.n, self.t = n, [0] * (n + 1)
    def add(self, i, d):
        i += 1
        while i <= self.n: self.t[i] += d; i += i & -i
    def kth(self, k):                        # smallest index with prefix count >= k
        pos, step = 0, 1 << self.n.bit_length()
        while step:
            nxt = pos + step
            if nxt <= self.n and self.t[nxt] < k: pos = nxt; k -= self.t[nxt]
            step >>= 1
        return pos

def fenwick_window_quantile(xs, k, q):
    vals = sorted(set(xs)); rank = {v: i for i, v in enumerate(vals)}
    fw, out = Fenwick(len(vals)), []
    r = max(1, math.ceil(Fraction(str(q)) * k))   # exact ceil(q*k)
    for i, x in enumerate(xs):
        fw.add(rank[x], 1)
        if i >= k: fw.add(rank[xs[i - k]], -1)
        if i >= k - 1: out.append(vals[fw.kth(r)])
    return out

This version compresses ranks offline from the whole input. In a live stream, fix the domain up front instead, for example one bucket per millisecond up to a cap plus an overflow bucket, which turns the structure into an exactly maintained histogram. The Fenwick tree article explains the descent in kth.

Measured: which method wins where

All three exact methods were run on 200,000 uniform random integers from 1 to 1,000, computing the moving median at three window sizes. The machine was a laptop running CPython 3.13, single run, so read the ratios rather than the absolute numbers:

The shapes match the analysis. The sorted list is fastest for small windows and degrades linearly in k once the memmove dominates. Two heaps are nearly flat, growing with log k. The Fenwick tree does not depend on k at all, only on V. In a compiled language the constant factors shift and heaps or an order-statistic tree win sooner, but the crossover logic is the same.

Time windows, sketches and aggregation

Time windows. Keep a FIFO queue of (timestamp, value) pairs. On each arrival, add the new value to the structure, then pop and remove every queued value older than T. Everything above still works, but k now varies, so recompute the target rank ceil(q k) on every query rather than caching it. The same expiry pattern drives the sliding window minimum. That case is easier, because a monotonic deque can discard dominated values, which no percentile structure can.

When to approximate. Exact methods store the whole window. For a p99 over an hour of traffic, or across thousands of series, store a sketch instead. A bucketed histogram with logarithmic buckets bounds the relative error. A t-digest is accurate at the tails and mergeable. Ring buffers of per-minute sketches give a sliding hour by merging 60 of them.

Percentiles do not compose. The average of per-host p99s is not the fleet p99, and the mean of per-minute medians is not the hourly median. Merge histograms or sketches, or recompute from raw values. Never aggregate percentiles directly.

Failure modes

• Unbounded lazy entries. The two-heap version without a rebuild threshold leaks memory on trending input. Monitor the physical heap size against k.
• NaN in the stream. NaN compares false with everything, so it corrupts heap and bisect invariants silently. Filter or map it before insertion.
• Mismatched definitions. The dashboard uses interpolation and the alert uses nearest rank, so the alert fires at a value the dashboard never shows.
• Tail percentile on a tiny window. A p99 over 50 samples is just the maximum. Size the window to about 10 / (1 - q) or report a lower percentile.
• Floating-point rank errors. A computed q times k that lands just above an integer rounds up to the wrong rank. Use integer or rational arithmetic for ranks.

Trade-offs

What to do next

1. Write down your percentile definition (nearest rank or interpolation, even-window median rule) and use it in every dashboard and alert.
2. Implement a brute-force reference and a randomised test that compares your structure against it on small windows with many duplicates.
3. Pick a structure with the decision path above, then benchmark it at your real k and value range rather than trusting the asymptotics.
4. If you use two heaps, add the rebuild threshold and a metric for physical heap size.
5. Size tail windows to about 10 / (1 - q) samples, and switch to mergeable sketches before aggregating percentiles across hosts or time.
6. Read median of medians for the one-shot selection problem behind all of this.