My log analyzer reported different p99 request durations depending on how many worker processes parsed the files. Same input, same code, different percentiles — and the numbers looked plausible every time. The cause was a three-line merge that replayed one reservoir sample into another.
TL;DR — Merging two reservoir samples with for v in other.buf: self.add(v) treats each buffered value as a single observation, so the receiving side crushes the incoming side no matter which way you merge. A correct merge fills the buffer by drawing from each side with probability proportional to its true population count n, then sets n = n_a + n_b.
The setup
The program is a log analyzer: a Python script that reads web server access logs and reports traffic statistics, including request-duration percentiles (p50, p95, p99). Millions of durations do not fit in memory, so it uses reservoir sampling: keep a buffer of at most k = 2,048 values, count every observation in n, and once the buffer is full let each new value overwrite a random slot with probability k/n. The buffer stays a uniform sample of everything seen, and sorting it approximates the percentiles.
With --workers N, each worker process builds its own reservoir over its share of the log files, and the parent merges them. The merge looked like this:
def merge(self, other: "Reservoir") -> None:
for v in other.buf:
self.add(v)
What I expected
Replaying the buffer through add() looks right because add() is exactly how the values entered in the first place. I expected the merged reservoir to be a fair sample of both populations, n to end up at the combined total, and A.merge(B) to be statistically the same as B.merge(A).
What actually happens
I fed two artificial populations through the real class: A = 1,000,000 copies of the value 100 (fast requests) and B = 60,000 copies of the value 10,000 (slow requests). Exact percentiles of the pooled 1,060,000 values: p50 = 100, p95 = 10,000, p99 = 10,000.
Merged percentiles [p50, p95, p99], identical in all 9 runs each way:
A.merge(B): [100, 100, 100] n = 1,002,048
B.merge(A): [10000, 10000, 10000] n = 62,048
exact: [100, 10000, 10000] n = 1,060,000
The two orders disagree on every percentile, and neither n is the true count. A.merge(B) erased the slow requests: of B’s 2,048 buffered values, only 1 to 7 survived into the merged buffer across nine runs, where a fair sample of the pooled data holds about 116. B.merge(A) kept its slow values and reported 10,000 as the median of a population that is 94% fast.
In the analyzer this merge ran once per worker result, delivered by multiprocessing.Pool.imap_unordered — so the merge order, and with it the reported percentiles, depended on --workers and on which worker happened to finish first. Like my optimal-solver tests that all passed for a solver that was not optimal, nothing crashed and every individual number looked believable.
The fix
Weight each side by its own n. Fill the output buffer slot by slot, taking the next value from self with probability n_self / (n_self + n_other), otherwise from the other side; shuffling each buffer once and popping gives sampling without replacement:
import random
def merge(self, other: "Reservoir") -> None:
rng = random.Random()
a, b = list(self.buf), list(other.buf)
rng.shuffle(a)
rng.shuffle(b)
total = self.n + other.n
out = []
for _ in range(min(self.k, len(a) + len(b))):
take_self = rng.random() < self.n / total
if (take_self and a) or not b:
out.append(a.pop())
else:
out.append(b.pop())
self.buf = out
self.n = total
With this merge, the same experiment reported [100, 10000, 10000] with n = 1,060,000 in 8 of 9 runs. The ninth reported p95 = 100 — sampling noise, not a bug: the slow population is 5.66% of the total, so the true p95 sits close to the fast/slow boundary and a 2,048-value sample lands on the wrong side of it occasionally.
Why it works
Each buffer is already a uniform sample of its own population. A uniform sample of the union follows from letting every output slot come from side A with probability n_a / (n_a + n_b) — the probability that a uniformly chosen element of the union belongs to A.
The replay version discards that weighting. Each of B’s 2,048 buffered values stands in for 60,000 / 2,048 ≈ 29 real observations, but add() offers it as one observation to a reservoir whose counter is already past 1,000,000, so it is kept with probability ≈ 2048 / 1,002,048 ≈ 0.2%. Expected survivors: 2,048 × 0.2% ≈ 4, matching the 1 to 7 I observed. And because add() increments n once per buffered value instead of adding other.n, the merged count lands at 1,002,048 instead of 1,060,000, which poisons any later merge or weighting that trusts n.
The bug is easy to miss in tests. While a reservoir holds fewer than k values, its buffer is its entire population, and replaying it through add() is exactly correct. Small test inputs never fill the buffer, so the merge only starts lying past 2,048 values per worker — that is, on real data.
What I did not test
- Only a two-way merge of two constant-valued populations. Real duration mixes will shift less dramatically, but the order dependence is structural, not an artifact of the extreme inputs.
- The fix draws with a fixed ratio instead of an exact hypergeometric draw (the ratio should shift slightly as values are consumed). With k = 2,048 and populations in the millions, the difference stayed below what nine runs could detect.
- Merging more than two reservoirs pairwise, and merges where both buffers are still partially filled.
- CPython 3.14 only.
Facts
context: parallel log analyzer computing request-duration percentiles via reservoir sampling (k=2048), one reservoir per worker, merged in the parent
problem: merge implemented as replaying other.buf through add() — order-dependent, erases the smaller side, and undercounts n, so percentiles change with worker count
solution: fill the merged buffer by drawing from each side with probability n_a/(n_a+n_b) without replacement, then set n = n_a + n_b
verified_on: 2026-08-17
applies_to: [any reservoir-sampling implementation merged across shards or processes, Python 3.14]
does_not_apply_to: [reservoirs that never fill (n <= k), where replaying the buffer is exact]