Learning Hub

How to Merge Latency Percentiles Across Servers in Python With a DDSketch and Stop Averaging p99

Averaging p99 across servers gives the wrong answer, so build a 100-line mergeable DDSketch in Python that keeps every percentile within 1 percent and merges by addition.

November 1, 2026 42 Min Read
26

Dashboards love one number: the p99 latency, the response time that 99 percent of requests beat. Reading it is easy. Computing it correctly gets hard the moment more than one server is involved. Fifty web servers each see a slice of the traffic, each can report its own p99, and the dashboard wants a single number for the whole fleet. Averaging the fifty p99s feels natural, and it is wrong, sometimes by a wide margin. The same trap shows up over time: a database that stores one p99 per minute cannot give you an honest p99 for the last hour.

Table Of Content

In this tutorial you will watch that failure happen with real numbers, then build the fix. The fix is a DDSketch: a small summary of a stream of numbers that keeps every percentile within a relative error you choose (1 percent here) and that merges with other sketches by plain addition. You will write one from scratch in about 100 lines of standard-library Python, test it, check it against Datadog’s reference package, and then recognize the same idea inside OpenTelemetry’s exponential histograms. DDSketch comes from the paper DDSketch: A Fast and Fully-Mergeable Quantile Sketch with Relative-Error Guarantees by Charles Masson, Jee E. Rim and Homin K. Lee of Datadog (VLDB 2019).

This is not a niche problem. The Prometheus guide to histograms and summaries lists, as a limit of pre-calculated quantiles, that “you cannot aggregate quantiles (e.g. to calculate the total 90th percentile latency for a service backed by multiple replicated workers)”. OpenTelemetry’s metrics data model says something similar about its legacy summary type: summary points “cannot always be merged in a meaningful way”. Both projects point toward the same remedy, which is to ship a distribution instead of a percentile. A DDSketch is one way to do that, and it is small enough to understand completely.

By the end you will be able to:

  • show, with numbers, why averaging percentiles gives the wrong answer;
  • explain what a fixed-bucket histogram costs you in accuracy;
  • build a mergeable sketch whose every percentile is within 1 percent of the truth, and prove that with a test;
  • merge sketches from many servers and from many minutes, and ship them as JSON safely;
  • spot the traps (mixed accuracy settings, bad values, wide value ranges, untrusted payloads) before they cost you a dashboard.

The terms used in this tutorial

Six terms appear all the time, so here they are once:

  • Quantile and percentile. A quantile answers “what value do a fraction q of my observations fall at or below?”, where q is between 0 and 1. The median is the 0.5 quantile. A percentile is the same idea out of 100, so p99 is the 0.99 quantile. If your p99 is 470 ms, then 99 of every 100 requests finished in 470 ms or less and 1 of every 100 took longer.
  • Exact versus estimated. An exact quantile needs every value kept and sorted. An estimated quantile comes from a summary and is allowed to be slightly off.
  • Relative error. The size of the error as a fraction of the true value. “Within 1 percent” means that when the true p99 is 400 ms, the estimate is somewhere between 396 and 404 ms. The paper writes this error as α (alpha), and the code calls it alpha.
  • Bucket and histogram. A bucket is a counter for the values that fall in one range, for example 100 ms to 250 ms. A histogram is a set of buckets.
  • Sketch. A small summary of a stream of values, built in one pass, that answers questions approximately.
  • Mergeable. Two sketches built from different data can be combined into one sketch that is as good as a single sketch that saw all the data.

Why latency needs a relative error

Many older quantile sketches promise a small error in rank: ask for the 99th percentile and you get a value somewhere between, say, the 98.5th and the 99.5th percentile. The DDSketch paper points out why that is a weak promise for response times, which are heavily skewed. On its example of 2 million web requests, such a window for the p99 spanned “anywhere from 2 to 20 seconds”. A promise of “within 1 percent of the true value” is far more useful, because it stays meaningful whether the value is 5 milliseconds or 5 seconds. That is the promise DDSketch makes.

Prerequisites

  • Python 3.13 or newer. I used Python 3.13.14 on Windows 11. Linux and macOS behave the same. Where I write python, use python3 if that is your command name. The sketch itself needs only the standard library.
  • pytest 9.2.0 for the tests, plus two packages used only in the last two steps: ddsketch 3.0.1 (Datadog’s reference implementation) and opentelemetry-sdk 1.45.2.
  • Basic Python: functions, classes, dictionaries and f-strings. Everything else is explained where it first appears.
  • About an hour to read and run everything. The whole lab runs in under ten seconds on my machine.

Create an empty folder called percentile_lab and work inside it. Every file goes in that folder, and every command assumes you run it from there. A fresh virtual environment keeps things tidy:

python -m venv .venv
# Windows: .venv\Scripts\activate    Linux or macOS: source .venv/bin/activate
pip install pytest==9.2.0 ddsketch==3.0.1 opentelemetry-sdk==1.45.2

Every step uses seeded random numbers, so on Python 3.13 you should see the same digits that appear below. Other Python versions may differ in the last digits. Only timings vary from machine to machine.

Step 1: Make some request times and compute exact percentiles

Before building anything clever, you need data and a way to get the right answer, so that every later estimate has something honest to be compared with. Create labkit.py. It holds four small helpers that every step shares.

"""labkit.py: shared helpers for the percentile lab (latency generator, exact quantile, relative error)."""
import math
import random


def make_latencies(n, seed, median=0.080, sigma=0.55, slow_share=0.01, slow_factor=12.0):
    """Return n request times in seconds: a log-normal body plus a small share of slow requests."""
    rng = random.Random(seed)
    mu = math.log(median)
    times = []
    for _ in range(n):
        t = rng.lognormvariate(mu, sigma)
        if rng.random() < slow_share:
            t *= slow_factor * rng.uniform(0.5, 2.0)
        times.append(t)
    return times


def make_fleet():
    """Three servers with very different traffic and tails. Steps 2 and 6 use this same fleet."""
    return {
        "web-1": make_latencies(95_000, seed=11, median=0.080, slow_share=0.01),
        "web-2": make_latencies(4_000, seed=12, median=0.085, slow_share=0.08),   # a sick canary
        "web-3": make_latencies(1_000, seed=13, median=0.300, slow_share=0.02),   # old hardware
    }


def exact_quantile(sorted_values, q):
    """The 'lower' quantile used in the DDSketch paper: the item of rank 1 + floor(q * (n - 1))."""
    return sorted_values[int(q * (len(sorted_values) - 1))]


def rel_err(estimate, truth):
    """Relative error |estimate - truth| / truth."""
    return abs(estimate - truth) / truth

What each helper does and why:

  • make_latencies invents request times in seconds. Real response times are positive and skewed: most requests cluster around a typical value and a few are many times slower (cold caches, retries, garbage collection pauses). A log-normal distribution has that shape, and the slow_share option adds a small crowd of extra-slow requests on top. The fixed seed makes every run repeat exactly.
  • make_fleet builds three servers with different traffic and different tails: a healthy majority, a sick canary (8 percent slow requests), and an old machine. I made the fleet unfair on purpose, because real fleets are.
  • exact_quantile returns the item of rank 1 + floor(q × (n − 1)) in a sorted list. That is the “lower” quantile definition the DDSketch paper uses. Other definitions interpolate between neighbors, and on large data the difference is tiny (the descriptive statistics tutorial on this site shows how much it matters on small data). We use the paper’s definition so that the guarantee in Step 5 is stated about the same number.
  • rel_err is the relative error, |estimate - truth| / truth.

Now create step01_baseline.py. It keeps all 100,000 request times of one server, sorts them, and reads off exact percentiles:

"""step01_baseline.py: one server, keep every value, read exact percentiles from the sorted list."""
import statistics
import sys

from labkit import exact_quantile, make_latencies

times = make_latencies(100_000, seed=1)
ordered = sorted(times)

print(f"requests: {len(times):,}")
print(f"mean:  {statistics.fmean(times) * 1000:7.1f} ms")
for name, q in [("p50", 0.50), ("p90", 0.90), ("p99", 0.99), ("p99.9", 0.999), ("max", 1.0)]:
    print(f"{name:<6}{exact_quantile(ordered, q) * 1000:7.1f} ms")

per_float = (sys.getsizeof(times) + sum(sys.getsizeof(t) for t in times)) / len(times)
print(f"memory: {per_float:.0f} bytes per request, so {per_float * 1_000_000 / 1e6:.0f} MB per million requests")

Run it:

python step01_baseline.py
requests: 100,000
mean:    106.5 ms
p50      80.6 ms
p90     166.3 ms
p99     467.7 ms
p99.9  2696.3 ms
max   10774.0 ms
memory: 32 bytes per request, so 32 MB per million requests

