Probability for Machine Learning

Most models output a probability or a score rather than a verdict: this email is 0.93 spam, this patient is low-risk, this forecast is within a few percent. Probability is the language those outputs are written in.

It is also where the expensive mistakes happen. Trusting a test that is 99% accurate, treating repeated measurements as if they were separate people, reading a conditional probability as if it were a plain one — these look like modeling errors and they are all probability errors.

A probability is a proportion

Start with the least mysterious version. If you do something many times and count how often a particular result comes up, the fraction you get estimates the probability. Saying a coin has probability 0.5 of landing heads is saying that if you flip it enough times, about half will be heads.

How many times is enough depends on how close you need to get, and the numbers below show the scale of it.

import numpy as np

rng = np.random.default_rng(0)
flips = rng.random(100_000) < 0.5           # True = heads

for n in (10, 100, 1_000, 10_000, 100_000):
    print(f"  after {n:7,} flips: heads {flips[:n].mean():.4f}")
# after      10 flips: heads 0.3000
# after     100 flips: heads 0.4400
# after   1,000 flips: heads 0.4730
# after  10,000 flips: heads 0.4990
# after 100,000 flips: heads 0.5010

Ten flips gave 30% heads. That is not a broken coin — it is what ten flips look like. In this run the proportion is still 0.4730 at a thousand flips and 0.5010 at a hundred thousand, so how many trials are enough depends on how close you need to be.

Random variables and where their values land

A random variable is just a number attached to an uncertain outcome. The number of heads in ten flips is one. Tomorrow’s temperature is one. A customer’s spend next month is one. By convention it gets a capital letter, \(X\), and a particular value it might take gets a small one, \(x\).

A distribution says how the values are spread out. When the variable counts things — heads, clicks, defects — you can ask directly for the probability of each value, written \(P(X = x)\), and those probabilities add up to 1. That is the easy case.

When the variable is continuous and has a density — a temperature, a price — the probability of any single exact value is zero. (Having infinitely many possible values is not by itself the reason: a count with no upper bound also has infinitely many values, and each still carries positive probability.) What is meaningful is the probability of landing in a range, which is the area under the density \(p(x)\) over that range. For a thin slice, height times width approximates that area, and for a uniform density it is exact.

A density can exceed 1 because it is a height, not a probability.

import numpy as np

rng = np.random.default_rng(1)
x = rng.uniform(0.0, 0.1, size=200_000)     # spread evenly over a width of 0.1

print("largest value seen :", round(float(x.max()), 4))
print("share landing in [0.02, 0.03]:", float(((x >= 0.02) & (x < 0.03)).mean()))
print("width of that slice          :", 0.01)
print("height needed to make the area come out right:",
      round(float(((x >= 0.02) & (x < 0.03)).mean()) / 0.01, 1))
# largest value seen : 0.1
# share landing in [0.02, 0.03]: 0.100775
# width of that slice          : 0.01
# height needed to make the area come out right: 10.1

The values are spread evenly across a window only 0.1 wide. A slice of width 0.01 catches about a tenth of them, so the probability of that slice is 0.10 — an ordinary probability, safely below 1. But to make height times width equal 0.10 when the width is 0.01, the height has to be 10.

So a density of 10 is not an error and does not mean anything is 1000% likely. This is why a likelihood reported by a continuous model can be a large positive number, and why comparing likelihoods between models that predict on different scales — one in dollars, one in log-dollars — is meaningless unless you account for the change of units.

Expectation is the long-run average

The expectation \(\mathbb{E}[X]\) is the average value you would get if you could repeat the situation forever. For a die it is 3.5 — a number the die can never actually show, which is a useful reminder that an expectation is a summary and need not be a value the variable can actually take.

For variables with finite expectations, linearity does not require independence: \(\mathbb{E}[aX + bY] = a\mathbb{E}[X] + b\mathbb{E}[Y]\) holds whether or not the two have anything to do with each other. Variance behaves differently, which is the subject of the next section.

What expectations do not do is survive being pushed through a curved function. \(\mathbb{E}[f(X)]\) is generally not \(f(\mathbb{E}[X])\), and the gap can be substantial. A prediction on the transformed scale may not estimate the quantity needed on the original scale.

