Distributions — Naming the Common Shapes of Uncertainty

You have 50,000 recorded checkout times from your web shop. Someone asks how many customers take longer than twelve seconds. Easy — count the rows, divide by 50,000, done. Then they ask a second question: what fraction take longer than two minutes?

You look. Your longest observed checkout is 103 seconds. So your data says the answer is zero, which is obviously wrong — you just have not seen one yet. And your histogram is no help, because the far right of it is built from four or five lonely observations that could each have gone a different way. The bit of the distribution you care about most is the bit your data describes worst.

There is a second problem. To answer any question at all you had to keep all 50,000 numbers. To answer it on a phone, or inside a training loop that runs a million times, that is not viable. And you cannot do algebra with a histogram: there is no way to say "if I add an independent delay with these properties, what does the total look like?"

The fix is to stop describing the data point by point and instead say which shape it has. Checkout times are roughly exponential with a rate of about 0.1 per second. That is one number. From it you get any probability you like, including for the two-minute case you never observed, plus closed-form answers about sums and averages. The cost is that if you pick the wrong shape, everything downstream is confidently wrong. So the skill is knowing the shapes, and knowing how to tell when one does not fit.

What a named family buys you, and what it costs

That last row is not a footnote. Assume request latency is normal with mean 200 ms and standard deviation 30 ms and the formula hands you a 99th percentile of $200 + 2.326 \times 30 \approx 270$ ms. Real latency is right-skewed with a long tail, and the true p99 might be 900 ms. Your alerting thresholds, your capacity plan and your service-level objective are all built on a number that is off by a factor of three — not because the maths was wrong, but because the family was.

A named distribution replaces thousands of data points with two numbers. That is a bargain when the shape is right and a disaster when it is not, so always check the fit before you trust the tail.

The Normal distribution, and why it is everywhere

That looks forbidding until you read it in pieces, from the inside out:

• $(x - \mu)$ is the distance from the centre. Only distance appears, never $x$ on its own, which is why the curve is perfectly symmetric about $\mu$.
• Squaring it makes left and right equivalent and punishes far-away points much harder than nearby ones.
• Dividing by $2\sigma^2$ measures that squared distance in units of variance. A large $\sigma$ makes the exponent grow slowly, so the curve is wide and flat.
• The minus sign inside $e^{\,\cdot\,}$ turns "far from the centre" into "very little density", and the exponential makes that fall-off ferociously fast.
• $\frac{1}{\sigma\sqrt{2\pi}}$ does no conceptual work at all. It is the number that makes the total area come out to exactly 1.

Two parameters, $\mu$ and $\sigma$, and they mean exactly what you want them to mean: where the bump sits and how wide it is.

The 68-95-99.7 rule

For any normal distribution, roughly 68% of the mass lies within one standard deviation of the mean, 95% within two, and 99.7% within three. With latency $\mu = 200$ ms, $\sigma = 30$ ms: 68% of requests land between 170 and 230 ms, 95% between 140 and 260 ms, and 99.7% between 110 and 290 ms. Anything outside that last range happens about three times in a thousand — half of that on each side, so a request slower than 290 ms is roughly a one-in-750 event.

This is the fastest sanity check in statistics. If a "normal" dataset has 5% of its points beyond three standard deviations, it is not normal — it has heavy tails, and every tail probability you compute from a normal fit will be far too small.

Z-scores

The z-score restates a value as "how many standard deviations from the mean". A 275 ms request gives $z = (275 - 200)/30 = 2.5$. That single number is comparable across any two quantities measured in any units, which is exactly why standardising features before training is standard practice — it puts every column on the same footing. And because every normal turns into the same standard normal after this transformation, one table of probabilities serves all of them: $P(Z \le 2.5) = 0.9938$, so about 0.62% of requests are slower than 275 ms.

The Central Limit Theorem is why the normal keeps appearing