Look at what the numbers say. The mean is 106.5 ms, which sits between the median (p50, 80.6 ms) and the p90 (166.3 ms). It describes nobody’s experience: half the requests were faster than 80.6 ms, and one in a hundred took 467.7 ms or longer. The tail stretches further: p99.9 is about 2.7 seconds and the slowest request took 10.8 seconds. This is why teams watch percentiles, not averages.

The last line shows the cost of exactness. Every request time takes about 32 bytes as a Python float in a list (24 for the float object and 8 for the pointer), so a million requests need about 32 MB. A service handling 5,000 requests per second would produce 18 million values an hour, roughly 576 MB at that rate (my arithmetic, not a measurement). And exact percentiles need the values from every server in one place before they can be sorted. Keep that in mind for Step 6.

Checkpoint: your output should match the numbers above. If your mean or percentiles differ, check that labkit.py was copied exactly, including the seed.

Step 2: See why averaging percentiles gives the wrong answer

Now the failure that motivates everything else. Three servers each compute their own p99, which is cheap and is what pre-calculated quantiles, such as a Prometheus summary, give you. The question is what single number to report for the fleet. Create step02_average_of_p99.py:

"""step02_average_of_p99.py: three servers each report a p99. Which combination gives the fleet p99?"""
from labkit import exact_quantile, make_fleet, rel_err

fleet = make_fleet()
p99 = {name: exact_quantile(sorted(times), 0.99) for name, times in fleet.items()}
sizes = {name: len(times) for name, times in fleet.items()}
for name in fleet:
    print(f"{name}: {sizes[name]:>6,} requests, p99 = {p99[name] * 1000:6.1f} ms")

everything = sorted(t for times in fleet.values() for t in times)
truth = exact_quantile(everything, 0.99)
print(f"\nfleet p99 from all {len(everything):,} raw values: {truth * 1000:6.1f} ms\n")

total = sum(sizes.values())
guesses = {
    "average of the p99s": sum(p99.values()) / len(p99),
    "request-weighted average": sum(p99[n] * sizes[n] for n in fleet) / total,
    "highest p99": max(p99.values()),
}
for label, value in guesses.items():
    print(f"{label:<26}{value * 1000:7.1f} ms  ({rel_err(value, truth):6.1%} off)")
python step02_average_of_p99.py
web-1: 95,000 requests, p99 =  482.9 ms
web-2:  4,000 requests, p99 = 2445.9 ms
web-3:  1,000 requests, p99 = 2259.7 ms

fleet p99 from all 100,000 raw values:  770.3 ms

average of the p99s        1729.5 ms  (124.5% off)
request-weighted average    579.1 ms  ( 24.8% off)
highest p99                2445.9 ms  (217.5% off)

The line that matters is the fleet p99 computed from all 100,000 raw values: 770.3 ms. That is the right answer, and it needed every raw value. Now compare the three tempting shortcuts:

  • The plain average of the three p99s gives 1,729.5 ms, more than double the truth. It treats web-3’s 1,000 requests as if they counted as much as web-1’s 95,000.
  • The request-weighted average fixes that and gets much closer, at 579.1 ms, but it is still 24.8 percent too low. It cannot see that web-2’s slow requests (8 percent of its traffic) pile up above the fleet’s p99.
  • The highest p99 is 217.5 percent too high. It is a worst case, not a percentile.

No arithmetic on the three p99 numbers can fix this. A percentile is a position in a sorted list. The p99 of web-1 says that 1 percent of web-1’s requests were slower than 482.9 ms. It says nothing about how slow those requests were, or how web-2’s slow requests compare. The fleet p99 depends on the whole shape of every server’s distribution, not just on one point of it. That is exactly the limit the Prometheus guide names when it says quantiles cannot be aggregated.

Checkpoint: the fleet p99 is 770.3 ms and all three shortcuts miss it. If you ever see a dashboard that averages per-host percentiles, you now know what it is hiding.

Step 3: See what a fixed-bucket histogram costs you

The standard fix, used by Prometheus for years, is to stop reporting percentiles and report counts per bucket instead. A classic Prometheus histogram has a list of upper edges chosen in advance, say 0.1, 0.25 and 0.5 seconds, and each server counts how many requests were at or below each edge. Counts add up across servers, so histograms merge. The percentile is estimated later from the merged counts by linear interpolation: find the bucket that holds the requested rank, then assume the requests inside it are spread evenly between its two edges.

That works, and it hides a cost that is easy to miss. Create step03_fixed_buckets.py. It implements the estimator and runs two experiments:

"""step03_fixed_buckets.py: a classic fixed-bucket histogram and the error that its bucket edges cause."""
import bisect

from labkit import exact_quantile, make_latencies, rel_err

INF = float("inf")


def classic_quantile(bounds, counts, q):
    """Prometheus-style histogram_quantile: find the bucket holding the rank, then interpolate inside it."""
    total = sum(counts)
    rank = q * total
    running = 0
    for i, c in enumerate(counts):
        if running + c >= rank:
            upper = bounds[i]
            if upper == INF:
                return bounds[i - 1]  # nothing to interpolate to: report the highest finite edge
            lower = bounds[i - 1] if i > 0 else 0.0
            return lower + (upper - lower) * (rank - running) / c
        running += c
    raise AssertionError("unreachable")


def fill(bounds, values):
    counts = [0] * len(bounds)
    for v in values:
        counts[bisect.bisect_left(bounds, v)] += 1  # buckets are 'less than or equal to the edge'
    return counts


# 1. The thought experiment from the Prometheus histogram guide: every request takes about 220 ms.
doc_bounds = [0.1, 0.2, 0.3, 0.45, INF]
for spike in (0.22, 0.32):
    counts = fill(doc_bounds, [spike] * 1000)
    print(f"all requests at {spike * 1000:.0f} ms -> estimated p95 = {classic_quantile(doc_bounds, counts, 0.95) * 1000:.1f} ms")

# 2. The default buckets of the Python Prometheus client, on our 100,000 requests, then on a 3x faster service.
DEFAULTS = [0.005, 0.01, 0.025, 0.05, 0.075, 0.1, 0.25, 0.5, 0.75, 1.0, 2.5, 5.0, 7.5, 10.0, INF]
times = make_latencies(100_000, seed=1)
for label, scale in [("today", 1.0), ("after a 3x speed-up", 1 / 3)]:
    values = [t * scale for t in times]
    ordered = sorted(values)
    counts = fill(DEFAULTS, values)
    print(f"\n{label}:")
    for name, q in [("p50", 0.5), ("p90", 0.9), ("p99", 0.99), ("p99.9", 0.999)]:
        truth = exact_quantile(ordered, q)
        guess = classic_quantile(DEFAULTS, counts, q)
        print(f"  {name:<6} true {truth * 1000:7.1f} ms   histogram {guess * 1000:7.1f} ms   {rel_err(guess, truth):6.1%} off")
python step03_fixed_buckets.py
all requests at 220 ms -> estimated p95 = 295.0 ms
all requests at 320 ms -> estimated p95 = 442.5 ms

today:
  p50    true    80.6 ms   histogram    81.4 ms     1.0% off
  p90    true   166.3 ms   histogram   216.7 ms    30.4% off
  p99    true   467.7 ms   histogram   495.4 ms     5.9% off
  p99.9  true  2696.3 ms   histogram  2838.5 ms     5.3% off

after a 3x speed-up:
  p50    true    26.9 ms   histogram    28.1 ms     4.8% off
  p90    true    55.4 ms   histogram    59.0 ms     6.5% off
  p99    true   155.9 ms   histogram   216.8 ms    39.1% off
  p99.9  true   898.8 ms   histogram   922.1 ms     2.6% off

Experiment 1 recreates a thought experiment from the Prometheus guide. Every request takes about 220 ms, so the true p95 is about 220 ms. With classic buckets at 0.1, 0.2, 0.3 and 0.45 seconds, all the requests land in the bucket from 200 to 300 ms, and the interpolated p95 comes out at 295.0 ms. Move every request to 320 ms (the guide adds a fixed 100 ms of routing delay) and the estimate is 442.5 ms, which the guide rounds to 443 ms. Those are the guide’s own numbers, reproduced. The guide’s comparison table says the same about classic histograms: “Error is limited by the width of the bucket the quantile is located in.”

Experiment 2 uses the default buckets of the Python Prometheus client. The client’s documentation says its default buckets “are intended to cover typical web/RPC request latency in seconds”, and the list is 0.005, 0.01, 0.025, 0.05, 0.075, 0.1, 0.25, 0.5, 0.75, 1, 2.5, 5, 7.5, 10 and an overflow bucket. On our 100,000 requests the p50 is off by 1.0 percent but the p90 is off by 30.4 percent, because the true p90 (166.3 ms) lies in the 150 ms wide bucket between 0.1 and 0.25 seconds. Then the service gets three times faster, which is good news, and the p99 estimate becomes 39.1 percent wrong because the p99 has moved into that same wide bucket. The service improved and the dashboard got less accurate.

