Reservoir Sampling and Priority Sampling on a Single Pass

A uniform sample of a stream is the simplest summary that supports arbitrary queries: any statistic that can be computed from a sample can be computed from a uniform one, and the error of the result follows from the sample size rather than from the query. The difficulty is that the sample must be uniform over a stream whose length is not known in advance and whose elements cannot be revisited. One pass, $k$ slots, no knowledge of $n$.

The exchange argument behind Algorithm R

Fill the first $k$ elements into the reservoir. For each subsequent element, the $i$-th in arrival order, keep it with probability $k/i$; if kept, replace a uniformly chosen slot.

The proof that the result is uniform is an induction on the number of elements seen. Assume after $i-1$ elements that each of them is in the reservoir with probability $k/(i-1)$, and each slot is filled. After element $i$ arrives, it is kept with probability $k/i$ and, if kept, occupies each slot with probability $1/k$, so each slot holds it with probability $(k/i)(1/k) = 1/i$, which is the required probability for the $i$-th element. For any earlier element $j$, the probability it survives is its old probability times the probability it is not evicted:

$$ \frac{k}{i-1}\left(1 - \frac{k}{i} \cdot \frac{1}{k}\right) = \frac{k}{i-1} \cdot \frac{i-1}{i} = \frac{k}{i} $$

The induction closes, so every element is in the reservoir with probability $k/n$ at the end. No assumption is made about $n$, and the algorithm never stores more than $k$ elements.

The cost is one random number per element, which at high arrival rates is the dominant cost. Vitter’s Algorithm Z skips this: instead of testing every element, it draws the number of elements to skip before the next replacement from a distribution derived from the acceptance probabilities, using rejection sampling to avoid expensive transcendental evaluations. The skip distribution is exact, so the result is still a uniform sample, but the expected number of random draws is proportional to the number of replacements rather than the number of elements — roughly $k \ln(n/k)$ draws instead of $n$. For a stream of a billion elements with a thousand-slot reservoir, that is the difference between a billion draws and about fourteen thousand.

import random

def reservoir(stream, k, rng=random):
    res = []
    for i, item in enumerate(stream):
        if i < k:
            res.append(item)
        else:
            j = rng.randrange(i + 1)     # uniform in [0, i]
            if j < k:
                res[j] = item
    return res

The loop above is Algorithm R written in the compact form where the decision and the slot choice are one draw: element $i$ replaces slot $j$ exactly when $j < k$, which happens with probability $k/(i+1)$, and the slot is uniform. This is the version to write when clarity matters; the skipping version is the one to write when the stream is fast.

Sampling with weights is a different problem

Many streams carry a weight per element — a flow’s byte count, an event’s monetary value — and the goal is a sample that supports estimation of the total. Uniform sampling of elements does not work: an unweighted sample of $k$ flows estimates the count of flows, not their total bytes, and the variance of the total estimate is dominated by whichever heavy flow happened to be included or excluded.

Priority sampling assigns each element $i$ with weight $w_i$ a priority $p_i = u_i / w_i$ where $u_i$ is uniform in $(0,1]$, and keeps the $k$ elements with the smallest priorities. The threshold $\tau$ is the $(k+1)$-th smallest priority. Each kept element is assigned the estimate $\hat{w}i = \max(w_i, 1/\tau)$, and the total is estimated as $\sum{i \in S} \hat{w}_i$.

Two properties make this work. First, the estimator is unbiased for the total: $\mathbb{E}[\sum_{i \in S} \hat{w}_i] = \sum_i w_i$, for any weights. Second, the estimates for individual elements are independent given the threshold, which is what allows confidence intervals to be computed from the sample alone. The second property is unusual; most sampling schemes produce correlated estimates, and the independence is the reason priority sampling is used for traffic measurement rather than simpler schemes.

The estimator’s variance is bounded, up to a constant, by $\sum_i w_i^2 / k$, so it degrades as the weights concentrate. An element heavier than the threshold is kept with certainty, so the heavy end of the distribution is measured exactly and contributes no variance at all. The noise comes from the elements below the threshold: each is kept with probability proportional to its weight, and every kept element is rescaled by $1/\tau$, so a sample that happens to include an unusually heavy tail element overestimates the total and one that omits it underestimates. Over the random priorities those two cases cancel in expectation, which is exactly the sense in which the estimator is unbiased.

