How to Calculate Descriptive Statistics in Python to See Past a Misleading Average

Imagine you are on call and a dashboard tells you the checkout API’s average response time last hour was 362.6 milliseconds. Is that good or bad? You cannot actually tell from that one number, and by the end of this tutorial you will understand exactly why. You are going to build a small, fully tested Python toolkit for describing a dataset the way engineers actually need it described: not just “the average,” but the typical value, the spread, the outliers, and the relationships hiding inside it. Every number in this article comes from code you can run yourself, not a hand-picked example.

This is a beginner-friendly tutorial. You do not need a statistics background. Every term (mean, variance, percentile, correlation) is defined in plain language the first time it appears, and every code block was actually executed to produce the output shown underneath it.

What You Will Build

You will generate a small, realistic dataset of API response times, then use it to learn and verify each of the following:

• Why the mean (the everyday “average”) can be dragged far away from what a typical request actually experienced, and why the median often tells a truer story.
• The difference between population variance and sample variance, a distinction that trips up experienced engineers, not just beginners.
• How percentiles like p95 and p99 are actually calculated, and why two “correct” methods can disagree by over 100 milliseconds on the same data.
• A simple, automatic way to flag outliers using the interquartile range (IQR), plus a real gotcha in how that method behaves.
• How to measure whether two things are related (like payload size and latency) with correlation, and how to prove correlation is not causation using your own data.
• How to package all of this into one reusable, pytest-verified describe() function.

Prerequisites

• Python 3.10 or later. This tutorial was written and verified against Python 3.13. Everything in the main tutorial uses only Python’s standard library statistics module, so you do not need to install anything to follow along. (Python 3.10 is the minimum because the correlation(), covariance(), and linear_regression() functions used in Step 6 were added in that release.)
• NumPy is used briefly in Steps 3 and 4 to cross-check two results against a second, independent implementation (pip install numpy). This is optional and skippable; the core lesson does not depend on it.
• pytest for the verification section (pip install pytest).
• Basic familiarity with Python: functions, lists, and importing your own modules. No prior statistics knowledge is assumed; every concept is explained before it is used.

Understanding Descriptive Statistics in Plain Language

Descriptive statistics just means describing a dataset with a small number of representative values instead of forcing someone to read every raw data point. It splits into three families, and this tutorial covers one section on each:

• Central tendency: what is a “typical” value in this data? (mean, median, mode)
• Variability: how spread out is the data around that typical value? (variance, standard deviation, percentiles, outliers)
• Correlation: how do two different measurements move together? (covariance, correlation, linear regression)

Python’s standard library ships a module called statistics (added in Python 3.4) specifically for this. Its own documentation is upfront that it is “aimed at the level of graphing and scientific calculators,” not a replacement for NumPy or pandas on large datasets. That is exactly the right tool for this tutorial: everything here works on ordinary Python lists, with no extra dependencies, and the same ideas apply whether you have 60 data points or 60 million.

Step 1: Build a Dataset Worth Describing

Statistics is easiest to learn on data you understand, so instead of an abstract list of numbers, you will simulate response times (in milliseconds) from a fictional checkout-api endpoint. Most requests hit a fast cache-hit path. A handful hit a slower cache-miss path that has to query a database. Two requests are genuine timeouts, the kind an upstream payment processor occasionally produces. Create a file called dataset.py:

"""Shared dataset module used by every step in this tutorial."""
import random

random.seed(42)

_fast = [round(random.gauss(70, 15), 1) for _ in range(50)]
_fast = [max(20.0, v) for v in _fast]
_slow = [round(random.gauss(950, 120), 1) for _ in range(8)]
_timeouts = [5210.4, 5202.9]

LATENCY_MS = _fast + _slow + _timeouts
random.shuffle(LATENCY_MS)

# A second dataset pairing request payload size (KB) with latency (ms),
# used in the correlation step. Larger payloads take a bit longer to
# parse and transmit, plus normal jitter.
random.seed(7)
PAYLOAD_KB = [round(random.uniform(1, 40), 1) for _ in range(50)]
LATENCY_FOR_PAYLOAD_MS = [
    round(35 + 2.1 * kb + random.gauss(0, 12), 1) for kb in PAYLOAD_KB
]