Two lessons follow. First, a fixed-bucket histogram’s accuracy depends on edges you chose before you saw the data, and data changes. Second, merging classic histograms only works when every server uses identical edges; the Prometheus guide’s comparison table says classic histograms can be aggregated “provided there are no changes in bucket boundaries”. What we want is buckets whose relative width is the same everywhere, so that the error is the same at 5 ms and at 5 seconds. That is the next step.

Checkpoint: you should see 295.0 and 442.5 ms in the first experiment, and an error of 30.4 percent at p90 followed by 39.1 percent at p99 after the speed-up.

Step 4: Build the sketch from buckets that grow by a fixed ratio

Here is the whole idea of DDSketch in one sentence: instead of buckets of equal width (0 to 100 ms, 100 to 200 ms, and so on), use buckets whose width grows with the value, so that every bucket is the same percentage wide. A bucket near 5 ms is about a tenth of a millisecond wide. A bucket near 60 seconds is over a second wide. If every bucket is about 2 percent wide, then reporting the right point inside a bucket is wrong by at most 1 percent, no matter where the value lives.

The few formulas you need

Pick the accuracy you want, α (alpha). For 1 percent, α = 0.01. From it, compute the growth ratio γ (gamma):

  • γ = (1 + α) / (1 − α). For α = 0.01 that is 1.020202, so each bucket’s upper edge is about 2 percent above its lower edge.
  • Bucket i holds the values x with γi−1 < x ≤ γi. To find i for a value, take the logarithm of the value and divide by the logarithm of γ, then round up: i = ceil(log(x) / log(γ)). The index can be zero or negative: values smaller than 1 land in bucket 0 or below.
  • The answer for bucket i is 2γi / (γ + 1). This is not an edge. It is the point whose relative distance to both edges is the same.

A line of algebra shows why that point keeps the error at or below α. Write the upper edge as U = γi and the lower edge as L = U / γ, and let E = 2U / (γ + 1) be the answer. A value sitting exactly at the upper edge is off by (U − E) / U = 1 − 2 / (γ + 1) = (γ − 1) / (γ + 1). A value sitting exactly at the lower edge is off by (E − L) / L = 2γ / (γ + 1) − 1 = (γ − 1) / (γ + 1). Replace γ by (1 + α) / (1 − α) and (γ − 1) / (γ + 1) simplifies to exactly α. Values in the middle of the bucket are off by less. The paper states the same result as its Lemma 2.

The first half of sketch.py

Create sketch.py and add this first part. The class stores nothing but a dictionary from bucket index to count, plus two more counters:

"""sketch.py: a small mergeable DDSketch for request times, after Masson, Rim and Lee (VLDB 2019)."""
import bisect
import heapq
import math

MIN_INDEXABLE = 1e-9  # seconds: anything at or below this goes into the zero bucket


class DDSketch:
    """Counts positive numbers in log-spaced buckets, so every quantile is within `alpha` of the truth."""

    def __init__(self, alpha=0.01, max_buckets=2048):
        if not 0 < alpha < 1:
            raise ValueError("alpha must be between 0 and 1")
        self.alpha = alpha
        self.gamma = (1 + alpha) / (1 - alpha)
        self._ln_gamma = math.log(self.gamma)
        self._max_index = int(700 / self._ln_gamma)  # keeps gamma ** index inside the float range
        self.max_buckets = max_buckets
        self.buckets = {}      # bucket index -> how many values landed in it
        self.zero_count = 0    # values at or below MIN_INDEXABLE
        self.count = 0

    def _index(self, x):
        return math.ceil(math.log(x) / self._ln_gamma)

    def _value(self, i):
        return 2 * self.gamma ** i / (self.gamma + 1)

    def _collapse(self):
        """Fold the lowest bucket into its neighbour until at most max_buckets remain (Algorithm 3 in the paper)."""
        while len(self.buckets) > self.max_buckets:
            lowest, second = heapq.nsmallest(2, self.buckets)
            self.buckets[second] += self.buckets.pop(lowest)

    def add(self, x):
        if not math.isfinite(x) or x < 0:
            raise ValueError(f"cannot add {x!r}: need a finite number that is 0 or more")
        if x <= MIN_INDEXABLE:
            self.zero_count += 1
        else:
            i = self._index(x)
            if abs(i) > self._max_index:
                raise ValueError(f"cannot add {x!r}: outside the range this sketch can index")
            self.buckets[i] = self.buckets.get(i, 0) + 1
            self._collapse()
        self.count += 1

Piece by piece:

  • MIN_INDEXABLE exists because the logarithm of zero is undefined. Anything at or below one nanosecond goes into a separate zero bucket, a plain counter. Request times of exactly 0.0 happen when a clock is coarse, and the sketch must not crash on them. The unit is seconds, which is also the unit Prometheus asks for (more on that in the traps step).
  • __init__ checks that alpha is between 0 and 1, computes γ, and precomputes log(γ) so that each insert needs only one logarithm. It also computes _max_index. A float cannot hold numbers above roughly 1.8e308, and gamma ** i would raise an OverflowError for a huge index, so we refuse indexes beyond what 700 / log(γ) allows. That turns a later crash into a clear error now.
  • buckets, zero_count and count are the entire state: bucket index to count, the zero bucket, and the total number of values.
  • _index and _value are the two formulas from the list above, written as code.
  • _collapse keeps the sketch from growing without limit: if there are more than max_buckets buckets, it folds the lowest bucket into its neighbor. For now you can ignore it. It does nothing until the sketch has more than 2,048 buckets, and Step 8 explains it properly.
  • add rejects values that are not finite or are negative (a negative duration means something is wrong with the timer, and a sketch that silently swallowed it would hide the bug), puts tiny values in the zero bucket, and otherwise increments the right bucket.

Now see the buckets at work. Create step04_buckets.py. It shows where six request times land and how far each bucket’s answer is from the real value:

"""step04_buckets.py: where do a few request times land, and how far is each bucket's answer from the truth?"""
from sketch import DDSketch

sketch = DDSketch(alpha=0.01)
print(f"alpha = {sketch.alpha}, gamma = {sketch.gamma:.6f}\n")
print(f"{'time':>9}  {'bucket':>6}  {'the bucket covers':>25}  {'answer':>10}  {'error':>6}")
for ms in (0.2, 5, 80, 220, 1500, 60000):
    x = ms / 1000
    i = sketch._index(x)
    low, high = sketch.gamma ** (i - 1), sketch.gamma ** i
    answer = sketch._value(i)
    print(f"{ms:>6} ms  {i:>6}  {low * 1000:>10.3f} to {high * 1000:>9.3f} ms  {answer * 1000:>7.3f} ms  {abs(answer - x) / x:6.2%}")
python step04_buckets.py
alpha = 0.01, gamma = 1.020202

     time  bucket          the bucket covers      answer   error
   0.2 ms    -425       0.199 to     0.203 ms    0.201 ms   0.69%
     5 ms    -264       4.991 to     5.092 ms    5.041 ms   0.81%
    80 ms    -126      78.860 to    80.453 ms   79.648 ms   0.44%
   220 ms     -75     218.701 to   223.119 ms  220.888 ms   0.40%
  1500 ms      21    1491.845 to  1521.983 ms  1506.763 ms   0.45%
 60000 ms     205   59153.515 to 60348.535 ms  59745.050 ms   0.42%

Read one row. The 80 ms request lands in bucket −126, which covers 78.860 to 80.453 ms. The sketch would answer 79.648 ms for anything in that bucket, which is 0.44 percent away from 80 ms. Now compare the widths. The bucket for 0.2 ms spans 0.199 to 0.203 ms (four thousandths of a millisecond), while the bucket for 60 seconds spans 59,153 to 60,348 ms (almost 1.2 seconds). Both are about 2 percent wide, and every error in the last column is below 1 percent.

Checkpoint: gamma prints as 1.020202 and every error is under 1.00 percent.

Step 5: Ask for percentiles by walking the buckets

A sketch cannot sort the values, because it threw them away. It can still count. The quantile q refers to the item at position floor(q × (n − 1)) in the sorted data (counting from 0). To find the bucket holding that item, walk the buckets from smallest to largest while keeping a running total, and stop at the first bucket whose running total is greater than q × (n − 1). That bucket must contain the item. Report its answer. This is Algorithm 2 in the paper.

Add this second part to sketch.py, right after add:

    # sketch.py (continued)
    def quantiles(self, qs):
        """Estimate several quantiles with one pass over the sorted buckets."""
        if self.count == 0:
            raise ValueError("the sketch is empty")
        keys = sorted(self.buckets)
        running, cumulative = self.zero_count, []
        for i in keys:
            running += self.buckets[i]
            cumulative.append(running)
        answers = []
        for q in qs:
            if not 0 <= q <= 1:
                raise ValueError("q must be between 0 and 1")
            rank = q * (self.count - 1)
            if self.zero_count > rank:
                answers.append(0.0)
            else:
                answers.append(self._value(keys[bisect.bisect_right(cumulative, rank)]))
        return answers

    def quantile(self, q):
        return self.quantiles([q])[0]