Take $X_1, X_2, \dots, X_n$ drawn independently from any distribution with a finite mean $\mu$ and a finite variance $\sigma^2$. As $n$ grows, the distribution of their average $\bar{X}_n$ approaches a normal distribution with mean $\mu$ and variance $\sigma^2/n$.

Read the conditions carefully, because they carry the content. The source distribution can be wildly skewed, discrete, bimodal — it does not matter. What must hold is independence and finite variance. Distributions with infinite variance do not obey it, and heavy-tailed real-world quantities behave badly for exactly this reason.

This is why so many measured quantities look normal: they are sums of many small independent contributions. It is also why the average of a sample is normal even when the sample is not, which is the foundation every confidence interval and t-test in the world stands on. Note the $\sigma^2/n$: quadruple your sample size and you halve the standard error, which is the arithmetic behind "you need four times the data to halve the error bar".

The Central Limit Theorem is a statement about averages, not about data. Your raw measurements need not be normal for the mean of them to be.

Counting things: Binomial and Poisson

Binomial: successes out of a fixed number of tries

Each of the three factors does one job. $p^k$ is the probability of the $k$ successes, $(1-p)^{n-k}$ the probability of the failures, and $\binom{n}{k}$ counts how many different orderings produce that same total. Use it when you have a fixed number of independent yes/no trials, each with the same probability.

Twenty visitors land on a checkout page and each converts with probability 0.15. Expected conversions: $np = 20 \times 0.15 = 3$. Variance: $np(1-p) = 20 \times 0.15 \times 0.85 = 2.55$, so a standard deviation of about 1.60. The probability of seeing exactly five:

Roughly a 10% chance. Note what that means for your dashboard: five conversions instead of the expected three is completely unremarkable, and treating it as a win would be reading noise.

Poisson: counts in a window, with no ceiling

Use it when events happen independently at a steady average rate and there is no natural maximum. Errors per hour, arrivals per minute, typos per page. The single parameter $\lambda$ is the average count per window, and remarkably it is also the variance.

Your service averages three errors an hour, so $\lambda = 3$. The chance of a clean hour is $P(X = 0) = e^{-3} = 0.0498$ — about one hour in twenty. The chance of exactly five errors:

Put that next to the binomial answer above: 0.101 against 0.103. That is not a coincidence. When $n$ is large and $p$ is small, the binomial converges to a Poisson with $\lambda = np$. Here $np = 3$ in both cases. The Poisson is the right choice when you cannot even define $n$ — you do not know how many people could have hit your server this hour, only how many did.

The mean-equals-variance property is also a diagnostic. If your observed counts have a variance much larger than their mean, the data is overdispersed and a Poisson model will understate the risk of a bad hour. That usually means the rate is not constant, and you need a mixture or a negative binomial instead.

Waiting for things: the Exponential

Poisson counts events per window; the exponential describes the gap between events. Same underlying process, opposite question. With errors arriving at $\lambda = 3$ per hour, the average gap is $1/3$ hour, or 20 minutes. The probability of going a whole hour without one is $e^{-3 \times 1} = 0.0498$ — the same 5%, arrived at from the other direction.

Memorylessness, plainly

The exponential has a property that feels wrong the first time:

In words: having already waited tells you nothing about how much longer you will wait. At the start, the chance of waiting more than 20 minutes for the next error is $e^{-3 \times (1/3)} = e^{-1} = 0.368$. After 40 quiet minutes, the chance of waiting another 20 is still 0.368. The process does not accumulate pressure and it is not "due".

This is the right model for genuinely random arrivals — cosmic rays, independent user requests, radioactive decay. It is the wrong model for anything that ages. A hard drive that has run for five years is more likely to fail this month than a new one; a bus that is ten minutes late is more likely to arrive in the next minute. Using an exponential for those understates the risk exactly when it matters.

Uniform and Beta: distributions over a bounded range

Uniform