The random.seed(42) call makes this deterministic: every time you run this file (or anything that imports it), you get the exact same 60 numbers, so your output will match this tutorial’s exactly. Run it to confirm:

from dataset import LATENCY_MS
print(f"n = {len(LATENCY_MS)}")
print(LATENCY_MS)

n = 60
[73.7, 67.8, 47.5, 48.2, 55.0, 68.4, 60.6, 56.3, 77.2, 75.0, 65.3, 5202.9, 75.6, 81.7, 1120.6, 67.4, 79.7, 51.7, 60.9, 71.7, 66.0, 87.5, 79.8, 818.1, 82.7, 70.6, 78.5, 30.7, 1230.3, 1009.9, 77.4, 54.8, 66.7, 73.7, 71.7, 5210.4, 851.7, 47.3, 65.6, 75.2, 80.0, 71.7, 68.1, 1042.2, 66.4, 893.6, 68.3, 920.0, 80.5, 59.2, 78.0, 82.3, 58.9, 83.1, 79.6, 68.7, 89.7, 83.1, 53.3, 73.5]

60 latency values, most clustered in the 50 to 90 millisecond range, with a visible smattering of much larger numbers (818, 1120, 5202, 5210) hidden in the middle. Keep this dataset in mind: it represents one hour of real-looking traffic, and everything below describes it.

Step 2: Mean vs Median: Why One Slow Request Can Lie to You

The mean

The mean is what most people call “the average”: add everything up, divide by how many values there are. Python’s statistics module gives you three ways to get it:

import statistics
from dataset import LATENCY_MS

n = len(LATENCY_MS)
naive_mean = sum(LATENCY_MS) / n
lib_mean = statistics.mean(LATENCY_MS)
lib_fmean = statistics.fmean(LATENCY_MS)

print(f"n = {n}")
print(f"naive mean (sum/len)   = {naive_mean:.2f} ms")
print(f"statistics.mean()      = {lib_mean:.2f} ms")
print(f"statistics.fmean()     = {lib_fmean:.2f} ms")

n = 60
naive mean (sum/len)   = 362.60 ms
statistics.mean()      = 362.60 ms
statistics.fmean()     = 362.60 ms

All three agree, which makes sense since they compute the same thing three ways. (Use statistics.fmean() when you specifically want a fast, float-only mean; use statistics.mean() when your data might be integers, fractions, or Decimal values you want preserved exactly.) But look at that number again: 362.60 ms. Does that sound like a typical request to you, given that the raw list above is mostly numbers in the 50s, 60s, 70s, and 80s?

The median

The median is the middle value once the data is sorted (or the average of the two middle values, if there is an even count). Add one line:

lib_median = statistics.median(LATENCY_MS)
print(f"statistics.median()    = {lib_median:.2f} ms")

statistics.median()    = 73.70 ms

73.70 ms. That is nearly 5 times smaller than the mean, and it matches what the raw list actually looks like. The official Python documentation for mean() puts this plainly: “the mean is strongly affected by outliers and is not necessarily a typical example of the data points. For a more robust… measure of central tendency, see median().”

Now prove it. Drop just the 2 timeout requests (3.3% of the traffic) and recompute both:

without_timeouts = [v for v in LATENCY_MS if v < 3000]
mean_without = statistics.mean(without_timeouts)
median_without = statistics.median(without_timeouts)

print(f"n without timeouts = {len(without_timeouts)}")
print(f"mean without timeouts   = {mean_without:.2f} ms "
      f"(dropped {lib_mean - mean_without:.2f} ms)")
print(f"median without timeouts = {median_without:.2f} ms "
      f"(dropped {lib_median - median_without:.2f} ms)")

n without timeouts = 58
mean without timeouts   = 195.56 ms (dropped 167.04 ms)
median without timeouts = 73.60 ms (dropped 0.10 ms)

Removing 2 out of 60 requests dropped the mean by 167 ms but moved the median by one tenth of a millisecond. This is the single most important lesson in this tutorial: a mean can be dominated by a tiny fraction of extreme values, while a median mostly ignores them. If your monitoring dashboard only shows "average latency," a handful of slow requests can make a perfectly healthy service look terrible, or (just as dangerous) a genuinely degraded service can hide behind a mean if most requests are still fast.