How it works:

  • quantiles takes a list of q values so that the bucket indexes are sorted and the running totals are built once, not once per question.
  • rank = q * (self.count - 1) is the position described above. If the zero bucket alone already holds more than that many values, the answer is 0.0.
  • Otherwise bisect.bisect_right(cumulative, rank) finds the first bucket whose running total is greater than rank. It is a binary search, which gives the same answer as the walk in the paper with less work. The result is turned into a number by _value.
  • quantile is a convenience wrapper for one q.

Now test the promise. Create step05_guarantee.py. It first replays the two spikes from Step 3, then feeds one million request times to a sketch and compares it with the exact answer at six quantiles and at 10,001 evenly spaced ones:

"""step05_guarantee.py: feed one million request times to a sketch and compare its quantiles with the exact ones."""
from labkit import exact_quantile, make_latencies, rel_err
from sketch import DDSketch

# First the spike test from step 3: every request takes the same time, and we ask for the p95.
for spike in (0.22, 0.32):
    tiny = DDSketch(alpha=0.01)
    for _ in range(1000):
        tiny.add(spike)
    print(f"all requests at {spike * 1000:.0f} ms -> sketch p95 = {tiny.quantile(0.95) * 1000:.1f} ms")
print()

times = make_latencies(1_000_000, seed=5)
ordered = sorted(times)

sketch = DDSketch(alpha=0.01)
for t in times:
    sketch.add(t)

print(f"values added: {sketch.count:,}, buckets used: {len(sketch.buckets)}\n")
print(f"{'quantile':>9}  {'exact':>10}  {'sketch':>10}  {'error':>6}")
for q in (0.5, 0.9, 0.99, 0.999, 0.9999, 1.0):
    truth, guess = exact_quantile(ordered, q), sketch.quantile(q)
    print(f"{q:>9}  {truth * 1000:>7.2f} ms  {guess * 1000:>7.2f} ms  {rel_err(guess, truth):6.3%}")

grid = [i / 10_000 for i in range(10_001)]
worst = max(rel_err(guess, exact_quantile(ordered, q)) for q, guess in zip(grid, sketch.quantiles(grid)))
print(f"\nworst error over {len(grid):,} evenly spaced quantiles: {worst:.4%} (promised: at most {sketch.alpha:.0%})")
python step05_guarantee.py
all requests at 220 ms -> sketch p95 = 220.9 ms
all requests at 320 ms -> sketch p95 = 323.0 ms

values added: 1,000,000, buckets used: 364

 quantile       exact      sketch   error
      0.5    80.44 ms    79.65 ms  0.984%
      0.9   166.62 ms   166.94 ms  0.193%
     0.99   452.75 ms   453.81 ms  0.233%
    0.999  2608.01 ms  2585.66 ms  0.857%
   0.9999  5060.32 ms  5103.90 ms  0.861%
      1.0  13792.80 ms  13874.29 ms  0.591%

worst error over 10,001 evenly spaced quantiles: 0.9999% (promised: at most 1%)

The first two lines answer the thought experiment from Step 3. The classic histogram said 295.0 and 442.5 ms for spikes at 220 and 320 ms. The sketch says 220.9 and 323.0 ms. For comparison, the Prometheus guide says its native histogram with a bucket factor of 1.1 estimates 228 ms for the first spike and 323 ms for the second.

The million-value run shows the real point. The sketch summarized one million requests in 364 buckets. Its p50, p90, p99 and p99.9 are all within 1 percent of the exact values, and the worst error over 10,001 evenly spaced quantiles is 0.9999 percent. That is just under the promise. It is the guarantee at work: the error can come close to 1 percent, but the paper’s proof says it cannot exceed it. The p50 row shows it, with 0.984 percent. Adding values is also fast. On my machine a loop of plain Python add calls ran at about 4 million per second (your number will differ).

Checkpoint: “worst error” is below 1 percent. If it is not, look at _index and _value first: swapping ceil for floor is a classic bug, and Step 10 plants exactly that bug to prove the tests would catch it.

Step 6: Merge sketches from many servers

This is the step that Step 2 asked for. Because the bucket edges depend only on γ and never on the data, two sketches built with the same alpha have identical buckets: bucket 57 means the same range in both. Merging is therefore just adding the counts of equal indexes. The paper puts it this way: any two sketches using the same γ “can be merged by simply summing up the buckets that share an index”. Add the third part of sketch.py:

    # sketch.py (continued)
    def merge(self, other):
        if other.alpha != self.alpha:
            raise ValueError(f"cannot merge sketches with different alpha: {self.alpha} and {other.alpha}")
        for i, c in other.buckets.items():
            self.buckets[i] = self.buckets.get(i, 0) + c
        self.zero_count += other.zero_count
        self.count += other.count
        self._collapse()

The first line of merge refuses sketches with a different alpha. Their buckets would have different edges, so adding counts by index would silently produce nonsense. The rest adds the bucket counts, the zero bucket and the total, and then calls _collapse so that the merged sketch obeys the size limit too.

Create step06_merge.py. It sketches each of the three servers from Step 2 on its own, merges the three sketches, and compares the result with the exact fleet percentiles. It also checks two properties of merging:

"""step06_merge.py: sketch each server on its own, merge the sketches, compare with the exact fleet quantiles."""
import itertools

from labkit import exact_quantile, make_fleet, rel_err
from sketch import DDSketch

fleet = make_fleet()
sketches = {}
for name, times in fleet.items():
    sketches[name] = DDSketch(0.01)
    for t in times:
        sketches[name].add(t)

merged = DDSketch(0.01)
for sketch in sketches.values():
    merged.merge(sketch)

everything = [t for times in fleet.values() for t in times]
ordered = sorted(everything)
for q in (0.5, 0.9, 0.99, 0.999):
    truth, guess = exact_quantile(ordered, q), merged.quantile(q)
    print(f"p{q * 100:g}: exact {truth * 1000:7.1f} ms, merged sketches {guess * 1000:7.1f} ms, {rel_err(guess, truth):.3%} off")

single = DDSketch(0.01)
for t in everything:
    single.add(t)
print("\nmerged buckets equal the buckets of one sketch fed every value:", merged.buckets == single.buckets)

results = set()
orders = list(itertools.permutations(sketches.values()))
for order in orders:
    again = DDSketch(0.01)
    for sketch in order:
        again.merge(sketch)
    results.add(tuple(sorted(again.buckets.items())))
print(f"{len(orders)} merge orders gave {len(results)} distinct result")
python step06_merge.py
p50: exact    81.6 ms, merged sketches    81.3 ms, 0.410% off
p90: exact   173.6 ms, merged sketches   173.8 ms, 0.094% off
p99: exact   770.3 ms, merged sketches   763.3 ms, 0.908% off
p99.9: exact  2982.6 ms, merged sketches  2974.2 ms, 0.279% off

merged buckets equal the buckets of one sketch fed every value: True
6 merge orders gave 1 distinct result

Compare the p99 line with Step 2. The exact fleet p99 is 770.3 ms. The merged sketches say 763.3 ms, which is 0.908 percent off. The plain average of the three p99s was 124.5 percent off, and the request-weighted average was 24.8 percent off. All four percentiles in the table are within 1 percent.

The last two lines are the surprising ones. The merged sketch’s buckets are identical to the buckets of a single sketch fed every value, not merely close. And all six orders in which three sketches can be merged give the same result. In practice that means you can merge sketches in any order, in a tree (each region merges its servers, then the regions merge), or one at a time as they arrive. Each server ships a few hundred bucket counts instead of its raw request times, and the fleet percentile is as accurate as if one machine had seen everything.

Checkpoint: “True” on the identical-buckets line, and “6 merge orders gave 1 distinct result”.

Step 7: Ship sketches as JSON and answer any time window

To send a sketch to a monitoring system, or to store one per minute, you need to turn it into plain data and back. Add the last part of sketch.py:

    # sketch.py (continued)
    def to_dict(self):
        return {"alpha": self.alpha, "zero": self.zero_count,
                "buckets": {str(i): c for i, c in self.buckets.items()}}

    @classmethod
    def from_dict(cls, payload, max_buckets=2048):
        """Rebuild a sketch from JSON data. A sketch from the network is untrusted input, so check it."""
        sketch = cls(payload["alpha"], max_buckets)
        zero = payload["zero"]
        if type(zero) is not int or zero < 0:
            raise ValueError("zero count must be an integer of 0 or more")
        if len(payload["buckets"]) > 100_000:
            raise ValueError("too many buckets")
        for key, c in payload["buckets"].items():
            i = int(key)
            if type(c) is not int or c < 1 or abs(i) > sketch._max_index:
                raise ValueError(f"bad bucket {key!r}: {c!r}")
            sketch.buckets[i] = c
        sketch.zero_count = zero
        sketch.count = zero + sum(sketch.buckets.values())
        sketch._collapse()
        return sketch

