How to Keep a Fair Random Sample of a Log Stream in Python With Reservoir Sampling and Merge Samples From Many Workers

Sooner or later every system runs into the same awkward question. You have far more data than you can keep, you can read it only once, and you still need a small, honest sample of it. A security team wants 1,000 authentication events to review out of tonight’s ten million. A data engineer wants a test fixture that looks like production. A database wants to estimate how the values in a column are spread without reading every row. In each case the sample has to be fair, meaning every item had the same chance of being picked, and it has to be built while the data streams past, because there is no second pass and nobody knows in advance how long the stream will be.

The classic answer is reservoir sampling. Wikipedia defines it as “a family of randomized online algorithms for choosing a simple random sample, without replacement, of k items from a population of unknown size n in a single pass over the items.” PostgreSQL uses it when ANALYZE samples a table’s rows. The comment in analyze.c says: “The first targrows sample rows are simply copied into the reservoir. Then we start replacing tuples in the sample until we reach the end of the relation.” It adds that “At all times the reservoir is a true random sample of the tuples we’ve passed over so far, so when we fall off the end of the relation we’re done.” Python’s standard library has no reservoir sampler. random.sample insists that “The population must be a sequence,” so you would have to hold the whole stream in memory before you could sample it.

This tutorial builds reservoir sampling in Python from the ground up and stress-tests every piece. You will start with a synthetic stream of one million authentication events that is rigged to expose weak samples, write the simplest reservoir algorithm, learn how to test a random algorithm with statistics, break it on purpose, make it fast, add weights, merge samples from several workers, keep rare events from being crowded out, and finish with a command-line tool and a test suite. By the end you will be able to:

• explain why “the first N lines”, “every Nth line” and “keep each line with probability p” are poor samples, and measure how they fail;
• write Algorithm R and explain, with a pencil, why every item ends up with the same chance;
• test a random algorithm with z scores, and see why a test on single items is not enough;
• speed the sampler up with Algorithm L and avoid two floating-point traps that no test run will ever hit by chance;
• sample in proportion to a weight with the A-Res algorithm and check it against exact probabilities;
• merge the samples of several workers into one fair sample, and see why “pool and draw again” is wrong;
• keep one reservoir per category so a rare event survives, and measure the memory you save.

Everything below was run on Python 3.13.14 and pytest 9.2.0 on Windows 11 on November 3, 2026, using only the standard library plus pytest. Every output block is copied from those runs. Timings differ from machine to machine, so compare the ratios, not the digits.

The terms used in this tutorial

• Stream. Items that arrive one at a time and can be read once: log lines, events from a queue, rows from a database cursor. You do not know in advance how many there will be.
• Sample. A small group of items picked from the stream to stand in for all of it. k is the sample size and n is how many items the stream held in the end.
• Uniform (fair) sample. A sample in which every item of the stream had the same chance, k / n, of being included. Statisticians call this a simple random sample.
• Reservoir. The list of k slots that holds the sample while the stream runs. At every moment it is a fair sample of everything seen so far.
• Algorithm R and Algorithm L. Two ways of maintaining a reservoir. R looks at every item. L skips ahead and looks only at the items that will actually be kept.
• Weighted sample. A sample in which heavier items are more likely to be picked, for example in proportion to the bytes a request sent.
• Standard error and z score. If a frequency should be p and you measure it over many runs, it wobbles around p by a typical amount called the standard error. The z score is the gap between your measurement and p, counted in standard errors. A z around 2 is ordinary noise; a z above 5 means something is really wrong.
• Seed. A number that starts a random generator, so the same seed replays the same “random” choices.

Prerequisites

• Python 3.13 (the version I tested; the code uses only the standard library). Where I write python, use python3 if that is your command name.
• pytest 9.2.0 for the tests. Nothing else is installed.
• Any operating system. I used Windows 11; the commands for Linux and macOS are shown where they differ.
• Comfort with Python loops, lists, functions and generators. No statistics background is assumed; the two statistics ideas you need are explained when they first appear.
• About an hour. The scripts themselves run in seconds, apart from two that take 10 to 15 seconds.

Step 0: Create the project folder

Put everything in one folder with its own virtual environment, so nothing touches your system Python.

mkdir reservoir-lab
cd reservoir-lab
python -m venv .venv
.venv\Scripts\activate
python -m pip install pytest==9.2.0
python -m pytest --version

On Linux or macOS, replace the activation line with source .venv/bin/activate. The last command should print pytest 9.2.0. Every file below goes into this folder, and every command is run from it.

Step 1: Meet the stream and watch three tempting samples fail

The idea

To judge a sampling method you need data whose true shape you already know. The file below generates a synthetic authentication log with three traps built in, all of them things real logs do:

• A rhythm. Every 100th event is a healthcheck from a monitoring tool.
• A burst. Between 70 and 72 percent of the way through the stream, 90 percent of the events are of a kind called stuffing, as in a credential-stuffing attack.
• A rare event. key_export events appear only at positions you choose. Step 7 uses them.

Everything is seeded, so you will get exactly the stream I did.

The code

Save this as events.py.

"""events.py: a synthetic authentication log with three traps built in (a rhythm, a burst, a rare event)."""
import random

PERIOD = 100  # every 100th event is a health check


def auth_events(n, seed=7, rare_positions=()):
    """Yield n (position, kind) pairs.

    kinds: login_ok, login_fail, mfa_prompt, password_reset,
           healthcheck (every 100th event), stuffing (a burst over 70 to 72 percent of the stream),
           key_export (only at the positions in rare_positions)
    """
    rng = random.Random(seed)
    rand = rng.random
    burst_start, burst_end = n * 70 // 100, n * 72 // 100
    rare = frozenset(rare_positions)
    for pos in range(n):
        if pos in rare:
            yield pos, "key_export"
        elif pos % PERIOD == 0:
            yield pos, "healthcheck"
        elif burst_start <= pos < burst_end and rand() < 0.9:
            yield pos, "stuffing"
        else:
            r = rand()
            if r < 0.80:
                yield pos, "login_ok"
            elif r < 0.95:
                yield pos, "login_fail"
            elif r < 0.995:
                yield pos, "mfa_prompt"
            else:
                yield pos, "password_reset"

Now a script that counts the true shape of a million-event stream and then takes three kinds of “obvious” sample of 1,000 events. Save it as step01_naive.py.

"""step01_naive.py: three tempting ways to sample a log stream, and what each one misses."""
import random
from collections import Counter
from itertools import islice

from events import auth_events

N, K = 1_000_000, 1_000

truth = Counter(kind for _, kind in auth_events(N))
print(f"the whole stream: {N:,} events")
for kind, count in truth.most_common():
    print(f"  {kind:<15}{count:>8,}  {count / N:6.2%}")
fair = {kind: truth[kind] / N * K for kind in ("stuffing", "healthcheck")}
print(f"a fair sample of {K:,} events should hold about {fair['stuffing']:.1f} stuffing and {fair['healthcheck']:.1f} healthcheck events")
print()


def show(name, sample):
    kinds = Counter(kind for _, kind in sample)
    print(f"{name:<38}{len(sample):>5} events   stuffing {kinds['stuffing']:>3}   healthcheck {kinds['healthcheck']:>4}")


show("first 1,000 events", list(islice(auth_events(N), K)))
show("every 1,000th event", list(islice(auth_events(N), 0, None, 1000)))
show("every 1,001st event", list(islice(auth_events(N), 0, None, 1001)))
for seed in (1, 2, 3):
    rng = random.Random(seed)
    sample = [event for event in auth_events(N) if rng.random() < K / N]
    show(f"keep each event with p = 0.001 (seed {seed})", sample)

Run it

python step01_naive.py

the whole stream: 1,000,000 events
  login_ok        777,722  77.77%
  login_fail      145,869  14.59%
  mfa_prompt       43,893   4.39%
  stuffing         17,722   1.77%
  healthcheck      10,000   1.00%
  password_reset    4,794   0.48%
a fair sample of 1,000 events should hold about 17.7 stuffing and 10.0 healthcheck events

first 1,000 events                     1000 events   stuffing   0   healthcheck   10
every 1,000th event                    1000 events   stuffing   0   healthcheck 1000
every 1,001st event                    1000 events   stuffing  17   healthcheck   10
keep each event with p = 0.001 (seed 1)  982 events   stuffing  22   healthcheck    7
keep each event with p = 0.001 (seed 2) 1009 events   stuffing  25   healthcheck    8
keep each event with p = 0.001 (seed 3)  949 events   stuffing  23   healthcheck   12

What the output tells you

The top half is the truth. Stuffing makes up 1.77 percent of the stream and health checks exactly 1.00 percent, so a fair sample of 1,000 events should hold about 17.7 stuffing events and 10.0 health checks. Now look at what each shortcut produced:

• The first 1,000 events contain no stuffing at all. The burst starts at event 700,000, so a sample taken from the head of the stream cannot know it exists.
• Every 1,000th event gave 1,000 events and every one is a health check. The step, 1,000, is a multiple of the 100-event rhythm, so every pick lands on the same point of the cycle. This is called aliasing.
• Every 1,001st event looks perfect: 17 stuffing events and 10 health checks. That is luck of arithmetic, because 1,001 shares no factor with 100. Nothing in your code warns you, and tomorrow’s log may have a different rhythm. Step 3 shows a deeper flaw that even a lucky step cannot escape.
• Keeping each event with probability 0.001 gave 982, 1,009 and 949 events. The size of the sample is itself random (a standard deviation of about 32 events, by the binomial formula), and you must choose the probability before you know how long the stream is. Too high and memory grows without limit; too low and you may get almost nothing.

So the target is a method that keeps exactly k items, gives every item of the stream the same chance, works in one pass without knowing n, and uses memory for k items only.

Step 2: Build Algorithm R, the simplest reservoir

The idea

Keep a list of k slots, the reservoir. Copy the first k items into it. After that, when the i-th item arrives, pick a random whole number j from 0 to i − 1. If j is smaller than k, the new item replaces whatever sits in slot j. Otherwise the item is discarded. In other words, item number i is accepted with probability k / i, and when it is accepted it evicts a randomly chosen member. This is Algorithm R. Wikipedia says it was “discovered independently by Alan G. Waterman and McLeod & Bellhouse.”

Why every item ends up with the same chance

