· Development

Picking k random items from a stream of unknown length: reservoir sampling (Algorithm R)

Reservoir sampling picks k random items in one pass over data of unknown size. Why each item enters with probability k/i, proofs of k/n and uniform k-subsets, and a common off-by-one.

DevelopmentAlgorithms

Say you want 100 random lines from a log file several gigabytes long, and you don’t know how many lines it has. Counting first and then picking means reading the file twice. If the input is standard input or a network stream, you can’t read it twice at all. Loading everything into memory and shuffling gets more expensive as the input grows.

Reservoir sampling is built for this. It reads the input once and needs only an array of k slots. Whenever the input ends, the k items in that array are a random sample.

Algorithm R

NIST’s Dictionary of Algorithms and Data Structures (DADS) defines it like this:

  1. Put the first k items into an array of size k (the reservoir).
  2. For the j-th item (j > k), choose a random integer M from 1 to j, inclusive.
  3. If M ≤ k, replace slot M of the array with the j-th item. Otherwise discard it.

The j-th item enters the reservoir with probability k/j, and if it does enter, the slot it overwrites is uniform over the k slots. One draw of M decides both. DADS describes this as an in-memory variant of Knuth’s Algorithm 3.4.2R and notes that Knuth credits Alan G. Waterman. This post calls it Algorithm R for short.

In JavaScript, with zero-based indexes, we draw m from 0 up to but not including seen instead of M, and replace when m < k:

// Uniform integer in [0, n), for 1 ≤ n ≤ 2^32
function randomBelow(n) {
  const limit = Math.floor(2 ** 32 / n) * n;
  const buf = new Uint32Array(1);
  do {
    crypto.getRandomValues(buf);
  } while (buf[0] >= limit);
  return buf[0] % n;
}

function reservoirSample(iterable, k) {
  if (!Number.isInteger(k) || k < 0) {
    throw new RangeError("k must be a non-negative integer");
  }
  const reservoir = [];
  if (k === 0) return reservoir;

  let seen = 0;
  for (const item of iterable) {
    seen++;
    if (seen <= k) {
      reservoir.push(item);
      continue;
    }
    if (seen > 2 ** 32) throw new RangeError("input exceeds randomBelow's range");
    const m = randomBelow(seen); // 0 ≤ m < seen
    if (m < k) reservoir[m] = item;
  }
  return reservoir;
}

randomBelow throws away any 32-bit value at or above floor(2^32 / n) × n and draws again. The number of values left is a multiple of n, so % n maps the same number of them to each result from 0 to n−1. It can’t handle n above 2^32; if your input might be longer than that, you need a wider random integer.

Each item costs constant work, so time is O(n) and extra memory is O(k). One random integer is drawn for each item after the k-th, so for k ≥ 1 that’s max(n − k, 0) draws in total. An input of k items or fewer needs no draws, and with k = 0 the function returns before reading anything, so it draws none either.

Why every item stays with probability k/n

Suppose the input ends after n items, with 1 ≤ k ≤ n. The calculation that follows is this post’s own. (With k = 0 the function returns an empty array, so when n > 0 every item stays with probability 0; when n = 0 there are no items to assign a probability to.)

First, the chance that a given reservoir slot gets overwritten while processing item j (j > k): M has to equal that slot’s number, so it’s 1/j. Whatever sits in that slot survives the step with probability 1 − 1/j = (j − 1)/j.

One of the first k items goes in with probability 1, then has to survive every step from k+1 through n:

(k/(k+1)) × ((k+1)/(k+2)) × … × ((n−1)/n) = k/n

Each numerator cancels the previous denominator, leaving k/n.

The i-th item (i > k) first enters with probability k/i, then has to survive steps i+1 through n. By the same cancellation that’s i/n:

(k/i) × (i/n) = k/n

Every item ends up in the final reservoir with the same probability, k/n. This holds wherever the input stops, which is why n never has to be known in advance.

Equal per-item odds aren’t the whole story

Each item having probability k/n is not enough to call a sample uniform. Take four items, split them into the pairs {1, 2} and {3, 4}, and flip a coin to pick one pair. Every item is chosen with probability 1/2, yet a set like {1, 3} can never come out.