to_dict is straightforward. JSON object keys must be strings, so the bucket indexes become strings. from_dict is where the care goes, because a sketch that arrives over a network is input from outside your trust boundary. It checks that the zero count is a non-negative integer, that there are at most 100,000 buckets, and that every count is a positive integer (type(c) is int also rejects True and floats) with a bucket index inside the allowed range. It computes the total from the buckets instead of trusting a separate count field, then applies the size cap.

Each check blocks something real. I tried both on a sketch without validation: a bucket index of 999,999 loads without complaint and then crashes the first quantile query with OverflowError, and a negative count loads silently and produces nonsense quantiles with no error at all. Step 10 turns these cases into tests.

Now use it the way a metrics database would. Create step07_windows.py. It simulates one hour of traffic with a seven-minute incident, stores one JSON sketch per minute, and answers “what was the p99 of the last 15 minutes?” by merging 15 stored sketches:

"""step07_windows.py: store one JSON sketch per minute, then answer 'p99 of the last 15 minutes' by merging."""
import json

from labkit import exact_quantile, make_latencies, rel_err
from sketch import DDSketch

# One hour of traffic, 1,500 requests a minute. Minutes 38 to 44 are an incident: slower and with more slow requests.
minutes = []
for m in range(60):
    incident = 38 <= m <= 44
    minutes.append(make_latencies(1_500, seed=100 + m, median=0.240 if incident else 0.080,
                                  slow_share=0.05 if incident else 0.01))

stored = []   # what a metrics database would keep: one JSON document per minute
for times in minutes:
    sketch = DDSketch(0.01)
    for t in times:
        sketch.add(t)
    stored.append(json.dumps(sketch.to_dict()))

raw_bytes = sum(len(json.dumps(times)) for times in minutes)
sketch_bytes = sum(len(doc) for doc in stored)
print(f"one hour as raw values in JSON: {raw_bytes / 1e6:.2f} MB; as 60 sketches in JSON: {sketch_bytes / 1e6:.2f} MB")

whole_hour = DDSketch(0.01)
for doc in stored:
    whole_hour.merge(DDSketch.from_dict(json.loads(doc)))
one_minute = DDSketch.from_dict(json.loads(stored[0]))
print(f"buckets in one minute's sketch: {len(one_minute.buckets)}; in the merged hour ({whole_hour.count:,} requests): {len(whole_hour.buckets)}")

window = range(31, 46)   # the last 15 minutes ending at minute 45
merged = DDSketch(0.01)
for m in window:
    merged.merge(DDSketch.from_dict(json.loads(stored[m])))

raw = sorted(t for m in window for t in minutes[m])
truth = exact_quantile(raw, 0.99)
per_minute = [exact_quantile(sorted(minutes[m]), 0.99) for m in window]
average = sum(per_minute) / len(per_minute)
print(f"\np99 of the last 15 minutes ({len(raw):,} requests)")
print(f"  exact, from the raw values      {truth * 1000:7.1f} ms")
print(f"  merged sketches                 {merged.quantile(0.99) * 1000:7.1f} ms  ({rel_err(merged.quantile(0.99), truth):.2%} off)")
print(f"  average of 15 per-minute p99s   {average * 1000:7.1f} ms  ({rel_err(average, truth):.2%} off)")
python step07_windows.py
one hour as raw values in JSON: 1.88 MB; as 60 sketches in JSON: 0.11 MB
buckets in one minute's sketch: 159; in the merged hour (90,000 requests): 359

p99 of the last 15 minutes (22,500 requests)
  exact, from the raw values       4033.7 ms
  merged sketches                  4014.8 ms  (0.47% off)
  average of 15 per-minute p99s    3093.1 ms  (23.32% off)

Three results stand out:

  • Size. One hour of raw request times as JSON is 1.88 MB in this toy. Sixty sketches take 0.11 MB. A sketch’s size depends on the spread of the values, not on how many there were: one minute’s sketch has 159 buckets, and the merged hour (90,000 requests) only has 359. So the gap widens as traffic grows.
  • The window answer. The exact p99 over the 22,500 requests of the last 15 minutes is 4,033.7 ms. The merged sketches say 4,014.8 ms, which is 0.47 percent off.
  • The trap again. Averaging the 15 per-minute p99s gives 3,093.1 ms, which is 23.32 percent too low. The seven incident minutes hold far slower requests than the eight quiet ones, yet an average of per-minute p99s gives every minute the same vote, so the quiet minutes drag the answer down.

OpenTelemetry’s data model describes this operation by name. Its list of data transformations includes temporal reaggregation: “Metrics that are collected at a high-frequency can be re-aggregated into longer intervals”. Merging per-minute sketches into a 15-minute sketch is exactly that.

Checkpoint: the merged answer is within 1 percent of the exact one, and the average of per-minute p99s is not.

Step 8: Keep the size bounded

A sketch has one bucket for every 2 percent step of value that it has seen, so its size grows with the spread of the data, not with the number of values. For request times that spread is small. But a bug can make it large: a unit mix-up, or a timer that sometimes returns garbage, can scatter values across many orders of magnitude, and each new 2 percent step adds a bucket. The paper’s remedy, called collapsing, is to cap the number of buckets at m. When an insert or a merge pushes the sketch over the cap, the lowest bucket is folded into its neighbor (Algorithm 3 in the paper). That is what _collapse in Step 4 does: heapq.nsmallest(2, self.buckets) picks the two lowest indexes, and self.buckets[second] += self.buckets.pop(lowest) moves the counts up. A while loop repeats it until the sketch fits.

The price is that quantiles that fall in a folded bucket become wrong, so the sketch gives up accuracy at the low end first. For latency that is the right end to give up, because nobody pages anyone for a request that was too fast. Create step08_collapse.py. It uses a deliberately absurd data set, request times spread evenly in log space from 1 microsecond to 1,000 seconds, and then counts the buckets different accuracy settings need:

"""step08_collapse.py: cap the number of buckets and see which quantiles pay for it."""
import random

from labkit import exact_quantile, rel_err
from sketch import DDSketch

rng = random.Random(8)
values = [10 ** rng.uniform(-6, 3) for _ in range(200_000)]   # 1 microsecond to 1,000 seconds: absurdly wide
ordered = sorted(values)

for cap in (2048, 300):
    sketch = DDSketch(0.01, max_buckets=cap)
    for v in values:
        sketch.add(v)
    print(f"max_buckets={cap}: {len(sketch.buckets)} buckets in use")
    for q in (0.01, 0.5, 0.9, 0.99):
        truth, guess = exact_quantile(ordered, q), sketch.quantile(q)
        within = "yes" if rel_err(guess, truth) <= sketch.alpha else "no"
        print(f"  q={q:<5} exact {truth:>10.4g} s   sketch {guess:>10.4g} s   within 1%? {within}")

probe = DDSketch(0.01)
low, high = probe._index(80e-6), probe._index(365 * 24 * 3600)
print(f"\n80 microseconds to 1 year needs {high - low + 1} buckets at alpha = 1%; the default allows 2048")

print("\nbuckets needed to cover 1 millisecond to 100 seconds:")
for alpha in (0.05, 0.01, 0.001):
    probe = DDSketch(alpha)
    print(f"  alpha = {alpha:<6} {probe._index(100) - probe._index(0.001) + 1:>5} buckets")
python step08_collapse.py
max_buckets=2048: 1037 buckets in use
  q=0.01  exact  1.237e-06 s   sketch  1.228e-06 s   within 1%? yes
  q=0.5   exact    0.02982 s   sketch    0.02989 s   within 1%? yes
  q=0.9   exact      124.9 s   sketch      125.2 s   within 1%? yes
  q=0.99  exact      811.2 s   sketch      804.5 s   within 1%? yes
max_buckets=300: 300 buckets in use
  q=0.01  exact  1.237e-06 s   sketch      2.534 s   within 1%? no
  q=0.5   exact    0.02982 s   sketch      2.534 s   within 1%? no
  q=0.9   exact      124.9 s   sketch      125.2 s   within 1%? yes
  q=0.99  exact      811.2 s   sketch      804.5 s   within 1%? yes

80 microseconds to 1 year needs 1336 buckets at alpha = 1%; the default allows 2048

buckets needed to cover 1 millisecond to 100 seconds:
  alpha = 0.05     117 buckets
  alpha = 0.01     577 buckets
  alpha = 0.001   5757 buckets