Take item number i, with i larger than k, in a stream that finally holds n items. It enters the reservoir with probability k / i. After that it has to survive every later step j, from i + 1 up to n. At step j the newcomer is accepted with probability k / j and then picks one of the k slots, so a particular member is evicted with probability (k / j) × (1 / k) = 1 / j, and survives with probability (j − 1) / j. Multiply the survival chances for all the later steps and nearly everything cancels: (i / (i + 1)) × ((i + 1) / (i + 2)) and so on up to ((n − 1) / n) equals i / n. The item’s total chance is therefore (k / i) × (i / n) = k / n. The i cancels, and that is the whole trick. The first k items start inside the reservoir and face the same product starting at step k + 1, which gives k / n as well. Wikipedia gives the formal proof by induction; the argument here is the same idea in a form you can check with a pencil.

The code

Save this as reservoir.py. The imports and the END marker at the top are used by later steps; adding them now keeps the file tidy. reservoir_r is Algorithm R.

"""reservoir.py: reservoir sampling for streams of unknown length (Python 3.13, standard library only)."""
import heapq
import itertools
import math
import random

END = object()  # what next() returns when the stream has run out


def reservoir_r(stream, k, rng=None):
    """Algorithm R: a uniform random sample of up to k items from a stream, in one pass."""
    if k < 0:
        raise ValueError("k must be 0 or more")
    rng = rng or random.Random()
    reservoir = []
    for i, item in enumerate(stream, start=1):  # i counts this item too
        if i <= k:
            reservoir.append(item)  # fill phase: keep the first k items
        else:
            j = rng.randrange(i)  # a slot number from 0 to i - 1
            if j < k:  # true with probability k / i
                reservoir[j] = item  # the new item takes over slot j
    return reservoir

Three details are worth a second look. enumerate(stream, start=1) makes i count the current item, which the algorithm needs. rng.randrange(i) returns a whole number from 0 up to and including i − 1, so the test j < k is true with probability exactly k / i. And the function takes the random generator as an argument, which lets us feed it a seeded generator now and a rigged one in Step 4.

A way to test a random algorithm

You cannot test a random algorithm with one run. Instead you run it many times and count how often each item lands in the sample; a fair sampler includes every item with probability k / n. The helper below does that counting. It passes range(n) as the stream, so each item is its own number, and it gives the sampler a seeded generator so the experiment is repeatable. For each item it computes the z score: the gap between the observed frequency and k / n, divided by the standard error, which is the square root of p × (1 − p) / runs. With 20 items the biggest z is typically around 2; a value above 5 means a real bias. The imports at the top include Fraction and combinations for later steps.

Save this as labkit.py.

"""labkit.py: helpers for testing randomized code. A sampler is right when many runs give the right frequencies."""
import math
import random
from collections import namedtuple
from fractions import Fraction
from itertools import combinations

Result = namedtuple("Result", "low high worst_z freqs")


def inclusion_z(sampler, n, k, runs, seed):
    """Run sampler(range(n), k, rng) many times and count how often each item lands in the sample.

    A fair sampler includes every item with probability k / n. The z score says how many standard errors
    an observed frequency is from that. Returns the lowest and highest frequency and the worst |z|.
    """
    rng = random.Random(seed)
    counts = [0] * n
    for _ in range(runs):
        for item in sampler(range(n), k, rng):
            counts[item] += 1
    p = k / n
    se = math.sqrt(p * (1 - p) / runs)
    freqs = [c / runs for c in counts]
    return Result(min(freqs), max(freqs), max(abs(f - p) / se for f in freqs), freqs)

Now a script that traces Algorithm R on eight letters, with a reservoir of three slots, and then runs the statistical test. Save it as step02_algorithm_r.py.

"""step02_algorithm_r.py: watch Algorithm R work on 8 items, then test it on many runs."""
import random

from labkit import inclusion_z
from reservoir import reservoir_r

# Part 1: a trace. Eight items go by and the reservoir has three slots.
rng = random.Random(4)
k, reservoir = 3, []
for i, item in enumerate("ABCDEFGH", start=1):
    if i <= k:
        reservoir.append(item)
        action = "fill phase: keep it"
    else:
        j = rng.randrange(i)
        if j < k:
            action = f"drew j = {j}, below {k}: replaces {reservoir[j]} in slot {j}"
            reservoir[j] = item
        else:
            action = f"drew j = {j}, not below {k}: discarded"
    print(f"item {i} ({item}): {action:<44} reservoir = {reservoir}")

# Part 2: the test. 20 items, 5 slots, 100,000 runs. Every item should appear in 25 percent of the samples.
n, k, runs = 20, 5, 100_000
result = inclusion_z(reservoir_r, n, k, runs, seed=11)
print(f"\nn = {n}, k = {k}, {runs:,} runs; every item should appear with probability {k / n:.4f}")
print(f"lowest frequency {result.low:.4f}, highest frequency {result.high:.4f}, worst |z| = {result.worst_z:.2f}")

Run it

python step02_algorithm_r.py

item 1 (A): fill phase: keep it                          reservoir = ['A']
item 2 (B): fill phase: keep it                          reservoir = ['A', 'B']
item 3 (C): fill phase: keep it                          reservoir = ['A', 'B', 'C']
item 4 (D): drew j = 1, below 3: replaces B in slot 1    reservoir = ['A', 'D', 'C']
item 5 (E): drew j = 2, below 3: replaces C in slot 2    reservoir = ['A', 'D', 'E']
item 6 (F): drew j = 0, below 3: replaces A in slot 0    reservoir = ['F', 'D', 'E']
item 7 (G): drew j = 5, not below 3: discarded           reservoir = ['F', 'D', 'E']
item 8 (H): drew j = 6, not below 3: discarded           reservoir = ['F', 'D', 'E']

n = 20, k = 5, 100,000 runs; every item should appear with probability 0.2500
lowest frequency 0.2473, highest frequency 0.2519, worst |z| = 1.95

What the output tells you

Read the trace line by line. Items A, B and C fill the three slots. Item D arrives with i = 4 and draws j = 1, which is below 3, so D takes slot 1 and B is evicted. Item E draws 2 and evicts C. Item F draws 0 and evicts A. Item G draws 5 and item H draws 6, neither below 3, so both are discarded. The chance of each acceptance was 3 / i, so later items are accepted less and less often, yet the argument above says every letter had the same final chance. The test confirms it on a bigger scale: over 100,000 runs with 20 items and 5 slots, the least frequent item appeared in 24.73 percent of the samples, the most frequent in 25.19 percent, and the worst z score was 1.95. That is ordinary noise around 25 percent.

Step 3: Break it on purpose, and learn what a test can and cannot see

The idea

A statistical test earns its keep when the code is subtly wrong. Here are four samplers that look right. The first three are plausible mistakes in Algorithm R. Wikipedia’s description draws j from 1 to i, while Python’s randrange(i) draws from 0 to i − 1, so a translation that writes randrange(i + 1) or randint(0, i) draws one value too many. The fourth, systematic sampling (pick a random starting offset, then take every (n / k)th item), is not a reservoir at all: it is “every 1,001st line” with a random starting offset added to make each item’s chance fair. It also needs the whole list in memory.

To catch the fourth, the test needs one more tool: instead of asking how often each item is sampled, ask how often each pair of items is sampled together. In a fair sample of k out of n, any two items appear together with probability k(k − 1) / (n(n − 1)). Add this function to the end of labkit.py.

# labkit.py (continued)
def pair_z(sampler, n, k, runs, seed):
    """The same test for pairs: any two items should be in the sample together with probability k(k-1) / (n(n-1))."""
    rng = random.Random(seed)
    counts = dict.fromkeys(combinations(range(n), 2), 0)
    for _ in range(runs):
        for pair in combinations(sorted(sampler(range(n), k, rng)), 2):
            counts[pair] += 1
    p = k * (k - 1) / (n * (n - 1))
    se = math.sqrt(p * (1 - p) / runs)
    return max(abs(c / runs - p) / se for c in counts.values())

The code

Save this as step03_bugs.py.

"""step03_bugs.py: samplers that look right. Only a statistical test tells them apart."""
from labkit import inclusion_z, pair_z
from reservoir import reservoir_r


def off_by_one(stream, k, rng):
    reservoir = []
    for i, item in enumerate(stream, start=1):
        if i <= k:
            reservoir.append(item)
        else:
            j = rng.randrange(i + 1)  # BUG: one slot too many, so late items are accepted too rarely
            if j < k:
                reservoir[j] = item
    return reservoir


def forgot_current_item(stream, k, rng):
    reservoir = []
    for i, item in enumerate(stream, start=1):
        if i <= k:
            reservoir.append(item)
        else:
            j = rng.randrange(i - 1)  # BUG: the current item was left out of the count
            if j < k:
                reservoir[j] = item
    return reservoir


def evict_oldest(stream, k, rng):
    reservoir = []
    for i, item in enumerate(stream, start=1):
        if i <= k:
            reservoir.append(item)
        elif rng.randrange(i) < k:
            reservoir.pop(0)  # BUG: the oldest item leaves, not a randomly chosen one
            reservoir.append(item)
    return reservoir


def systematic(stream, k, rng):
    """Not a reservoir at all: a random start, then every (n / k)-th item. It needs the whole list."""
    items = list(stream)
    step = len(items) // k
    return items[rng.randrange(step)::step]


n, k, runs = 20, 5, 100_000
print(f"n = {n}, k = {k}, {runs:,} runs; a fair item frequency is {k / n:.4f}, a fair pair frequency is {k * (k - 1) / (n * (n - 1)):.4f}")
print(f"{'sampler':<22}{'lowest':>8}{'highest':>9}{'item |z|':>10}{'pair |z|':>10}  verdict")
for name, sampler in [("reservoir_r", reservoir_r), ("off_by_one", off_by_one), ("forgot_current_item", forgot_current_item),
                      ("evict_oldest", evict_oldest), ("systematic", systematic)]:
    items = inclusion_z(sampler, n, k, runs, seed=11)
    pairs = pair_z(sampler, n, k, runs, seed=12)
    verdict = "PASS" if max(items.worst_z, pairs) < 5 else "FAIL"
    print(f"{name:<22}{items.low:>8.4f}{items.high:>9.4f}{items.worst_z:>10.1f}{pairs:>10.1f}  {verdict}")

Run it

python step03_bugs.py

