The distributions you will actually meet

A distribution describes how a quantity tends to come out when it varies: request latencies, support tickets per hour, measurement error. Knowing which one fits a quantity tells you what counts as normal, what is genuinely unusual, and which summary numbers will mislead you. This section is a field guide to the few you actually meet, with the code to sample and fit each, and a warning about the heavy-tailed ones, where the average is not a useful number at all.

Two numbers summarise any of them#

Whatever the shape, two numbers describe it. The mean is the centre. The variance is the spread: take each value's distance from the mean, square it so that being below counts the same as being above, and average those. Squaring leaves the units wrong — a variance of latencies is in milliseconds squared — so the number usually quoted is its square root, the standard deviation, which is back in the original units and can be compared with the mean.

Both are expectations, which is the word for a probability-weighted average: the mean is the expectation of the value, and the variance is the expectation of the squared distance from the mean. When you have a distribution written down, that weighting is by each outcome's probability. When you have a pile of data instead, every observation is equally likely, so the weighting collapses into the ordinary average and the two definitions become the same arithmetic.

import numpy as np

latencies = np.array([90, 95, 100, 105, 110])
mean = latencies.mean()
variance = ((latencies - mean) ** 2).mean()
print(mean, variance, np.sqrt(variance))
# 100.0 50.0 7.0710678118654755

# Same mean, far more spread. The mean alone tells you nothing about this.
spiky = np.array([10, 20, 100, 180, 190])
print(spiky.mean(), spiky.std())
# 100.0 76.15773105863909

A field guide#

Quantity                                   Distribution to try first
Did the user convert? (yes/no)             Bernoulli; conversions out of n: binomial
Support tickets per hour                   Poisson
Time until the next request                exponential
Measurement error, heights, averages       normal
Request latency, file sizes, incomes       log-normal or other heavy-tailed
Which of k classes                         categorical

The rest of this section takes that table one row at a time. For each distribution: the signal in the problem statement that should make you reach for it, the two numbers that describe it, and a check in code. Recognition is the skill being trained here — the maths of each is a lookup, but noticing that a quantity is a count per period rather than a proportion is not.

Yes or no: Bernoulli and binomial#

The recognition signal is a single trial with two outcomes. Did the user convert, did the request succeed, did the eval case pass. That is a Bernoulli variable: it takes the value 1 with probability p and 0 otherwise. Its mean is p, because the average of a pile of ones and zeros is the share of ones. Its variance is p times one minus p, which is largest at p equals a half and shrinks to nothing at either end — a system that succeeds 99 per cent of the time varies far less from sample to sample than one that succeeds half the time.

The binomial is the same trial repeated n times, counting successes. Its mean is n times p and its variance is n times p times one minus p, which is just n independent Bernoulli variances added. Any pass rate on an eval set of fixed size is binomial, which is why the margin-of-error formula for a proportion works and where it comes from.

rng = np.random.default_rng(0)

one = rng.binomial(1, 0.3, 200_000)              # Bernoulli(0.3)
print(one.mean(), one.var())                     # p, then p(1 - p) = 0.21
# 0.299145 0.209657268975

k = rng.binomial(20, 0.3, 200_000)               # Binomial(20, 0.3)
print(k.mean(), k.var())                         # np = 6, then np(1 - p) = 4.2
# 6.00427 4.2180517671

Counts per period: Poisson#

The recognition signal is a count of events in a window, where there is no natural maximum. Support tickets per hour, errors per minute, arrivals per second. The binomial needs an n; the Poisson does not, which is exactly what makes it the right shape for counts. It has one parameter, the average rate, and its defining oddity is that its variance equals its mean. If you measure a count whose variance is much larger than its mean, the events are not arriving independently — they are bursty, and a Poisson model will badly understate your worst hours.

The Poisson and the exponential are the same story told twice. If events arrive independently at a steady rate, the count per window is Poisson and the gaps between consecutive events are exponential. You can derive one from the other in code rather than on paper: generate exponential gaps, add them up, and count how many land in the first hour.

lam = 12.0                                             # 12 tickets an hour on average
gaps = rng.exponential(1 / lam, size=(50_000, 40))     # exponential time between arrivals
arrival = np.cumsum(gaps, axis=1)
per_hour = (arrival <= 1.0).sum(axis=1)

print(per_hour.mean(), per_hour.var())     # a Poisson has variance equal to its mean
# 11.9993 11.954779509999998

from scipy import stats
print((per_hour == 20).mean(), stats.poisson(mu=12).pmf(20))
# 0.00976 0.00968203216821916

Waiting times: exponential#