With the default cap of 2,048 buckets, the nine decades of values need 1,037 buckets and every quantile is within 1 percent. With the cap set to 300, the p90 and p99 are still within 1 percent, but the p1 and the median both come back as 2.534 seconds, because everything below that was folded into one bucket. The paper gives the condition under which a quantile is safe (its Proposition 4): the largest value must be no more than γm−1 times the quantile’s value. With m = 300 that factor is about 395.5, and 1,000 seconds divided by 395.5 is 2.5 seconds, which matches the output.

The next line checks a claim in the paper. For α = 0.01, it says a sketch of size 2,048 “can handle values from 80 microseconds to 1 year”. The reference package’s README repeats it. The count above shows why it holds: that whole range needs 1,336 buckets, which is under 2,048.

The last table is the one to remember when you choose α. Covering 1 millisecond to 100 seconds takes 117 buckets at α = 0.05, 577 at α = 0.01 and 5,757 at α = 0.001. Ten times more accuracy costs about ten times more buckets, so pick the loosest α that your decisions can tolerate.

Checkpoint: “within 1%? no” appears only for the cap of 300 and only for q = 0.01 and 0.5.

Step 9: Know the traps before they bite

Five mistakes are worth knowing before you rely on a sketch. Create step09_pitfalls.py. It triggers each one on purpose so you can read the error message:

"""step09_pitfalls.py: the mistakes that bite first, and the error messages that catch them."""
import time

from sketch import DDSketch


def show(label, action):
    try:
        action()
        print(f"{label}: no error")
    except ValueError as error:
        print(f"{label}: ValueError: {error}")


show("merge alpha 0.01 with 0.02", lambda: DDSketch(0.01).merge(DDSketch(0.02)))
show("add(nan)", lambda: DDSketch().add(float("nan")))
show("add(-0.003)", lambda: DDSketch().add(-0.003))
show("quantile of an empty sketch", lambda: DDSketch().quantile(0.99))
show("sketch received with a huge index", lambda: DDSketch.from_dict({"alpha": 0.01, "zero": 0, "buckets": {"999999": 5}}))

zeros = DDSketch()
for _ in range(90):
    zeros.add(0.0)
for _ in range(10):
    zeros.add(0.5)
print(f"\n90 zero-length and 10 half-second requests: p50 = {zeros.quantile(0.5)}, p95 = {zeros.quantile(0.95):.4f}")

print()
for name in ("time", "perf_counter"):
    info = time.get_clock_info(name)
    print(f"time.{name}: monotonic = {info.monotonic}, adjustable = {info.adjustable}")
python step09_pitfalls.py
merge alpha 0.01 with 0.02: ValueError: cannot merge sketches with different alpha: 0.01 and 0.02
add(nan): ValueError: cannot add nan: need a finite number that is 0 or more
add(-0.003): ValueError: cannot add -0.003: need a finite number that is 0 or more
quantile of an empty sketch: ValueError: the sketch is empty
sketch received with a huge index: ValueError: bad bucket '999999': 5

90 zero-length and 10 half-second requests: p50 = 0.0, p95 = 0.5015

time.time: monotonic = False, adjustable = True
time.perf_counter: monotonic = True, adjustable = False
  • Merging sketches with different accuracy. The first line shows the refusal. The fix is organizational: choose one alpha as a shared constant for every producer, and let the payload carry it (ours does) so that a stray sketch is caught at merge time rather than blended in.
  • Bad values. NaN and negative numbers are rejected. A negative request time almost always means the timing code used the wrong clock. The Python documentation of time.get_clock_info explains the two flags printed at the bottom of the output: adjustable is “True if the clock can be set to jump forward or backward in time”, and monotonic is “True if the clock cannot go backward”. time.time() is adjustable and not monotonic, so an adjustment between two readings can produce a negative duration. time.perf_counter() is the clock to time requests with, as in the snippet below. Decide what your request path does with a ValueError (count it and carry on, usually) so that a metrics bug never fails a customer’s request.
  • Zero-length requests. A coarse clock reports a duration of exactly 0.0 for a very fast request. The zero bucket handles them: 90 zeros and 10 half-second requests give a median of 0.0 and a p95 of 0.5015, which is 0.3 percent above the true 0.5.
  • Untrusted sketches. The hostile bucket index of 999,999 is refused at load time, as designed in Step 7.
  • Mixed units. The code cannot see this one, so it needs a rule. Use seconds everywhere. Prometheus says the same in its guide (“you should always use base units in Prometheus”). A sketch fed milliseconds still works on its own, but merging it with a sketch of seconds produces numbers that look plausible and mean nothing.
# example only (not a lab file): time a request and record it
import time

start = time.perf_counter()
handle_request()
try:
    sketch.add(time.perf_counter() - start)
except ValueError:
    errors_total += 1   # count the problem, never fail the request because of metrics

One more trap is not in the script because the fix is a single line: add updates a dictionary in several separate steps, so if several threads call it on the same sketch, guard it with a threading.Lock.

Step 10: Test it, then break it on purpose

You have seen the guarantee hold on one data set. A test suite should prove it on many, and fail loudly when someone changes the code wrongly. Create test_sketch.py:

"""test_sketch.py: pytest checks for sketch.py. Run with: python -m pytest -q test_sketch.py"""
import json
import math
import random

import pytest

from labkit import exact_quantile, make_latencies, rel_err
from sketch import DDSketch

QS = [i / 1000 for i in range(1001)] + [0.9999, 0.99999]
BIG = 100_000   # a bucket cap high enough that nothing collapses in these tests


def build(values, alpha=0.01):
    sketch = DDSketch(alpha, max_buckets=BIG)
    for v in values:
        sketch.add(v)
    return sketch


@pytest.mark.parametrize("alpha", [0.05, 0.01, 0.001])
@pytest.mark.parametrize("seed", [1, 2, 3])
def test_every_quantile_is_within_alpha(seed, alpha):
    values = make_latencies(20_000, seed)
    ordered = sorted(values)
    for q, guess in zip(QS, build(values, alpha).quantiles(QS)):
        assert rel_err(guess, exact_quantile(ordered, q)) <= alpha * (1 + 1e-9)


def test_merge_equals_one_big_sketch_in_any_order():
    values = make_latencies(30_000, 4)
    whole = build(values)
    parts = [build(values[i::5]) for i in range(5)]
    random.Random(0).shuffle(parts)
    merged = DDSketch(0.01, max_buckets=BIG)
    for part in parts:
        merged.merge(part)
    assert merged.buckets == whole.buckets
    assert merged.count == whole.count


def test_merging_different_alphas_is_refused():
    with pytest.raises(ValueError, match="different alpha"):
        DDSketch(0.01).merge(DDSketch(0.02))


@pytest.mark.parametrize("bad", [math.nan, math.inf, -math.inf, -1.0, 1e306])
def test_bad_values_are_rejected(bad):
    with pytest.raises(ValueError):
        DDSketch().add(bad)


def test_zero_bucket_and_extremes():
    sketch = DDSketch()
    for _ in range(10):
        sketch.add(0.0)
    sketch.add(2.0)
    assert sketch.quantile(0.5) == 0.0
    assert rel_err(sketch.quantile(1.0), 2.0) <= 0.01


def test_merge_carries_the_zero_bucket_across():
    a, b = DDSketch(), DDSketch()
    a.add(1.0)
    for _ in range(5):
        b.add(0.0)
    a.merge(b)
    assert a.zero_count == 5 and a.count == 6
    assert a.quantile(0.0) == 0.0


def test_empty_sketch_and_bad_q():
    with pytest.raises(ValueError, match="empty"):
        DDSketch().quantile(0.5)
    sketch = build([1.0])
    with pytest.raises(ValueError):
        sketch.quantile(1.5)


def test_collapse_keeps_the_high_quantiles_accurate():
    rng = random.Random(5)
    values = [10 ** rng.uniform(-6, 3) for _ in range(20_000)]
    ordered = sorted(values)
    sketch = DDSketch(0.01, max_buckets=300)
    for v in values:
        sketch.add(v)
    assert len(sketch.buckets) <= 300
    assert sketch.count == len(values)
    for q in (0.9, 0.99, 0.999):
        assert rel_err(sketch.quantile(q), exact_quantile(ordered, q)) <= 0.01 * (1 + 1e-9)
    assert rel_err(sketch.quantile(0.01), exact_quantile(ordered, 0.01)) > 0.01   # the price of collapsing


def test_merge_respects_the_bucket_cap():
    a, b = DDSketch(0.01, max_buckets=50), DDSketch(0.01, max_buckets=50)
    for i in range(1, 50):
        a.add(1.02 ** i)
        b.add(1.02 ** (i + 50))
    a.merge(b)
    assert len(a.buckets) <= 50 and a.count == 98


