Statistics for IT — Math for IT

TL;DR

The mean is pulled up by a few extreme values; the median, the middle value, is not. Latency data almost always has such a long tail.
The p95 is the value that about 95% of observations stay at or below. There are several definitions, and on small samples they give clearly different answers; Python's statistics.quantiles default can even return a value above the maximum.
p95 and p99 latency matter more than the mean because they describe what the slowest users actually experience, and one page load often waits for many requests.
The standard deviation measures how far values typically lie from the mean; use the sample version (divide by n − 1) for a sample.
A strong correlation shows that two quantities move together, not that one causes the other.

Describing data with a few numbers

Monitoring tools, ticket systems and smart meters produce far more numbers than anyone can read. Descriptive statistics summarize a set of numbers in a few values: where the middle is, how spread out the values are, and how extreme the worst cases are. Every dashboard you will ever look at is built from these summaries, so it pays to know exactly what each one means and what it hides.

The running example on this page is response times (latencies) of Kestrel Energy's customer portal, in milliseconds (ms). Python's statistics module computes everything shown here; the SQL aggregate functions such as AVG do the same inside a database.

Mean, median and mode

Ten requests to the portal took these times:

112  98  105  101  97  120  103  99  950  104

The mean

The mean (arithmetic mean, the everyday "average") is the sum of all values divided by how many there are. For values x_1, x_2, \ldots, x_n:

\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i

Read it as: "add up all x_i for i from 1 to n, then divide by n". The bar over the x is the usual symbol for a mean. For the ten requests:

Sum: 112 + 98 + 105 + 101 + 97 + 120 + 103 + 99 + 950 + 104 = 1,889.
Divide by 10: 188.9 ms.

Nine of the ten requests took between 97 and 120 ms, yet the mean is 188.9 ms, slower than every one of them. The single 950 ms request has dragged it up.

The median

The median is the middle value after sorting: half of the values are at or below it, half at or above. With an even number of values there are two middle values, and the median is their mean.

Sort: 97, 98, 99, 101, 103, 104, 105, 112, 120, 950.
Ten values, so the middle ones are the 5th and 6th: 103 and 104.
Median: (103 + 104) / 2 = 103.5 ms.

The median barely notices the 950 ms request: if it had taken 95,000 ms the median would still be 103.5. That property is called being robust. The mode, the most frequent value, is a third kind of middle; it is mainly useful for categories ("which error code occurs most?").

import statistics as st

latencies = [112, 98, 105, 101, 97, 120, 103, 99, 950, 104]   # ms
print(st.mean(latencies), st.median(latencies))
print(st.mode([404, 500, 404, 503]), st.multimode([1, 1, 2, 2]))

188.9 103.5
404 [1, 2]

Outliers and skew

An outlier is an observation that lies far away from the rest of the data, such as the 950 ms request. An outlier can be a measurement error, but in IT it is often the most interesting value in the set: a timeout, a lock wait, a garbage-collection pause, a failing disk. Never delete outliers just to make a chart look calmer; find out what caused them.

A common rule of thumb for flagging outliers (from the statistician John Tukey) uses the quartiles, the values that cut the sorted data into four equal parts. The first quartile Q1 has a quarter of the data below it, the third quartile Q3 three quarters. Their difference is the interquartile range, IQR = Q3 − Q1, the width of the middle half. Values more than 1.5 IQR below Q1 or above Q3 are flagged:

import statistics as st

data = [120, 85, 97, 110, 300, 101, 99, 105, 92, 88]     # ms
q1, q2, q3 = st.quantiles(data, n=4, method="inclusive")
iqr = q3 - q1
low, high = q1 - 1.5 * iqr, q3 + 1.5 * iqr
print(q1, q2, q3, iqr)
print(low, high)
print([x for x in data if x < low or x > high])

93.25 100.0 108.75 15.5
70.0 132.0
[300]

Applying Tukey's fences per group, and what to do with each flag, is shown under finding outliers.