Three numbers are enough to see it. Take 1, 10 and 100, transform them with a logarithm, average the results, then transform back.

import numpy as np

values = np.array([1.0, 10.0, 100.0])
print("average of the values       :", values.mean())
print("average of their logs       :", np.log(values).mean().round(6))
print("that average, exponentiated :", np.exp(np.log(values).mean()).round(6))
# average of the values       : 37.0
# average of their logs       : 2.302585
# that average, exponentiated : 10.0

Averaging the values gives 37. Transforming, averaging and transforming back gives 10. Neither is a mistake — they are answers to different questions, and the round trip through the logarithm does not return you to the average you started from. Everything in the income example below is this, at scale.

Incomes are the standard example. They are skewed, so people model the logarithm, which is better behaved, and then exponentiate the prediction to get back to money.

import numpy as np

rng = np.random.default_rng(2)
income = rng.lognormal(mean=10.0, sigma=0.8, size=200_000)

log_income = np.log(income)
print(f"true average income          {income.mean():12,.0f}")
print(f"exp(average of log income)   {np.exp(log_income.mean()):12,.0f}")
print(f"median income                {np.median(income):12,.0f}")
# true average income                30,282
# exp(average of log income)         21,993
# median income                      22,060

Exponentiating the average of the logs gives 21,993 when the true average is 30,282 — low by 27%. And look at the third line: the number it produced is essentially the median, 22,060. Exponentiating the mean of the logs gives the geometric mean, which in this log-normal example coincides with the population median. On other positive distributions it does not: the geometric mean of 1, 1 and 100 is about 4.64 while the median is 1.

The median and the mean part company whenever the distribution is skewed, and income is skewed by construction. If you are forecasting total revenue, the mean is the quantity you need. In this example the back-transformed figure sits about 27% below it, and a model that fits beautifully on the log scale will keep producing that shortfall.

Variance, and why it does not simply add

Variance measures spread — how far values typically sit from their average. Its square root, the standard deviation, is in the same units as the data and is usually the one worth reporting.

For variables with finite second moments, the variance of a sum also depends on their covariance: \(\operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y) + 2\operatorname{Cov}(X, Y)\), where the last term is the covariance — a measure of whether the two move together.

import numpy as np

rng = np.random.default_rng(3)
n = 200_000
a = rng.normal(size=n)

for rho in (0.0, 0.5, 0.9, -0.9):
    b = rho * a + np.sqrt(1 - rho**2) * rng.normal(size=n)
    print(f"  corr {rho:+.1f}   Var(a) {a.var():.3f}   Var(b) {b.var():.3f}"
          f"   Var(a+b) {(a + b).var():.3f}   sum of the two {a.var() + b.var():.3f}")
# corr +0.0   Var(a) 0.998   Var(b) 1.003   Var(a+b) 2.007   sum of the two 2.002
# corr +0.5   Var(a) 0.998   Var(b) 1.001   Var(a+b) 2.997   sum of the two 1.999
# corr +0.9   Var(a) 0.998   Var(b) 0.999   Var(a+b) 3.795   sum of the two 1.997
# corr -0.9   Var(a) 0.998   Var(b) 1.001   Var(a+b) 0.199   sum of the two 1.999

Every row has the same two individual variances, near 1.0 each. Only the relationship between the variables changes. When they are unrelated the variance of the sum is 2.007, matching the naive addition. At correlation 0.9 it is 3.795 — nearly double what adding would give. At correlation -0.9 it collapses to 0.199, a tenth of the naive sum.

That last row is why averaging models with less-correlated errors can reduce variance. Averaging models reduces variance exactly to the degree their errors fail to move together; averaging ten copies of the same model reduces nothing at all, because the correlation is 1.

The first row is the one your tooling assumes. Cross-validation folds that share observations produce correlated scores — the middle rows — so the ordinary standard error across folds understates the real uncertainty. Repeated measurements per patient and multiple rows per user do the same thing.

One caution before moving on: zero covariance does not mean independent. If \(X\) is symmetric around zero and \(Y = X^2\), then \(Y\) is completely determined by \(X\) and their covariance is exactly zero. Covariance only sees straight-line association, which is why a correlation matrix is a weak screen for redundant features.

Conditional probability, and what Naive Bayes assumes