n = 20, k = 5, 100,000 runs; a fair item frequency is 0.2500, a fair pair frequency is 0.0526
sampler                 lowest  highest  item |z|  pair |z|  verdict
reservoir_r             0.2473   0.2519       1.9       2.4  PASS
off_by_one              0.2357   0.2874      27.3      28.3  FAIL
forgot_current_item     0.2086   0.2648      30.2      26.3  FAIL
evict_oldest            0.0001   0.3970     182.5     146.3  FAIL
systematic              0.2492   0.2515       1.1     281.8  FAIL

What the output tells you

• reservoir_r passes both tests, with z scores of 1.9 and 2.4.
• off_by_one fails clearly. The lowest and highest frequencies are 0.2357 and 0.2874 instead of 0.2500. I worked out what this bug does: the first five items end up with probability 6 / 21 = 0.2857 and every later item with 5 / 21 = 0.2381, and the measured extremes agree within noise. One wrong character, and the early items end up 20 percent more likely than the later ones.
• forgot_current_item and evict_oldest fail too. Evicting the oldest item is the worst: its least likely item appears in only 0.01 percent of the samples.
• systematic is the lesson. Its item z score is 1.1, better than the correct reservoir’s, yet its pair z score is 281.8. With n = 20 and k = 5 the step is 4, so there are only four possible samples, one per starting offset, out of the 15,504 different five-item subsets that exist. Two neighbouring items can never be sampled together. Every item has the right chance, but the sample as a whole is not fair.

The takeaway: “every item has the same chance” is necessary but not sufficient. A fair sampler can produce every subset with equal probability, which is a much stronger promise. The pair test checks a slice of that promise cheaply, and it is why a test that looks only at single items would have passed the fourth sampler.

Step 4: Make it fast with Algorithm L

The idea

Algorithm R draws a random number for every item. That is wasteful: by item 1,000,000 a reservoir of 1,000 accepts a newcomer only once in a thousand tries, and the other 999 random numbers are thrown away. Wikipedia notes that Algorithm R “has linear runtime, since it generates a random number for every element in the input stream,” and that Jeffrey Vitter showed faster algorithms are possible “when the algorithm is allowed to skip past sections of the input stream without examining them.” Kim-Hung Li described a simple one, Algorithm L.

Here is the idea in plain words. Imagine every item draws a random score between 0 and 1 and the sample is the k items with the lowest scores; that is a fair sample. Let w be the highest score currently in the sample. A new item gets in only if its score is below w, which happens with probability w. So the number of items skipped before the next one gets in follows a geometric distribution, and one random number is enough to draw it: floor(log(u) / log(1 - w)). When an item does get in, it takes the place of the member holding the highest score, which by symmetry is equally likely to be any of the k members, so the code simply picks a random slot; and w shrinks by a random factor. The algorithm never looks at the skipped items at all. Wikipedia sums up the cost: it “computes three random numbers for each item that becomes part of the reservoir, and does not spend any time on items that do not.”

The code

Add this to the end of reservoir.py. The helper open01 and the line that creates w look fussier than the textbook pseudocode; the “Two float traps” section below explains why.

# reservoir.py (continued)
def open01(rng):
    """A random float in (0, 1]. random() itself is in [0, 1), and log(0) is an error."""
    return 1.0 - rng.random()


def reservoir_l(stream, k, rng=None):
    """Algorithm L: the same kind of sample as Algorithm R, but it skips items instead of looking at each one."""
    if k < 0:
        raise ValueError("k must be 0 or more")
    rng = rng or random.Random()
    it = iter(stream)
    reservoir = list(itertools.islice(it, k))
    if k == 0 or len(reservoir) < k:
        return reservoir  # nothing to sample, or the stream was shorter than k
    w = min(math.exp(math.log(open01(rng)) / k), 1.0 - 2.0**-53)  # keep w strictly below 1
    while True:
        ratio = math.log(open01(rng)) / math.log1p(-w)  # how many items to skip, as a float
        skip = math.floor(ratio)
        item = next(itertools.islice(it, skip, skip + 1), END)
        if item is END:
            return reservoir
        reservoir[rng.randrange(k)] = item
        w *= math.exp(math.log(open01(rng)) / k)

Line by line: list(itertools.islice(it, k)) fills the reservoir with the first k items. The loop then draws a skip count, throws away that many items with itertools.islice(it, skip, skip + 1) and takes the next one (next(..., END) returns the marker END when the stream has run out). It puts the item into a random slot and shrinks w. Walking over items you do not use is cheap in Python compared with drawing a random number and running a few bytecodes per item.

To measure the savings we need a random generator that counts. We also need its opposite, a generator that always returns the smallest value random() may return, which we will use to probe the float traps. Add both to the end of labkit.py.

# labkit.py (continued)
class CountingRandom:
    """A stand-in for random.Random that counts every random number it hands out."""

    def __init__(self, seed):
        self._rng = random.Random(seed)
        self.draws = 0

    def random(self):
        self.draws += 1
        return self._rng.random()

    def randrange(self, n):
        self.draws += 1
        return self._rng.randrange(n)


class ZeroRandom:
    """The unluckiest generator possible: random() always returns 0.0, which random.random() is allowed to return."""

    def random(self):
        return 0.0

    def randrange(self, n):
        return 0

Now the experiment. Save this as step04_algorithm_l.py. Parts 1 to 3 count random numbers, time both algorithms and test fairness. Part 4 rewrites Algorithm L without the safeguards, so you can watch it fail.

"""step04_algorithm_l.py: Algorithm L keeps the same kind of sample while asking for far fewer random numbers."""
import itertools
import math
import random
import time

from labkit import CountingRandom, ZeroRandom, inclusion_z, pair_z
from reservoir import END, reservoir_l, reservoir_r

# Part 1: how many random numbers does each algorithm ask for?
n, k = 1_000_000, 1_000
for name, sampler in (("Algorithm R", reservoir_r), ("Algorithm L", reservoir_l)):
    rng = CountingRandom(5)
    sample = sampler(iter(range(n)), k, rng)
    print(f"{name}: kept {len(sample):,} items and drew {rng.draws:,} random numbers")
print(f"theory for L: about 3 * k * ln(n / k) = {3 * k * math.log(n / k):,.0f} random numbers")

# Part 2: how long does each take on 5 million items?
print()
n = 5_000_000
for name, sampler in (("Algorithm R", reservoir_r), ("Algorithm L", reservoir_l)):
    start = time.perf_counter()
    sampler(iter(range(n)), k, random.Random(5))
    print(f"{name}: {n:,} items in {time.perf_counter() - start:.2f} s")

# Part 3: is L's sample as fair as R's?
print()
for name, sampler in (("Algorithm R", reservoir_r), ("Algorithm L", reservoir_l)):
    items = inclusion_z(sampler, 20, 5, 100_000, seed=11)
    pairs = pair_z(sampler, 20, 5, 100_000, seed=12)
    print(f"{name}: item frequencies {items.low:.4f} to {items.high:.4f}, worst item |z| {items.worst_z:.2f}, worst pair |z| {pairs:.2f}")


# Part 4: two float traps. The same algorithm, written without the two safeguards.
def l_without_safeguards(stream, k, rng, open_interval, clamp):
    it = iter(stream)
    reservoir = list(itertools.islice(it, k))
    draw = (lambda: 1.0 - rng.random()) if open_interval else rng.random
    w = math.exp(math.log(draw()) / k)
    if clamp:
        w = min(w, 1.0 - 2.0**-53)
    while True:
        skip = math.floor(math.log(draw()) / math.log1p(-w))
        item = next(itertools.islice(it, skip, skip + 1), END)
        if item is END:
            return reservoir
        reservoir[rng.randrange(k)] = item
        w *= math.exp(math.log(draw()) / k)


print("\nrandom() always returns 0.0, the smallest value it is allowed to return:")
for label, call in (
    ("raw random()", lambda: l_without_safeguards(range(20), 5, ZeroRandom(), False, False)),
    ("1 - random(), no clamp on w", lambda: l_without_safeguards(range(20), 5, ZeroRandom(), True, False)),
    ("reservoir_l from reservoir.py", lambda: reservoir_l(iter(range(20)), 5, ZeroRandom())),
):
    try:
        print(f"  {label:<30} returned {call()}")
    except ValueError as err:
        print(f"  {label:<30} raised ValueError: {err}")

kk = 1000
m = 0
while math.exp(math.log(1.0 - m * 2.0**-53) / kk) == 1.0:
    m += 1
print(f"\nfor k = {kk}, w rounds to exactly 1.0 for the {m} smallest values random() can return: a chance of {m * 2.0**-53:.1e} per run")

Run it

python step04_algorithm_l.py

Algorithm R: kept 1,000 items and drew 999,000 random numbers
Algorithm L: kept 1,000 items and drew 20,609 random numbers
theory for L: about 3 * k * ln(n / k) = 20,723 random numbers

Algorithm R: 5,000,000 items in 0.82 s
Algorithm L: 5,000,000 items in 0.04 s

Algorithm R: item frequencies 0.2473 to 0.2519, worst item |z| 1.95, worst pair |z| 2.38
Algorithm L: item frequencies 0.2472 to 0.2525, worst item |z| 2.02, worst pair |z| 3.06

What the output tells you

• Random numbers. Out of a million items Algorithm R drew 999,000 random numbers, one for every item after the first thousand. Algorithm L drew 20,609, close to the theory of about 3 × k × ln(n / k) = 20,723. That is roughly 48 times fewer.
• Time. On five million items Algorithm R took 0.79 seconds and Algorithm L 0.04 seconds, about 20 times faster. This is the best case, because walking over a range is nearly free. When every item has to be read from a file the reading dominates: in Step 8 the same comparison on a 200,000-line log gives 0.060 seconds against 0.028 seconds, a gain of about 2 times.
• Fairness. The two algorithms pass the same tests with similar scores. Algorithm L’s worst item z is 2.02 and its worst pair z is 3.06, both inside the noise you saw for Algorithm R.

Two float traps

The second half of the output shows what part 4 of the script found.

random() always returns 0.0, the smallest value it is allowed to return:
  raw random()                   raised ValueError: math domain error
  1 - random(), no clamp on w    raised ValueError: math domain error
  reservoir_l from reservoir.py  returned [19, 1, 2, 3, 4]

for k = 1000, w rounds to exactly 1.0 for the 500 smallest values random() can return: a chance of 5.6e-14 per run

Trap 1: random() can return 0.0