Latency data is almost never symmetric. Most requests are fast and bunched together, a few are much slower, and nothing can be faster than zero. The values form a long tail to the right: the data is right-skewed. In right-skewed data the mean is usually pulled above the median, so a mean latency describes hardly any real request.

Percentiles

A percentile describes a position in the sorted data: the 95th percentile, written p95, is a value that about 95% of the observations are at or below. If the p95 latency is 300 ms, then roughly 1 request in 20 is slower than 300 ms. The median is the 50th percentile, and the quartiles are p25, p50 and p75. Monitoring tools show p50, p90, p95, p99 and p99.9 (often written p999).

"About" and "roughly" are there on purpose. With a finite set of values there is usually no value that has exactly 95% of the data below it, so a rule is needed, and there are several rules in use. They agree on large data sets but not on small ones.

The nearest-rank method

The simplest definition: sort the n values and take the value at rank

r = \lceil p \cdot n \rceil

counting from 1, where p is the percentile as a fraction (0.95 for p95) and \lceil \; \rceil means round up to a whole number. The result is always one of the measured values, and it is the smallest value with at least a share p of the data at or below it. Worked example with the numbers 1 to 20:

p \cdot n = 0.95 \cdot 20 = 19, already whole, so r = 19.
The 19th value in sorted order is 19, so p95 = 19.
Check: 19 of the 20 values (95%) are at or below 19.

Several pages on this site use nearest rank, for example analyzing structured logs in Python and the Kanban service level expectation.

Methods that interpolate

Other methods treat the sorted data as points on a line and interpolate between the two neighboring values, so the answer can lie between measured values. They differ in where they put the position:

Inclusive (linear) method: the minimum is the 0th percentile and the maximum the 100th, and the position runs evenly in between: position 1 + p \cdot (n - 1), counting from 1. For 1 to 20 the position of p95 is 1 + 0.95 \cdot 19 = 19.05, so the answer is 5% of the way from 19 to 20: 19.05. This method never goes outside the data. Spreadsheets call it PERCENTILE.INC, and it is the default in many data tools.
Exclusive method: position p \cdot (n + 1). It assumes the true range of values extends beyond what you measured. For 1 to 20 the position is 0.95 \cdot 21 = 19.95, so the answer is 19.95. Spreadsheets call it PERCENTILE.EXC.

Three methods, three answers for the same 20 numbers: 19, 19.05 and 19.95. Here the differences are small. On a small sample with a long tail they are not.

Percentiles in Python

statistics.quantiles(data, n=100) returns the 99 cut points between 100 equal groups, so index 94 (counting from 0) is p95. Its default method is "exclusive"; pass method="inclusive" for the other one. Python has no nearest-rank function, but it is two lines:

import math
import statistics as st

data = list(range(1, 21))      # 1, 2, ..., 20
print(st.quantiles(data, n=100)[94])                      # default: exclusive
print(st.quantiles(data, n=100, method="inclusive")[94])
rank = math.ceil(0.95 * len(data))                         # nearest rank
print(sorted(data)[rank - 1])
print(st.quantiles(data, n=4), st.quantiles(data, n=4, method="inclusive"))

19.95
19.05
19
[5.25, 10.5, 15.75] [5.75, 10.5, 15.25]

Now a realistic small sample: ten portal requests, the slowest of which took 410 ms.

import math
import statistics as st

# p95 latency: nearest rank versus interpolation on the same 10 values (ms)
latency = [87, 91, 95, 98, 102, 105, 110, 118, 140, 410]
rank = math.ceil(0.95 * len(latency))
print("nearest rank:", sorted(latency)[rank - 1])
print("exclusive:", st.quantiles(latency, n=100)[94])
print("inclusive:", st.quantiles(latency, n=100, method="inclusive")[94])
print("maximum:", max(latency))

nearest rank: 410
exclusive: 531.5
inclusive: 288.5
maximum: 410

