Probability, and reading an experiment honestly

Most of the probability worth knowing here serves one practical question: my new version scored better, should I believe it? This section covers what you need to answer that, distributions and expectation, conditional probability and Bayes, and how much a measurement drifts purely because of which cases happened to be sampled. It is the difference between reporting a result and reporting a fluctuation.

What you actually need#

The probability that matters in practice is not the theory. It is being able to answer one question: my new version scored better, is that real?

A random variable, in plain words#

A random variable is neither a variable nor random in the everyday sense. It is a rule that turns something you cannot predict into a number. Roll a die and read the top face: that is a random variable. Serve a web request and record how many milliseconds it took: that is a random variable. Send a prompt to a model, have a grader score the answer out of two, and the score is a random variable. The rule never changes. What changes is which outcome the world hands you this time.

A distribution is the whole picture of one such variable: every value it can take, and how often. Write the values in one column and their long-run shares in another, and you have written the distribution down. The shares add up to one, because some outcome always happens. Nothing else is required of it, and the word carries no more mystery than that.

Make it concrete. A grader awards 0, 1 or 2 points per answer, and across everything your system produces, 20 per cent score 0, 50 per cent score 1 and 30 per cent score 2. Those three numbers are the distribution of the score. Every other quantity in this section — the average score, how much it wobbles, how far a 50-case eval can drift from the truth — is arithmetic on those three numbers.

import numpy as np
rng = np.random.default_rng(0)

values = np.array([0, 1, 2])          # what the grader can award
shares = np.array([0.2, 0.5, 0.3])    # how often each score happens; they sum to 1

draws = rng.choice(values, size=200_000, p=shares)
print(np.bincount(draws) / len(draws))
# [0.200215 0.50064  0.299145]

Expectation and variance, by hand then by counting#

The expectation — also called the mean, or the expected value — is the long-run average: each value multiplied by its share, added up. Here that is 0 times 0.2, plus 1 times 0.5, plus 2 times 0.3, which is 0 plus 0.5 plus 0.6, so 1.1. Notice what 1.1 is not. It is not a score any single answer can receive. An expectation does not have to be an attainable value, which is why an average of anything counted needs reading with care.

The variance measures spread: it is the expectation of the squared distance from the mean. The distance is squared so that two below the mean counts the same as two above, and so that one large miss counts for more than several small ones. By hand: a 0 is 1.1 away, squared 1.21, weighted by its share 0.2, contributing 0.242. A 1 is 0.1 away, squared 0.01, weighted 0.5, contributing 0.005. A 2 is 0.9 away, squared 0.81, weighted 0.3, contributing 0.243. Adding the three gives 0.49.

Squaring left the units wrong — 0.49 is in points squared — so the number usually quoted is its square root, 0.7, the standard deviation. That one is in points and can be compared directly with the mean of 1.1. If you remember one thing from the pair, it is that the mean tells you where a quantity sits and the standard deviation tells you how much a single observation of it is worth trusting.

mean = (values * shares).sum()
variance = (shares * (values - mean) ** 2).sum()
print(mean, variance, variance ** 0.5)
# 1.1 0.49 0.7

print(draws.mean(), draws.var())
# 1.09893 0.4895728551000001

Distributions#

• Normal: heights, measurement error, averages of many things. Symmetric around a mean.
• Binomial: counts of successes in trials. Pass rates on an eval set are binomial.
• Long-tailed: latency, request sizes, income. The mean is misleading and the median plus a high percentile is honest.

Conditional probability#

The probability of something given that something else is true. The definition is one line: the probability of A given B is the probability that A and B both happen, divided by the probability of B. In words, and this is the version to actually hold in your head: throw away every case where B did not happen, then ask what share of the cases left over had A. Conditioning is a filter, not a formula.

The place it goes wrong is that two facts which sound the same are different filters. Flip two fair coins. Four outcomes are equally likely: heads-heads, heads-tails, tails-heads, tails-tails. Told that at least one coin came up heads, you throw away tails-tails and three outcomes survive, of which one is two heads, so the answer is one in three. Told instead that the first coin came up heads, you keep heads-heads and heads-tails, so the answer is one in two. Same coins, same question, different filter, and most people give one half to both.