Wikipedia’s pseudocode for Algorithm L says that “random() generates a uniform (0,1) random number.” Python’s random() does not: the documentation says it returns a float “in the range 0.0 <= X < 1.0,” so 0.0 is allowed, and math.log(0.0) raises ValueError: math domain error. The chance is tiny, about 1 in 9 quadrillion per call, which is exactly why it survives testing and then ruins a production run. PostgreSQL’s sampler handles it explicitly. Its comment in sampling.c notes that its random function returns a value that can be 0.0, “so we must reject 0.0,” and the code loops until it gets a nonzero value. Our fix is open01, which returns 1.0 - random(), a number in (0, 1] where the logarithm always exists.

Trap 2: w can round to exactly 1.0

The skip formula divides by math.log1p(-w), the logarithm of 1 minus w. If w is exactly 1.0 that is the logarithm of zero. The first value of w is exp(log(u) / k), and when u is very close to 1 the result rounds to 1.0. The last line of the output measures it: for k = 1000, w rounds to exactly 1.0 for the 500 smallest values random() can return, a chance of 5.6e-14 per run, roughly one run in 18 trillion. No test will ever hit that by luck, but the rigged generator hits it instantly, as the second line of the output shows. The fix is a clamp that keeps w just below 1.0. The skip formula uses log1p because, as the math documentation says, its result “is calculated in a way which is accurate for x near zero,” and late in a long stream w is tiny, so computing 1 − w directly would throw digits away.

Step 5: Sample by weight with A-Res

The idea

Sometimes events should not be equally likely. You might want requests sampled in proportion to the bytes they sent, or alerts in proportion to a risk score. Efraimidis and Spirakis published a reservoir algorithm for this in 2006 (the paper is listed in the sources); Wikipedia calls it A-Res. Give every item the key u ** (1 / weight), where u is a fresh random number in (0, 1), and keep the k items with the largest keys. A heavy item raises u to a tiny power, which pushes its key toward 1, the top of the ranking. A min-heap (Python’s heapq) holds the k best keys, so each new item is compared with the smallest kept key and replaces it only if it is larger.

What does “right” mean here? One common definition of weighted sampling without replacement is to draw one item with probability proportional to its weight, remove it, and draw the next the same way from what is left. A-Res reproduces exactly that distribution. For a small case you can enumerate every possible draw and compute the exact probability of each subset, which gives us something to test against.

The code

Add this to the end of reservoir.py. The key is written as log(u) / weight, the logarithm of u ** (1 / weight). The logarithm never reorders the keys, but it avoids a rounding problem you will measure below. The tuple (key, n, item) carries the arrival number n so that ties are broken without Python ever comparing two items.

# reservoir.py (continued)
def reservoir_weighted(stream, k, rng=None):
    """A-Res (Efraimidis and Spirakis): a weighted sample without replacement from (item, weight) pairs."""
    if k < 0:
        raise ValueError("k must be 0 or more")
    rng = rng or random.Random()
    heap = []  # a min-heap of (key, arrival number, item)
    for n, (item, weight) in enumerate(stream):
        if not (weight > 0 and math.isfinite(weight)):
            raise ValueError(f"weight must be a positive, finite number, got {weight!r}")
        key = math.log(open01(rng)) / weight  # the log of u ** (1 / weight): never above 0, bigger is better
        if len(heap) < k:
            heapq.heappush(heap, (key, n, item))
        elif k and key > heap[0][0]:
            heapq.heapreplace(heap, (key, n, item))
    return [item for _, _, item in sorted(heap, reverse=True)]

And add the exact-probability checker to the end of labkit.py.

# labkit.py (continued)
def exact_subset_probabilities(weights, k):
    """Exact chance of every size-k subset when k items are drawn one at a time,
    each draw proportional to the weights of the items still left (successive sampling)."""
    probs = {}

    def walk(chosen, prob, left):
        if len(chosen) == k:
            key = frozenset(chosen)
            probs[key] = probs.get(key, 0) + prob
            return
        for i, w in enumerate(weights):
            if i not in chosen:
                walk(chosen + (i,), prob * Fraction(w, left), left - w)

    walk((), Fraction(1), sum(weights))
    return probs


def weighted_z(sampler, weights, k, runs, seed):
    """Compare a weighted sampler with the exact probabilities. Returns (worst |z|, one row per subset)."""
    exact = exact_subset_probabilities(weights, k)
    rng = random.Random(seed)
    stream = list(enumerate(weights))
    seen = {}
    for _ in range(runs):
        subset = frozenset(sampler(stream, k, rng))
        seen[subset] = seen.get(subset, 0) + 1
    rows = []
    for subset, p in sorted(exact.items(), key=lambda kv: sorted(kv[0])):
        p = float(p)
        freq = seen.get(subset, 0) / runs
        rows.append((sorted(subset), p, freq, (freq - p) / math.sqrt(p * (1 - p) / runs)))
    return max(abs(row[3]) for row in rows), rows

Save this as step05_weighted.py. It also tests the tempting formula weight * random(), compares how often the heaviest and the lightest item get picked, and counts how many distinct key values the two key formulas produce for very large weights.

"""step05_weighted.py: weighted sampling without replacement (A-Res), checked against exact probabilities."""
import math
import random

from labkit import weighted_z
from reservoir import reservoir_weighted

weights = [1, 2, 3, 4, 10]
k, runs = 2, 200_000


def product_key(stream, k, rng):
    """The tempting formula: score each item as weight * random number and keep the k highest."""
    scored = sorted(stream, key=lambda pair: rng.random() * pair[1], reverse=True)
    return [item for item, _ in scored[:k]]


print(f"weights {weights}, k = {k}, {runs:,} runs; subsets are written as item numbers (0 is the weight-1 item)")
worst, rows = weighted_z(reservoir_weighted, weights, k, runs, seed=5)
print(f"{'subset':<10}{'exact':>9}{'A-Res':>9}{'z':>7}")
for subset, exact, freq, z in rows:
    print(f"{str(subset):<10}{exact:>9.4f}{freq:>9.4f}{z:>7.2f}")
print(f"worst |z| for A-Res: {worst:.2f}")
worst_wrong, wrong_rows = weighted_z(product_key, weights, k, runs, seed=5)
print(f"worst |z| for weight * random(): {worst_wrong:.1f}")


def inclusion_of(table, item, column):
    """Add up the chance (column 1 = exact, column 2 = observed) of every subset that contains the item."""
    return sum(row[column] for row in table if item in row[0])


print(f"how often the weight-10 item (item 4) is in the sample: exact {inclusion_of(rows, 4, 1):.4f}, "
      f"A-Res {inclusion_of(rows, 4, 2):.4f}, weight * random() {inclusion_of(wrong_rows, 4, 2):.4f}")
print(f"how often the weight-1 item (item 0) is in the sample: exact {inclusion_of(rows, 0, 1):.4f}, "
      f"A-Res {inclusion_of(rows, 0, 2):.4f}, weight * random() {inclusion_of(wrong_rows, 0, 2):.4f}")

print("\ndistinct key values among 10,000 random numbers (more is better; 1 means every item ties):")
rng = random.Random(3)
draws = [1.0 - rng.random() for _ in range(10_000)]
print(f"{'weight':>8}{'u ** (1/w)':>13}{'log(u) / w':>13}")
for w in (1e6, 1e12, 1e16, 1e18):
    plain = len({u ** (1 / w) for u in draws})
    logged = len({math.log(u) / w for u in draws})
    print(f"{w:>8.0e}{plain:>13}{logged:>13}")

Run it

python step05_weighted.py

weights [1, 2, 3, 4, 10], k = 2, 200,000 runs; subsets are written as item numbers (0 is the weight-1 item)
subset        exact    A-Res      z
[0, 1]       0.0108   0.0112   1.48
[0, 2]       0.0167   0.0167  -0.22
[0, 3]       0.0230   0.0231   0.16
[0, 4]       0.0763   0.0758  -0.89
[1, 2]       0.0343   0.0339  -0.94
[1, 3]       0.0472   0.0478   1.31
[1, 4]       0.1556   0.1552  -0.41
[2, 3]       0.0728   0.0724  -0.66
[2, 4]       0.2382   0.2398   1.62
[3, 4]       0.3250   0.3241  -0.83
worst |z| for A-Res: 1.62
worst |z| for weight * random(): 147.6
how often the weight-10 item (item 4) is in the sample: exact 0.7951, A-Res 0.7949, weight * random() 0.8568
how often the weight-1 item (item 0) is in the sample: exact 0.1269, A-Res 0.1267, weight * random() 0.0175

distinct key values among 10,000 random numbers (more is better; 1 means every item ties):
  weight   u ** (1/w)   log(u) / w
   1e+06        10000        10000
   1e+12         7828        10000
   1e+16            9        10000
   1e+18            1        10000

What the output tells you

• A-Res is exact. Across all ten possible two-item subsets the observed frequencies match the exact probabilities, with a worst z score of 1.62. The two heaviest items, 3 and 4 (weights 4 and 10), appear together 32.5 percent of the time, while the two lightest, 0 and 1 (weights 1 and 2), appear together 1.1 percent of the time.
• The tempting formula is wrong. Scoring each item as weight * random() gives a worst z score of 147.6. It favours heavy items too much: the weight-10 item lands in the sample 85.68 percent of the time instead of 79.51, and the weight-1 item only 1.75 percent of the time instead of 12.69. A-Res gets 79.49 and 12.67. The product formula looks plausible and does rank heavy items higher, which is why it is a common mistake.
• The log matters for huge weights. The key u ** (1 / w) is a number just below 1.0, and floats near 1.0 are spaced about 1.1e-16 apart. At weight 1e12 only 7,828 of 10,000 random numbers still give distinct keys, at 1e16 only 9, and at 1e18 all keys round to 1.0, so the sample would ignore the weights and be decided by arrival order. The key log(u) / w keeps all 10,000 distinct at every weight.

The function also refuses weights that are zero, negative, nan or infinite, because each would silently corrupt the ranking; the tests in Step 9 check that. Wikipedia also describes A-ExpJ, a faster variant that jumps over items the way Algorithm L does. This tutorial does not implement it.

Step 6: Merge reservoirs from many workers

The idea

One process rarely sees the whole stream. Logs are sharded by host, jobs by hour, work by partition. Each worker can keep its own reservoir, and at the end you want one fair sample of everything. There is a wrong way and a right way.

The wrong way is to pool the two reservoirs and draw k items again. Each worker contributes k items however much data it saw, so a worker that saw little is over-represented.