A conditional probability \(P(A \mid B)\) is the probability of \(A\) once you already know \(B\) happened — it restricts attention to the \(B\) cases and asks what fraction of those are also \(A\). Two events are independent when knowing one tells you nothing about the other, which is the same as \(P(A \cap B) = P(A)P(B)\).

There is a different notion of independence that matters more in practice. Events are conditionally independent given \(C\) when they are independent inside each group defined by \(C\), even though they may be strongly related overall. Neither kind of independence implies the other, and that gap is the whole basis of the Naive Bayes classifier.

The words “free” and “offer” appear together in email far more often than chance would allow, so they are plainly not independent. Naive Bayes never claims they are. It claims they are independent within spam and independent within non-spam — that once you know the class, one word tells you nothing more about the other.

import numpy as np

rng = np.random.default_rng(0)
n = 200_000
spam = rng.random(n) < 0.3

def word(p_spam, p_ham):                       # each word depends only on the class
    return rng.random(n) < np.where(spam, p_spam, p_ham)

w1, w2 = word(0.6, 0.1), word(0.5, 0.05)
corr = lambda a, b: np.corrcoef(a.astype(float), b.astype(float))[0, 1]

print(f"marginal corr(w1, w2)  {corr(w1, w2):+.4f}")
print(f"given spam             {corr(w1[spam], w2[spam]):+.4f}")
print(f"given not spam         {corr(w1[~spam], w2[~spam]):+.4f}")
# marginal corr(w1, w2)  +0.2782
# given spam             -0.0030
# given not spam         +0.0026

The two words were built so that each depends only on the class and on nothing else — conditional independence holds by construction. Yet measured across all the email at once they correlate at +0.2782. Split by class, the correlation vanishes: -0.0030 and +0.0026.

The apparent association was entirely the class creating it. Both words are common in spam and rare elsewhere, so they rise and fall together across the whole corpus while being unrelated inside either group. Naive Bayes is often called a crude approximation, and it is — but the assumption it makes is the second measurement, not the first.

Bayes’ theorem, by counting people

Bayes’ theorem reverses a conditional probability: it turns “how often does an ill person test positive” into “how often is a positive person ill”. Those two sound interchangeable and are not. Confusing them is what makes a 99% accurate test look trustworthy when it is not.

The formula is \(P(A \mid B) = P(B \mid A)P(A) / P(B)\). It can also be understood through expected counts, which is the route taken here. Take a test with 99% sensitivity and 99% specificity — it catches 99% of ill people and wrongly flags 1% of healthy ones — and apply it to 100,000 people where 1 in 1,000 is ill.

people = 100_000
prevalence = 0.001
sick = int(people * prevalence)
healthy = people - sick

true_positive = int(sick * 0.99)
false_positive = int(healthy * 0.01)

print(f"  {people:,} people, {prevalence:.1%} of them ill")
print(f"  ill      : {sick:6,}  ->  {true_positive:6,} test positive")
print(f"  healthy  : {healthy:6,}  ->  {false_positive:6,} test positive")
print(f"  positives total {true_positive + false_positive:,}"
      f"   of which actually ill {true_positive:,}")
print(f"  chance of being ill given a positive test: "
      f"{true_positive / (true_positive + false_positive):.4f}")
# 100,000 people, 0.1% of them ill
# ill      :    100  ->      99 test positive
# healthy  : 99,900  ->     999 test positive
# positives total 1,098   of which actually ill 99
# chance of being ill given a positive test: 0.0902

There is no probability theory in that block — it is arithmetic on head counts, and you can check every line by hand. Out of 100,000 people, 100 are ill and the test finds 99 of them. The other 99,900 are healthy, and 1% of them, 999 people, get a positive result anyway.

So 1,098 people are told they tested positive, and only 99 of them are ill. The chance that a positive result means illness is 99/1098, about 9%. Ten of every eleven alarms are false, from a test with 99% sensitivity and 99% specificity. That does not settle whether the test is worth running: screening programs are judged on the cost of the follow-up a positive triggers, the cost of the cases a negative misses, and what the alternative is.

The reason is that the 1% error is charged against a group 999 times larger than the ill group. A small error rate on a huge population easily produces more cases than a large success rate on a tiny one. That is what Bayes’ theorem encodes, and the piece practitioners drop is the size of the group the rate is applied to, which is exactly what the prior, or base rate, records.