def test_json_round_trip():
    sketch = build(make_latencies(5_000, 6))
    again = DDSketch.from_dict(json.loads(json.dumps(sketch.to_dict())), max_buckets=BIG)
    assert again.buckets == sketch.buckets
    assert again.count == sketch.count
    assert again.quantile(0.99) == sketch.quantile(0.99)


@pytest.mark.parametrize("payload", [
    {"alpha": 2, "zero": 0, "buckets": {}},
    {"alpha": 0.01, "zero": -1, "buckets": {}},
    {"alpha": 0.01, "zero": 0, "buckets": {"5": 0}},
    {"alpha": 0.01, "zero": 0, "buckets": {"999999": 1}},
    {"alpha": 0.01, "zero": 0, "buckets": {"1": 1.5}},
    {"alpha": 0.01, "zero": 0, "buckets": {"1": True}},
])
def test_hostile_payloads_are_refused(payload):
    with pytest.raises(ValueError):
        DDSketch.from_dict(payload)
python -m pytest -q test_sketch.py
............................ [100%]
28 passed in 0.22s

The most important test is test_every_quantile_is_within_alpha. It runs three accuracy settings (5 percent, 1 percent and 0.1 percent) on three seeds, and on each one compares 1,003 quantiles with the exact answer. A tiny tolerance of one part in a billion lets floating-point rounding through. Next in importance are the merge test (any order gives the same buckets as one big sketch) and the hostile-payload test, which turns the cases from Step 7 into permanent checks. test_collapse_keeps_the_high_quantiles_accurate even asserts that the low quantile is wrong, to document the price of collapsing.

Passing tests prove little unless they can also fail. Mutation checking tests the tests: plant one bug at a time in a copy of the code and see whether the suite notices. Create mutations.py:

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

HERE = Path(__file__).parent
MUTANTS = [
    ("floor instead of ceil in the bucket index", "return math.ceil(math.log(x) / self._ln_gamma)",
     "return math.floor(math.log(x) / self._ln_gamma)"),
    ("answer with the bucket's upper edge", "return 2 * self.gamma ** i / (self.gamma + 1)",
     "return self.gamma ** i"),
    ("rank off by one in the quantile walk", "bisect.bisect_right(cumulative, rank)",
     "bisect.bisect_left(cumulative, rank)"),
    ("merge forgets the zero bucket", "        self.zero_count += other.zero_count\n", ""),
    ("merge skips the alpha check", "if other.alpha != self.alpha:", "if False:"),
    ("collapse throws counts away", "self.buckets[second] += self.buckets.pop(lowest)",
     "self.buckets.pop(lowest)"),
    ("from_dict trusts the counts", "if type(c) is not int or c < 1 or abs(i) > sketch._max_index:",
     "if False:"),
]

source = (HERE / "sketch.py").read_text(encoding="utf-8")
caught = 0
for label, old, new in MUTANTS:
    assert source.count(old) == 1, f"pattern for {label!r} must match once"
    with tempfile.TemporaryDirectory() as tmp:
        for name in ("labkit.py", "test_sketch.py"):
            shutil.copy(HERE / name, tmp)
        (Path(tmp) / "sketch.py").write_text(source.replace(old, new), encoding="utf-8", newline="\n")
        done = subprocess.run([sys.executable, "-B", "-m", "pytest", "-q", "-x", "-p", "no:cacheprovider", "test_sketch.py"],
                              cwd=tmp, capture_output=True, text=True, encoding="utf-8")
    status = "caught" if done.returncode != 0 else "SURVIVED"
    caught += done.returncode != 0
    print(f"{status:<9} {label}")
print(f"\n{caught} of {len(MUTANTS)} planted bugs were caught")
python mutations.py
caught    floor instead of ceil in the bucket index
caught    answer with the bucket's upper edge
caught    rank off by one in the quantile walk
caught    merge forgets the zero bucket
caught    merge skips the alpha check
caught    collapse throws counts away
caught    from_dict trusts the counts

7 of 7 planted bugs were caught

All seven planted bugs are caught. It did not start that way. On my first run, six were caught and “merge forgets the zero bucket” survived, which meant the original tests never merged a sketch that contained zeros. That was a real hole, so I added test_merge_carries_the_zero_bucket_across, and the check now reports 7 of 7. This is the habit worth copying: a surviving mutant is a bug report about your tests.

Step 11: Check your sketch against the reference package

Your own implementation agrees with the paper on every test you wrote, but those tests were also written by you. An independent check is worth more. Datadog publishes the reference implementation as the ddsketch package (source: DataDog/sketches-py). Its README states the same guarantee: if the true q-quantile is x, the package returns a y such that |x − y| / x is below the relative-error parameter. The README also says that its collapsing variants default to 2,048 buckets, which cover 80 microseconds to 1 year when the data is time in seconds. Create step11_reference.py. It feeds the same one million values to both implementations and compares 10,001 quantiles:

"""step11_reference.py: run our sketch and Datadog's ddsketch package side by side on the same values."""
import importlib.metadata

from ddsketch import DDSketch as Reference

from labkit import make_latencies, rel_err
from sketch import DDSketch

times = make_latencies(1_000_000, seed=5)
ours, reference = DDSketch(0.01), Reference(relative_accuracy=0.01)
for t in times:
    ours.add(t)
    reference.add(t)

print("ddsketch package", importlib.metadata.version("ddsketch"))
grid = [i / 10_000 for i in range(10_001)]
differences = [rel_err(a, reference.get_quantile_value(q)) for q, a in zip(grid, ours.quantiles(grid))]
print(f"quantiles compared: {len(grid):,}")
print(f"largest relative difference between the two answers: {max(differences):.1e}")
print(f"quantiles where the answers differ by more than 1e-9: {sum(1 for d in differences if d > 1e-9)}")
python step11_reference.py
ddsketch package 3.0.1
quantiles compared: 10,001
largest relative difference between the two answers: 2.4e-14
quantiles where the answers differ by more than 1e-9: 0

The two answers agree to about 2.4 parts in a hundred trillion at worst, which is floating-point rounding (the package computes the bucket’s answer with slightly different arithmetic). No quantile differs by more than one part in a billion. Same buckets, same walk, same answers: your 100 lines do what the reference package does for these inputs. The package is the better choice in production: it is maintained by the authors’ company, it handles negative values, and it offers collapsing variants with a fixed bucket budget.

Step 12: Recognize the same idea in OpenTelemetry and Prometheus

You do not always have to build this yourself, and it helps to know where the same idea already lives. OpenTelemetry’s exponential histogram is a bucketed distribution with log-spaced buckets, close kin to a DDSketch. The data model specification describes it like this:

  • The resolution is set by an integer scale, and bucket edges sit at powers of a base: base = 2**(2**(-scale)). Scale 3 gives a base of 1.09051, scale 5 gives 1.02190 and scale 6 gives 1.01089.
  • The bucket with a given index “represents values in the population that are greater than base**index and less than or equal to base**(index+1)”, which is exactly the bucket rule from Step 4 with γ equal to the base.
  • Every base is a power-of-two root, and the base at one scale is the square of the base one scale down, so the bucket edges at a lower scale are a subset of the edges at a higher scale. The spec calls this perfect subsetting: buckets at one scale map exactly into buckets at a lesser scale, “which allows consumers to lower the resolution of a histogram (i.e., downscale) without introducing error”.

The Python SDK can produce these histograms, and the metrics SDK specification says its default of 160 buckets was chosen “to cover a long-tail latency distribution from 1ms to 100s with less than 5% relative error”. Create step12_otel.py. It records request times into an SDK exponential histogram, reads the data point back, and estimates percentiles from its buckets with the same walk as DDSketch.quantiles. It also tests the subsetting claim with plain arithmetic:

"""step12_otel.py: read percentiles back out of an OpenTelemetry exponential histogram (opentelemetry-sdk)."""
import math
import random

from opentelemetry.sdk.metrics import MeterProvider
from opentelemetry.sdk.metrics.export import InMemoryMetricReader
from opentelemetry.sdk.metrics.view import ExponentialBucketHistogramAggregation, View

from labkit import exact_quantile, make_latencies, rel_err


def record(values, max_size=160):
    reader = InMemoryMetricReader()
    view = View(instrument_name="latency", aggregation=ExponentialBucketHistogramAggregation(max_size=max_size))
    meter = MeterProvider(metric_readers=[reader], views=[view]).get_meter("lab")
    histogram = meter.create_histogram("latency", unit="s")
    for v in values:
        histogram.record(v)
    return reader.get_metrics_data().resource_metrics[0].scope_metrics[0].metrics[0].data.data_points[0]


