Your stream produces 500,000 events per second. You need a statistically valid sample. You cannot store the stream first.
This is the constraint that defines stream algorithms. The data does not fit. You get one pass. You must answer questions about the whole while holding only a piece of it.
Reservoir sampling is the canonical answer. It maintains a fixed-size sample from a stream of unknown length. At any moment, every element seen so far has an equal probability of being in the sample. No buffering. No two passes.
The setting: measurements you cannot store
A metrics service receives 500,000 HTTP request events per second. Each event carries a latency measurement. The team wants P99 latency — the value that 99% of requests are faster than.
Storing every event for an exact P99 would require petabytes per day. Computing P99 from a sample is acceptable. The question is not whether to sample. The question is how to sample without bias when you have never seen the full stream.
Stream algorithms process each element exactly once and maintain bounded state. They cannot revisit previous elements. Three constraints define them:
- One pass: each element processed once, in order.
- Bounded memory: state does not grow with stream length.
- Approximate results: exact answers often require unbounded memory; stream algorithms return approximate answers with error bounds.
Any algorithm that maintains exact state for a stream of length N requires O(N) memory. Stream algorithms break this by accepting bounded error in exchange for O(1) or O(k) memory.
Reservoir sampling: the algorithm
Given a stream of unknown length, maintain a sample of exactly k elements. Each element must have an equal probability of being in the sample at any time.
Algorithm R, from Vitter, 1985:
- Fill the reservoir with the first k elements.
- For each subsequent element i (i > k), generate a random integer j in [0, i).
- If j < k, replace reservoir[j] with element i. Otherwise, discard element i.
import random
def reservoir_sample(stream, k: int) -> list:
reservoir = []
for i, element in enumerate(stream):
if i < k:
reservoir.append(element)
else:
j = random.randint(0, i)
if j < k:
reservoir[j] = element
return reservoir
The correctness argument is short. At step i, each previously seen element sits in the reservoir with probability k/i. Element i replaces a reservoir element with probability k/(i+1). By induction, after N elements, each element occupies the reservoir with probability k/N. That is uniform sampling.
The sample gives you percentiles. Sort the reservoir. Read the 99th percentile at index 0.99 * k. For k=10,000 latency measurements, the P99 estimate has roughly 1% error. That is usually acceptable for monitoring. Storing 10,000 values instead of hundreds of millions decides the design.
Weighted reservoirs
If elements have different weights — some measurements matter more than others — Algorithm A-Chao generalises the trick. Generate a random key per element proportional to its weight. Keep the k elements with the highest keys in a min-heap. Elements with higher weight have higher expected keys and land in the reservoir more often.
import heapq, math, random
def weighted_reservoir_sample(stream, k, weight_fn=lambda x: 1.0):
heap = []
for element in stream:
w = weight_fn(element)
if w <= 0:
continue
key = random.random() ** (1.0 / w)
if len(heap) < k:
heapq.heappush(heap, (key, element))
elif key > heap[0][0]:
heapq.heapreplace(heap, (key, element))
return [e for _, e in heap]
Same one-pass discipline. Same bounded memory. Different distribution over the sample.
Welford's algorithm: running mean and variance
Reservoir sampling holds a sample. Some statistics do not need one. Mean and variance can be maintained exactly with O(1) state.
Welford's online algorithm updates both incrementally:
class RunningStats:
def __init__(self):
self.n = 0
self.mean = 0.0
self._M2 = 0.0
def update(self, value: float) -> None:
self.n += 1
delta = value - self.mean
self.mean += delta / self.n
delta2 = value - self.mean
self._M2 += delta * delta2
@property
def variance(self) -> float:
return self._M2 / (self.n - 1) if self.n > 1 else 0.0
Welford's algorithm is numerically stable. It does not accumulate floating-point errors the way the naive formula does — compute sum and sum-of-squares, then subtract. That subtraction loses precision on large streams. Production monitoring systems like Prometheus and Grafana use running statistics rather than storing raw measurements.
The running median: two heaps in balance
Some questions need order statistics, not sums. What is the median request latency, updated every second, on a stream that never ends? What is the median transaction value in the last hour, at 50,000 transactions per second?
Sorting the full window per query is O(N log N). Doing that every second on a million-value window is a small fortune in CPU.
The exact answer is O(log N) per insertion with two heaps. A max-heap holds the lower half — its root is the largest small value. A min-heap holds the upper half — its root is the smallest large value. The median sits at one of the roots, or is the average of both.
import heapq
class RunningMedian:
def __init__(self):
self._lower = [] # max-heap (negated)
self._upper = [] # min-heap
def add(self, x: float) -> None:
if not self._lower or x <= -self._lower[0]:
heapq.heappush(self._lower, -x)
else:
heapq.heappush(self._upper, x)
if len(self._lower) > len(self._upper) + 1:
heapq.heappush(self._upper, -heapq.heappop(self._lower))
elif len(self._upper) > len(self._lower):
heapq.heappush(self._lower, -heapq.heappop(self._upper))
def median(self) -> float:
if len(self._lower) > len(self._upper):
return -self._lower[0]
return (-self._lower[0] + self._upper[0]) / 2
The two-heap trick generalises. K-th smallest of a stream is a max-heap of size k, rejecting new values larger than the root. LLM decoders that use top-K sampling maintain a bounded heap of candidate tokens for exactly this reason. Bounded-heap extreme tracking is asymptotically cheaper than repeated sort.
The tradeoffs
AT2 — Latency vs Throughput. Reservoir sampling collapses per-event work to O(1): one random draw, one comparison, one possible write. That is what lets the system absorb 500,000 events per second. Exact computation of P99 needs a full sort at query time — high latency per query, and query latency that grows with stream length. Sampling flattens both.
AT9 — Correctness vs Performance. Exact P99 stores every measurement. Reservoir sampling with k samples estimates P99 with error roughly O(1/√k). At k=10,000, that is about 1%. For a monitoring dashboard, 1% error is invisible. For a billing calculation, it is not.
AT4 — Precomputation vs On-Demand. Welford's algorithm precomputes incrementally. Each new value takes O(1) to incorporate. On-demand percentile computation requires O(N log N) at query time. For statistics queried often, the incremental path always wins.
Where it fails
FM9 — Silent Data Corruption. A uniform reservoir sample answers questions about the stream distribution. But the numerical machinery around it can still lie. Welford's algorithm loses precision on pathological inputs — millions of values very close to the mean, or values with extreme variance. For financial calculations, floating-point running statistics are the wrong tool. Use exact integer arithmetic with sufficient precision.
The subtler corruption is bias hiding in a valid-looking sample. If your stream is non-stationary — a traffic spike causes elevated latency for 30 seconds — a reservoir sample taken over 24 hours dilutes the spike into invisibility. The dashboard reports a healthy P99. The 30 seconds when users saw failures do not show up. The sample was uniform. The uniformity was the problem.
The fix is a time-windowed sample: keep only the last 10 minutes. Now recent behaviour dominates, at the cost of losing long-run baselines. That is another AT4 decision, made explicit.
What this means for your system right now
Every stream-processing decision resolves the same shape. You are choosing between exact answers over the whole stream (memory you do not have) and approximate answers with bounded state (error you must own).
What error can you tolerate? For P99 latency on a dashboard, 1% is fine. For a fraud threshold that decides whether to block a transaction, 1% may not be. Match the sample size to the decision the number drives.
What time horizon do you care about? Reservoir sampling is unbiased across the whole stream. That is exactly wrong when only recent behaviour matters. Use a time-windowed sample when the distribution moves.
What breaks first when the stream grows? If per-event work is O(1) — a random draw, a heap push — the system scales with throughput. If per-query work is O(N log N) over a growing window, the dashboard slows down as the stream grows. The per-event cost is the one that matters at 500,000 events per second.
The signal that tells you this applies to your system
You cannot store the stream but you need statistics about it. Percentiles, means, variances, top-K, medians. Your queries do not need exact answers — a 1% error would not change any decision. The moment those two conditions hold together, stream algorithms are already the right tool. The only question left is which one.
The full framework treatment — compression blocks, three-level exercises, and the complete AT/FM mapping — is in Book 2: Algorithm Engineering, Chapter 16 (Reservoir Sampling and Stream Algorithms). Free chapter available at computingseries.com/books/book2.