Sweeping the base rate while holding the test fixed shows how much of the answer it controls.

sensitivity, specificity = 0.99, 0.99

for prevalence in (0.5, 0.01, 0.001):
    tp = prevalence * sensitivity
    fp = (1 - prevalence) * (1 - specificity)
    print(f"prevalence {prevalence:6.3f}"
          f"   P(disease | positive) {tp / (tp + fp):.4f}")
# prevalence  0.500   P(disease | positive) 0.9900
# prevalence  0.010   P(disease | positive) 0.5000
# prevalence  0.001   P(disease | positive) 0.0902

The test never changed across those three rows. It gives a 99% answer when half the population is affected, a coin flip at 1%, and 9% at 0.1%. Nothing about the classifier is broken; the base rate is doing all the work.

This is the same arithmetic that makes fraud detection, rare-disease screening and predictive maintenance hard, and it is why accuracy alone can be misleading on imbalanced data — a model that always answers “no” is 99.9% accurate at 0.1% prevalence. It is also why a model trained where positives are common and deployed where they are rare returns miscalibrated probabilities even though nothing else changed.

Which distribution models what

A distribution is a claim about how your data was produced. Choosing one is not decoration — a distributional model can motivate a corresponding loss function, because most of the common regression and classification losses are the negative log-likelihood of some distribution. Hinge loss and 0-1 loss are the notable exceptions. Use the table below as a reference for where each distribution tends to be used; you do not need to memorize every row before continuing, and the ones this series needs get their own treatment later.

Note the pair in the middle. Sums of many small effects tend toward a Normal; products of many small effects tend toward a log-normal, because taking logs turns a product into a sum. That is why incomes, session lengths and response times are right-skewed, and why the log trap from earlier keeps appearing.

The practical use is choosing a loss. Squared error is the maximum-likelihood loss under a Normal with constant variance, though it remains a reasonable choice for predicting a conditional mean without that assumption. If the target is a count, Poisson deviance is a natural candidate — but counts are often overdispersed or zero-inflated, and Poisson assumes neither. If the target is strictly positive and right-skewed, modeling its logarithm and assuming normality on that scale gives a log-normal model, and it changes what “average error” means.

Why sample averages are trustworthy — and when they are not

Return to the coin flips at the start. The proportion settled down as the count grew, and two results explain why. The law of large numbers says a sample mean converges to the true mean. For independent and identically distributed draws with finite variance, \(\operatorname{SD}(\bar{X}) = \sigma/\sqrt{n}\) exactly, at every \(n\); no theorem about large samples is needed for that. The central limit theorem is about shape: once standardized, the distribution of the sample mean approaches a Normal whatever the shape of the thing you sampled.

This approximation underlies many of the confidence intervals in common use. The experiment below draws \(n\) values from an exponential distribution — strongly skewed, nothing like a bell curve — records their average, and repeats that 20,000 times. Two different distributions are in play, and it is worth keeping them apart: the exponential is the distribution of the individual values, while the 20,000 averages have a distribution of their own, and that second one is what the numbers below describe. The block uses SciPy for its skewness function, which is the first thing in this series that NumPy alone will not give you.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
for n in (1, 5, 30, 200):
    means = rng.exponential(1.0, size=(20_000, n)).mean(1)
    print(f"n={n:4d}  sd {means.std():.4f}  theory {1 / np.sqrt(n):.4f}"
          f"  skew {stats.skew(means):+.3f}")
# n=   1  sd 0.9956  theory 1.0000  skew +1.962
# n=   5  sd 0.4463  theory 0.4472  skew +0.883
# n=  30  sd 0.1819  theory 0.1826  skew +0.381
# n= 200  sd 0.0709  theory 0.0707  skew +0.161

Two things happen at different speeds, and the difference is useful. The standard deviation matches the predicted \(\sigma/\sqrt{n}\) almost exactly from \(n = 5\) onward — 0.4463 against 0.4472. The skewness, which measures how lopsided the shape still is, falls much more slowly: 1.96, 0.88, 0.38, 0.16.