The right way needs one extra number from each worker: how many items it saw, which the code below calls seen. Of the k slots in the final sample, some belong to worker A and the rest to worker B. To decide how many, draw k positions without replacement from an imaginary combined stream of seen A plus seen B positions. Each draw lands in A with probability (positions left in A) / (positions left in A and B). Then take that many random items from A’s reservoir and the rest from B’s. A random subset of a fair sample is itself a fair sample of the same stream, so the result is fair. This is a short derivation, not a quotation from a paper; the experiment below is what justifies trusting it.

The code

Add this to the end of reservoir.py. Reservoir applies the same rule as Algorithm R, but keeps its state in an object so it can be fed one item at a time, inspected, and merged.

# reservoir.py (continued)
class Reservoir:
    """A reservoir you can feed one item at a time, inspect, and merge with another reservoir."""

    def __init__(self, k, rng=None):
        if k < 1:
            raise ValueError("k must be 1 or more")
        self.k = k
        self.seen = 0  # how many items have gone by so far
        self.items = []
        self.rng = rng or random.Random()

    def add(self, item):
        self.seen += 1
        if len(self.items) < self.k:
            self.items.append(item)
        else:
            j = self.rng.randrange(self.seen)
            if j < self.k:
                self.items[j] = item

    def merge(self, other):
        """Return a new Reservoir that is a uniform sample of both streams taken together."""
        if other.k != self.k:
            raise ValueError("both reservoirs must have the same k")
        total = self.seen + other.seen
        take = min(self.k, total)
        from_self, left_self, left_other = 0, self.seen, other.seen
        for _ in range(take):  # draw `take` positions without replacement from the combined stream
            if self.rng.randrange(left_self + left_other) < left_self:
                from_self, left_self = from_self + 1, left_self - 1
            else:
                left_other -= 1
        merged = Reservoir(self.k, self.rng)
        merged.seen = total
        merged.items = self.rng.sample(self.items, from_self) + self.rng.sample(other.items, take - from_self)
        return merged

Save this as step06_merge.py. It splits a stream of 1,000 items between two workers, one that sees 100 items and one that sees 900, and merges both ways. Then it tests four workers with unequal chunks, merged one after another.

"""step06_merge.py: merging reservoirs from several workers, the wrong way and the right way."""
from functools import reduce

from labkit import inclusion_z
from reservoir import Reservoir


def right_merge(a, b, rng):
    return a.merge(b).items


def naive_merge(a, b, rng):
    return rng.sample(a.items + b.items, a.k)  # pool both samples, then draw again


def two_workers(merge, split):
    """Items below `split` go to worker A, the rest to worker B; then the two reservoirs are merged."""

    def sampler(stream, k, rng):
        a, b = Reservoir(k, rng), Reservoir(k, rng)
        for item in stream:
            (a if item < split else b).add(item)
        return merge(a, b, rng)

    return sampler


k, split, n, runs = 10, 100, 1000, 20_000
print(f"worker A saw {split} items, worker B saw {n - split}; k = {k}; every item should appear with probability {k / n:.4f}")
for name, merge in (("pool and draw again", naive_merge), ("Reservoir.merge", right_merge)):
    result = inclusion_z(two_workers(merge, split), n, k, runs, seed=7)
    a_freq = sum(result.freqs[:split]) / split
    b_freq = sum(result.freqs[split:]) / (n - split)
    print(f"{name:<22} item in A: {a_freq:.4f}   item in B: {b_freq:.4f}   worst |z| {result.worst_z:6.1f}")


def chunked(sizes):
    """Four workers with unequal chunks of the stream, merged one after another."""
    bounds = [sum(sizes[:i]) for i in range(len(sizes) + 1)]

    def sampler(stream, k, rng):
        parts = [Reservoir(k, rng) for _ in sizes]
        for item in stream:
            parts[max(i for i, start in enumerate(bounds[:-1]) if item >= start)].add(item)
        return reduce(lambda x, y: x.merge(y), parts).items

    return sampler


sizes, k, runs = (20, 30, 50, 100), 5, 50_000
n = sum(sizes)
items = inclusion_z(chunked(sizes), n, k, runs, seed=8)
print(f"\nfour workers with chunks {sizes}, k = {k}, {runs:,} runs; fair item frequency {k / n:.4f}")
print(f"item frequencies {items.low:.4f} to {items.high:.4f}, worst item |z| {items.worst_z:.2f}")

a, b = Reservoir(10), Reservoir(10)
for x in range(3):
    a.add(x)
for x in range(100, 102):
    b.add(x)
merged = a.merge(b)
print(f"\nmerging reservoirs that saw 3 and 2 items (k = 10): seen = {merged.seen}, items = {sorted(merged.items)}")

Run it

python step06_merge.py

worker A saw 100 items, worker B saw 900; k = 10; every item should appear with probability 0.0100
pool and draw again    item in A: 0.0500   item in B: 0.0056   worst |z|   62.6
Reservoir.merge        item in A: 0.0100   item in B: 0.0100   worst |z|    3.4

four workers with chunks (20, 30, 50, 100), k = 5, 50,000 runs; fair item frequency 0.0250
item frequencies 0.0233 to 0.0266, worst item |z| 2.41

merging reservoirs that saw 3 and 2 items (k = 10): seen = 5, items = [0, 1, 2, 100, 101]

What the output tells you

• Pooling is badly biased. Every item should appear with probability 0.0100. With “pool and draw again”, an item from the small worker appears with probability 0.0500 and an item from the big worker with 0.0056. The worker that saw 10 percent of the stream supplies half of the final sample (z score 62.6).
• Reservoir.merge is fair. Both groups come out at 0.0100 and the worst z score is 3.4, ordinary for the biggest of 1,000 comparisons.
• Chains of merges work. With four workers that saw 20, 30, 50 and 100 items, the item frequencies run from 0.0233 to 0.0266 around the fair 0.0250, worst z 2.41. This test looks at single items only, which is the check I could afford at this size.
• Small streams are kept whole. Merging reservoirs that saw 3 and 2 items with k = 10 gives a merged seen of 5 and all five items.

One precaution: give each worker its own generator, with a different seed or none. The argument above assumes the workers’ random choices are independent of each other.

Step 7: Keep rare events with one reservoir per category, and check the memory

The idea

A uniform sample is a fine tool for estimating shares and a poor tool for finding needles. If an event happens 5 times in 2 million, a fair sample of 1,000 will almost never contain it, no matter how good the algorithm is. When you need every category represented, keep one reservoir per category. The class below does that. A fair sample is still fair within each category; across categories it deliberately is not.

The code

Add this to the end of reservoir.py.

# reservoir.py (continued)
class Stratified:
    """One independent Reservoir per category, so a rare category cannot be crowded out."""

    def __init__(self, k, rng=None):
        self.k = k
        self.rng = rng or random.Random()
        self.groups = {}

    def add(self, category, item):
        group = self.groups.get(category)
        if group is None:
            group = self.groups[category] = Reservoir(self.k, self.rng)
        group.add(item)

Save this as step07_stratified.py. It streams 2 million events, five of them rare key_export events, through a uniform reservoir of 1,000 and a per-kind reservoir of 100 at the same time. Then it computes how likely a uniform sample is to catch a rare event, and measures the memory the sampler needs.

"""step07_stratified.py: a uniform sample cannot promise to hold a rare event; one reservoir per kind can."""
import math
import random
import tracemalloc
from collections import Counter

from events import auth_events
from reservoir import Reservoir, Stratified, reservoir_l

N = 2_000_000
RARE = (123_457, 456_789, 789_123, 1_234_567, 1_876_543)

rng = random.Random(21)
uniform, by_kind = Reservoir(1000, rng), Stratified(100, rng)
for pos, kind in auth_events(N, rare_positions=RARE):
    uniform.add((pos, kind))
    by_kind.add(kind, (pos, kind))

in_uniform = Counter(kind for _, kind in uniform.items)
print(f"{N:,} events, one pass; the uniform reservoir holds {len(uniform.items):,}, each per-kind reservoir holds up to {by_kind.k}")
print(f"{'kind':<16}{'seen':>10}{'uniform':>9}{'per kind':>10}")
for kind, group in sorted(by_kind.groups.items(), key=lambda kv: -kv[1].seen):
    print(f"{kind:<16}{group.seen:>10,}{in_uniform[kind]:>9}{len(group.items):>10}")

exact = 1 - math.comb(N - len(RARE), 1000) / math.comb(N, 1000)
runs = 20_000
hits = sum(bool(set(random.Random(s).sample(range(N), 1000)) & set(RARE)) for s in range(runs))
print(f"\nchance that a uniform sample of 1,000 holds at least one of the {len(RARE)} rare events: {exact:.2%}")
print(f"measured with {runs:,} uniform samples: {hits} of {runs:,} held one")

print("\npeak memory while building a sample of 1,000 events:")


def peak_mib(build):
    tracemalloc.start()
    build()
    peak = tracemalloc.get_traced_memory()[1]
    tracemalloc.stop()
    return peak / 2**20


for n in (50_000, 200_000, 800_000):
    keep_all = peak_mib(lambda: random.sample(list(auth_events(n)), 1000))
    streaming = peak_mib(lambda: reservoir_l(auth_events(n), 1000, random.Random(1)))
    print(f"{n:>9,} events: keep everything then random.sample {keep_all:7.2f} MiB   reservoir_l {streaming:5.2f} MiB")

Run it

python step07_stratified.py

2,000,000 events, one pass; the uniform reservoir holds 1,000, each per-kind reservoir holds up to 100
kind                  seen  uniform  per kind
login_ok         1,555,042      793       100
login_fail         292,250      130       100
mfa_prompt          87,374       43       100
stuffing            35,665       18       100
healthcheck         20,000        9       100
password_reset       9,664        7       100
key_export               5        0         5

chance that a uniform sample of 1,000 holds at least one of the 5 rare events: 0.25%
measured with 20,000 uniform samples: 43 of 20,000 held one

peak memory while building a sample of 1,000 events:
   50,000 events: keep everything then random.sample    4.68 MiB   reservoir_l  0.04 MiB
  200,000 events: keep everything then random.sample   18.28 MiB   reservoir_l  0.05 MiB
  800,000 events: keep everything then random.sample   73.46 MiB   reservoir_l  0.05 MiB

What the output tells you

