Probability for Machine Learning

A classifier may report a probability, a ranking model a score, and a forecast an estimated value or interval. Interpreting these outputs requires knowing what quantity each represents. This article develops probability through repeated trials, transformations, and conditional groups. Run the body code blocks in order; they use NumPy and, for one experiment, SciPy.

A test’s overall accuracy does not give the probability of illness after a positive result. Repeated measurements from one person need not carry as much information as measurements from separate people. These are different questions, and the examples below show the calculations needed for each.

Estimating a probability from a proportion

For independent repetitions of the same experiment with a fixed event probability, the observed fraction estimates that probability. A fair coin has probability 0.5 of heads; as the number of independent flips grows, the observed fraction converges to 0.5. This long-run interpretation does not promise a particular proportion in any finite run.

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

This fair-coin simulation produced three heads in ten flips. Such a result is possible under the stated probability; it does not by itself establish that a real coin is fair or biased. 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 value in its discrete support can carry 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 model spreads values uniformly across a window of width 0.1. Its density is exactly 10; the printed 10.1 estimates that density from the finite-sample fraction 0.100775. 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.

A density of 10 does not mean an event is 1000% likely. Densities, and likelihoods built from them, can exceed 1. Comparing a density for income with one for log income requires expressing both on the same scale, including the transformation’s Jacobian factor; their raw values are not directly comparable.

Expectation is the long-run average

When it is finite, the expectation \(\mathbb{E}[X]\) is the probability-weighted average of possible values. For a fair six-sided die it is \((1+2+3+4+5+6)/6=3.5\), even though the die cannot show 3.5. Under independent, identically distributed repetition with a finite absolute mean, the sample average converges to this expectation. Some distributions do not have a finite expectation.

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.

A nonlinear transformation generally changes which average is being calculated. \(\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.

For positive, right-skewed quantities such as income, modeling logarithms can make some patterns easier to represent. The quantity obtained by exponentiating a fitted log-scale mean still needs to be distinguished from an original-scale mean.

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"sample 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}")
# sample 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 sample average is 30,282 — low by 27%. And look at the third line: the number it produced is essentially the median, 22,060. These calculations use the generated sample. At the population level, exponentiating the expected log income gives the geometric mean, which equals the median for this log-normal distribution. On other positive distributions it does not: the geometric mean of 1, 1 and 100 is about 4.64 while the median is 1.

For this log-normal population, \(\log X\) has mean 10 and variance \(0.8^2\). The population geometric mean and median are \(e^{10}\), while the arithmetic mean is \(e^{10+0.8^2/2}\). Exponentiating only the log mean therefore misses the variance adjustment. The 27% gap above is a sample comparison illustrating this difference. Expected total revenue is the sum of expected individual revenues; predicting those means requires the appropriate original-scale quantities, including any conditioning features.

Variance, and why it does not simply add

Variance is the expected squared deviation from the mean: \(\operatorname{Var}(X)=\mathbb{E}[(X-\mathbb{E}[X])^2]\), when finite. For equally likely values 1 and 3, the mean is 2, the variance is \((1+1)/2=1\), and the standard deviation is 1. Standard deviation is the square root of variance, restoring the original measurement units.

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 \(\operatorname{Cov}(X,Y)\) 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.

For a fixed collection of finite-variance prediction errors, averaging gives variance \(M^{-2}\sum_{i,j}\operatorname{Cov}(e_i,e_j)\). Lower covariance helps when the individual variances are held fixed. For example, two errors each with variance 1 give an average-error variance of \((1+\rho)/2\). Averaging identical copies changes nothing. Models with very different error variances also need their individual contributions considered; correlation alone does not determine whether an ensemble improves on a particular member.

The independent-sample standard error omits covariance terms. In ordinary K-fold cross-validation the validation subsets are disjoint, but training sets overlap, so fold scores need not be independent. Treating them as independent replicates can understate uncertainty; neither the direction nor the size of the error follows from the number of folds alone. Repeated measurements within patients or users require their own dependence model.

One caution before moving on: zero covariance does not mean independent. For a standard Normal \(X\) 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\) given an event \(B\) with \(P(B)>0\) — 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. Naive Bayes uses conditional independence of features given the class.

Words such as “free” and “offer” can appear together more often because both are common in spam. A Bernoulli Naive Bayes model assumes their presence indicators are independent within each class. The following simulation makes that assumption true by construction; it does not establish it for real email.

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. Within each class, the measured correlations are near zero: -0.0030 and +0.0026. Independence comes from the generating process here, not from a near-zero sample correlation alone.