So the spread of the sample mean matches \(\sigma/\sqrt{n}\) from the start, while the shape takes much longer to become Normal. Note what this experiment does and does not show: it measures the spread across 20,000 repeated samples, so it says nothing about how reliably a standard error estimated from one sample of 30 would perform.

The theorem needs two things, and both fail regularly in practice. It needs finite variance, and some heavy-tailed quantities do not have it at all — for those the theorem simply does not apply. For others the variance is finite but the distribution is strongly skewed or heavy in the tails, and it is that shape — not the size of the variance — that makes the Normal approximation slow to arrive. Rescaling a distribution changes \(\sigma\) without changing the standardized shape at all. It also needs independent draws. Under dependence the independent-sample standard error formula no longer applies, and depending on the dependence structure the rate itself can change; the exercises measure one such case.

Eight ways this goes wrong

The opening said that most probability errors in practice look like modeling errors. Here is the list, in the order the article met them.

When interpreting a model’s output, ask what the reported quantity estimates and which assumptions it relies on. A rate needs a population attached to it. A density needs an interval. Applying a nonlinear transformation to a mean generally differs from averaging the transformed values. A standard error needs independence, or a correction for its absence.

Exercises

1. The screening trade-off. For a 99%/99% test, report the positive predictive value, the negative predictive value, and the number of false alarms per true case, at prevalences 0.5, 0.01, and 0.001. Explain which of the three numbers stays reassuring.

Expected observation: a negative predictive value that looks excellent at every prevalence, including the one where a positive result leaves substantial uncertainty.

sens, spec = 0.99, 0.99
for prev in (0.5, 0.01, 0.001):
    tp = prev * sens
    fp = (1 - prev) * (1 - spec)
    ppv = tp / (tp + fp)
    npv = ((1 - prev) * spec) / ((1 - prev) * spec + prev * (1 - sens))
    print(f"prev {prev:6.3f}  PPV {ppv:.4f}  NPV {npv:.6f}"
          f"  false alarms per true case {fp / tp:.1f}")
# prev  0.500  PPV 0.9900  NPV 0.990000  false alarms per true case 0.0
# prev  0.010  PPV 0.5000  NPV 0.999898  false alarms per true case 1.0
# prev  0.001  PPV 0.0902  NPV 0.999990  false alarms per true case 10.1

2. Explaining away. Generate two independent causes and condition on a common effect. Report the correlation between the causes before and after conditioning, and name the structure.

Expected observation: two variables that are independent by construction and strongly correlated once you select on their consequence.

import numpy as np

rng = np.random.default_rng(0)
n = 500_000
talent = rng.normal(size=n)                    # independent by construction
luck = rng.normal(size=n)
admitted = (talent + luck) > 2.0               # a common effect: the collider

corr = lambda a, b: np.corrcoef(a, b)[0, 1]
print(f"marginal            {corr(talent, luck):+.4f}")
print(f"given admitted      {corr(talent[admitted], luck[admitted]):+.4f}")
print(f"given not admitted  {corr(talent[~admitted], luck[~admitted]):+.4f}")
print(f"admitted fraction   {admitted.mean():.4f}")
# marginal            -0.0010
# given admitted      -0.7285
# given not admitted  -0.1444
# admitted fraction   0.0787

3. Effective sample size under correlated observations. Measure the variance of a sample mean when the observations share a common component, for equicorrelation 0, 0.3, and 0.9 at \(n = 25\). Report the effective sample size in each case.

Expected observation: twenty-five observations that behave like one, from a correlation people routinely ignore.

import numpy as np

rng = np.random.default_rng(1)
n, reps = 25, 20_000
for rho in (0.0, 0.3, 0.9):
    shared = np.sqrt(rho) * rng.normal(size=(reps, 1))
    own = np.sqrt(1 - rho) * rng.normal(size=(reps, n))
    empirical = (shared + own).mean(1).var()
    theory = (1 + (n - 1) * rho) / n
    print(f"rho {rho:.1f}  var(mean) {empirical:.4f}  theory {theory:.4f}"
          f"  effective n {1 / empirical:.1f}")
# rho 0.0  var(mean) 0.0404  theory 0.0400  effective n 24.8
# rho 0.3  var(mean) 0.3286  theory 0.3280  effective n 3.0
# rho 0.9  var(mean) 0.8881  theory 0.9040  effective n 1.1