The default method reports a p95 of 531.5 ms although no request took longer than 410 ms. What happened: the exclusive position is 0.95 \cdot 11 = 10.45, beyond the 10th and last value. Python does not stop at the maximum; it continues the straight line through the last two values (140 and 410) another 0.45 of a step: 410 + 0.45 \cdot (410 - 140) = 531.5. That is a deliberate design choice (the method assumes the data could go higher than what was measured), but it is rarely what an SLO (service level objective) report should say. Other tools that use the exclusive positions usually refuse such a percentile or return the maximum instead.

With the exclusive method a p95 stays inside the data only when 0.95 \cdot (n + 1) \le n, that is when there are at least 19 values. More generally: a p95 from fewer than about 20 values, or a p99 from fewer than about 100, is mostly a statement about one or two individual requests.

Why p95 and p99 matter more than the mean

Here are 100 portal requests from a quiet morning, drawn as a histogram: each bar counts the requests in a 25 ms range.

A typical right-skewed latency distribution: the mean lies to the right of the median, and the p95 and p99 sit far out in the tail.

The mean of 106.7 ms suggests a fast portal, and for most requests it is. But the slowest 5% of requests took 178 ms or more, and the slowest 1% at least 303 ms. Those slow requests are not spread evenly over imaginary "average users"; they are real people who waited, and they are often the heaviest users (big accounts, long histories) or the ones hitting a struggling server. The mean hides them completely: you could double the time of the slowest 5% and the mean would rise only modestly.

The effect is stronger because one action rarely means one request. A portal page that loads its data with 50 backend calls is only as fast as the slowest of them. If each call has a 1% chance of being slower than the p99, the chance that a page load hits at least one slow call is (by the complement rule from Probability, assuming independent calls) 1 - 0.99^{50} \approx 0.39:

# Chance that a page load with this many backend calls
# includes at least one call slower than the p99
for calls in (1, 10, 50, 100):
    print(f"{calls:>3} call(s): {1 - 0.99 ** calls:.1%}")

  1 call(s): 1.0%
 10 call(s): 9.6%
 50 call(s): 39.5%
100 call(s): 63.4%

So with 50 calls, the "1 in 100" slow case shows up in nearly 4 out of 10 page loads. That is why SLOs are usually written on percentiles ("99% of requests within 300 ms") rather than on the mean, and why golden-signal dashboards show p95 and p99 next to the median. See alerting on service level objectives for how such targets drive alerts.

Spread: variance and standard deviation

Two services can have the same mean response time while one is steady and the other jumps about. A measure of spread captures that difference. The idea: look at how far each value lies from the mean (its deviation), and summarize those distances.

Worked example: eight tickets were resolved in 2, 4, 4, 4, 5, 5, 7 and 9 hours.

Mean: (2 + 4 + 4 + 4 + 5 + 5 + 7 + 9) / 8 = 40 / 8 = 5 hours.
Deviations from the mean: −3, −1, −1, −1, 0, 0, 2, 4. They always add up to 0, so their plain average is useless.
Square them to make them all positive: 9, 1, 1, 1, 0, 0, 4, 16. Their sum is 32.
Average the squares: 32 / 8 = 4. This is the variance.
Take the square root to get back to hours: \sqrt{4} = 2 hours. This is the standard deviation.

The variance is the mean of the squared deviations from the mean, and the standard deviation is its square root, a "typical distance from the mean" in the same unit as the data. The Greek letter \sigma (sigma) is the usual symbol:

\sigma^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2 \qquad \sigma = \sqrt{\sigma^2}

Sample or population?

The formula above divides by n: it describes a complete population, all the values there are. Usually you have a sample, such as some of last month's tickets, and want to estimate the spread of the whole process. A sample's values lie, on average, a little closer to their own mean than to the true mean, so dividing by n slightly underestimates the spread. The sample variance divides by n - 1 instead: 32 / 7 \approx 4.571, giving a sample standard deviation of about 2.14 hours.

import statistics as st

hours = [2, 4, 4, 4, 5, 5, 7, 9]
print(st.pvariance(hours), st.pstdev(hours))      # population: divide by n
print(st.variance(hours), round(st.stdev(hours), 3))   # sample: divide by n - 1

latencies = [112, 98, 105, 101, 97, 120, 103, 99, 950, 104]
print(round(st.stdev(latencies), 1), round(st.pstdev(latencies), 1))