What a sample can and cannot answer

A uniform sample of size $k$ from $n$ elements estimates a proportion $p = c/n$ with standard error approximately $\sqrt{p(1 - p)/k}$, and the count estimate $n\hat{p}$ carries $n$ times that. For a proportion near one half and $k = 1000$, the error on the proportion is about 1.6 percent. For a rare event with $c = 10$, the sample contains about $10k/n$ of them, which for large $n$ is zero — and no estimator can recover a count from an empty sample.

This is the fundamental limit of sampling-based summaries and the reason sketches exist. Sampling preserves the distribution of the stream; a sketch preserves a statistic of it. When the query is known in advance and is a linear statistic, a sketch is dramatically more space-efficient. When the query is arbitrary — “what fraction of requests came from mobile clients” — sampling is the only structure that answers it without prior knowledge, and the answer’s quality depends on the sample size relative to the query’s selectivity.

The same limit explains why sampling is not mergeable in the naive sense. Two uniform samples of two halves of a stream cannot be combined into a uniform sample of the whole without knowing the halves’ lengths, and even with the lengths, combining requires a weighted merge that is uniform only if the sampling fractions are known exactly. Priority sampling merges correctly: take the union of the two samples, keep the $k$ smallest priorities, and recompute the threshold. The merged sample is a priority sample of the concatenation, which is the property that makes it usable in a distributed setting.

Stratification and the assumptions it needs

A sample can be made more precise by stratifying: partition the stream by a key, sample each partition separately, and combine the estimates with the known stratum sizes. If the strata are more homogeneous than the whole, the variance falls.

The catch is that stratification requires knowing the stratum sizes, and in a stream those are not known until the stream ends. The standard approach is to maintain a per-stratum reservoir and combine at query time using the observed per-stratum counts, which are themselves estimates from the sample. The combination is unbiased only if the stratum sizes used in the weighting are the true ones, and using sample-derived weights introduces a small bias that is usually negligible but is not zero.

Stratification also multiplies the memory: $s$ strata with $k$ slots each cost $sk$. A common compromise allocates slots proportionally to expected stratum size, which requires a prior distribution over strata and degrades when the prior is wrong.

Where sampling fails quietly

Three failure modes are easy to miss.

The stream is not a population. Sampling assumes the elements are exchangeable — that the sample’s distribution matches the stream’s. A stream sorted by key, or one where the arrival order correlates with the value being measured, breaks this. Sampling the first $k$ elements of a stream sorted by timestamp gives a sample of the earliest events, not a sample of the stream. Any sampling scheme that uses arrival position as a randomiser is only valid if arrival order is independent of the measured value.

The randomiser is weak. Algorithm R consumes one random number per element, and at high rates implementations substitute a cheap hash of the element or a linear congruential generator with a short period. A generator with period $2^{32}$ used on a stream of $10^{10}$ elements repeats, and the repetitions correlate the inclusion decisions across the stream. The sample then over- or under-represents whatever the generator’s cycle happens to favour.

The sample is not the summary anyone queries. A reservoir sample supports “what fraction of items satisfy predicate P” but not “how many distinct items satisfy P”, because distinct counting needs the identities that were dropped. For distinct counts under a predicate, the right structure is a per-predicate HyperLogLog, not a sample.

Sampling versus sketching, decided by the query

The choice is not a matter of taste. Three questions decide it.

If the query is a sum or count over the whole stream, a linear sketch is smaller by orders of magnitude. The Tug-of-War construction in frequency moments estimates the sum of squares in a few hundred counters; a sample that estimates the same quantity to the same precision would need counters proportional to the stream’s skew.

If the query is a top-$k$ list, a counter-based summary is the right structure, as described in heavy hitters. A sample of size $k$ will typically miss the heavy hitters entirely, because a heavy hitter is one item among many and a uniform sample of $k$ items rarely contains a specific one.

If the query is arbitrary and unknown in advance, sampling is the only option, and the sample size is the only tuning parameter. The sample is the most general summary that fits in fixed space, and its generality is paid for in precision on any specific question.

Priority sampling sits between the two: it answers subset-sum queries — totals over an arbitrary subset of keys, where the subset is chosen after the sketch is built — with a guarantee that a pure sketch cannot provide and a size that a sample cannot match for the same precision. That combination is why it appears in traffic measurement, where the query set is chosen by an operator long after the packets have been summarised.