• The uniform sample has no rare events. It holds 793 login_ok events against about 778 expected from their 77.8 percent share, but 0 of the 5 key_export events. The exact chance that a uniform sample of 1,000 out of 2,000,000 holds at least one of five particular events is 1 − C(N − 5, 1000) / C(N, 1000), which is 0.25 percent. Drawing 20,000 uniform samples with random.sample(range(N), 1000), which is quick on a range and just as fair, 43 held one. The per-kind reservoirs hold all 5.
• Do not read shares off a stratified sample. The per-kind reservoirs hold 605 events in total, and 5 of them are key_export, a share of 0.83 percent, when the true share is 0.00025 percent. To estimate overall shares, weight each kept item by (events seen) / (events kept) for its kind. By that rule each kept login_ok event stands for 1,555,042 / 100 = 15,550 events and each key_export event for exactly 1.
• Memory is flat. Keeping everything and then calling random.sample needs 4.68, 18.28 and 73.46 MiB for 50,000, 200,000 and 800,000 events, growing with the stream. reservoir_l needs 0.04 to 0.05 MiB at every size, however long the stream is. (I measured with Python’s tracemalloc; our tracemalloc walkthrough explains how.)
• Bound the number of categories. Per-category reservoirs use memory for categories × k items. Do not use an unbounded field such as an IP address or user ID as the category without a cap.

Step 8: Wrap it in a command-line tool

The idea

Now turn reservoir_l into something you can point at any big file or pipe. The tool reads lines from one or more files or from standard input, keeps k of them, and prints them, optionally with line numbers and in the original order. To try it you also need a log file, so a small script writes 200,000 synthetic lines to auth.log using the same generator as before.

The code

Save this as make_log.py.

"""make_log.py: write a synthetic authentication log (200,000 lines by default) to auth.log."""
import random
import sys
from datetime import datetime, timedelta, timezone

from events import auth_events

n = int(sys.argv[1]) if len(sys.argv) > 1 else 200_000
rng = random.Random(99)
start = datetime(2026, 11, 3, 0, 0, 0, tzinfo=timezone.utc)
with open("auth.log", "w", encoding="utf-8", newline="\n") as out:
    for pos, kind in auth_events(n, seed=7):
        stamp = (start + timedelta(seconds=pos * 0.4)).strftime("%Y-%m-%dT%H:%M:%SZ")
        out.write(f"{stamp} {kind} user=u{rng.randrange(5000):04d} ip=203.0.113.{rng.randrange(1, 255)}\n")
print(f"wrote auth.log with {n:,} lines")

Save this as sample_lines.py.

"""sample_lines.py: print k uniformly random lines from files or standard input, using memory for k lines only.

usage: python sample_lines.py -k 5 [--seed 7] [--keep-order] [--numbers] FILE [FILE ...]    (- means standard input)
"""
import argparse
import itertools
import random
import sys

from reservoir import reservoir_l


def numbered_lines(paths):
    """Yield (line number, text) for every line of every input, numbering on across files."""
    counter = itertools.count(1)
    for path in paths:
        handle = sys.stdin if path == "-" else open(path, encoding="utf-8", errors="replace", newline="")
        with handle:
            for text in handle:
                yield next(counter), text.rstrip("\r\n")


def main(argv=None):
    parser = argparse.ArgumentParser(description="Print k uniformly random lines from the input, in one pass.")
    parser.add_argument("paths", nargs="+", help="files to read; use - for standard input")
    parser.add_argument("-k", type=int, required=True, help="how many lines to keep")
    parser.add_argument("--seed", type=int, help="seed for a repeatable sample")
    parser.add_argument("--keep-order", action="store_true", help="print the lines in the order they appear in the input")
    parser.add_argument("--numbers", action="store_true", help="put the line number and a tab in front of each line")
    args = parser.parse_args(argv)
    if args.k < 0:
        parser.error("-k must be 0 or more")
    sys.stdin.reconfigure(encoding="utf-8", errors="replace")
    sys.stdout.reconfigure(encoding="utf-8", newline="\n")
    sample = reservoir_l(numbered_lines(args.paths), args.k, random.Random(args.seed))
    if args.keep_order:
        sample.sort()
    for number, text in sample:
        print(f"{number}\t{text}" if args.numbers else text)
    return 0


if __name__ == "__main__":
    sys.exit(main())

Each sampled item is a pair of line number and text. Keeping the number costs nothing and lets --keep-order sort the sample back into stream order, which a reservoir does not preserve on its own. - means standard input, and the reading code numbers lines on across several files.

Run it

python make_log.py

wrote auth.log with 200,000 lines

python sample_lines.py -k 5 --seed 7 --numbers auth.log

114174	2026-11-03T12:41:09Z login_ok user=u2715 ip=203.0.113.4
155420	2026-11-03T17:16:07Z login_ok user=u4133 ip=203.0.113.15
13939	2026-11-03T01:32:55Z mfa_prompt user=u2085 ip=203.0.113.188
33434	2026-11-03T03:42:53Z login_ok user=u3661 ip=203.0.113.89
51704	2026-11-03T05:44:41Z login_ok user=u1592 ip=203.0.113.106

The numbers on the left are line numbers in auth.log. They are spread across the whole log, as a fair sample should be. Now another seed with the lines in their original order:

python sample_lines.py -k 5 --seed 8 --numbers --keep-order auth.log

9417	2026-11-03T01:02:46Z login_ok user=u3320 ip=203.0.113.154
96776	2026-11-03T10:45:10Z login_ok user=u4276 ip=203.0.113.33
109692	2026-11-03T12:11:16Z mfa_prompt user=u1371 ip=203.0.113.113
137099	2026-11-03T15:13:59Z login_ok user=u1751 ip=203.0.113.171
195994	2026-11-03T21:46:37Z login_ok user=u0254 ip=203.0.113.222

And standard input, using cat on Linux and macOS or Git Bash (use type auth.log | in cmd):

cat auth.log | python sample_lines.py -k 3 --seed 1 -

2026-11-03T11:10:46Z login_ok user=u4572 ip=203.0.113.74
2026-11-03T00:42:30Z login_ok user=u0323 ip=203.0.113.211
2026-11-03T11:06:32Z login_ok user=u2621 ip=203.0.113.156

Finally, a check script that runs the tool the way a user would and compares a 1,000-line sample with the whole file. It also times Algorithms R and L over the file. Save it as step08_cli_check.py.

"""step08_cli_check.py: run sample_lines.py on auth.log and compare the sample with the whole file."""
import random
import subprocess
import sys
import time
from collections import Counter

from reservoir import reservoir_l, reservoir_r
from sample_lines import numbered_lines

K = 1000


def run(*args):
    done = subprocess.run([sys.executable, "sample_lines.py", *args, "auth.log"], capture_output=True, text=True, encoding="utf-8", check=True)
    return done.stdout.splitlines()


def kinds(lines):
    return Counter(line.split(" ")[1] for line in lines)


whole = open("auth.log", encoding="utf-8").read().splitlines()
sample = run("-k", str(K), "--seed", "7")
all_kinds, got = kinds(whole), kinds(sample)
print(f"{len(whole):,} lines in auth.log; the sample holds {len(sample):,} lines")
print(f"{'kind':<16}{'whole file':>11}{'sample':>9}{'expected':>10}")
for kind, count in all_kinds.most_common():
    print(f"{kind:<16}{count / len(whole):>11.2%}{got[kind] / len(sample):>9.2%}{count / len(whole) * K:>10.1f}")

again = run("-k", str(K), "--seed", "7")
other = run("-k", str(K), "--seed", "8")
ordered = run("-k", str(K), "--seed", "7", "--keep-order", "--numbers")
numbers = [int(line.split("\t")[0]) for line in ordered]
print(f"\nsame seed gives the same sample: {sample == again}")
print(f"another seed gives a different sample: {sample != other}")
print(f"--keep-order lists the line numbers in increasing order: {numbers == sorted(numbers)}")
print(f"every sampled line is a line of the file: {set(sample) <= set(whole)}")


def seconds(sampler):
    start = time.perf_counter()
    sampler(numbered_lines(["auth.log"]), K, random.Random(7))
    return time.perf_counter() - start


print()
for name, sampler in (("Algorithm R", reservoir_r), ("Algorithm L", reservoir_l)):
    print(f"{name} over the {len(whole):,} lines of auth.log: {min(seconds(sampler) for _ in range(3)):.3f} s (best of 3)")

python step08_cli_check.py

200,000 lines in auth.log; the sample holds 1,000 lines
kind             whole file   sample  expected
login_ok             77.72%   76.50%     777.2
login_fail           14.62%   15.40%     146.2
mfa_prompt            4.41%    4.20%      44.1
stuffing              1.78%    2.20%      17.8
healthcheck           1.00%    1.20%      10.0
password_reset        0.47%    0.50%       4.7

same seed gives the same sample: True
another seed gives a different sample: True
--keep-order lists the line numbers in increasing order: True
every sampled line is a line of the file: True

Algorithm R over the 200,000 lines of auth.log: 0.060 s (best of 3)
Algorithm L over the 200,000 lines of auth.log: 0.028 s (best of 3)

What the output tells you

The sample’s shares are close to the file’s: 76.50 percent login_ok against 77.72 percent, 2.20 percent stuffing against 1.78 percent. With only 1,000 lines a few tenths of a percent of wobble is expected, and it shrinks as k grows. The same seed gives the same sample, a different seed gives a different one, --keep-order lists the line numbers in increasing order, and every sampled line really is a line of the file. The last two lines show Algorithm L about twice as fast as R when the items come from a file, as promised earlier.

Two cautions about seeds. First, a seed makes a sample repeatable, but the Python documentation promises less than you might hope: “Most of the random module’s algorithms and seeding functions are subject to change across Python versions, but two aspects are guaranteed not to change,” and the second is that “The generator’s random() method will continue to produce the same sequence when the compatible seeder is given the same seed.” Our tool also uses randrange, so if you must replay a sample exactly, record the seed and the Python version. Second, a seed makes the sample predictable. The same documentation warns that “The pseudo-random generators of this module should not be used for security purposes.” If an audit sample must not be guessable by someone who might want to hide from it, pass random.SystemRandom() as the generator. Its documentation says it “Does not rely on software state, and sequences are not reproducible.” I checked that reservoir_l, reservoir_r, reservoir_weighted and Reservoir.merge all accept it.

Step 9: Turn the checks into tests, and prove the tests can fail

The tests

Everything so far was a script you read once. Tests keep the guarantees true after you edit the code. Save this as test_reservoir.py. It checks the boring things (a short stream is kept whole, k = 0 returns nothing, a negative k is refused, there are no duplicates) and the statistical things (fair items, fair pairs, exact weighted probabilities, fair merges, the rare-kind guarantee). It also feeds Algorithm L the rigged generator from Step 4 and bounds the number of random numbers it draws.