DADS states that with Algorithm R every possible selection of k items is equally likely. Induction confirms it; again, the calculation is this post’s. One step below divides by k, so the proof assumes 1 ≤ k ≤ n. For k = 0 the only possible subset is the empty one, and the code above always returns an empty array, so there is nothing to prove. Assume that after t items (t ≥ k), each of the C(t, k) k-subsets of the first t items has probability 1/C(t, k), and process item t+1. The steps below use C(t+1, k) = C(t, k) × (t+1)/(t+1−k).

  • A subset S without the new item arises only if the previous reservoir was S and the new item was discarded: (1/C(t, k)) × (t+1−k)/(t+1) = 1/C(t+1, k).
  • A subset S with the new item: call the other k−1 items T. The previous reservoir must have been T plus one more item x, and there are t − (k−1) candidates for x. The new item must enter (k/(t+1)) and land on x’s slot (1/k), for a total of (t−k+1) × (1/C(t, k)) × (1/(t+1)) = 1/C(t+1, k).

Right after the first k items there is only one possible subset, so the assumption holds, and each step preserves it. When the input ends, all C(n, k) subsets are equally likely.

Common mistakes and edge cases

Off by one in the range. In a loop counting from zero, index i is item number i+1. Calling randomBelow(i) instead of randomBelow(i + 1) makes the entry probability k/i instead of k/(i+1). Assume 1 ≤ k < n. At i = k, the (k+1)-th item, randomBelow(k) always returns something below k, so that item always gets in. Under the same assumption, carrying the calculation through as in the previous section (this post’s own) gives the first k items a final probability of (k−1)/(n−1) and every other item k/(n−1). The first k items are under-sampled. With k = 0 the faulty code calls randomBelow(0) on the very first item. That is outside the range stated in the helper’s comment (1 ≤ n), so the call is simply invalid and there are no probabilities to talk about. With n ≤ k the faulty call is never reached, and the formulas above don’t apply. The code above sidesteps the confusion by incrementing seen first and calling randomBelow(seen) with a one-based count.

Fewer than k items. The function above returns everything it saw, an array shorter than k. With a stream of unknown length you can’t rule this out beforehand, so the caller has to decide whether it’s an error or an acceptable short result. For comparison, Python’s random.sample(population, k) raises ValueError if k is larger than the population.

k of zero or less. For k = 0 the code returns an empty array without reading the input, and it rejects negative or non-integer values with RangeError. Other policies are fine, but pick one on purpose. Leave a negative k unchecked and seen <= k is never true, so you quietly get an empty result.

Order of the result. Which items end up in the reservoir is uniform; their order inside the array is not random. With k or fewer items the output is in input order, and any of the first k items that survive stay in their original slots. If order matters downstream, shuffle the result with Fisher–Yates.

Quality of the random numbers. The proof assumes M is exactly uniform over 1 to j. Multiplying a random float by j and rounding down needs a closer look at precision as j grows. Python’s documentation says its default generator, the Mersenne Twister, is deterministic and completely unsuitable for cryptographic purposes, and points to the secrets module for security uses. If someone predicting your sample would be a problem, look at the random source before the algorithm.

When you know the length, or need it faster

If all the items are already in memory and you know how many there are, you don’t need a reservoir. Run Fisher–Yates for k steps, or in Python use random.sample. Note that random.sample takes a sequence as its population, and since Python 3.11 it no longer converts sets to lists automatically. To use it on a stream you’d have to turn the whole thing into a list first, which defeats the one-pass setting.

At the other extreme, when n is much larger than k, most items cost a random draw only to be discarded. DADS points to Jeffrey Scott Vitter’s 1985 paper “Random Sampling with a Reservoir” (ACM TOMS 11(1)) and explains that its Algorithms X, Y and Z use far fewer random numbers by choosing how many items to skip instead of deciding item by item. That’s the direction to look in if random-number generation is your bottleneck.

Sources

All checked on 2026-10-13.

Back to the blog