The recognition signal is time until the next thing happens, when the thing has no memory of how long you have already waited. Its mean is one over the rate, so at twelve arrivals an hour the average gap is five minutes, and its standard deviation equals its mean, which is why waiting times feel so unpredictable.

Memorylessness is the property worth checking, because it is the one that violates intuition. Having already waited three minutes tells you nothing about how much longer you will wait: the distribution of remaining time is identical to the distribution you faced at the start. This is why a retry after a timeout on an exponential service does not improve your position, and it is exactly the assumption that fails for real queues, where a request that has already been slow is usually stuck on something and will keep being slow.

w = rng.exponential(2.0, 500_000)        # mean wait of 2

print((w > 2).mean())                    # 0.367862   chance of waiting past 2
print((w[w > 3] > 5).mean())             # 0.36882405745062835
# given you have already waited 3, the chance of another 2 is the same number

Measurements: normal#

The recognition signal is a quantity built from many small independent contributions: heights, measurement error, and above all averages of anything at all, by the central limit theorem. It takes two parameters, the mean and the standard deviation, and those two describe it completely — there is no third number to learn. The figures worth memorising are that 68 per cent of the mass lies within one standard deviation of the mean, 95 per cent within two, and 99.7 per cent within three, which is where the 1.96 in every confidence interval comes from.

z = rng.normal(0, 1, 500_000)
for k_ in (1, 2, 3):
    print(k_, (np.abs(z) < k_).mean())
# 1 0.682944
# 2 0.954696
# 3 0.997334

Which of k: categorical, and the softmax#

The recognition signal is a choice among k options with no ordering, and it is the distribution you will meet most often in this field, because it is what a language model emits at every step. A categorical distribution is just k probabilities that sum to one. Models do not produce those directly; they produce k unbounded scores, and the softmax turns scores into a categorical distribution by exponentiating and dividing by the total.

Sampling from it is then one call. The temperature knob divides the scores before the exponential, so a temperature below one exaggerates the gaps and concentrates the mass on the leader, while a temperature above one flattens everything towards uniform. Watching the sampled frequencies match the computed probabilities is the check that the sampler is doing what you think.

logits = np.array([4.0, 3.0, 1.0, 0.5, 0.0])       # raw model scores over 5 tokens

def softmax(z, T=1.0):
    e = np.exp(z / T - (z / T).max())
    return e / e.sum()

for T in (0.5, 1.0, 2.0):
    p = softmax(logits, T)
    draws = rng.choice(len(logits), size=100_000, p=p)
    print(T, np.round(p, 3), np.round(np.bincount(draws, minlength=5) / 100_000, 3))

# 0.5 [0.878 0.119 0.002 0.001 0.   ] [0.878 0.118 0.002 0.001 0.   ]
# 1.0 [0.682 0.251 0.034 0.021 0.012] [0.682 0.251 0.034 0.021 0.013]
# 2.0 [0.468 0.284 0.104 0.081 0.063] [0.467 0.285 0.104 0.082 0.063]

With SciPy#

SciPy gives every one of these the same four methods, which is worth knowing because it means learning one distribution object teaches you all of them. The pmf or pdf is the height at a value, the cdf is the mass at or below a value, ppf is the inverse of the cdf and is how you get a percentile, and rvs draws samples. Nearly every question you have about a distribution is one of those four.

print(stats.binom(n=1000, p=0.03).interval(0.95))   # (20.0, 41.0) conversions to expect
print(stats.poisson(mu=12).pmf(20))                 # 0.00968203216821916
print(stats.norm(loc=0, scale=1).cdf(1.96))         # 0.9750021048517795

Heavy tails break averages#

Latency is not normally distributed. Most requests are fast and a few are extremely slow. The mean is dragged up by the slow ones and describes almost nobody's experience. That is why services report the median and high percentiles such as p95 and p99.

lat = rng.lognormal(mean=5, sigma=1, size=100_000)       # a realistic latency shape
print(f"mean {lat.mean():.0f}  median {np.median(lat):.0f}  p99 {np.percentile(lat, 99):.0f}")
# mean 245  median 147  p99 1527

# What a normal model would have predicted for the p99:
print(round(float(lat.mean() + 2.326 * lat.std())))
# 996

Fitting the right family instead takes one line and gives back the parameters you generated with. stats.lognorm.fit(lat, floc=0) returns a shape of 1.003 and a log scale of 4.999, against the 1.0 and 5.0 used to make the data. The reason to fit rather than to eyeball is that the fitted object then answers questions the sample cannot: the chance of exceeding a timeout you have never yet hit, or the percentile a new SLA would land on.