"""test_reservoir.py: tests for reservoir.py and sample_lines.py. Run with: python -m pytest -q"""
import math
import pathlib
import random
import subprocess
import sys
from functools import reduce

import pytest

from labkit import CountingRandom, ZeroRandom, inclusion_z, pair_z, weighted_z
from reservoir import Reservoir, Stratified, reservoir_l, reservoir_r, reservoir_weighted

SAMPLERS = [reservoir_r, reservoir_l]
LIMIT = 5.0  # a result more than 5 standard errors from the truth counts as a failure
SCRIPT = pathlib.Path(__file__).with_name("sample_lines.py")


@pytest.mark.parametrize("sampler", SAMPLERS)
def test_a_short_stream_is_kept_whole(sampler):
    assert sorted(sampler(iter(range(3)), 5, random.Random(1))) == [0, 1, 2]
    assert sampler(iter([]), 5, random.Random(1)) == []


@pytest.mark.parametrize("sampler", SAMPLERS)
def test_sample_has_k_distinct_items_from_the_stream(sampler):
    sample = sampler(iter(range(1000)), 25, random.Random(2))
    assert len(sample) == 25 and len(set(sample)) == 25 and set(sample) <= set(range(1000))


@pytest.mark.parametrize("sampler", SAMPLERS)
def test_k_zero_gives_nothing_and_negative_k_is_refused(sampler):
    assert sampler(iter(range(10)), 0, random.Random(1)) == []
    with pytest.raises(ValueError):
        sampler(iter(range(10)), -1, random.Random(1))


@pytest.mark.parametrize("sampler", SAMPLERS)
def test_every_item_is_equally_likely(sampler):
    assert inclusion_z(sampler, 12, 4, 30_000, seed=3).worst_z < LIMIT


@pytest.mark.parametrize("sampler", SAMPLERS)
def test_every_pair_is_equally_likely(sampler):
    assert pair_z(sampler, 10, 4, 30_000, seed=4) < LIMIT


def test_algorithm_l_asks_for_few_random_numbers():
    n, k = 100_000, 100
    rng = CountingRandom(5)
    reservoir_l(iter(range(n)), k, rng)
    assert rng.draws < 4 * k * (1 + math.log(n / k))


def test_algorithm_l_survives_a_generator_that_returns_zero():
    assert len(reservoir_l(iter(range(20)), 5, ZeroRandom())) == 5


def test_weighted_sample_matches_the_exact_probabilities():
    worst, _ = weighted_z(reservoir_weighted, [1, 2, 3, 4, 10], 2, 60_000, seed=5)
    assert worst < LIMIT


@pytest.mark.parametrize("bad", [0, -1, float("nan"), float("inf")])
def test_weighted_refuses_unusable_weights(bad):
    with pytest.raises(ValueError):
        reservoir_weighted([("a", 1), ("b", bad)], 1, random.Random(1))


def test_weighted_still_ranks_enormous_weights():
    rng = random.Random(6)
    stream = [("light", 10**18), ("heavy", 4 * 10**18)]
    runs = 20_000
    heavy = sum(reservoir_weighted(stream, 1, rng) == ["heavy"] for _ in range(runs))
    se = math.sqrt(0.8 * 0.2 / runs)
    assert abs(heavy / runs - 0.8) / se < LIMIT


def chunked(sizes):
    bounds = [sum(sizes[:i]) for i in range(len(sizes) + 1)]

    def sampler(stream, k, rng):
        parts = [Reservoir(k, rng) for _ in sizes]
        for item in stream:
            parts[max(i for i, start in enumerate(bounds[:-1]) if item >= start)].add(item)
        return reduce(lambda x, y: x.merge(y), parts).items

    return sampler


def test_merged_reservoirs_match_one_pass_over_the_whole_stream():
    assert inclusion_z(chunked((10, 30, 80)), 120, 6, 20_000, seed=7).worst_z < LIMIT


def test_merging_small_reservoirs_keeps_everything():
    a, b = Reservoir(10), Reservoir(10)
    for x in range(3):
        a.add(x)
    for x in range(100, 102):
        b.add(x)
    merged = a.merge(b)
    assert merged.seen == 5 and sorted(merged.items) == [0, 1, 2, 100, 101]


def test_merge_refuses_different_k():
    with pytest.raises(ValueError):
        Reservoir(3).merge(Reservoir(4))


def test_stratified_keeps_every_event_of_a_rare_kind():
    strat = Stratified(5, random.Random(1))
    for i in range(10_000):
        strat.add("common", i)
    for i in range(3):
        strat.add("rare", i)
    assert len(strat.groups["rare"].items) == 3 and len(strat.groups["common"].items) == 5


def run_cli(*args):
    return subprocess.run([sys.executable, str(SCRIPT), *args], capture_output=True, text=True, encoding="utf-8", cwd=SCRIPT.parent)


def test_command_line_tool_is_repeatable_and_can_keep_order(tmp_path):
    log = tmp_path / "log.txt"
    log.write_text("".join(f"line {n}\n" for n in range(1, 1001)), encoding="utf-8")
    first = run_cli("-k", "5", "--seed", "3", str(log))
    assert first.returncode == 0 and len(first.stdout.splitlines()) == 5
    assert run_cli("-k", "5", "--seed", "3", str(log)).stdout == first.stdout
    ordered = run_cli("-k", "20", "--seed", "3", "--numbers", "--keep-order", str(log)).stdout.splitlines()
    numbers = [int(line.split("\t")[0]) for line in ordered]
    assert numbers == sorted(numbers) and len(numbers) == 20


def test_command_line_tool_refuses_a_negative_k(tmp_path):
    log = tmp_path / "log.txt"
    log.write_text("a\n", encoding="utf-8")
    assert run_cli("-k", "-1", str(log)).returncode == 2

python -m pytest -q

........................                                                                                                                                           [100%]
24 passed in 1.70s

How strict should a statistical test be?

The tests use a limit of 5 on the z score. That choice is not arbitrary. A correct sampler still produces a z score above 2 fairly often when you compare many numbers, because you are looking at the biggest of many noisy values. This script runs the same correct sampler under 200 different seeds and counts false alarms at four limits. Save it as step09_seed_sweep.py.

"""step09_seed_sweep.py: how strict should a statistical test be? Run the same correct sampler with 200 different seeds."""
from labkit import inclusion_z
from reservoir import reservoir_r

seeds = range(1, 201)
worst = [inclusion_z(reservoir_r, 12, 4, 20_000, seed).worst_z for seed in seeds]
print(f"reservoir_r, n = 12, k = 4, 20,000 runs per seed, {len(worst)} seeds; worst |z| per seed ranged from {min(worst):.2f} to {max(worst):.2f}")
for limit in (2, 3, 4, 5):
    print(f"  false alarms with a limit of {limit}: {sum(z > limit for z in worst)} of {len(worst)} seeds")

python step09_seed_sweep.py

reservoir_r, n = 12, k = 4, 20,000 runs per seed, 200 seeds; worst |z| per seed ranged from 0.89 to 3.37
  false alarms with a limit of 2: 84 of 200 seeds
  false alarms with a limit of 3: 4 of 200 seeds
  false alarms with a limit of 4: 0 of 200 seeds
  false alarms with a limit of 5: 0 of 200 seeds

A limit of 2 raised a false alarm for 84 of the 200 seeds, a limit of 3 for 4 seeds, and limits of 4 and 5 for none. The more numbers you compare, the larger the biggest z naturally gets, so the limit has to grow with the number of comparisons. A limit of 5 is safe for the tests here. A bug large enough to matter produced z scores from 27 to 182 in Step 3, far beyond it. (A limit this strict would miss very small biases; if you need to detect those, raise the number of runs rather than lowering the limit.)

Can the tests actually fail?

A test suite that has never failed proves little. The script below copies the lab into a temporary folder, plants one bug at a time (a changed operator, a missing clamp, a wrong sort) and checks that the tests notice. Save it as mutations.py.

"""mutations.py: plant one bug at a time in a copy of the lab and check that the tests notice."""
import os
import pathlib
import shutil
import subprocess
import sys
import tempfile

HERE = pathlib.Path(__file__).parent
MUTANTS = [
    ("R accepts with probability k / (i + 1)", "reservoir.py", "j = rng.randrange(i)", "j = rng.randrange(i + 1)"),
    ("R accepts when j equals k", "reservoir.py", "if j < k:", "if j <= k:"),
    ("R fills only k - 1 slots first", "reservoir.py", "if i <= k:", "if i < k:"),
    ("L without the clamp on w", "reservoir.py", "w = min(math.exp(math.log(open01(rng)) / k), 1.0 - 2.0**-53)", "w = math.exp(math.log(open01(rng)) / k)"),
    ("L uses raw random(), which can be 0.0", "reservoir.py", "return 1.0 - rng.random()", "return rng.random()"),
    ("L skips one item too many", "reservoir.py", "itertools.islice(it, skip, skip + 1)", "itertools.islice(it, skip + 1, skip + 2)"),
    ("L always replaces slot 0", "reservoir.py", "reservoir[rng.randrange(k)] = item", "reservoir[0] = item"),
    ("A-Res multiplies by the weight", "reservoir.py", "key = math.log(open01(rng)) / weight", "key = math.log(open01(rng)) * weight"),
    ("A-Res keeps the lowest keys", "reservoir.py", "elif k and key > heap[0][0]:", "elif k and key < heap[0][0]:"),
    ("Reservoir.add draws from k, not seen", "reservoir.py", "j = self.rng.randrange(self.seen)", "j = self.rng.randrange(self.k)"),
    ("merge flips a fair coin", "reservoir.py", "if self.rng.randrange(left_self + left_other) < left_self:", "if self.rng.random() < 0.5:"),
    ("Stratified forgets its groups", "reservoir.py", "group = self.groups.get(category)", "group = self.groups.get(None)"),
    ("--keep-order reverses instead of sorting", "sample_lines.py", "sample.sort()", "sample.reverse()"),
]