Continuous uniform on $[a, b]$ has $f(x) = 1/(b-a)$, mean $(a+b)/2$ and variance $(b-a)^2/12$. Discrete uniform on $\{1, \dots, n\}$ gives every value probability $1/n$. Its flatness is its meaning: every value in range is equally plausible, so it is the natural "I know only the bounds" assumption.

You meet it constantly in ML plumbing. Glorot and He weight initialisation draw from a uniform whose width is set by the layer size. Random hyperparameter search samples uniformly over a range (or log-uniformly, when the sensible range spans orders of magnitude). And every other sampler is built on it: draw $u \sim U(0,1)$, apply the inverse CDF $F^{-1}(u)$, and out comes a sample from any distribution you like.

Beta: a distribution over a probability

The Beta lives entirely on $[0, 1]$, which makes it the natural way to express uncertainty about a probability itself. Beta$(\alpha, \beta)$ has mean $\alpha/(\alpha+\beta)$, and Beta$(1,1)$ is flat — total ignorance.

Its magic property is conjugacy with the binomial: if your belief before seeing data is Beta$(\alpha, \beta)$ and you then observe $s$ successes and $f$ failures, your belief afterwards is

You just add the counts. No integration, no sampling. Start from Beta$(1,1)$ and run 20 visitors through checkout, of whom 7 convert: the posterior is Beta$(1+7,\, 1+13) =$ Beta$(8, 14)$, with mean $8/22 = 0.364$. Compare that to the raw estimate $7/20 = 0.35$.

The difference is small here, and enormous with little data. Suppose a variant gets 2 conversions from 3 visitors. Raw estimate: 0.667, and a naive dashboard now declares it the best-performing variant on the site. The Beta posterior is Beta$(3, 2)$ with mean $3/5 = 0.6$ and a distribution so wide it overlaps almost everything. The prior acts as a couple of imaginary prior observations that stop three data points from shouting.

This is exactly the machinery behind Thompson sampling. Keep a Beta posterior per variant, draw one random sample from each, and send the next user to whichever variant drew highest. Variants with little data have wide posteriors and get explored; variants that are genuinely better win more draws over time. You get the explore-exploit trade-off for free, out of arithmetic on two counters.

Two distributions that describe estimates rather than data

The chi-square with $k$ degrees of freedom is what you get by squaring $k$ independent standard normals and adding them:

Because it is a sum of squares it lives on $[0, \infty)$ and is right-skewed. Nobody models raw data with it. It appears whenever a test statistic is built from squared deviations: goodness-of-fit tests, tests of independence in contingency tables, and the sampling distribution of a sample variance.

Student's t exists because of one honest admission. To standardise a sample mean you want $(\bar{x} - \mu)/(\sigma/\sqrt{n})$, but you do not know $\sigma$ — you estimate it as $s$ from the same small sample. That estimate is itself uncertain, and sometimes too small, which makes the ratio blow up more often than a normal would. The t distribution accounts for that with heavier tails, controlled by degrees of freedom $\nu = n-1$.

The practical size of the effect: for a 95% two-sided interval, the normal uses a multiplier of 1.96. With $\nu = 4$ (a sample of five) the t needs 2.776 — a 42% wider interval for the same confidence. By $\nu = 30$ the multiplier is 2.04 and the difference has almost vanished. That is the whole story of t versus normal: it is the price of not knowing $\sigma$, and it is only expensive when $n$ is small.

Choosing a family without guessing

(Beta variance is $\alpha\beta / [(\alpha+\beta)^2(\alpha+\beta+1)]$ — rarely needed by hand.)

In practice you can get to the right family with four questions, asked in this order:

1. Is the outcome a count? If there is a fixed number of trials, binomial. If events accumulate in a window with no ceiling, Poisson.
2. Is it a waiting time or a duration? Exponential if the process has no memory; something with an ageing hazard otherwise.
3. Is it a proportion, bounded between 0 and 1? Beta.
4. Is it unbounded in both directions, roughly symmetric, plausibly a sum of many small effects? Normal. If it is unbounded but heavily skewed, consider modelling its logarithm instead.

