Why quantiles matter beyond the mean

Most A/B test analyses center on average metric changes, which map cleanly to business value and come with convenient mathematical properties for uncertainty quantification. But averages hide how different user segments respond. A change that improves app load time may disproportionately help users on slower devices; an engagement initiative might only lift highly active users while leaving the least engaged untouched. Quantiles — and specifically differences in quantiles — reveal such distributional shifts. The catch is that quantiles lack the neat analytic properties of means, so quantifying uncertainty in quantile changes requires more robust methods. Bootstrap inference is one such approach, though its computational cost has historically limited its use with large samples.

How Poisson bootstrap changes the picture

Bootstrap inference, introduced by Bradley Efron in 1979, offers a non-parametric way to estimate the sampling distribution of any statistic by resampling from the observed data with replacement. The logic: treat the sample as if it were the population, draw repeated resamples, record the statistic of interest each time, and use the collection of estimates to build confidence intervals.

Consider a sample of five observations: 4, 8, 11, 13, 15. A bootstrap resample might draw 11, 11, 13, 13, 15, whose median is 13. Repeating this thousands of times yields a distribution of bootstrap medians, from which a (1-alpha)100% confidence interval is obtained by taking the alpha/2 and (1-alpha/2) quantiles of that distribution. The theoretical validity of bootstrap for quantiles is well established (see Ghosh et al., 1984; Falk and Reiss, 1989).

The catch is computational. Standard bootstrap requires repeated resampling, which becomes intractable for the millions or hundreds of millions of observations common in tech-industry A/B tests. Several workarounds exist. The Poisson bootstrap, proposed by Hanley and MacGibbon (2006) and refined by Chamandy et al. (2012), replaces the exact resampling with Poisson-distributed weights per observation, enabling scalable implementations via map-reduce (Dean and Ghemawat, 2008). The “little bag of bootstrap” (Kleiner et al., 2014) splits the data into subsets, applies Poisson bootstrap to each, and combines results. Some companies, like Netflix, compress data before bootstrapping. Most others rely on normal approximations, such as results from Liu et al. (2019), but those parametric shortcuts depend heavily on assumptions about the data-generating process.

Unlocking scalable quantile bootstrap

Despite these advances, bootstrap inference for quantiles at scale remains computationally demanding because the standard approach treats bootstrap as a black box: resample, compute the quantile, repeat. The key insight is that the structure of the Poisson bootstrap, combined with properties of quantile estimators, can be exploited to cut the complexity dramatically.

In the Poisson bootstrap, each observation receives a Poisson-distributed weight. For quantile estimation, this weighting has a special property: the quantile of the weighted sample depends on the cumulative weight distribution, not on any particular ordering tied to the resampling process. This means the entire bootstrap distribution of a quantile — or a difference in quantiles — can be derived analytically from the sorted observations and the Poisson weights, without explicitly generating thousands of resampled datasets.

The upshot is a two-pronged efficiency gain. First, no explicit resampling is needed; the bootstrap replicate estimates emerge from a single weighted pass over the sorted data. Second, because the Poisson weights for independent observations are independent but their sums follow tractable distributions, quantile confidence intervals can be computed directly. At Spotify, this approach has made it routine to run bootstrap confidence intervals for difference-in-quantiles in A/B tests involving hundreds of millions of observations — a scale that previously required complex distributed implementations or parametric shortcuts.

Turning the Bootstrap Inside Out

The cost of bootstrapping quantiles at scale comes from the resampling step itself. Rebuilding resampled datasets and locating the relevant order statistic in each one quickly becomes prohibitive as data volumes grow. The key insight behind a faster alternative is that, for a given quantile, the bootstrap can be re-expressed as a problem about which index in the original sorted data gets selected per bootstrap round — not about the data values at all.

The estimator that unlocks this view is a slightly crude quantile estimator. For a sample of size N and a quantile of interest q, the estimator is the order statistic at position (N+1)q. When (N+1)q is not an integer, one of the two neighboring order statistics is chosen at random with equal probability. In a sample of 100 values, for example, the median is either the 50th or the 51st ordered observation, each with 50% probability.

With this estimator, a Poisson bootstrap sample needs just one piece of information from the original data: which index in the ordered sample lands at the target quantile. Simulating the bootstrap then reduces to studying the distribution of that index across repeated rounds. Empirically, with a sample of size 200 and q = 0.2, most ordered positions are never selected as the bootstrap quantile across 1 million bootstrap samples, and the distribution of the selected indexes closely tracks a binomial with parameters N+1 and q. For the difference-in-quantiles case, this yields a stark simplification: sort the data once, then draw indexes from that binomial distribution and read off the corresponding order statistics to form each bootstrap estimate. Details of the complexity reduction are in the paper by Schultzberg and Ankargren (2022).