caught = 0
for name, filename, old, new in MUTANTS:
    with tempfile.TemporaryDirectory() as tmp:
        work = pathlib.Path(tmp)
        for source in HERE.glob("*.py"):
            shutil.copy(source, work / source.name)
        text = (work / filename).read_text(encoding="utf-8")
        assert text.count(old) == 1, f"{name}: expected one match, found {text.count(old)}"
        (work / filename).write_text(text.replace(old, new), encoding="utf-8", newline="\n")
        env = dict(os.environ, PYTHONDONTWRITEBYTECODE="1", COLUMNS="170")
        done = subprocess.run([sys.executable, "-m", "pytest", "-q", "-x", "-p", "no:cacheprovider"], cwd=work, env=env,
                              capture_output=True, text=True, encoding="utf-8", errors="replace")
        failed = [line for line in done.stdout.splitlines() if line.startswith("FAILED")]
        status = "caught" if done.returncode != 0 else "SURVIVED"
        caught += done.returncode != 0
        first = failed[0].split(" - ")[0].replace("FAILED ", "") if failed else ""
        print(f"{status:<9}{name:<44}{first}")
print(f"{caught} of {len(MUTANTS)} planted bugs were caught")

python mutations.py

caught   R accepts with probability k / (i + 1)      test_reservoir.py::test_every_item_is_equally_likely[reservoir_r]
caught   R accepts when j equals k                   test_reservoir.py::test_sample_has_k_distinct_items_from_the_stream[reservoir_r]
caught   R fills only k - 1 slots first              test_reservoir.py::test_sample_has_k_distinct_items_from_the_stream[reservoir_r]
caught   L without the clamp on w                    test_reservoir.py::test_algorithm_l_survives_a_generator_that_returns_zero
caught   L uses raw random(), which can be 0.0       test_reservoir.py::test_algorithm_l_survives_a_generator_that_returns_zero
caught   L skips one item too many                   test_reservoir.py::test_every_item_is_equally_likely[reservoir_l]
caught   L always replaces slot 0                    test_reservoir.py::test_every_item_is_equally_likely[reservoir_l]
caught   A-Res multiplies by the weight              test_reservoir.py::test_weighted_sample_matches_the_exact_probabilities
caught   A-Res keeps the lowest keys                 test_reservoir.py::test_weighted_sample_matches_the_exact_probabilities
caught   Reservoir.add draws from k, not seen        test_reservoir.py::test_merged_reservoirs_match_one_pass_over_the_whole_stream
caught   merge flips a fair coin                     test_reservoir.py::test_merged_reservoirs_match_one_pass_over_the_whole_stream
caught   Stratified forgets its groups               test_reservoir.py::test_stratified_keeps_every_event_of_a_rare_kind
caught   --keep-order reverses instead of sorting    test_reservoir.py::test_command_line_tool_is_repeatable_and_can_keep_order
13 of 13 planted bugs were caught

All 13 planted bugs were caught. Note the first two rows of the Algorithm L group: both float traps from Step 4 are caught by the rigged-generator test, since no honest random run would ever reach them.

Step 10: Run everything with one command

Finally, one script that runs every step, the tests and the mutation check, and saves each output into an out folder. Save it as run_all.py.

"""run_all.py: run every step, the tests and the mutation check; save the output of each into out/NAME.txt."""
import os
import pathlib
import subprocess
import sys
import time

HERE = pathlib.Path(__file__).parent
(HERE / "out").mkdir(exist_ok=True)
env = dict(os.environ, PYTHONUTF8="1", PYTHONIOENCODING="utf-8", COLUMNS="170")

JOBS = [
    ("step01_naive", ["step01_naive.py"]),
    ("step02_algorithm_r", ["step02_algorithm_r.py"]),
    ("step03_bugs", ["step03_bugs.py"]),
    ("step04_algorithm_l", ["step04_algorithm_l.py"]),
    ("step05_weighted", ["step05_weighted.py"]),
    ("step06_merge", ["step06_merge.py"]),
    ("step07_stratified", ["step07_stratified.py"]),
    ("make_log", ["make_log.py"]),
    ("cli_sample", ["sample_lines.py", "-k", "5", "--seed", "7", "--numbers", "auth.log"]),
    ("cli_ordered", ["sample_lines.py", "-k", "5", "--seed", "8", "--numbers", "--keep-order", "auth.log"]),
    ("cli_stdin", ["sample_lines.py", "-k", "3", "--seed", "1", "-"]),
    ("step08_cli_check", ["step08_cli_check.py"]),
    ("step09_seed_sweep", ["step09_seed_sweep.py"]),
    ("pytest", ["-m", "pytest", "-q", "-p", "no:cacheprovider"]),
    ("mutations", ["mutations.py"]),
]
only = sys.argv[1:]
for stem, args in JOBS:
    if only and stem not in only:
        continue
    start = time.perf_counter()
    feed = (HERE / "auth.log").read_text(encoding="utf-8") if stem == "cli_stdin" else None  # standard input for the "-" case
    done = subprocess.run([sys.executable, *args], cwd=HERE, env=env, input=feed, capture_output=True, text=True, encoding="utf-8", errors="replace")
    text = done.stdout + (("\n" + done.stderr) if done.stderr.strip() else "")
    (HERE / "out" / f"{stem}.txt").write_text(text, encoding="utf-8", newline="\n")
    print(f"{stem:<22} exit {done.returncode}  {time.perf_counter() - start:6.1f} s  {len(text.splitlines())} lines")

python run_all.py

step01_naive           exit 0     0.7 s  15 lines
step02_algorithm_r     exit 0     0.4 s  11 lines
step03_bugs            exit 0     2.5 s  7 lines
step04_algorithm_l     exit 0     2.4 s  16 lines
step05_weighted        exit 0     0.4 s  21 lines
step06_merge           exit 0    15.0 s  8 lines
step07_stratified      exit 0     6.5 s  17 lines
make_log               exit 0     0.6 s  1 lines
cli_sample             exit 0     0.1 s  5 lines
cli_ordered            exit 0     0.1 s  5 lines
cli_stdin              exit 0     0.1 s  3 lines
step08_cli_check       exit 0     0.8 s  16 lines
step09_seed_sweep      exit 0     6.8 s  5 lines
pytest                 exit 0     2.0 s  2 lines
mutations              exit 0    11.7 s  14 lines

Every job should print exit 0. The slowest are Step 6 and the mutation check, at about 15 and 12 seconds on my machine.

Common mistakes and how to spot them

• Counting from 0 when the textbook counts from 1. randrange(i + 1) or randint(0, i) in place of randrange(i) made the first five items 20 percent more likely than the later ones (Step 3). Only a statistical test notices.
• Calling random.sample on a stream. It needs a sequence, so you hold the whole stream in memory first: 73 MiB for 800,000 small events against 0.05 MiB (Step 7). If you do know n in advance, random.sample(range(n), k) is the right tool; the Python documentation calls this “especially fast and space efficient for sampling from a large population.” PostgreSQL does exactly that kind of thing for blocks: its comment says “Since we know the total number of blocks in advance, we can use the straightforward Algorithm S from Knuth 3.4.2, rather than Vitter’s algorithm.”
• Taking a logarithm that can hit zero. Use 1.0 - random() so the argument is never 0.0, and keep w below 1 so that the logarithm of 1 minus w exists (Step 4).
• Merging by pooling. Always carry the count of items seen and merge with it (Step 6).
• Reading overall shares off a stratified sample. Weight each category by seen over kept (Step 7).
• Testing with one run, or with a limit that false-alarms. Run thousands of times and choose a limit using a seed sweep (Step 9).
• Treating a seed as a promise. It repeats a run on the same Python version; it is not secret and not guaranteed across versions (Step 8).
• Expecting a uniform sample to find rare events. It estimates; it does not search. Use per-category reservoirs or a different method for needles (Step 7).

Confirm it all works end to end

From a clean folder with the files above, these four commands should all succeed: python run_all.py prints exit 0 for all 15 jobs; python -m pytest -q reports 24 passed; python mutations.py ends with 13 of 13 planted bugs were caught; and python sample_lines.py -k 5 --seed 7 --numbers auth.log reproduces the five lines of Step 8 on the same Python version. If your timings differ, compare the ratios. If a statistical line differs slightly on another Python version, that is the seed caveat from Step 8, not a bug.

Where to go next

• Vitter’s Algorithm Z. This is Vitter’s own optimal algorithm; Wikipedia notes that Li’s Algorithm L achieves the same asymptotic runtime with a simpler method. PostgreSQL’s sampling.c says its routines “embody Algorithm Z from ‘Random sampling with a reservoir’ by Jeffrey S. Vitter,” using the simpler Algorithm X until the stream is long enough. Reading that file is a good exercise once Algorithm L feels familiar.
• A-ExpJ. Add jumps to the weighted sampler, following the pseudocode in the Wikipedia article.
• How big should k be? PostgreSQL’s analyze.c derives its sample size from a published bound and sets the minimum rows to 300 times the statistics target; the documentation gives the default target as 100, which works out to 30,000 rows (my arithmetic). Work out what error you can accept, then size k for it.
• More sketches. Reservoir sampling belongs to a family of “keep a small summary of a stream” tools. We cover others in our tutorials on the Bloom filter, cuckoo filter, HyperLogLog, Count-Min sketch, MinHash and DDSketch. The last one uses the same “merge the summaries from each worker” idea as Step 6.
• Test more cleverly. Our Hypothesis tutorial shows how to generate unusual inputs automatically, which is a natural way to hunt for edge cases in a sampler.
• Sample an audit trail. If you keep a tamper-evident log, as in our audit log tutorial, a reservoir with random.SystemRandom() gives reviewers a fair, unpredictable sample of it. Measure summary statistics with the help of our descriptive statistics tutorial.

Sources

• Wikipedia, Reservoir sampling: definition, Algorithm R (credited to Waterman and to McLeod and Bellhouse), Algorithm L, A-Res and A-ExpJ.
• PostgreSQL source, src/backend/commands/analyze.c and src/backend/utils/misc/sampling.c, read on November 3, 2026, and the PostgreSQL documentation for ANALYZE and default_statistics_target.
• Python documentation: random, math and itertools.
• Jeffrey S. Vitter, “Random sampling with a reservoir,” ACM Transactions on Mathematical Software 11(1), 1985, doi:10.1145/3147.3165.
• Kim-Hung Li, “Reservoir-Sampling Algorithms of Time Complexity O(n(1+log(N/n))),” ACM Transactions on Mathematical Software 20(4), 1994, doi:10.1145/198429.198435.
• Pavlos S. Efraimidis and Paul G. Spirakis, “Weighted random sampling with a reservoir,” Information Processing Letters 97(5), 2006, doi:10.1016/j.ipl.2005.11.003.