The mode

The mode is simply the most frequently occurring value. It matters more for categorical data than for continuous measurements like latency. Try it on a window of HTTP status codes from the same traffic:

status_codes = [200] * 41 + [404] * 3 + [500] * 2 + [200] * 4 + [503] * 1 + [404, 404]

print(f"status codes sample: {status_codes[:10]}... (n={len(status_codes)})")
print(f"statistics.mode()      = {statistics.mode(status_codes)}")
print(f"statistics.multimode() = {statistics.multimode(status_codes)}")

tied = [200, 200, 404, 404, 500]
print(f"tied dataset: {tied}")
print(f"statistics.mode(tied)      = {statistics.mode(tied)}")
print(f"statistics.multimode(tied) = {statistics.multimode(tied)}")

status codes sample: [200, 200, 200, 200, 200, 200, 200, 200, 200, 200]... (n=53)
statistics.mode()      = 200
statistics.multimode() = [200]
tied dataset: [200, 200, 404, 404, 500]
statistics.mode(tied)      = 200
statistics.multimode(tied) = [200, 404]

A gotcha worth knowing: ties

Before Python 3.8, calling mode() on data with more than one most-common value raised an error. That changed: the official documentation notes it was "changed in version 3.8: Now handles multimodal datasets by returning the first mode encountered. Formerly, it raised StatisticsError when more than one mode was found." You can see that above: on the tied dataset, mode() quietly returns 200 (the first one encountered) with no warning that 404 was tied with it. If ties matter for your use case, use multimode() instead, which returns every tied value as a list.

Step 3: Variance and Standard Deviation: Measuring the Spread

Mean and median tell you where the "center" of your data is. Variance and standard deviation tell you how spread out the data is around that center. A small variance means most values sit close to the mean; a large variance means they are scattered widely.

Population variance vs sample variance

This is where Python's statistics module has a gotcha that catches experienced engineers, not just beginners: there are two different variance functions, and they give you two different numbers on the exact same data.

without_timeouts = [v for v in LATENCY_MS if v < 3000]

pop_var = statistics.pvariance(without_timeouts)
sample_var = statistics.variance(without_timeouts)
pop_std = statistics.pstdev(without_timeouts)
sample_std = statistics.stdev(without_timeouts)

n = len(without_timeouts)
print(f"n = {n}")
print(f"population variance (pvariance, divides by N)   = {pop_var:.2f}")
print(f"sample variance     (variance,  divides by N-1) = {sample_var:.2f}")
print(f"population std dev  (pstdev)                     = {pop_std:.2f} ms")
print(f"sample std dev      (stdev)                      = {sample_std:.2f} ms")

n = 58
population variance (pvariance, divides by N)   = 102463.64
sample variance     (variance,  divides by N-1) = 104261.25
population std dev  (pstdev)                     = 320.10 ms
sample std dev      (stdev)                      = 322.90 ms

Same 58 numbers, two different variances. Here is why. Python's statistics module documentation for variance() says: "return the sample variance of data... use this function when your data is a sample from a population." pvariance() says: "return the population variance of data." The difference is one word, but it changes the math: sample variance divides by N - 1 (called Bessel's correction) instead of N, which slightly inflates the result to correct for the fact that a sample tends to underestimate how spread out the true population really is.

The rule of thumb: if your data is the entire thing you care about (these literally are the only 58 requests that happened during this exact hour), use pvariance(). If your data is a sample meant to represent something bigger (these 58 requests are meant to represent typical traffic all month), use variance(). Getting this backwards will not crash your program; it will just quietly give you a number that is close, but not exactly the number you meant to compute.

The NumPy gotcha: default ddof

If you also work with NumPy, there is a second version of this same trap waiting for you. NumPy's var() function has a parameter called ddof ("Delta Degrees of Freedom"). Its official documentation states plainly: "the divisor used in the calculation is N - ddof, where N represents the number of elements. By default ddof is zero." A default of ddof=0 means the divisor is N - 0 = N, which is population variance, not sample variance. Verify it directly against the two stdlib results above:

import numpy as np