In the one-sample case, this result converges with classical order-statistic-based confidence intervals for quantiles described by Gibbons and Chakraborti (2010), but it derives them from the bootstrap perspective. The extra value is that knowing the full index distribution is what makes two-sample comparisons possible. The code below shows how the binomial approximation turns a one-sample quantile confidence interval into a handful of lines.

import numpy as np
from scipy.stats import binom

alpha=.05
quantile_of_interest=0.5
sample_size=10000
number_of_bootstrap_samples=1000000
outcome_sorted = np.sort(np.random.normal(1,1,sample_size))

ci_indexes = binom.ppf([alpha/2,1-alpha/2],sample_size+1, quantile_of_interest)
bootstrap_confidence_interval = outcome_sorted[[int(np.floor(ci_indexes[0])), int(np.ceil(ci_indexes[1]))]]

f"The sample median is {np.quantile(outcome_sorted, quantile_of_interest)}, the {(1-alpha)*100}%\
confidence interval is given by ({bootstrap_confidence_interval})."

The one-sample computation is essentially a single sort, because the boundaries of the confidence interval correspond to the alpha/2 and 1-alpha/2 tails of the binomial index distribution. For the two-sample case, which is the primary contribution, a difference still has to be computed per bootstrap iteration. That implementation is similarly compact.

import numpy as np
from numpy.random import normal, binomial

alpha=.05
quantile_of_interest=0.5
sample_size=10000
number_of_bootstrap_samples=1000000
outcome_control_sorted = np.sort(normal(1,1,sample_size))
outcome_treatment_sorted = np.sort(normal(1.2,1,sample_size))

bootstrap_difference_distribution = outcome_treatment_sorted[binomial(sample_size+1, quantile_of_interest,
number_of_bootstrap_samples)] - outcome_control_sorted[binomial(sample_size+1,
                        quantile_of_interest, number_of_bootstrap_samples)]
bootstrap_confidence_interval = np.quantile(bootstrap_difference_distribution,
[alpha/2 , 1-alpha/2])

f"The sample difference-in-medians is \
{np.quantile(outcome_treatment_sorted, quantile_of_interest)-np.quantile(outcome_control_sorted, quantile_of_interest)},\
the {(1-alpha)*100}% confidence interval for the difference-in-medians is given by ({bootstrap_confidence_interval})."

Simulation confirms that the binomial-approximated intervals behave like standard bootstrap intervals and have correct statistical properties. The code below checks the false positive rate of a one-sided confidence interval for the median.

import numpy as npfrom scipy.stats import binom

alpha=.05quantile_of_interest=0.5sample_size=10000number_of_bootstrap_samples=1000000replications = 10000

bootstrap_confidence_intervals = []ci_index = int(np.floor(binom.ppf(alpha,sample_size+1,quantile_of_interest)))for i in range(replications):outcome_sorted = np.sort(np.random.normal(0,1,sample_size))bootstrap_confidence_intervals.append(outcome_sorted[ci_index])

f"The empirical false positive rate of the test using the bootstrap confidence interval is {np.mean([1 if i>0  else 0 for i in bootstrap_confidence_intervals])*100}%, the intended false positive rate is {alpha*100}%"

The practical impact of removing resampling is considerable. In Julia benchmarks for a two-sample median confidence interval with a total of 2,000 observations and 10,000 bootstrap samples, a standard Poisson bootstrap implementation had a median runtime of 1,821 milliseconds and consumed 2.4 GB of memory. The binomial-based version finished in 2.2 milliseconds using 407 KiB. A faster Julia script for comparative simulations for the difference-in-quantiles case is available in this GitHub repository.

A SQL-Friendly Procedure

The binomial approximation also exposes an implementation that fits naturally into a data warehouse pipeline. The confidence interval can be computed in two steps:

  1. In Python, calculate the order-statistic indexes corresponding to the alpha/2 upper and lower tails of the Binom(N+1, q) distribution.
  2. In SQL, return the order statistics from the outcome table at those indexes.

This keeps the logic simple and pushes only two values into the query, avoiding materializing any resampled data.

Practical Inference for Quantile Metrics

Quantile metrics like latency percentiles are standard in A/B testing, but their inference has historically been bottlenecked by computationally heavy resampling. The index-based approach replaces that machinery with a binomial approximation, making difference-in-quantiles bootstrap confidence intervals tractable on data sizes encountered in production experimentation platforms. This opens up assumption-free inference for quantile metrics in dashboards and experiments at scale.