4 2.0
4.571428571428571 2.138
267.5 253.8

The last line shows that the standard deviation, like the mean, is very sensitive to outliers: one 950 ms request turns a set of values that mostly differ by a few milliseconds into a "typical distance" of over 250 ms. The number of standard deviations a value lies from the mean is its z-score, a simple basis for anomaly alerts; see the z-score idea.

How far an estimate from a sample can be trusted (confidence intervals, sampling bias) is on Statistics for Analysts.

Distributions

A distribution describes how often each value (or range of values) occurs; the histogram above is a picture of one. Many textbook methods assume the normal distribution, the symmetric bell curve. For normally distributed data a handy rule of thumb holds: about 68% of values lie within one standard deviation of the mean, about 95% within two and about 99.7% within three.

Heights of people and small measurement errors are roughly normal. Latencies, file sizes, ticket resolution times and download counts are not: they are right-skewed with long tails. Applying the 68–95–99.7 rule or other normal-distribution formulas to them gives wrong answers. The example below treats the 100 requests from the histogram as if they were normal, with the same mean and standard deviation, and compares the percentiles that model predicts with the real ones. For skewed data, describe the middle with the median and the tail with percentiles, and look at the histogram before choosing a method.

import math
import random
import statistics as st

random.seed(18)
latency = [round(40 + random.lognormvariate(4.0, 0.55)) for _ in range(100)]
ranked = sorted(latency)

def nearest_rank(p):
    return ranked[math.ceil(p * len(ranked)) - 1]

mean, sd = st.mean(latency), st.stdev(latency)
print(f"median {st.median(latency)}  mean {mean:.1f}  stdev {sd:.1f}")
print(f"p95 {nearest_rank(0.95)}  p99 {nearest_rank(0.99)}  max {max(latency)}")
print(f"normal model: p95 {st.NormalDist(mean, sd).inv_cdf(0.95):.1f}  p99 {st.NormalDist(mean, sd).inv_cdf(0.99):.1f}")
print(sum(1 for x in latency if x > nearest_rank(0.95)), "requests slower than the p95")

median 96.0  mean 106.7  stdev 43.7
p95 178  p99 303  max 326
normal model: p95 178.5  p99 208.3
5 requests slower than the p95

The normal model happens to get the p95 about right for this data, but it puts the p99 at 208 ms, while the real p99 is 303 ms. A bell curve badly underestimates how slow the tail gets, and the tail is exactly what an SLO is about.

Correlation is not causation

Do two quantities move together? Kestrel compared the outside temperature with the gas usage of a group of customers over five months:

Average temperature (°C)	4	7	10	14	18
Gas usage (m³ per month)	520	470	430	360	300

Correlation in statistics measures how closely two quantities follow a straight-line relationship. The usual measure, Pearson's correlation coefficient r, runs from −1 to 1: close to 1 means "when one goes up, the other goes up", close to −1 means "when one goes up, the other goes down", and near 0 means no straight-line relationship. Here r \approx -0.999: warmer months, less gas, almost perfectly. A fitted straight line (linear regression, see classic machine learning) says usage drops by about 15.7 m³ per degree:

import statistics as st

temps = [4, 7, 10, 14, 18]          # outside temperature, degrees C
usage = [520, 470, 430, 360, 300]   # gas usage, m3 per month (fictional)
print(round(st.correlation(temps, usage), 3))
fit = st.linear_regression(temps, usage)
print(round(fit.slope, 2), round(fit.intercept, 2))

-0.999
-15.73 582.74

In this case the cause is clear (cold weather makes people heat their homes), but correlation alone never proves it. Two quantities can move together because a third factor drives both, or by coincidence. A real IT example: servers with more installed monitoring agents have more incidents. The agents do not cause incidents; the most critical servers get both the most monitoring and the most load. Before acting on a correlation, ask what else changed at the same time, and test the cause directly where you can.

Confounding variables, spurious correlations and Simpson's paradox, worked through on real data, are on Exploratory Data Analysis.