In this simulation, class membership creates the marginal association: both indicators are more likely in spam. For real text, within-class dependence can remain. Naive Bayes approximates the within-class joint distribution by a product of individual feature distributions.

Bayes’ theorem, by counting people

Bayes’ theorem relates two directions of conditioning. A test’s sensitivity describes positive results among ill people; its positive predictive value describes illness among people with positive results. The population’s base rate connects those two quantities.

For events with positive conditioning probabilities, 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

The block expresses the probability calculation as expected counts that can be checked by hand. Among 100,000 people with the stated prevalence, 100 are expected to be ill and 99 of those to test positive. Among the other 99,900, the expected number of false positives is 999. Actual screening counts would fluctuate; the integer expectations happen to be exact for these chosen rates.

These expectations give 1,098 positive results, including 99 true positives. 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

Sensitivity and specificity are held fixed across these rows. The positive predictive value changes from 99% to 50% to about 9% as prevalence falls. This isolates the effect of the base rate under the assumed fixed test properties.

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. If class proportions change while class-conditional feature distributions stay fixed, posterior probabilities generally change too. Unadjusted probabilities calibrated for the old population may then be miscalibrated in the new one; the effect depends on how informative the features are.

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.

Distribution	Models	Where it appears
Bernoulli / Binomial	one binary outcome / success counts in fixed trials	binary classification, conversion rates
Categorical / Multinomial	one of \(k\) outcomes / counts across those outcomes	multiclass output, bag-of-words
Poisson	event counts in fixed exposure	arrivals, defects, Poisson regression
Normal	continuous values; also approximates suitably standardized sums	residuals, the CLT, Gaussian models
Log-normal	positive values whose logarithms are Normal	income, session length, response times
Exponential	waiting time at constant hazard	time to failure, survival baselines
Beta	a probability that is itself uncertain	priors on rates, calibration curves

Under suitable assumptions, including independence and finite variance in the standard CLT, standardized sums approach a Normal distribution. Positive products can be related to sums by taking logs, but an approximate log-normal model requires suitable behavior of those log terms. This is a possible modeling rationale, not proof that income or response-time data follow a log-normal distribution.

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

For independent, identically distributed draws with a finite absolute mean, the law of large numbers says the sample mean converges to the population mean. If the variance is also finite, \(\operatorname{SD}(\bar X)=\sigma/\sqrt n\) exactly for every \(n\). If \(0<\sigma^2<\infty\), the standard central limit theorem says \(\sqrt n(\bar X-\mu)/\sigma\) converges in distribution to a standard Normal. The first result concerns convergence of the average, the exact formula concerns spread, and the CLT concerns the standardized shape.

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 supplies a ready-made skewness function; the statistic could also be computed from NumPy operations.

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

The simulated standard deviations are close to their exact theoretical values at all four sample sizes. The remaining differences come from estimating spread with 20,000 repetitions. Skewness, which measures asymmetry, decreases from 1.96 to 0.88, 0.38, and 0.16 as the sample size increases.

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 i.i.d. finite-variance version of the CLT does not cover infinite-variance distributions. Among finite-variance distributions, strong skewness or heavy tails can make the Normal approximation poor at practical sample sizes. Rescaling changes \(\sigma\) but not standardized shape. Dependent observations require other assumptions and a variance calculation that includes covariance; some forms of weak dependence still permit a CLT and a \(\sqrt n\) rate.

Eight ways this goes wrong

The examples distinguish the quantities below. Each interpretation depends on what was measured and how the observations were generated.

The mistake	What is actually true
Reading a rate off a few dozen trials	Ten flips gave 30% heads on a fair coin. How many trials you need depends on the accuracy you want.
Treating a density above 1 as a bug	A density is a height. An interval probability is the area under it, which height times width only approximates.
Exponentiating a log-scale mean as an estimate of the original-scale mean	Exponentiating the mean of logs gives the geometric mean; for conditional predictions, the corresponding conditional log mean is required. It equals the median under a log-normal model. The sample example was 27% below its arithmetic mean.
Adding variances of things that move together	The covariance term is real: at correlation 0.9 the sum was almost double.
Screening features by correlation alone	Zero covariance does not imply independence; the symmetric finite-moment example above uses \(X\) and \(X^2\).
Saying Naive Bayes assumes words are unrelated	It assumes they are unrelated within each class, which is a within-class assumption rather than a marginal one.
Quoting a test’s accuracy as its reliability	At 0.1% prevalence a 99%/99% test is right 9% of the time it fires.
Putting error bars on correlated folds	The independent-sample formula assumes independence; the correction depends on the dependence structure.

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.