def quantile_from_buckets(point, q):
    """Same walk as DDSketch.quantile: bucket `index` holds (base**index, base**(index + 1)]."""
    base = 2 ** (2 ** -point.scale)
    rank = q * (point.count - 1)
    seen = point.zero_count
    if seen > rank:
        return 0.0
    for k, c in enumerate(point.positive.bucket_counts):
        seen += c
        if seen > rank:
            low = base ** (point.positive.offset + k)
            high = low * base
            return 2 * low * high / (low + high)   # the harmonic middle, the same answer as 2 * gamma**i / (gamma + 1)


def report(label, values, max_size=160):
    point = record(values, max_size)
    base = 2 ** (2 ** -point.scale)
    ordered = sorted(values)
    used = sum(1 for c in point.positive.bucket_counts if c)
    print(f"{label}: scale {point.scale}, base {base:.5f}, bound {(base - 1) / (base + 1):.2%}, {used} buckets in use")
    for q in (0.5, 0.9, 0.99, 0.999):
        truth, guess = exact_quantile(ordered, q), quantile_from_buckets(point, q)
        print(f"  q={q:<5} exact {truth * 1000:9.2f} ms   from buckets {guess * 1000:9.2f} ms   {rel_err(guess, truth):6.2%} off")


report("100,000 requests, default 160 buckets", make_latencies(100_000, seed=1))
rng = random.Random(8)
wide = [10 ** rng.uniform(-3, 1) for _ in range(100_000)]   # 1 millisecond to 10 seconds, evenly spread in log space
report("\nwide range, default 160 buckets", wide)
report("\nwide range, 1,000 buckets allowed", wide, max_size=1000)


def index(x, scale):
    return math.ceil(math.log2(x) * 2 ** scale) - 1


xs = [10 ** rng.uniform(-6, 3) for _ in range(100_000)]
same = sum(1 for x in xs if index(x, 3) == index(x, 4) >> 1)
print(f"\nscale 4 index shifted right by one equals the scale 3 index for {same:,} of {len(xs):,} values")
python step12_otel.py
100,000 requests, default 160 buckets: scale 3, base 1.09051, bound 4.33%, 80 buckets in use
  q=0.5   exact     80.60 ms   from buckets     77.54 ms    3.79% off
  q=0.9   exact    166.28 ms   from buckets    169.12 ms    1.71% off
  q=0.99  exact    467.71 ms   from buckets    478.35 ms    2.28% off
  q=0.999 exact   2696.35 ms   from buckets   2705.97 ms    0.36% off

wide range, default 160 buckets: scale 3, base 1.09051, bound 4.33%, 107 buckets in use
  q=0.5   exact     97.35 ms   from buckets    100.56 ms    3.29% off
  q=0.9   exact   3965.95 ms   from buckets   3826.82 ms    3.51% off
  q=0.99  exact   9118.16 ms   from buckets   9101.77 ms    0.18% off
  q=0.999 exact   9882.63 ms   from buckets   9925.55 ms    0.43% off

wide range, 1,000 buckets allowed: scale 6, base 1.01089, bound 0.54%, 851 buckets in use
  q=0.5   exact     97.35 ms   from buckets     96.91 ms    0.46% off
  q=0.9   exact   3965.95 ms   from buckets   3978.34 ms    0.31% off
  q=0.99  exact   9118.16 ms   from buckets   9159.64 ms    0.45% off
  q=0.999 exact   9882.63 ms   from buckets   9881.06 ms    0.02% off

scale 4 index shifted right by one equals the scale 3 index for 100,000 of 100,000 values

What the output shows:

  • With the default of 160 buckets the SDK chose scale 3, a base of 1.09051. A DDSketch with the same buckets would promise (base − 1) / (base + 1) = 4.33 percent, which is under the spec’s 5 percent. The measured errors, at most 3.79 percent on the first data set and 3.51 percent on the wide one, stay inside that bound.
  • Accuracy is a trade between range and bucket count, as in Step 8. The SDK keeps the bucket count fixed and lowers the scale when the data spreads out. Allow 1,000 buckets for the wide data set and it picks scale 6, with a bound of 0.54 percent and measured errors below 0.5 percent.
  • The last line confirms perfect subsetting: shifting a scale 4 bucket index right by one gives the scale 3 index for all 100,000 values. In other words, lowering the scale by one is just adding aligned pairs of neighboring buckets.

That last point is the difference in design. DDSketch lets you choose any α, so two sketches with different α have unrelated bucket edges and cannot be merged exactly. OpenTelemetry fixes the base as a power-of-two root, so two histograms of different scales can be merged: lower the finer one to the coarser scale, then add counts. Prometheus does the same. Its histogram guide says that with the standard exponential schema, reconciling different bucket layouts works “essentially by falling back to the lowest common resolution among all involved histograms”, and that OpenTelemetry’s exponential histograms “are always converted into the usual native histograms (with standard exponential buckets)”. If your tooling already speaks exponential histograms, use them. Build your own sketch when you control both ends, such as in a batch job, a test harness or a custom pipeline, or when you want to understand exactly what a dashboard’s percentile can and cannot tell you.

Run everything and confirm it works

One script runs all the steps, the tests and the mutation check, and saves each output to an out folder. Create run_all.py:

"""run_all.py: run every step, the tests and the mutation check; save each output to out/NAME.txt."""
import subprocess
import sys
from pathlib import Path

HERE = Path(__file__).parent
OUT = HERE / "out"
OUT.mkdir(exist_ok=True)
STEPS = ["step01_baseline", "step02_average_of_p99", "step03_fixed_buckets", "step04_buckets", "step05_guarantee",
         "step06_merge", "step07_windows", "step08_collapse", "step09_pitfalls", "step11_reference", "step12_otel"]


def run(name, args):
    done = subprocess.run([sys.executable, "-B", *args], cwd=HERE, capture_output=True, text=True, encoding="utf-8")
    text = (done.stdout + done.stderr).rstrip("\n") + "\n"
    with open(OUT / f"{name}.txt", "w", newline="\n", encoding="utf-8") as f:
        f.write(text)
    print(f"== {name} (exit {done.returncode})", flush=True)
    print(text, flush=True)


for step in STEPS:
    run(step, [f"{step}.py"])
run("pytest", ["-m", "pytest", "-q", "-p", "no:cacheprovider", "test_sketch.py"])
run("mutations", ["mutations.py"])
python run_all.py

On my machine it finished in about ten seconds. Confirm these lines in the output, which together prove the whole tutorial end to end:

  • Step 2: fleet p99 of 770.3 ms, with the average of the p99s 124.5 percent off.
  • Step 5: “worst error over 10,001 evenly spaced quantiles: 0.9999% (promised: at most 1%)”.
  • Step 6: “True” for identical buckets and “6 merge orders gave 1 distinct result”.
  • Step 7: merged sketches 0.47 percent off on the 15-minute window, against 23.32 percent for averaged p99s.
  • Step 10: “28 passed” and “7 of 7 planted bugs were caught”.
  • Step 11: no quantile differing from the reference package by more than 1e-9.

Mistakes to avoid

  • Averaging percentiles across hosts, regions or minutes. Ship sketches or bucket counts and merge them.
  • Choosing bucket edges in advance for a fixed-bucket histogram and trusting the percentile it prints. The error is up to the width of the bucket the percentile falls in.
  • Merging sketches with different alpha or different units. Fix both as shared constants.
  • Setting alpha tighter than you need. Ten times more accuracy costs about ten times more buckets.
  • Timing requests with time.time(). Use time.perf_counter().
  • Loading a sketch without validating it. A payload from the network can carry a bucket index that crashes your queries or a count that quietly corrupts them.
  • Reading “within 1 percent” as “1 percent of requests”. The promise is about the value returned for the percentile you asked about, not about how many requests were slower than it.

Where to go next

To use this in a service, wrap request handling in a timer like the snippet in Step 9, keep one sketch per time slice, rotate it each minute, and ship the old one as JSON to wherever you store metrics. If you already run OpenTelemetry, the FastAPI instrumentation tutorial and the Prometheus histograms to OTLP tutorial on this site show the pipeline side, and Step 12 above shows what its exponential histograms carry.

For measuring a service under overload, pair this with the backpressure and load shedding tutorial and the LLM latency benchmark tutorial, both of which compute percentiles from a list of raw request times, which is the part a mergeable sketch would replace. DDSketch also belongs to a family of small summaries that this site has covered one by one: the Bloom filter, the HyperLogLog counter, the Count-Min sketch, MinHash and the cuckoo filter. Each one answers a different question about a stream (membership, distinct count, frequency, similarity, and here, percentiles) in memory that does not grow with the stream. If you want to go deeper, read the paper’s Section 3, which proves how few buckets heavy-tailed data needs, and try extending DDSketch with a second sketch for negative numbers, as Section 2.2 of the paper describes.

Sources

Tags:

DDSketchObservabilityOpenTelemetryPercentilesPythonStreaming Algorithms

Share

No Comment! Be the first one.

Leave a Reply

Your email address will not be published. Required fields are marked *

Related Posts