Fitting and checking in SciPy

Every distribution in scipy.stats exposes the same four methods, and learning them once covers all of them.

1from scipy import stats2import numpy as np34lat = stats.norm(loc=200, scale=30)   # a "frozen" distribution56lat.pdf(275)        # 0.00058 -- a DENSITY, not a probability7lat.cdf(275)        # 0.99379 -- P(X <= 275)8lat.sf(275)         # 0.00621 -- the survival function, 1 - cdf9lat.ppf(0.95)       # 249.35  -- the p95: inverse of the cdf10lat.rvs(size=5, random_state=0)       # five random draws1112stats.binom(n=20, p=0.15).pmf(5)      # 0.1028  discrete: .pmf not .pdf13stats.poisson(mu=3).pmf(5)            # 0.1008  the same answer, near enough14stats.expon(scale=1/3).sf(1.0)        # 0.0498  P(wait > 1 hour)15stats.beta(a=8, b=14).mean()          # 0.3636  the posterior above

Fitting is one call, and so is the sanity check that stops you trusting a bad fit. A Kolmogorov-Smirnov test compares your data's empirical CDF against the fitted one; a small p-value means the family is being rejected.

1rng = np.random.default_rng(0)2data = rng.exponential(scale=4.0, size=2000)   # right-skewed, definitely not normal34mu, sigma = stats.norm.fit(data)5print(stats.kstest(data, stats.norm(mu, sigma).cdf).pvalue)6# effectively zero -- the normal family is decisively rejected78loc, scale = stats.expon.fit(data, floc=0)9print(stats.kstest(data, stats.expon(loc, scale).cdf).pvalue)10# comfortably large -- no evidence against the exponential fit1112# The normal fit still "works": it happily returns a p99...13print(stats.norm(mu, sigma).ppf(0.99))    # ~13.614print(stats.expon(loc, scale).ppf(0.99))  # ~18.6  the honest answer

Both fits happily produce a p99. Only one of them is right, and nothing but the goodness-of-fit check tells you which.

So always pair a fit with a check — a KS test, or simply plotting the fitted density over the histogram and looking hard at the tail.

Where these show up in the models you build

Linear regression is a normal assumption in disguise. Least squares is the maximum-likelihood fit precisely when the residuals are normal with constant variance. That assumption is what makes the reported standard errors and p-values meaningful. So plot your residuals: if they are skewed, or fan out as the prediction grows, your coefficients may still be usable but your uncertainty estimates are not.

Naive Bayes is a pile of likelihoods. Gaussian naive Bayes fits one $\mu$ and one $\sigma$ per feature per class and multiplies the resulting densities. Multinomial naive Bayes swaps in count-based likelihoods for text. Same algorithm, different distributional assumption — and choosing the wrong one is the usual reason a naive Bayes model underperforms.

Initialisation is variance engineering. He and Glorot initialisation draw weights from a normal or uniform with variance scaled to the layer's width, chosen so activation variance neither explodes nor collapses as signals pass through many layers. Get it wrong and a deep network either saturates or produces nothing but noise.

Several standard layers are just sampling. Dropout multiplies activations by a Bernoulli mask. A variational autoencoder's encoder outputs a $\mu$ and a $\log \sigma^2$ per latent dimension and samples a Gaussian from them. Data augmentation applies uniformly random crops and flips.

Count targets need count models. Predicting clicks, purchases or defects with ordinary least squares lets the model output negative counts and assume constant variance when variance actually grows with the mean. Poisson regression fixes both. And when the observed variance runs well ahead of the mean, that is your signal that a single Poisson rate does not describe your data and something richer is required.

The habit worth building is small. Before you fit anything, plot the target and ask which of these shapes it looks like. After you fit, plot what the model implies and compare it to what you actually saw. Most modelling failures are not clever failures — they are a right-skewed quantity that someone assumed was symmetric.