arr = np.array(without_timeouts)
np_var_default = np.var(arr)          # ddof=0 -> population
np_var_sample = np.var(arr, ddof=1)   # ddof=1 -> sample (Bessel's correction)

print(f"numpy np.var(arr)          (ddof=0, population) = {np_var_default:.2f}")
print(f"numpy np.var(arr, ddof=1)  (sample)              = {np_var_sample:.2f}")
print(f"matches statistics.pvariance()? "
      f"{abs(np_var_default - pop_var) < 0.005}")
print(f"matches statistics.variance()?  "
      f"{abs(np_var_sample - sample_var) < 0.005}")
print(f"does default np.var() match statistics.variance() (sample)? "
      f"{abs(np_var_default - sample_var) < 0.005}")

numpy np.var(arr)          (ddof=0, population) = 102463.64
numpy np.var(arr, ddof=1)  (sample)              = 104261.25
matches statistics.pvariance()? True
matches statistics.variance()?  True
does default np.var() match statistics.variance() (sample)? False

If you assume np.var(arr) gives you the same "sample variance" that statistics.variance() gives you by default, it does not. You need ddof=1 for that. This is a genuinely common source of quiet, small, hard-to-notice discrepancies when a team mixes stdlib and NumPy code, or compares numbers computed in Python against numbers computed in a spreadsheet (Excel's VAR.S is sample variance, matching ddof=1; VAR.P is population variance, matching the NumPy default).

Step 4: Percentiles: The Metric Engineers Actually Care About

A percentile answers a question like "what value did 95% of requests fall at or below?" That single number, usually written p95 or p99, is what most real latency dashboards report, precisely because it captures the tail behavior that a mean or even a median can hide.

A naive percentile function, and where it breaks

It is tempting to write your own percentile function by sorting the data and indexing into it:

without_timeouts = sorted(v for v in LATENCY_MS if v < 3000)
n = len(without_timeouts)


def naive_percentile(data_sorted, p):
    """A common first attempt: nearest-rank percentile, no interpolation."""
    idx = int(p / 100 * len(data_sorted))
    return data_sorted[idx]


print(f"n = {n}")
for p in (50, 90, 95, 99):
    print(f"  p{p} = {naive_percentile(without_timeouts, p)}")

print("Now try p100 (the max value):")
try:
    naive_percentile(without_timeouts, 100)
except IndexError as exc:
    print(f"  IndexError: {exc!r}")
    print(f"  Why: int(100/100 * {n}) = {int(100 / 100 * n)}, "
          f"but valid indices only go up to {n - 1}.")

n = 58
  p50 = 73.7
  p90 = 893.6
  p95 = 1042.2
  p99 = 1230.3
Now try p100 (the max value):
  IndexError: IndexError('list index out of range')
  Why: int(100/100 * 58) = 58, but valid indices only go up to 57.

That is a real bug, not a contrived one: at exactly p100 the formula computes an index equal to the length of the list, which is always one past the last valid index. This is exactly the kind of off-by-one mistake that a naive percentile implementation runs into, and it is why the standard library has a dedicated function for this instead.

statistics.quantiles(): n and method, explained

Python's quantiles() function divides your data into equal-probability groups. Its documentation describes the signature as quantiles(data, *, n=4, method='exclusive') and explains: "divide data into n continuous intervals with equal probability. Returns a list of n - 1 cut points separating the intervals... set n to 100 for percentiles which gives the 99 cuts points that separate data into 100 equal sized groups." So for percentiles, you pass n=100, and the p-th percentile is at index p - 1 in the returned list (percentile 1 is index 0, percentile 99 is index 98).

q_exclusive = statistics.quantiles(without_timeouts, n=100, method="exclusive")
q_inclusive = statistics.quantiles(without_timeouts, n=100, method="inclusive")

print(f"{'p':>4} {'naive':>10} {'exclusive':>10} {'inclusive':>10}")
for p in (50, 90, 95, 99):
    print(f"{p:>4} {naive_percentile(without_timeouts, p):>10.1f} "
          f"{q_exclusive[p - 1]:>10.1f} {q_inclusive[p - 1]:>10.1f}")

   p      naive  exclusive  inclusive
  50       73.7       73.6       73.6
  90      893.6      896.2      864.3
  95     1042.2     1046.1     1014.7
  99     1230.3     1275.3     1167.8

A gotcha worth knowing: method choice changes your reported number

Look at p99: the naive method says 1230.3, the default exclusive method says 1275.3, and inclusive says 1167.8. That is a spread of over 100 milliseconds on the exact same 58 data points, purely from choosing a different (equally legitimate) interpolation method. If your team ever compares a p99 number from one tool against a p99 number from another tool and they do not match, this is very often why: they are using different quantile algorithms, not measuring different things.

Cross-checking with NumPy

arr = np.array(without_timeouts)
for p in (50, 90, 95, 99):
    print(f"  p{p} = {np.percentile(arr, p):.1f}")

  p50 = 73.6
  p90 = 864.3
  p95 = 1014.7
  p99 = 1167.8

NumPy's default np.percentile() (linear interpolation) matches statistics.quantiles(..., method="inclusive") exactly on this data. That is a useful fact to know if you are trying to reconcile numbers between the two libraries.

Step 5: Flagging Outliers Automatically With the IQR Method

You have been manually filtering out the 2 timeouts with v < 3000 so far. That works when you already know what an outlier looks like, but a real monitoring pipeline needs to detect outliers automatically. A common, simple technique is the interquartile range (IQR) method: compute the 25th percentile (Q1) and 75th percentile (Q3), take the range between them (the IQR), and flag anything more than 1.5 times that range beyond either edge.

data_sorted = sorted(LATENCY_MS)
n = len(data_sorted)

q1, q2, q3 = statistics.quantiles(data_sorted, n=4, method="inclusive")
iqr = q3 - q1
lower_fence = q1 - 1.5 * iqr
upper_fence = q3 + 1.5 * iqr

print(f"n = {n}")
print(f"Q1 = {q1:.1f} ms, Q2 (median) = {q2:.1f} ms, Q3 = {q3:.1f} ms")
print(f"IQR = Q3 - Q1 = {iqr:.1f} ms")
print(f"lower fence = Q1 - 1.5*IQR = {lower_fence:.1f} ms")
print(f"upper fence = Q3 + 1.5*IQR = {upper_fence:.1f} ms")

outliers = [v for v in data_sorted if v < lower_fence or v > upper_fence]
print(f"flagged as outliers ({len(outliers)} of {n}): {outliers}")

n = 60
Q1 = 65.9 ms, Q2 (median) = 73.7 ms, Q3 = 82.4 ms
IQR = Q3 - Q1 = 16.5 ms
lower fence = Q1 - 1.5*IQR = 41.2 ms
upper fence = Q3 + 1.5*IQR = 107.1 ms
flagged as outliers (11 of 60): [30.7, 818.1, 851.7, 893.6, 920.0, 1009.9, 1042.2, 1120.6, 1230.3, 5202.9, 5210.4]

A gotcha worth knowing: the IQR fence flags both directions

Look at the first value in that outlier list: 30.7. That is not a slow request, it is an unusually fast one, faster than 41.2 ms, the lower fence. The IQR method has no idea that "fast is good"; it only knows that a value sits statistically far from the bulk of the data, in either direction. When you build outlier detection for something like latency, decide up front whether you actually want to flag low-direction outliers too, or whether you should only apply the upper fence.

Also worth confirming: does this method produce false alarms on genuinely uniform data? Run the same fence math on just the 50 fast-path requests, with no slow or timeout traffic mixed in:

fast_only = sorted(v for v in LATENCY_MS if v < 200)
q1f, _, q3f = statistics.quantiles(fast_only, n=4, method="inclusive")
upper_fence_fast = q3f + 1.5 * (q3f - q1f)
would_flag = [v for v in fast_only if v > upper_fence_fast]
print(f"points flagged as outliers within the fast-only subset: {would_flag}")

points flagged as outliers within the fast-only subset: []

Zero false positives on genuinely uniform data. That is a good sign the method is behaving sensibly, not just flagging things at random.

Step 6: Correlation: Does Payload Size Actually Predict Latency?

Correlation measures how strongly two variables move together, on a scale from -1 (perfectly opposite) to 1 (perfectly together), with 0 meaning no linear relationship at all. Use the second dataset from Step 1, pairing request payload size (in KB) with latency:

from dataset import PAYLOAD_KB, LATENCY_FOR_PAYLOAD_MS

r = statistics.correlation(PAYLOAD_KB, LATENCY_FOR_PAYLOAD_MS)
cov = statistics.covariance(PAYLOAD_KB, LATENCY_FOR_PAYLOAD_MS)
slope, intercept = statistics.linear_regression(PAYLOAD_KB, LATENCY_FOR_PAYLOAD_MS)

print(f"n = {len(PAYLOAD_KB)}")
print(f"covariance(payload_kb, latency_ms) = {cov:.2f}")
print(f"correlation(payload_kb, latency_ms) = {r:.4f}")
print(f"linear_regression: latency ~= {slope:.2f} * payload_kb + {intercept:.2f}")

n = 50
covariance(payload_kb, latency_ms) = 217.90
correlation(payload_kb, latency_ms) = 0.8433
linear_regression: latency ~= 1.95 * payload_kb + 37.27

A correlation of 0.8433 is a strong positive relationship: bigger payloads really do tend to come with higher latency in this dataset. statistics.correlation() also has a built-in guardrail worth knowing about: it refuses to silently compute a wrong answer if your two lists do not line up.

try:
    statistics.correlation(PAYLOAD_KB, LATENCY_FOR_PAYLOAD_MS[:-1])
except statistics.StatisticsError as exc:
    print(f"correlation() on mismatched lengths raises: {exc!r}")

correlation() on mismatched lengths raises: StatisticsError('correlation requires that both inputs have same number of data points')

Correlation is not causation: the shuffle test

You have probably heard "correlation does not imply causation" as an abstract warning. Here is a way to see it concretely, on your own data. Take the exact same latency values (so the mean and standard deviation of the latency data do not change at all) and shuffle which payload they are paired with:

import random

shuffled_latency = LATENCY_FOR_PAYLOAD_MS[:]
random.seed(99)
random.shuffle(shuffled_latency)
r_shuffled = statistics.correlation(PAYLOAD_KB, shuffled_latency)

print(f"mean(latency) unchanged: {statistics.mean(shuffled_latency):.2f} ms "
      f"(original {statistics.mean(LATENCY_FOR_PAYLOAD_MS):.2f} ms)")
print(f"correlation after shuffling = {r_shuffled:.4f}")

mean(latency) unchanged: 71.32 ms (original 71.32 ms)
correlation after shuffling = -0.2326

Every individual latency value is identical to before (the mean proves it), but the correlation collapsed from a strong 0.8433 down to a weak -0.2326 the moment the pairing between payload size and latency was broken. Correlation is a property of the relationship between two specific paired values, not a property of either list on its own. That is exactly why correlation between two real-world measurements can hint at a connection, but can never prove one by itself: you would need to actually control the payload size experimentally to know it is causing the latency change, rather than both being driven by some third factor.

Step 7: Put It All Together: A Reusable describe() Function

Now combine everything above into one function you can call on any dataset. Save this as describe.py:

"""A reusable descriptive-statistics report."""
import statistics
from dataclasses import dataclass


@dataclass
class Summary:
    n: int
    mean: float
    median: float
    pvariance: float
    variance: float
    pstdev: float
    stdev: float
    q1: float
    q3: float
    iqr: float
    p95: float
    p99: float
    outliers: list


def describe(data):
    if len(data) < 2:
        raise ValueError("describe() needs at least 2 data points")

    data_sorted = sorted(data)
    q1, _, q3 = statistics.quantiles(data_sorted, n=4, method="inclusive")
    iqr = q3 - q1
    lower_fence = q1 - 1.5 * iqr
    upper_fence = q3 + 1.5 * iqr
    outliers = [v for v in data_sorted if v < lower_fence or v > upper_fence]

    q_inclusive = statistics.quantiles(data_sorted, n=100, method="inclusive")

    return Summary(
        n=len(data),
        mean=statistics.mean(data),
        median=statistics.median(data),
        pvariance=statistics.pvariance(data),
        variance=statistics.variance(data),
        pstdev=statistics.pstdev(data),
        stdev=statistics.stdev(data),
        q1=q1,
        q3=q3,
        iqr=iqr,
        p95=q_inclusive[94],
        p99=q_inclusive[98],
        outliers=outliers,
    )


def format_report(summary: Summary, unit: str = "ms") -> str:
    lines = [
        f"n            = {summary.n}",
        f"mean         = {summary.mean:.1f} {unit}",
        f"median       = {summary.median:.1f} {unit}",
        f"stdev (samp) = {summary.stdev:.1f} {unit}",
        f"pstdev (pop) = {summary.pstdev:.1f} {unit}",
        f"Q1 / Q3      = {summary.q1:.1f} / {summary.q3:.1f} {unit}",
        f"p95 / p99    = {summary.p95:.1f} / {summary.p99:.1f} {unit}",
        f"outliers     = {len(summary.outliers)} "
        f"({100 * len(summary.outliers) / summary.n:.1f}%)",
    ]
    return "\n".join(lines)


if __name__ == "__main__":
    from dataset import LATENCY_MS

    summary = describe(LATENCY_MS)
    print(format_report(summary))

Run it against the full, original 60-value dataset, timeouts included:

n            = 60
mean         = 362.6 ms
median       = 73.7 ms
stdev (samp) = 961.0 ms
pstdev (pop) = 953.0 ms
Q1 / Q3      = 65.9 / 82.4 ms
p95 / p99    = 1126.1 / 5206.0 ms
outliers     = 11 (18.3%)

This single report tells the whole story that the mean alone hid. The median (73.7 ms) says a typical request is fast. The p95 (1126.1 ms) shows that the slow cache-miss path is a real, measurable problem affecting the tail. And the p99 (5206.0 ms) is essentially reporting the timeout value itself, because with only 60 samples, "the 99th percentile" and "the second-worst thing that happened" end up being nearly the same number. That last point matters operationally: p99 on a small sample size can be extremely noisy and dominated by one or two events, which is worth remembering before paging someone over a p99 spike computed from a handful of requests.

Common Mistakes and Gotchas: A Recap

• Reporting only the mean. A single slow-but-rare event can move the mean far more than it should, while barely touching the median. Report both, or report percentiles.
• Mixing up population and sample variance. statistics.variance() divides by N - 1; statistics.pvariance() divides by N. Use sample variance when your data represents a sample of a larger population; use population variance when your data is the entire population you care about.
• Assuming NumPy's default matches the stdlib's default. np.var() defaults to ddof=0 (population variance). If you want it to match statistics.variance(), pass ddof=1 explicitly.
• Comparing percentiles across tools without checking the method. "Exclusive" and "inclusive" quantile methods, or different interpolation rules entirely, can disagree by a meaningful amount on the same data, especially at extreme percentiles like p99 on a small sample.
• Assuming IQR outlier detection only flags the "bad" direction. It flags anything statistically unusual in either direction; an unusually fast value gets flagged exactly like an unusually slow one.
• Treating a strong correlation as proof of causation. A high correlation coefficient tells you two things moved together in your specific dataset; it does not tell you one caused the other.

Step 8: Verify Everything With pytest

Manually re-running scripts and eyeballing the output is fine while you are exploring, but a reusable function deserves tests that lock its behavior in place. Save this as test_describe.py in the same folder as describe.py:

"""pytest suite for the descriptive-statistics helpers."""
import statistics

import pytest

from describe import describe


def test_describe_rejects_too_few_points():
    with pytest.raises(ValueError):
        describe([1.0])


def test_mean_matches_hand_calculation():
    data = [10.0, 20.0, 30.0, 40.0]
    summary = describe(data)
    assert summary.mean == pytest.approx(25.0)


def test_median_of_even_length_averages_middle_two():
    data = [10.0, 20.0, 30.0, 40.0]
    summary = describe(data)
    assert summary.median == pytest.approx(25.0)


def test_population_vs_sample_variance_differ_by_bessel_correction():
    data = [10.0, 12.0, 23.0, 23.0, 16.0, 23.0, 21.0, 16.0]
    summary = describe(data)
    n = len(data)
    assert summary.variance == pytest.approx(
        summary.pvariance * n / (n - 1)
    )


def test_stdev_is_sqrt_of_variance():
    data = [10.0, 12.0, 23.0, 23.0, 16.0, 23.0, 21.0, 16.0]
    summary = describe(data)
    assert summary.stdev == pytest.approx(summary.variance ** 0.5)


def test_extreme_outlier_moves_mean_far_more_than_median():
    baseline = [70.0] * 20
    with_outlier = baseline + [5000.0]
    base_summary = describe(baseline)
    outlier_summary = describe(with_outlier)

    mean_shift = outlier_summary.mean - base_summary.mean
    median_shift = outlier_summary.median - base_summary.median

    assert mean_shift > 200.0
    assert abs(median_shift) < 1.0


def test_iqr_fence_flags_the_injected_outlier():
    baseline = [70.0, 71.0, 69.0, 72.0, 70.5, 69.5, 71.5, 70.0] * 3
    with_outlier = baseline + [500.0]
    summary = describe(with_outlier)
    assert 500.0 in summary.outliers
    assert len(summary.outliers) == 1


def test_p95_and_p99_are_never_less_than_median():
    data = [10.0, 12.0, 23.0, 23.0, 16.0, 23.0, 21.0, 16.0, 45.0, 5.0]
    summary = describe(data)
    assert summary.p95 >= summary.median
    assert summary.p99 >= summary.p95


def test_correlation_requires_equal_length_inputs():
    with pytest.raises(statistics.StatisticsError):
        statistics.correlation([1, 2, 3], [1, 2])

Run it:

$ pytest test_describe.py -v
test_describe.py::test_describe_rejects_too_few_points PASSED
test_describe.py::test_mean_matches_hand_calculation PASSED
test_describe.py::test_median_of_even_length_averages_middle_two PASSED
test_describe.py::test_population_vs_sample_variance_differ_by_bessel_correction PASSED
test_describe.py::test_stdev_is_sqrt_of_variance PASSED
test_describe.py::test_extreme_outlier_moves_mean_far_more_than_median PASSED
test_describe.py::test_iqr_fence_flags_the_injected_outlier PASSED
test_describe.py::test_p95_and_p99_are_never_less_than_median PASSED
test_describe.py::test_correlation_requires_equal_length_inputs PASSED

============================== 9 passed in 0.05s ==============================

Notice what test_extreme_outlier_moves_mean_far_more_than_median actually checks: it injects one single extreme value into an otherwise uniform dataset and asserts, numerically, that the mean moves by more than 200 ms while the median moves by less than 1 ms. That test encodes the exact lesson from Step 2 as a permanent, automated check. If someone later "simplifies" your reporting code to only use the mean, this kind of test is how you would catch the regression in behavior, not just in syntax.

How to Confirm It All Works End to End

Put the four files (dataset.py, describe.py, and optionally the step scripts) in one folder and run through this checklist:

1. Run python dataset.py style imports and confirm you see the same 60-value list shown in Step 1. If you changed the random seed, your numbers will differ, and that is fine, the shape of the lesson stays the same.
2. Run python describe.py and confirm the report matches Step 7's output.
3. Run pytest test_describe.py -v and confirm all 9 tests pass.
4. As a sanity check on your own understanding, try this by hand: take describe.LATENCY_MS, remove the two timeout values, and predict whether the mean or the median will move more before you rerun describe(). If your prediction matches Step 2's result, you have internalized the core lesson.

Next Steps

Once mean, median, variance, percentiles, and correlation feel natural, a few directions to keep going:

• Apply describe() to a real dataset from your own systems, such as request durations pulled from your logs or a metrics backend, instead of simulated data.
• Read How to Correctly Benchmark LLM Inference Latency With Python and Ollama to see percentiles applied to a real, end-to-end benchmarking harness, including the concurrency effects that can distort them.
• Read How to Turn Slow SQL Queries Into Actionable Reliability Metrics With OpenTelemetry to see traffic-weighted ranking and anomaly detection built on top of the same descriptive-statistics foundation.
• If you need to compare two sets of judgments (not just describe one dataset), How to Calibrate an LLM-as-Judge Against Human Ratings With Cohen's Kappa covers a related but different statistical tool: agreement measurement.