flips = rng.integers(0, 2, size=(200_000, 2))     # 1 is heads
both = flips.sum(axis=1) == 2
at_least_one = flips.sum(axis=1) >= 1
first_is_head = flips[:, 0] == 1

print(both[at_least_one].mean())     # 0.3345792524485309   -> one third
print(both[first_is_head].mean())    # 0.5008677438659486   -> one half

Independence is the case where the filter changes nothing. Two events are independent when the probability of both is the product of the two separate probabilities, which is the same as saying that learning one of them leaves your belief about the other exactly where it was. The first coin and the second coin are independent. The first coin and 'at least one head' are not, because knowing the first is heads guarantees the second fact.

second_is_head = flips[:, 1] == 1

print(first_is_head.mean() * second_is_head.mean(),
      (first_is_head & second_is_head).mean())
# 0.2507678055 0.251085     -> the product matches, so independent

print(first_is_head.mean() * at_least_one.mean(),
      (first_is_head & at_least_one).mean())
# 0.376200585 0.5013        -> nowhere near, so not independent

One more, because interviewers reach for it and because it breaks the most common wrong mental model. Uncorrelated does not mean independent. Correlation only detects a straight-line relationship. Take x spread evenly between minus one and one, and let y be x squared. Knowing x tells you y exactly — there is no relationship more dependent than that — and the correlation is zero.

x = rng.uniform(-1, 1, 200_000)
y = x ** 2

print(np.corrcoef(x, y)[0, 1])       # 0.0022397797187538398  -> effectively zero
print((y > 0.25).mean())             # 0.499755   unconditionally, a coin flip
print((y[x > 0.5] > 0.25).mean())    # 1.0        given x > 0.5, a certainty

Bayes' theorem#

Bayes' theorem reverses the direction of a conditional probability. You usually know how the evidence behaves given the state of the world — how often the test fires when someone is ill — and you want the state of the world given the evidence, which is the opposite. The rule is: the probability of A given B equals the probability of B given A, times the probability of A, divided by the probability of B. Read as a procedure it is even shorter. Take what you believed before the evidence, multiply by how likely that evidence would be if the belief were true, do the same for every competing belief, and rescale so the results add to one.

The example every interviewer uses is the medical test, so work it properly once. A test is 99 per cent accurate for a condition that affects 1 person in 1000. Spell out what 99 per cent accurate means, because the phrase hides two numbers: the test catches 99 per cent of the people who have the condition, and it correctly clears 99 per cent of the people who do not, so it wrongly flags 1 per cent of them. You test positive. What is the chance you have it?

A test is 99% accurate. The condition affects 1 in 1000.

Out of 100,000 people:
  100 have it        ->  99 test positive        (99% of 100)
  99,900 do not      -> 999 test positive        (1% of 99,900)

  1,098 positives in total, of which 99 are real.
  99 / 1,098 = 0.090

A positive result means about a 9% chance of actually having it.

The symbolic version gives the identical number. The probability of being ill given a positive is 0.99 times 0.001, divided by that same product plus 0.01 times 0.999. The numerator is 0.00099 and the denominator is 0.00099 plus 0.00999, which is 0.01098, so the answer is 0.0902. It is the counts argument with the 100,000 divided out.

def posterior(prior, sensitivity, specificity):
    hit = sensitivity * prior                       # ill, and flagged
    false_alarm = (1 - specificity) * (1 - prior)   # healthy, and flagged anyway
    return hit / (hit + false_alarm)

for prior in (0.001, 0.01, 0.10):
    print(prior, round(posterior(prior, 0.99, 0.99), 4))
# 0.001 0.0902
# 0.01 0.5
# 0.1 0.9167

Check it by simulation rather than believing the algebra, which takes one more block. Two million people, one in a thousand ill, everyone tested.

n = 2_000_000
has_it = rng.random(n) < 0.001
tests_positive = np.where(has_it, rng.random(n) < 0.99, rng.random(n) < 0.01)

print(tests_positive.sum(), has_it[tests_positive].mean())
# 21848 0.09044306114976199

There is a form of the rule worth memorising because it makes repeated evidence trivial. Work in odds rather than probabilities: prior odds times the likelihood ratio gives posterior odds, and the likelihood ratio for a positive is sensitivity divided by the false positive rate, here 0.99 over 0.01, which is 99. Prior odds of 1 to 999, multiplied by 99, gives 0.099 to 1, which is the 9 per cent again. Test a second time with an independent test and multiply by 99 once more: odds of 9.81 to 1, or 91 per cent.

odds = 0.001 / 0.999            # prior odds: 1 case per 999 healthy people
for test in (1, 2, 3):
    odds = odds * 99            # each independent positive multiplies by the ratio
    print(test, round(odds, 3), round(odds / (1 + odds), 4))
# 1 0.099 0.0902
# 2 9.811 0.9075
# 3 971.27 0.999

None of this is medical trivia. It is the arithmetic of every rare-event classifier you will ship. A fraud model at a 1 per cent base rate with a 1 per cent false positive rate sends an analyst one real case for every one false alarm, and that is the good version; at a 0.1 per cent base rate the analyst sees ten false alarms per real case and stops reading the queue. The quantity that decides whether an alerting system is usable is precision at the operating point, and precision is this calculation.

Sampling variation#

Measuring on 50 examples gives a number with real uncertainty in it. Two systems scoring 82 percent and 86 percent on 50 examples may well be identical, and the difference is noise.

import math

def margin(p, n):
    """Rough 95% margin of error for a proportion."""
    return 1.96 * math.sqrt(p * (1 - p) / n)

print(margin(0.85, 50), margin(0.85, 500), margin(0.85, 5000))
# 0.09897514839594836 0.031298690068435775 0.009897514839594836

That formula is worth understanding rather than memorising, and it comes straight out of the previous two sections. A pass on one eval case is a random variable worth 1 or 0; if the true pass rate is p its variance is p times one minus p. Average n of them and the variance divides by n, so the standard deviation of the average is the square root of p times one minus p over n. The 1.96 turns a standard deviation into a 95 per cent range, for reasons the central limit theorem supplies later in this topic.

Now put a number on the mistake in the headline. Take two systems that are genuinely identical, both passing 85 per cent of the time, and score each on its own 50 cases. How often does the pair differ by four points or more — enough for someone to write that one is better?

a = rng.binomial(50, 0.85, 200_000) / 50      # 200,000 evals of system A
b = rng.binomial(50, 0.85, 200_000) / 50      # and of an identical system B

print((np.abs(a - b) >= 0.04).mean())    # 0.624755
print((np.abs(a - b) >= 0.10).mean())    # 0.13899

This is the most useful thing on this page. Before claiming an improvement, check whether your eval set is large enough for the difference to mean anything.

Paired comparison is stronger#

The 62 per cent above is pessimistic on purpose: it assumed the two systems were scored on different samples. In practice you score both on the same examples, and that is a much better experiment, because the hard cases are hard for both and the easy cases are easy for both. All the shared difficulty cancels. So throw away the aggregate scores and look only at the cases where the two systems disagree.

Say A scores 41 out of 50 and B scores 43, the familiar 82 against 86. On 45 of the 50 cases they gave the same verdict. Of the five where they differed, B won three and A won two. If the two systems were equally good, each disagreement would go either way like a coin flip, so the question is whether 3 out of 5 is a surprising number of coin flips. It is not: the two-sided probability of a split at least that lopsided is 1.0, meaning every possible outcome is at least as extreme as this one.

from scipy import stats

print(stats.binomtest(3, 5, 0.5).pvalue)   # 1.0      the actual result: 3 wins of 5
print(stats.binomtest(5, 5, 0.5).pvalue)   # 0.0625   a clean sweep of all 5
print(stats.binomtest(6, 6, 0.5).pvalue)   # 0.03125  six disagreements, all one way

That reframes what to ask for. The number that decides whether a comparison can conclude anything is not the size of the eval set, it is the number of cases the two systems actually disagree on. Two prompts that differ on five of fifty cases cannot be separated no matter how the five fall. Either find more cases where they behave differently — which usually means harder cases — or accept that the change is too small for this eval to see.