InterviewPrepKit

Home / Learn / Math

01 — Probability

Machine learning runs on a small core of probability:

Every result below is paired with the assumption it rests on and the specific way that assumption is violated in production. A formula applied outside its assumptions produces confident wrong answers.

By the end you should be able to update a belief from a base rate, derive the cross-entropy loss from maximum likelihood, state precisely what the central limit theorem does and does not promise, and solve the four puzzles that come up in interviews.

Probability sits underneath every loss function. A loss function is the formula that scores how wrong a prediction is. Training is the search for the parameter settings that make that score smallest.

Two concrete examples of “underneath”:

An interviewer asking a probability question is rarely asking about dice. They are asking whether you can name the assumption your model relies on.

What goes in, what comes out, and the words used throughout

Two things need pinning down first: what a probability model requires and what it returns, and the vocabulary the later sections use without redefining.

The contract

A probability model is a machine for turning a description of an uncertain situation into numbers.

The block below is that contract. The left column is what you have to supply. The right column is what you are then allowed to ask for. Every rule in this chapter computes one of the four things on the right.

INPUT                                    OUTPUT
  a list of everything that could         P(A), a number between 0 and 1,
    happen  (the sample space)              for any event A you can describe
  a rule assigning a probability to       E[X], one number: the long-run
    each basic outcome                      average of a quantity X
  (optionally) something you have         Var(X), one number: how much X
    observed                                spreads around that average
                                          P(H | E), an updated probability
                                            once evidence E is in hand

So: a description of what can happen goes in, and a number in [0, 1] comes out for every question you can pose. Read [0, 1] aloud as “the closed interval from zero to one” — the square brackets mean zero and one are both allowed values.

Everything in this chapter is either a rule for computing one of those outputs, or a warning about when the rule does not apply.

Six words that carry the whole chapter

Read these once now. They are used without further comment everywhere below.

Discrete, continuous, and i.i.d.

Two shapes of distribution appear throughout, and they differ in how you get a probability out of them.

A discrete random variable takes separated values you can count: 0, 1, 2, … purchases. Its probability mass function (PMF) gives the probability of each value outright — P(X = 2) = 0.18 is a direct lookup.

A continuous random variable takes any value in a range: a waiting time of 3.417 seconds. Here any single exact value has probability zero, so there is nothing to look up. Instead a probability density function (PDF) gives probability per unit of the axis, and you integrate it over an interval to get a probability. You ask for P(3 < T < 4), never for P(T = 3.417).

Both shapes share a cumulative distribution function (CDF), written F(t) = P(X <= t) and read “the probability that X is at most t”. It is a curve that starts at 0, never decreases, and ends at 1.

One more piece of vocabulary, and it is the most important one in the chapter. i.i.d. stands for independent and identically distributed: each draw comes from the same distribution, and no draw tells you anything about another. It is the single assumption most often stated and least often checked.

Notation, and how to say it out loud

You do not need to memorize this. Come back to it whenever a symbol below is unfamiliar. Each line is one symbol and the English sentence it stands for.

NOTATION              READ ALOUD AS
P(A)                  "the probability of A"
P(not A)              "the probability that A does not happen"
P(A and B)            "the probability that A and B both happen"
P(A or B)             "the probability that A happens, or B happens, or both"
P(A | B)              "the probability of A given B" — B is known to have happened
X ~ D                 "X is distributed as D" — X is a draw from distribution D
E[X]                  "the expectation of X" — its probability-weighted average
Var(X), sd(X)         "the variance of X", "the standard deviation of X"
Cov(X, Y)             "the covariance of X and Y"
Corr(X, Y)            "the correlation of X and Y"
sum_i a_i             "the sum over i of a-sub-i" — add the terms up
prod_i a_i            "the product over i of a-sub-i" — multiply the terms together
C(n, k)               "n choose k" — how many k-element subsets n things have
1[y = k]              "the indicator that y equals k" — 1 when true, 0 when false
argmax_t f(t)         "the argument that maximizes f" — the t at which f is largest
log, exp              natural logarithm, and its inverse e^x, with e = 2.71828...
theta                 "theta", the Greek letter conventionally used for parameters
sigma, mu, rho        "sigma", "mu", "rho" — standard deviation, mean, correlation

1. The setup, in the only form that matters

With the vocabulary fixed, the machinery is small. A probability model is three things:

  1. A sample space S listing everything that could happen.
  2. A set of events, each event being a subset of S.
  3. A probability measure P — a function that hands every event a number.

P is not free to hand out any numbers it likes. It obeys exactly three rules, called the axioms:

Everything else in probability is bookkeeping derived from those three. Two derived rules do most of the work:

inclusion-exclusion:   P(A or B) = P(A) + P(B) - P(A and B)
complement:            P(A) = 1 - P(not A)

The first line repairs the additivity axiom when events overlap. Adding P(A) and P(B) counts the outcomes lying in both twice, so you subtract that shared piece once.

The second line follows from the last two axioms together: A and not A are disjoint and between them cover S, so their probabilities add to P(S) = 1.

The complement rule is the move you will use most

“At least one” is almost always harder to count than “none.” Counting “at least one” means summing over exactly-one, exactly-two, exactly-three, and so on. Counting “none” is a single product.

So flip it: P(at least one) = 1 - P(none).

Three standard problems are all that one trick:

What has to be true, and how it gets misused. Additivity requires the events to be disjoint; inclusion-exclusion is what you owe when they are not. The classic error is quoting a union of overlapping risks as a sum — “a 3% chance of a cache miss plus a 4% chance of a timeout is 7% of failing” — which double-counts every request that does both and can even produce a probability above 1.

The complement rule requires that not A be taken inside the same sample space you started with. The complement of “at least one match” is “no matches at all”, not “exactly one match”. The complement of “the model beats baseline” includes ties.

2. Conditional probability and independence

The axioms describe events in isolation. Probability becomes useful when events inform each other, when knowing one changes the odds of another.

Conditional probability is the probability of one event once you know another has happened:

P(A | B) = P(A and B) / P(B)          defined when P(B) > 0

Read P(A | B) aloud as “the probability of A given B”. The vertical bar means “given” — everything to its right is known to have happened.

Read the formula itself as a renormalization. Discard every outcome outside B. What is left no longer sums to 1, because you threw some probability away. Dividing by P(B) rescales it back to 1.

Here is that on a single die. Let A be “the roll is a 2” and B be “the roll is even”. Then P(A and B) = 1/6 and P(B) = 3/6, so P(A | B) = (1/6) / (3/6) = 1/3. Once you know the roll is even, only three outcomes survive and 2 is one of them.

The chain rule

Multiply both sides of the definition by P(B) and the division disappears. Running the same argument with the roles of A and B swapped gives the second form on the first line:

P(A and B) = P(A | B) · P(B) = P(B | A) · P(A)
P(x_1, ..., x_n) = P(x_1) · P(x_2 | x_1) · P(x_3 | x_1, x_2) · ...

The second line generalizes the first from 2 events to n, by applying it repeatedly. The probability of a whole sequence is: the probability of the first item, times the probability of the second given the first, times the probability of the third given the first two, and so on.

That is exactly what an autoregressive language model computes. Autoregressive means it emits one token at a time — a token being a chunk of text roughly the size of a short word — with each token conditioned on everything it has already emitted.

The model’s joint probability over a token sequence is factorized left to right, one conditional per position. There is no approximation anywhere in it. It is the chain rule applied n times.

Independence

Independence is the statement P(A and B) = P(A)·P(B). Equivalently, P(A | B) = P(A): learning B tells you nothing about A.

Two rolls of a die are independent — P(both 6) = (1/6)(1/6) = 1/36. “Roll is a 2” and “roll is even” from the example above are not independent, because P(A | B) = 1/3 while P(A) = 1/6. Knowing B doubled the probability.

What has to be true, and how it gets misused. Conditioning needs P(B) > 0. You cannot condition on an event of probability zero without extra machinery, which is why “given that the latency was exactly 200.000 ms” is not a well-posed condition for a continuous variable and has to be widened to an interval.

The chain rule, by contrast, assumes nothing at all. It is an identity, true for any dependence structure whatsoever, which is why a language model factorizing a sentence left to right loses no generality.

Independence is the assumption people supply for free and almost never earn: two requests from the same user, two rows from the same session, two trees fit on the same data. The classic misuse is multiplying probabilities — of failures, of features, of test results — as though the factors were independent. That understates the probability of a joint failure by however much the events actually move together.

Conditional independence is not independence, in either direction

X and Y are conditionally independent given Z when P(X, Y | Z) = P(X | Z)·P(Y | Z). In words: once you know Z, learning X tells you nothing further about Y.

Candidates routinely treat this as a weaker or stronger form of plain independence. It is neither. The two conditions are logically unrelated, and there is a clean counterexample in each direction.

The diagram below shows those two counterexamples as causal pictures. Arrows point from cause to effect. Note that the two boxes have their arrows pointing opposite ways relative to Z — that is the whole difference.

flowchart LR
    subgraph FORK["Fork -- dependent, conditionally INdependent"]
        Z1(("Z<br/>true skill")) --> X1(("X<br/>test 1"))
        Z1 --> Y1(("Y<br/>test 2"))
    end
    subgraph COLL["Collider -- independent, conditionally DEpendent"]
        X2(("X<br/>coin A")) --> Z2(("Z<br/>XOR"))
        Y2(("Y<br/>coin B")) --> Z2
    end

    style Z1 fill:#2d6a4f,color:#fff
    style Z2 fill:#9d0208,color:#fff

On the left, one cause feeds two effects. That shape is called a fork. On the right, two causes feed one effect; that shape is called a collider, because the two arrows collide on it.

The fork: dependent, but conditionally independent

A candidate’s scores on test 1 and test 2 correlate, because both are driven by an unobserved true skill Z. Unobserved quantities like Z are called latent.

Now condition on Z — meaning, look only at candidates of equal skill. The correlation vanishes. All the shared information ran through Z, and there is none left once Z is held fixed.

This is the structure Naive Bayes assumes, and the next subsection takes it apart.

The collider: independent, but conditionally dependent

Take two fair coins X, Y with values in {0,1}, and let Z = X xor Y — the exclusive-or, which is 1 when exactly one of the two coins is 1.

The coins are independent by construction. But given Z = 0, X determines Y exactly: the only way the exclusive-or is 0 is for the two coins to agree.

Conditioning on a common effect creates dependence between its causes. This is known as “explaining away”: once you know the alarm rang, learning it was a burglar makes an earthquake less likely.

Why the collider matters: selection bias

Selection bias is the distortion you get from analyzing a sample that was filtered on an outcome, rather than drawn at random. It is a collider, and it is everywhere in real data.

Take X, Y independent Bernoulli(0.5) — each is 1 with probability one-half and 0 otherwise. Keep a row only when X + Y >= 1. That filter is the collider: you are conditioning on an effect of both variables.

The block below enumerates all four possible rows, then drops the one the filter rejects. Notice what happens to the two conditional probabilities on the last two lines — they should be equal if X and Y are independent, and they are not.

population, 4 equally likely cells:   (0,0) (0,1) (1,0) (1,1)
admitted (drop (0,0)):                      (0,1) (1,0) (1,1)

P(Y=1 | X=1, admitted) = 1/2     from cells (1,0), (1,1)
P(Y=1 | X=0, admitted) = 1       from cell  (0,1)

Walk the arithmetic. Among admitted rows with X = 1 there are two equally likely cells, (1,0) and (1,1), and only one of them has Y = 1 — so 1/2. Among admitted rows with X = 0 only one cell survives the filter at all, (0,1), and it has Y = 1 — so 1.

So X = 1 predicts a lower chance of Y = 1. Two independent variables now look negatively correlated, and no amount of data fixes it, because the bias is in the sampling rule rather than the sample size.

Every dataset filtered on an outcome has this structure: approved loans, admitted students, clicked impressions. It is the cleanest explanation of why “the model does badly on rejected applicants” is a data-generation problem, not a modeling one.

What has to be true, and how it gets misused. Conditional independence is a claim about a specific conditioning set, and it survives nothing about changing that set. Adding a collider to the set destroys independence that held without it. Removing a fork destroys independence that held with it.

The classic misuse is “control for everything you have”, which sweeps colliders into the regression alongside genuine confounders — variables that really do drive both things you are comparing — and so manufactures associations that are not there.

The second classic misuse is reading a filtered dataset as a population. Any correlation measured on approved loans is a correlation among the approved, and it does not carry over to applicants in general.

Why Naive Bayes survives an assumption that is always false

Naive Bayes is a classifier that assumes every feature is conditionally independent of every other, given the class label. Under that assumption P(y | x) is proportional to P(y) · prod_j P(x_j | y) — the class’s base rate, times one factor per feature.

That assumption is false everywhere. The words new and york are not conditionally independent given “spam.”

Take logarithms and the reason it survives anyway appears. Two words in the formula below need defining first:

log odds(y=1 | x)  =  log odds(y=1)  +  sum_j log [ P(x_j | y=1) / P(x_j | y=0) ]
                      ^^^^^^^^^^^^^     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                      prior term        one log-likelihood-ratio per feature

So the log-odds of the class start at the prior, and each feature then adds one term: the logarithm of how much more likely that feature is under class 1 than under class 0.

Correlated features contribute the same evidence several times, so the sum overshoots in magnitude. But it overshoots in the direction the true evidence pointed. That kind of distortion is monotone — it stretches the numbers without ever reordering them — and a monotone distortion preserves both the ranking and the argmax, the class with the largest score.

Naive bayes carries this to numbers. Three correlated tokens each with a likelihood ratio of 4 get multiplied as 4^3 = 64 instead of the true 4. Starting from prior odds of 0.25, the true posterior odds are 0.25 × 4 = 1.0, giving p = 0.500; Naive Bayes reports 0.25 × 64 = 16.0, giving p = 16/17 = 0.941. The probability is badly wrong. The classification at a 0.5 threshold is identical.

Use Naive Bayes scores for ranking and argmax, never for thresholding or expected-value arithmetic without calibrating first. Calibration is the post-hoc adjustment that makes a score usable as a probability, covered in Calibration what it means and when it matters.

What has to be true, and how it gets misused. The conditional-independence assumption is required for the probabilities to be right. It is not required for the ordering to be right. That is the whole content of the paragraph above.

The misuse follows directly: taking a Naive Bayes output of 0.94 as a 94% chance, then feeding it into an expected-loss calculation or a fixed 0.5 threshold tuned on a different feature set. The more correlated your features, the further that number is from a probability — and adding more correlated features makes it worse, not better.

3. Bayes’ theorem, and the question everyone gets wrong

Most of applied inference is one move: you know how likely the evidence is under each hypothesis, and you want the reverse — how likely each hypothesis is given the evidence. Bayes’ theorem is that reversal.

P(H | E) = P(E | H) · P(H) / P(E)      with  P(E) = sum_h P(E | h)·P(h)

odds form:   posterior odds = prior odds × likelihood ratio

Here H is a hypothesis (you have the disease) and E is evidence (the test came back positive).

Four names attach to the four pieces, and they recur everywhere in machine learning:

The theorem itself is one line of algebra. The chain rule from Conditional probability and independence says P(H and E) can be written two ways, P(H | E)·P(E) and P(E | H)·P(H). Set them equal and divide by P(E). That is all Bayes’ theorem is.

The odds form says the same thing with the division already done. Multiply your prior odds by the likelihood ratio — the ratio of how probable this evidence is under the hypothesis versus against it — and you have your posterior odds.

The medical test, computed

A disease affects 1 in 1,000. A test has 99% sensitivity and 99% specificity. You test positive. What is the probability you have the disease?

Two definitions first. Sensitivity is the probability the test says positive when the disease is present. Specificity is the probability it says negative when the disease is absent, so a specificity of 99% means a 1% false-positive rate.

Translate the question into symbols, with D for “has the disease” and + for “tests positive”:

P(D)         = 0.001      1 in 1,000 have it
P(+ | D)     = 0.99       sensitivity
P(+ | not D) = 0.01       1 - specificity, the false-positive rate

Route 1: straight substitution

The first step is P(+), the marginal probability of a positive result. Collect the two disjoint ways a positive can arise — a sick person testing positive, or a healthy person testing positive by error — and add them. Then Bayes’ theorem divides the first by the total.

P(+) = 0.99 × 0.001  +  0.01 × 0.999
     = 0.00099       +  0.00999          =  0.01098
        sick and +      healthy and +

P(D | +) = 0.00099 / 0.01098 = 0.09016   ->  about 9%

Route 2: the odds form

This gets there in one line and no long division.

prior odds       = 0.001 / 0.999 = 1/999
likelihood ratio = 0.99 / 0.01   = 99
posterior odds   = 99 / 999      = 0.0991
p = odds / (1 + odds) = 0.0991 / 1.0991 = 0.0902

The last line converts odds back to a probability. If you remember only one form of Bayes’ theorem, remember this one.

Route 3: count people

This is the route to say out loud in an interview, because it converts the whole problem into counting.

The tree below tracks 100,000 tested people down two branches — sick and healthy — and then splits each by test result. Look at the two boxes that lead into the green node: those are the positives, and they are wildly unequal in size.

flowchart TD
    P["100,000 people tested"] --> D["100 have the disease"]
    P --> H["99,900 do not"]
    D --> TP["99 test positive<br/>true positives"]
    D --> FN["1 tests negative"]
    H --> FP["999 test positive<br/>FALSE positives"]
    H --> TN["98,901 test negative"]
    TP --> R["1,098 positives total<br/>99 of them real<br/>99 / 1,098 = 9.0%"]
    FP --> R

    style FP fill:#9d0208,color:#fff
    style R fill:#2d6a4f,color:#fff

Walk the tree one step at a time:

  1. Of 100,000 people tested, at a base rate of 1 in 1,000, exactly 100 have the disease. The other 99,900 do not.
  2. Among the 100 sick, 99% sensitivity means 99 test positive — the true positives — and 1 tests negative.
  3. Among the 99,900 healthy, a 1% error rate means 0.01 × 99,900 = 999 test positive — the false positives — while 98,901 test negative.
  4. Adding the two positive branches: 99 + 999 = 1,098 positives in total, of which 99 are real.
  5. So a positive result is right 99 times out of 1,098, which is 9.0%.

The false positives outnumber the true positives 10 to 1, because the 1% error rate applies to a group 999 times larger. The arithmetic above backs up that sentence.

Why the intuition fails

Almost everyone answers 99%, and the failure is mechanical: they substitute P(+ | D) for P(D | +).

Those two are different conditionals over different denominators. P(+ | D) divides by the sick. P(D | +) divides by everyone who tested positive. Bayes’ theorem is the exchange rate between them — here a factor of P(D)/P(+) = 0.001/0.01098 = 0.091.

The term people drop is the base rate: the prevalence of the disease in the population being tested. Its leverage is enormous.

The table below holds the test fixed at 99% sensitivity and 99% specificity, and changes only the base rate. Read down the middle column and watch the answer swing from 1% to 92%.

Prevalence P(D)P(D | +)What changed
0.00010.98%Test unchanged, disease 10× rarer, answer 9× smaller
0.0019.02%The stated problem
0.0150.0%Prior odds 1/99 × LR 99 = 1
0.1091.7%Now the intuition is nearly right

The third column’s “LR” is the likelihood ratio, which is 99 for this test in every row. Identical test, four answers, and the only thing that moved was the prior.

The 0.01 row is the landmark worth memorizing: a 99/99 test breaks even at exactly 1% prevalence. Prior odds of 1/99 cancel the likelihood ratio of 99 to give posterior odds of exactly 1, which is a probability of one-half.

This governs every rare-event classifier you will ship: fraud at 0.1%, anomaly detection on healthy traffic, a safety filter on a benign corpus.

A 99%-accurate detector on a 0.1% base rate produces ten false alarms per true one, and the humans reviewing those alerts will stop trusting it within a week. The fix is not a better threshold. It is raising the base rate by pre-filtering the population you run the detector on.

What has to be true, and how it gets misused. Bayes’ theorem needs P(E) > 0, and needs the hypotheses summed over in P(E) to be mutually exclusive and exhaustive. Drop a hypothesis and every posterior is inflated.

The worked example adds two assumptions of its own: that the sensitivity and specificity measured in a trial transfer unchanged to the person in front of you, and that the base rate you plugged in is the rate in the population you are actually testing. Both fail routinely.

Screening a self-selected population that came in because of symptoms has a far higher base rate than the general population, so quoting the general-population 9% to a symptomatic patient is wrong in the pessimistic direction. Quoting a trial’s prevalence for a mass-screening program is wrong in the optimistic direction.

The third misuse is repeat testing. Two positives from the same test get treated as two independent likelihood ratios of 99, multiplying to 9,801. But a person whose biology fools the test once will fool it again: the errors are correlated, and the true combined evidence is far weaker.

4. Random variables, expectation, variance, covariance

Updating beliefs about events is half the job; the other half is summarizing a number attached to the outcome — a revenue, a latency, a count. Four summary quantities do that work, and the rules they obey split sharply into the ones that hold unconditionally and the ones that quietly require independence. That split decides whether an argument you are making is valid.

Expectation

Expectation is the probability-weighted average of a random variable: multiply each value by its probability, and add. Write it E[X] and read it “the expectation of X” or “the expected value of X”. It is the number the average of many draws settles down to.

Do it once on a fair die. Each of the six faces has probability 1/6:

E[X] = (1/6)(1) + (1/6)(2) + (1/6)(3) + (1/6)(4) + (1/6)(5) + (1/6)(6)
     = (1 + 2 + 3 + 4 + 5 + 6) / 6
     = 21 / 6
     = 3.5

Note that 3.5 is not a face the die can show. An expectation is a weighted average, not a prediction of any single draw.

The property that matters most is that expectation is linear, with no independence assumption whatsoever:

E[aX + bY + c] = a·E[X] + b·E[Y] + c            always
E[XY] = E[X]·E[Y]                               only if uncorrelated

That asymmetry between the two lines does more work than it looks like.

Linearity-without-independence is what makes indicator arguments work. To count something, write the count as a sum of indicators — variables that are 1 when an event happens and 0 otherwise. The expectation of an indicator is just the probability of its event, since E = 1·p + 0·(1-p) = p. Then add those probabilities up, regardless of dependence.

The coupon collector problem is that move. The birthday problem takes the other route, a complement-rule product, and the two are worth contrasting.

Variance, covariance, correlation

Expectation tells you where a random variable sits. The next three quantities tell you how much it moves.

Var(X)    = E[X^2] - (E[X])^2
Cov(X, Y) = E[XY] - E[X]·E[Y]
Corr(X,Y) = Cov(X,Y) / (sd(X)·sd(Y))            in [-1, 1], unit-free

Variance is the expected squared distance from the mean. Because the distances are squared, variance is in squared units. Its square root, the standard deviation written sd(X), returns to the original units, which is why you report standard deviations to people and carry variances through algebra.

The formula Var(X) = E[X^2] - (E[X])^2 is the computational shortcut. Do it on the same die. You already have E[X] = 3.5, so you need E[X^2] — the average of the squared faces:

E[X^2] = (1 + 4 + 9 + 16 + 25 + 36) / 6 = 91 / 6 = 15.1667

Var(X) = 15.1667 - (3.5)^2
       = 15.1667 - 12.25
       = 2.9167

sd(X)  = sqrt(2.9167) = 1.708

Sanity-check that against the definition directly: the squared distances from 3.5 are 6.25, 2.25, 0.25, 0.25, 2.25, 6.25, which average to 17.5 / 6 = 2.9167. Same number. The shortcut is just algebra, not a different quantity.

Covariance is the same construction run across two variables. It is positive when they tend to be large together, negative when one is large as the other is small, and zero when there is no linear tendency either way.

Covariance inherits the units of both variables — dollars times seconds — which makes its magnitude uninterpretable on its own. Dividing by both standard deviations cancels the units and produces the correlation, a unit-free number pinned between -1 and 1.

What has to be true, and how it gets misused. All four quantities assume the relevant averages are finite, which is not automatic. The Cauchy distribution at the end of The distributions and what each one models has no mean at all, and heavy-tailed business metrics behave like it at practical sample sizes.

Linearity of expectation is unconditional, so it is never the step that breaks. E[XY] = E[X]·E[Y] is the step that breaks.

The classic misuse is treating the expectation of a product, a ratio, or any other nonlinear function as the function of the expectations. E[1/X] is not 1/E[X]. The expected value of a revenue-per-click ratio is not the ratio of expected revenue to expected clicks. That last one gets shipped as a metric constantly.

Covariance measures the linear part and nothing else

Let X ~ Uniform(-1, 1), meaning X is equally likely to land anywhere between -1 and 1. Let Y = X^2.

Work the covariance term by term:

Y is completely determined by X — you can compute it exactly — and their correlation is still exactly zero.

The lesson: uncorrelated means no linear association, and nothing more. Independent implies uncorrelated, but not the reverse.

The consequence is common in practice. A correlation-based feature screen silently discards every non-monotone predictor: a U-shape, a threshold effect, an interaction. Correlation near zero, cut before a tree ever sees it. That is the mechanism behind the warning in Feature selection.

Variance of a sum, and the bagging floor derived

Variance does not simply add; it adds plus a covariance correction, and scaling a variable scales its variance by the square:

Var(X + Y) = Var(X) + Var(Y) + 2·Cov(X, Y)          Var(aX) = a^2·Var(X)

The second rule is what makes averaging useful. The first is what limits how useful.

Take B identically distributed variables — think of them as B decision trees fit on resampled copies of one dataset, which is the technique called bagging. Each has variance s^2, and any two of them have correlation rho. Average them.

Var( (1/B) sum_b T_b ) = (1/B^2) [ B·s^2 + B(B-1)·rho·s^2 ]
                       = s^2/B  +  (B-1)·rho·s^2 / B
                       = rho·s^2  +  (1 - rho)·s^2 / B

Line by line:

  1. The 1/B^2 out front is the scaling rule Var(aX) = a^2·Var(X) with a = 1/B. Inside the bracket, expanding the variance of a sum gives B variance terms (one per tree, each s^2) plus B(B-1) covariance terms (one per ordered pair), and each covariance is rho·s^2 because Corr = Cov/(sd·sd) and both standard deviations are s.
  2. Divide both bracket terms by B^2.
  3. Rearrange. Check it: rho·s^2 + s^2/B - rho·s^2/B equals s^2/B + rho·s^2·(B-1)/B, which is line 2.

The second term of the final line vanishes as B grows. The first term does not depend on B at all.

Two sanity checks confirm the algebra. Set rho = 0 and you get s^2/B, the familiar shrinkage of an average of independent draws. Set rho = 1 and you get s^2, since averaging identical copies buys nothing.

Put numbers on it. Suppose s^2 = 1 and the trees are mildly correlated at rho = 0.05:

B = 100    ->  0.05 + 0.95/100  = 0.0595
B = 1000   ->  0.05 + 0.95/1000 = 0.05095

Ten times as many trees moved the variance by less than 2%, because the 0.05 floor is fixed.

That floor is the entire theory of random forests. Averaging more trees cannot get past rho·s^2, so the only remaining lever is decorrelation: lowering rho by giving each tree a different random subset of features to split on, which is what the max_features setting does. Bagging a variance fix and why it needs unstable base learners uses this formula; this is where it comes from.

What has to be true, and how it gets misused. The derivation assumes the B variables share one variance s^2 and one common pairwise correlation rho. It does not assume independence — which is exactly why it applies to trees grown on overlapping bootstrap samples.

The classic misuse is the opposite simplification: quoting s^2/B and concluding that enough trees, enough random seeds, or enough repeated experiment runs drive variance to zero. They do not. The residual is rho·s^2.

The same error appears in monitoring. The standard error of a metric — the standard deviation of the metric’s own estimate, which is what a confidence interval is built from — gets computed as if n correlated events from the same users were n independent samples, producing intervals that are too narrow by exactly this factor.

5. Total expectation and total variance

Often the easiest way to compute an average or a variance is in stages: condition on something, work inside each group, then average over the groups. Two identities license that, and the bias-variance decomposition — the standard account of why models fail — falls out of the second one.

The pattern to notice in both is that each computes something about Y by first conditioning on a second variable X, and then averaging over X.

total expectation:  E[Y]   = E[ E[Y | X] ]
total variance:     Var(Y) = E[ Var(Y | X) ]  +  Var( E[Y | X] )
                             ^^^^^^^^^^^^^^     ^^^^^^^^^^^^^^^^
                             within-group        between-group

Read E[Y | X] as “the expectation of Y given X”. The thing to notice is that because X is itself random, E[Y | X] is also a random variable, not a fixed number — it is a different average for each value X takes.

Think of X as splitting your data into groups. Then:

Every variance decomposition you will meet is the second identity with a different choice of grouping variable.

Worked: revenue per user

This is the calculation those two lines exist for.

Set up the numbers. A user makes N purchases, where N has mean 3 and variance 9. Each purchase is worth V, with mean 40 and variance 400. The values are independent of N and of each other. Total revenue Y is the sum of N values of V.

The trick is to condition on N, because once you fix the number of purchases both inner quantities become easy:

Now substitute those into the two identities:

E[Y]   = E[ E[Y | N] ] = E[ 40·N ] = 40 × 3 = 120

Var(Y) = E[ Var(Y|N) ] + Var( E[Y|N] )
       = E[ 400·N ]    + Var( 40·N )
       = 400 × E[N]    + 40^2 × Var(N)
       = 400 × 3       + 1600 × 9
       = 1,200         + 14,400          = 15,600     sd = 124.9

The 1600 is 40^2, from the scaling rule Var(aX) = a^2·Var(X). The sd is sqrt(15,600) = 124.9.

Of the 15,600 total, 14,400 sits in the between-group term — that is 14,400 / 15,600 = 92%.

92% of the variance in revenue per user comes from how often people buy, not how much they spend per purchase.

That decomposition is useful three ways:

What has to be true, and how it gets misused. The identities themselves need only that the relevant means and variances are finite. There is no independence requirement in them at all.

The worked example is where the assumptions live: the purchase values must be independent of the purchase count and of each other, with the same mean and variance for each.

The classic misuse is applying this compound formula when value and frequency are correlated. Heavy buyers usually buy cheaper items per order, which the formula above assumes away, and the true variance will not match.

The second misuse is reading the 92% as a causal claim about which lever to pull. It is a statement about where variance sits today, not about what an intervention would move.

The bias-variance split, derived from the same identity

Set up the problem. Data is generated as y = f(x) + e, where f is the true underlying function and e is noise with E[e] = 0 and Var(e) = sigma^2. You fit a model f_hat on a random training set D. Fix one input x, and take the expectation over both the noise and the random draw of D.

E[(y - f_hat(x))^2]
  = E[(f(x) - f_hat(x))^2] + sigma^2                  e independent, mean 0

write  Ef := E_D[f_hat(x)]  and add-subtract it:

  = E[ ( (f(x) - Ef) + (Ef - f_hat(x)) )^2 ] + sigma^2
  = (f(x) - Ef)^2 + E[(Ef - f_hat(x))^2] + 0 + sigma^2       cross term is 0
  =    Bias^2     +      Variance        +     sigma^2

Two steps in there are where people stall. Both are worth doing slowly.

Step 1, peeling off the noise. Substitute y = f(x) + e and expand:

E[(f(x) + e - f_hat(x))^2]
  = E[(f(x) - f_hat(x))^2] + 2·E[ e·(f(x) - f_hat(x)) ] + E[e^2]

The middle term is zero because e is independent of the training set, so the expectation factors into 2·E[e]·E[f(x) - f_hat(x)], and E[e] = 0. The last term is E[e^2] = Var(e) = sigma^2, since e has mean zero. That leaves the first line of the block above.

Step 2, the add-subtract trick. Define Ef := E_D[f_hat(x)], the average prediction at x over all the training sets you might have drawn. Adding and subtracting Ef inside the square changes nothing, and then you expand the square of a sum:

( (f - Ef) + (Ef - f_hat) )^2
  = (f - Ef)^2  +  2·(f - Ef)·(Ef - f_hat)  +  (Ef - f_hat)^2

Take E_D of each of the three pieces:

So Bias is how far the average fit sits from the truth, and Variance is how much any one fit scatters around that average.

Bias^2 and Variance are the between-group and within-group terms of the total-variance identity, with “group” being the training set you happened to draw. The decomposition everything follows from turns this into an ensemble strategy.

What has to be true, and how it gets misused. This decomposition assumes squared-error loss, noise with mean zero and constant variance that is independent of the input, and a fixed true f. The expectation is over training draws, at one input x.

Change the loss and the clean three-way split does not survive. There is no equally tidy bias-variance identity for 0-1 classification error, which just counts wrong labels, nor for log loss, which is the cross-entropy of Entropy cross entropy kl. Quoting one anyway is the most common misuse.

The second misuse is calling the sigma^2 term “bias”. It is irreducible noise, and no model, no matter how large, removes it. If your reported error is already near sigma^2, a bigger model buys nothing.

6. The distributions, and what each one models

Six distributions cover nearly every situation an interview or a production system will throw at you, and they are not six separate facts to memorize: each arises as a limit or a variation of another, and each carries a physical assumption that has to hold before you are entitled to use it.

The diagram below is a family tree. There is one ancestor at the top, and every other distribution is reached by an arrow labelled with what you did to get there. The three green nodes are the three you reach by taking a limit, which is also why they are the three that carry an assumption that can fail.

flowchart TD
    B["Bernoulli(p)<br/>one yes/no trial"] --> BI["Binomial(n, p)<br/>successes in n trials"]
    BI -->|"n large, p small<br/>np = lambda fixed"| PO["Poisson(lambda)<br/>events per interval"]
    BI -->|"n large, np large<br/>CLT"| NO["Normal(mu, sigma^2)<br/>sums of many things"]
    PO -->|"time BETWEEN events"| EX["Exponential(lambda)<br/>memoryless waiting"]
    B -->|"trials until<br/>first success"| GE["Geometric(p)<br/>discrete waiting"]
    GE -.->|"discrete analog"| EX

    style PO fill:#2d6a4f,color:#fff
    style NO fill:#2d6a4f,color:#fff
    style EX fill:#2d6a4f,color:#fff

Walk the tree from the top:

Why the shading matters: Bernoulli, binomial and geometric are exact descriptions of a finite number of yes/no trials, and you can write their probabilities down by counting. Poisson, normal and exponential are what those become as the trials get numerous — rare ones give the Poisson, summed ones give the normal, and trials shrunk to instants give the exponential.

Bernoulli and binomial

Bernoulli(p) is one trial that succeeds with probability p. Its two summary numbers are E = p and Var = p(1-p).

That variance is maximized at p = 0.5 and collapses at the extremes:

p = 0.50  ->  Var = 0.50 × 0.50 = 0.25
p = 0.01  ->  Var = 0.01 × 0.99 = 0.0099        25× smaller

That drives A/B sample sizes for rare conversions: a metric with less variance needs fewer samples to measure to the same precision.

The same p(1-p) shows up in a second, unrelated-looking place. It is the second derivative of the log loss of Entropy cross entropy kl with respect to a model’s raw score. So a row whose predicted probability is near 0 or 1 contributes almost nothing to that curvature, while a row near 0.5 contributes the most. Boosted trees use exactly that number as a per-row weight when they score a split, where it is written h_iXgboost second order and a regularizer inside the split criterion is the only place you need it.

Binomial(n, p) counts successes in n independent trials, with E = np and Var = np(1-p). Both come straight from the variance-of-a-sum rule with every covariance zero: n copies of the same Bernoulli, added.

What has to be true, and how it gets misused. The binomial requires a fixed number of independent trials with the same success probability throughout. Both halves fail in production: conversion rates drift across the day, and two events from one user or one session are not independent trials.

The classic misuse is computing a confidence interval on a conversion rate from np(1-p) when the unit of randomization was the user and the unit of counting was the pageview. That understates the true variance — the same error of treating clustered rows as independent ones that Variance of a sum and the bagging floor derived flagged.

Poisson, derived as the rare-event limit of the binomial

The Poisson distribution counts how many events land in a fixed interval, when each individual opportunity for an event is rare but there are very many of them.

Derive it by sending n -> infinity and p -> 0 while holding the product np = lambda fixed, so the expected count stays put.

Start from the binomial PMF, substitute p = lambda/n, and sort the pieces into four factors. The three underlined factors each have a limit; the one that is not underlined is the answer.

P(X = k) = C(n,k) p^k (1-p)^(n-k)

  = (lambda^k / k!) · [ n(n-1)...(n-k+1) / n^k ] · (1 - lambda/n)^n · (1 - lambda/n)^(-k)
                       ^^^^^^^^^^^^^^^^^^^^^^^     ^^^^^^^^^^^^^^^^   ^^^^^^^^^^^^^^^^^^
                              -> 1                     -> e^-lambda        -> 1

  = e^-lambda · lambda^k / k!

The rearrangement on the second line is pure bookkeeping, and it is worth checking once. C(n,k) is n(n-1)...(n-k+1) / k! — a falling product of k terms, over k!. Substituting p = lambda/n turns p^k into lambda^k / n^k. And (1-p)^(n-k) splits into (1 - lambda/n)^n times (1 - lambda/n)^(-k), because subtracting k in the exponent is dividing by k copies. Pull lambda^k / k! out front, park the n^k under the falling product, and you have the four factors shown.

Now take the three underlined factors one at a time:

Multiply the survivors and you have e^-lambda · lambda^k / k!.

The convergence is fast enough to be useful at ordinary sizes. Check it at n = 1,000 and p = 0.002 — so lambda = np = 2 — asking for the probability of exactly k = 3 events:

binomial:  C(1000,3) · 0.002^3 · 0.998^997 = 166,167,000 · 8e-9 · 0.13587 = 0.18062
poisson:   e^-2 · 2^3 / 3!                 = 0.135335 · 8/6              = 0.18045

That is three-decimal agreement at n of only 1,000.

The Poisson has the striking property that its mean and its variance are the same number: E = Var = lambda. That equality is a testable claim, not a definition. Real counts are usually overdispersed, meaning their variance exceeds their mean, which is why insurance and click models reach for the negative binomial instead.

What has to be true, and how it gets misused. A Poisson count requires events that arrive independently at a constant rate, with no clustering and no bursts. Almost no real arrival process satisfies that: requests arrive in bursts, errors arrive in cascades, clicks arrive in sessions.

The classic misuse is sizing capacity or alert thresholds from a Poisson tail — “lambda = 100 per minute, so 130 is a 3-sigma event” — when the real traffic is bursty and 130 happens hourly.

The cheap diagnostic is the assumption itself. Compute the sample mean and the sample variance and compare them. If the variance is much bigger, the Poisson is the wrong model and its tail probabilities are far too small.

Exponential, and memorylessness proved

The exponential distribution answers the other question about a stream of events: not how many arrive, but how long you wait for the next one.

If events arrive as a Poisson process at rate lambda, the wait to the next one is Exponential(lambda), with

P(T > t) = e^(-lambda·t)          E[T] = 1/lambda

Read P(T > t) as “the probability that the wait T exceeds t”. At lambda = 2 events per second the mean wait is half a second, and the chance of waiting more than a full second is e^-2 = 0.135.

The proof

Its defining property drops straight out of that formula. Ask: given that you have already waited s, what is the chance you wait another t?

Apply the definition of conditional probability. The numerator is P(T > s+t and T > s), but waiting past s + t already implies waiting past s, so the “and” collapses to just P(T > s+t):

P(T > s + t | T > s) = e^(-lambda(s+t)) / e^(-lambda·s) = e^(-lambda·t) = P(T > t)

The s cancels because dividing exponentials subtracts their exponents. The distribution of the remaining wait does not depend on how long you have already waited.

That property is called memorylessness. The exponential is the only continuous distribution with it, and the geometric is its unique discrete counterpart.

Knowing when memorylessness is wrong

The general way to talk about this is the hazard rate: the instantaneous chance that the thing ends in the next moment, given that it has lasted this long. Memorylessness is exactly a constant hazard rate.

Sort real processes by whether the hazard rate is constant:

Where the hazard rate rises, a fixed timeout plus retry is the correct policy: the request that has run long is the one worth abandoning.

Assuming exponential where the hazard rate rises is how a retry policy amplifies an outage instead of surviving it, because the retries pile onto a system whose slow requests were never going to finish.

What has to be true, and how it gets misused. The exponential requires a constant hazard rate and — for the Poisson-process link — independent arrivals at a constant average rate.

The classic misuse is the timeout policy above. The second is estimating a mean wait as 1/lambda from a sample that was truncated by an existing timeout, which throws away exactly the long waits that carry the mean.

Normal, and the CLT stated precisely

The central limit theorem (CLT) is the reason the normal distribution is everywhere, and it says something narrower than most people remember. Let X_1, ..., X_n be independent, identically distributed, with mean mu and finite variance sigma^2. Write Xbar_n for the sample mean of the first n of them. Then

sqrt(n) · (Xbar_n - mu) / sigma   ->   Normal(0, 1)   in distribution

Read the left side as “the sample mean’s distance from the true mean, rescaled by sqrt(n) over the standard deviation”. The rescaling matters: without it the quantity Xbar_n - mu just collapses to zero as n grows, and there would be no shape left to describe.

Normal(0, 1) is the standard normal — mean 0, variance 1, the familiar bell curve.

“Converges in distribution” is the weakest of the several kinds of convergence, and the choice of word is load-bearing. It means the CDF of the left side approaches the CDF of a standard normal at every point. It says nothing about any individual draw.

Now the part that separates a real answer from a recited one. The table below lists four things people say the CLT means. All four are wrong, and the right-hand column says why.

Claim people makeWhy it is wrong
“The data becomes normal with enough samples”The X_i never change. Only the standardized sample mean converges.
“The sample mean is normal for n >= 30Convergence is asymptotic. n = 30 is a rule of thumb about the center, nothing more.
“So a p-value of 1e-6 from a t-test is trustworthy”Convergence is in distribution — pointwise on the CDF. Relative error in the far tail converges far more slowly.
“It applies to any data”Finite variance is required, and so is (roughly) independence.

Two terms from the third row. A p-value is the probability of seeing a result at least this extreme if nothing were going on. A t-test is the standard procedure that produces one for a difference in means. Both are covered in Hypothesis testing. The row’s point is that a p-value of 1e-6 lives in the extreme tail, where the CLT’s guarantee is weakest.

Berry-Esseen: how wrong the approximation is at n = 30

The tail warning has an actual number attached to it.

The Berry-Esseen theorem bounds the uniform CDF error — the largest gap between the true CDF and the normal approximation, anywhere on the axis:

max CDF error  <=  C · rho / (sigma^3 · sqrt(n))        with rho = E|X - mu|^3

Three things about that formula.

C ≈ 0.4748 is the sharpest known constant for identically distributed summands, due to Shevtsova (2011). It is not derivable here; it is quoted, unlike everything else in this chapter.

rho is the third absolute central moment — the average cubed distance from the mean, ignoring sign. Note the absolute-value bars. This is not the signed quantity that appears in the numerator of skewness, and the two differ. For Exponential(1), E[(X-1)^3] = 2 but E|X-1|^3 = 2.41. Using the signed one understates the bound.

The bound shrinks like 1/sqrt(n), which is slow. Quadrupling your sample size only halves the guaranteed error.

Evaluate it for Exponential(1) at n = 30. That distribution has sigma = 1, so sigma^3 = 1 and the denominator’s cube disappears:

0.4748 × 2.41 / sqrt(30) = 1.1443 / 5.477 = 0.209

A worst-case CDF error of 21 percentage points, at the sample size everyone treats as sufficient. The error is smaller near the center and worse in the tail, which is exactly why the bootstrap — the resampling procedure in Bootstrap — beats a normal approximation on skewed metrics.

The Cauchy distribution: when the CLT does not apply at all

The pathological case is worth naming. The Cauchy distribution has tails so heavy that its variance is infinite, which breaks the CLT’s finite-variance requirement outright.

The consequence is startling: the mean of n Cauchy draws is itself Cauchy with the same scale. There is no concentration at any n. Averaging more data buys literally nothing.

Heavy-tailed metrics such as revenue, session length, and token counts behave like this at practical sample sizes.

What has to be true, and how it gets misused. The CLT needs independence, identical distribution, and finite variance. Each fails in a specific, recognizable way:

The classic misuse is the n >= 30 folklore applied to a skewed revenue metric. The sample mean at n = 30 is nowhere near normal in the tail, and the resulting confidence interval and p-value are both wrong in the direction that flatters the result.

7. Log-probabilities, with the arithmetic

Every real system stores and combines probabilities as logarithms — not out of mathematical taste, but because floating-point numbers run out. Standard double-precision floating point, called float64, cannot represent arbitrarily small positive numbers.

Compare the two numbers below. The left is the smallest normal value float64 can hold. The right is what a modest sequence probability actually is. They are not close.

float64 smallest normal ~ 2.2e-308      a 1,000-token sequence at p = 1e-4:  1e-4000

Chain-rule a probability per token and the product shrinks by a factor of 1e-4 each step. Two things go wrong, in order:

  1. After 308/4 = 77 tokens the product drops below 2.2e-308 and leaves the normal floating-point range. You are now in the subnormals — special small values below the normal range that trade away precision to reach down to 5e-324. The number is still nonzero, but it is shedding bits of the mantissa, the significant digits. The answer is quietly wrong without looking broken.
  2. Around 81 tokens the subnormals are exhausted too and the product becomes exactly the number 0. Every ratio downstream is 0/0 and every log is -inf.

In log space the same quantity is 1,000 × log(1e-4) = 1,000 × (-9.2103) = -9,210.3, which float64 represents at full relative precision without effort.

Three rules follow.

Rule 1 — multiplication becomes addition. log(a·b) = log a + log b. Chain rules, Naive Bayes products, sequence likelihoods: all become sums.

Rule 2 — addition needs log-sum-exp. Adding probabilities is the awkward case, because logs do not distribute over sums. The identity is:

log(sum_i exp(x_i)) = m + log(sum_i exp(x_i - m))        with  m = max_i x_i

Subtracting the max makes the largest term exp(0) = 1 and puts every other term in (0, 1]. Nothing can overflow, and nothing that mattered can underflow.

Rule 3 — softmax must subtract the max. Softmax is the function that turns a vector of arbitrary real scores — called logits — into probabilities that sum to 1, by exponentiating each score and dividing by the total.

Logits [1000, 1001, 1002] overflow exp to inf, and inf/inf is nan. Subtract the max, 1002, and the answer is unchanged — the factor e^-m appears in every numerator and in the denominator, so it cancels:

shifted [-2, -1, 0]  ->  exp [0.1353, 0.3679, 1.0000], sum 1.5032
                     ->  softmax [0.0900, 0.2447, 0.6652] = softmax([0, 1, 2])

Check one entry by hand: 1.0000 / 1.5032 = 0.6652.

The two functions below implement rules 2 and 3. The three assert lines at the bottom are the point of the block — they check that log-sum-exp returns the shifted-and-corrected value, that softmax still sums to 1, and that the largest logit gets the 0.665 computed above, all on inputs that would overflow a naive implementation.

import math


def log_sum_exp(xs):
    """Numerically stable log(sum(exp(x)))."""
    m = max(xs)
    if m == -math.inf:
        return -math.inf
    return m + math.log(sum(math.exp(x - m) for x in xs))


def softmax(logits):
    m = max(logits)
    exps = [math.exp(z - m) for z in logits]
    total = sum(exps)
    return [e / total for e in exps]


assert abs(log_sum_exp([1000.0, 1001.0, 1002.0]) - (1002.0 + math.log(1.5032))) < 1e-4
assert abs(sum(softmax([1000.0, 1001.0, 1002.0])) - 1.0) < 1e-12
assert abs(softmax([1000.0, 1001.0, 1002.0])[2] - 0.66524) < 1e-4

This is why a large language model (LLM) returns logprobs — log-probabilities — rather than probabilities.

It is also why perplexity, the standard measure of how uncertain a language model is, is defined as exp(mean negative log-likelihood). The mean is taken in log space where the arithmetic is safe, and the exponential is applied exactly once, at the end.

What has to be true, and how it gets misused. None of this is statistical. It is arithmetic about float64, and the only assumption is that you never materialize the raw product.

Three ways people materialize it anyway:

The assertions in the block above exist to catch exactly these.

8. Entropy, cross-entropy, KL

Log-probabilities are also the native language of information theory, whose three central quantities turn out to be where every standard loss function comes from.

entropy         H(p)     = - sum_x p(x)·log p(x)
cross-entropy   H(p, q)  = - sum_x p(x)·log q(x)
KL divergence   KL(p||q) = sum_x p(x)·log( p(x)/q(x) )  =  H(p,q) - H(p)

Here p and q are two distributions over the same set of outcomes. KL stands for Kullback-Leibler, and KL(p||q) is read “the KL divergence from q to p”.

The coding interpretation makes all three concrete. Imagine you have to transmit a stream of draws, and you get to design a code first.

Cost is measured in bits or in nats, which are the same quantity computed with base-2 or natural logarithms respectively.

KL is always non-negative, and that is the whole reason cross-entropy works as a loss. The surcharge is zero exactly when q = p, so driving cross-entropy down drives your model toward the truth.

Entropy on two coins, with the arithmetic:

fair coin, p = 0.5:
  H = -0.5·log2(0.5) - 0.5·log2(0.5) = 0.5 + 0.5             = 1 bit

biased coin, p = 0.9:
  H = -0.9·log2(0.9) - 0.1·log2(0.1) = 0.137 + 0.332         = 0.469 bits

The biased coin costs less than half as much to transmit, because you are usually right guessing heads. Push it to p = 1 and there is nothing left to encode: H = 0.

Cross-entropy loss, derived from maximum likelihood

This is the connection most candidates cannot make, and it takes four steps. Maximum likelihood estimation (MLE) is the principle that you should pick the parameters under which the data you actually saw is most probable.

Step 1 — the likelihood. Your model outputs q_theta(y | x), the probability it assigns to label y for input x under parameters theta. Assuming the rows are independent given their inputs, the probability of the whole dataset is the product L(theta) = prod_i q_theta(y_i | x_i).

Step 2 — log and negate. Maximizing a product is the same as minimizing the negative log of it, since the logarithm is increasing and negation flips the direction, and the log turns the product into a sum:

theta* = argmax_theta  sum_i log q_theta(y_i | x_i)
       = argmin_theta  - (1/N) sum_i log q_theta(y_i | x_i)          <- the NLL

The quantity on the second line is the negative log-likelihood (NLL), averaged over the N rows.

Step 3 — recognize the sum. This is the step people cannot reproduce, and it is entirely bookkeeping.

Let p_hat(· | x_i) be the empirical distribution for row i: the distribution implied by the data alone. For a single labelled row that is a point mass — probability 1 on the observed label y_i, and 0 on every other class.

Now rewrite log q(y_i | x_i) as a sum over all K classes, using an indicator to select the one you want:

- (1/N) sum_i log q(y_i | x_i)  =  - (1/N) sum_i sum_k 1[y_i = k]·log q(k | x_i)
                                =  E_{x ~ data} [ H( p_hat(·|x), q(·|x) ) ]

Why the first line is an identity: 1[y_i = k] is 1 for exactly one value of k and 0 for the rest, so the inner sum has exactly one surviving term, log q(y_i | x_i). You added a sum that changes nothing.

Why the second line follows: 1[y_i = k] is p_hat(k | x_i), the point mass. Substitute it and the inner sum reads - sum_k p_hat(k|x_i)·log q(k|x_i), which matches the definition of H(p, q) from the top of this section term for term.

That is cross-entropy, exactly. Not analogous to it — the same expression.

Step 4 — connect to KL. Recall H(p,q) = H(p) + KL(p||q), which is the rearranged definition at the top of the section. Here p is p_hat, which is fixed by the data and contains no theta. Minimizing over theta cannot touch H(p_hat), so dropping it changes nothing:

argmin_theta H(p_hat, q_theta) = argmin_theta KL(p_hat || q_theta)

The diagram below chains those four steps left to right, with two special cases branching off the middle node.

flowchart LR
    A["Maximize<br/>likelihood"] --> B["Minimize<br/>negative log-likelihood"]
    B --> C["Minimize<br/>cross-entropy<br/>H(p_hat, q)"]
    C --> D["Minimize<br/>KL(p_hat || q)"]
    C --> E["binary case:<br/>-[y log p + (1-y) log(1-p)]"]
    C --> F["Gaussian noise:<br/>MSE"]

    style C fill:#2d6a4f,color:#fff
    style E fill:#40916c,color:#fff
    style F fill:#40916c,color:#fff

Maximum likelihood, minimum cross-entropy, and minimum KL to the empirical distribution are three names for one optimization.

The two branches off the middle node are the follow-up questions an interviewer asks next, and they now answer themselves.

Binary cross-entropy

In the binary case there are two classes and one number p to predict.

The Bernoulli likelihood for one row is p^y (1-p)^(1-y). Check both cases: when y = 1 the exponents are 1 and 0, giving p; when y = 0 they are 0 and 1, giving 1-p. That is what you wanted.

Its negative log is -[y·log p + (1-y)·log(1-p)], the standard binary cross-entropy loss. Nothing was invented.

Why MSE for regression

Assume the target is the true function plus Gaussian noise — Gaussian being the other name for the normal distribution — so y = f_theta(x) + e with e ~ Normal(0, sigma^2).

Write down the log of the normal density:

log q(y|x) = -0.5·log(2·pi·sigma^2) - (y - f_theta(x))^2 / (2·sigma^2)

Only the second term contains theta, and sigma^2 is a positive constant that cannot change which theta wins. So maximizing the likelihood is minimizing sum_i (y_i - f_theta(x_i))^2 — the sum of squared errors.

“MSE for regression, cross-entropy for classification” is one principle under two noise models.

Stating it that way tells you what to do when neither noise model holds. Poisson counts get a Poisson NLL. Heavy-tailed targets get a Laplace NLL, which works out to mean absolute error (MAE).

Label smoothing is the same idea approached from the other side. A one-hot target is the vector that is 1 at the true class and 0 everywhere else. Replace it with (1-eps)·onehot + eps/K over K classes. Because no target probability is 0 or 1 any more, no logit can run to infinity chasing one.

Cross-entropy and perplexity on one row

Take 3 classes, with the true class being class 2, and a model that outputs q = (0.1, 0.7, 0.2).

cross-entropy = -log(0.7)     = 0.3567 nats = 0.5146 bits
perplexity    = exp(0.3567)   = 1.4286      = 1 / 0.7

The target is one-hot, so only the 0.7 the model gave the correct class survives the sum — the other two terms are multiplied by zero.

Perplexity is exp(cross-entropy), and it is read as “effectively choosing uniformly among this many options.” For a one-hot target it works out to exactly 1/q(correct), as the last line shows.

A language model at 2.0 nats has perplexity exp(2.0) = 7.39: it is as uncertain as someone guessing among 7.39 equally likely tokens.

What has to be true, and how it gets misused. The derivation assumes rows are independent given their inputs, that the model family can represent something close to the truth, and that the labels are actual observed outcomes rather than someone else’s predictions.

The MSE branch adds two more: Gaussian noise, and constant noise variance across inputs, called homoscedasticity.

The classic misuse is MSE on a target whose spread grows with its level — revenue, latency, counts. There the constant-variance assumption is false and the fit is dominated by the large-valued rows. The fix is a log transform or the right NLL, not a bigger model.

The second classic misuse is treating cross-entropy as if it required labels to be one-hot. It does not. The same formula handles soft labels, label smoothing, and distillation — training a small model on a large model’s output probabilities rather than on hard labels — without any change at all.

KL is asymmetric, and the asymmetry is the whole story

KL(p||q) and KL(q||p) are different numbers, and which one you minimize changes what your model does.

Compute both on the same pair. Take p = (0.5, 0.5, 0.0) and q = (0.98, 0.01, 0.01), and note that p assigns zero to the third outcome while q does not.

KL(p||q) = 0.5·ln(0.5/0.98) + 0.5·ln(0.5/0.01) + 0
         = -0.336 + 1.956 = 1.620 nats

KL(q||p) = 0.98·ln(0.98/0.5) + 0.01·ln(0.01/0.5) + 0.01·ln(0.01/0)  =  INFINITE

In the first line the third term is 0·ln(0/0.01), which is taken as 0 — no mass, no cost.

In the second line the third term divides by zero, because q puts mass on an outcome that p rules out entirely. The general rule: KL(a||b) blows up wherever a > 0 and b = 0.

Which distribution sits in which slot therefore decides what your model is punished for. The table below names the two directions, what makes each one explode, and the behavior that produces.

DirectionBlows up whenMinimizer isWhere you meet it
Forward KL(p_data || q_model)model puts 0 where data has massmode-covering — spreads mass, hedgesmaximum likelihood, every standard classifier
Reverse KL(q_model || p)model puts mass where target has 0mode-seeking — locks onto one modevariational inference, RLHF’s KL penalty to a reference policy

A mode is a peak of the distribution. So mode-covering means smearing probability over every peak, to avoid ever assigning zero anywhere the data lives. Mode-seeking means committing to one peak and ignoring the rest.

The last cell names two places the reverse direction shows up. Variational inference is the technique that fits a simple distribution to approximate a complicated one. RLHF (reinforcement learning from human feedback) adds a KL penalty to keep a fine-tuned model from drifting away from the reference model it started as.

Maximum likelihood minimizes the forward KL, and that is why likelihood-trained generative models hedge. They are penalized infinitely for assigning zero probability anywhere the data appears, and not at all for assigning probability where it never does.

Two familiar symptoms of that: blurry samples from a variational autoencoder (VAE), a generative model trained with exactly this objective, and an LLM’s willingness to emit a plausible-but-wrong continuation. Both are the objective doing what it was told.

What has to be true, and how it gets misused. KL is only finite when the second distribution is positive everywhere the first one is, and it is defined only for distributions over the same outcome set.

The classic misuse is estimating a KL from two histograms with different bins, or with an empty bin in the denominator. Depending on the library that returns infinity or a silently clipped number. Psi and kl computed works through the drift-monitoring version of exactly this.

The second misuse is quoting “the KL” without saying the direction, which throws away the mode-covering versus mode-seeking distinction that is the entire point.

9. Worked interview problems

Everything above converges on four classic puzzles, the ones interviewers actually ask. Each also maps to a production situation. Work each to the arithmetic once and the production version becomes clear.

Birthday

How many people until two share a birthday with probability 1/2?

Use the complement rule: “no two people share” is a single product, while “at least two share” is not.

For n people over N = 365 days, each new person has to avoid the birthdays already taken:

P(no match) = (1 - 0/365)(1 - 1/365)(1 - 2/365) ··· (1 - (n-1)/365)

Person 1 is free. Person 2 must avoid 1 date, person 3 must avoid 2, and so on. At n = 23 that product is 0.4927, so P(match) = 1 - 0.4927 = 0.5073 — just over half.

The approximation that makes it explainable

You cannot multiply 23 terms in your head. Use 1 - x ≈ e^-x for small x, which turns the product into a single exponential whose exponent is the sum 0/365 + 1/365 + ... + (n-1)/365 = n(n-1)/(2·365):

P(no match) ≈ exp( - n(n-1) / (2·365) )

n = 23:  pairs = 253      253/365 = 0.693      exp(-0.693) = 0.500
n = 57:  pairs = 1,596  1,596/365 = 4.373      exp(-4.373) = 0.013   ->  98.7%

Both lines are the approximation, and the second shows how much it drifts. At n = 23 the approximation’s 0.500 sits 0.7 points under the exact 0.5073. At n = 57 its 98.7% sits 0.3 points under the exact product, which is 99.01%. Read these as the right order of magnitude and the right scaling law, never as an exact figure.

The answer is 23 because 23 people make 23 × 22 / 2 = 253 pairs, and the question was always about pairs.

The pair count n(n-1)/2 grows like n^2, which is why the threshold is so much smaller than intuition suggests. Solving exp(-n^2/(2N)) = 0.5 gives n ≈ 1.177·sqrt(N) = 1.177 × 19.1 = 22.5. The sqrt(N) scaling is the takeaway.

Where it shows up in machine learning

Feature hashing is the trick of mapping feature names to a fixed number of slots with a hash function, instead of keeping a dictionary. Hash 10,000 features into 2^20 = 1,048,576 buckets and the expected number of colliding pairs is C(10000,2) / 2^20 = 49,995,000 / 1,048,576 = 47.7.

Collisions are certain, not rare. That they cost little is a separate argument, made in ml/01; that they happen at all is pure birthday arithmetic.

The same applies to 64-bit request identifiers: generate 2^32 of them and there is a 39% chance of a duplicate.

What has to be true, and how it gets misused. The calculation assumes birthdays are uniform over 365 days and independent across people. Real birthdays are neither — they cluster by season and by weekday, and twins break independence outright. Both deviations make matches more likely, so 23 is the pessimistic end and the true threshold is slightly lower.

The misuse that matters is on the hashing side: assuming a hash function’s outputs are uniform when it is being fed adversarial or highly structured keys. That is how a hash table degrades to a linked list.

Coupon collector

n distinct coupon types, each draw uniform. How many draws to collect all n?

Decompose the total wait into stages, indexed by how many distinct types you already hold.

Suppose you hold i distinct types. A fresh draw is new with probability (n-i)/n, since n-i of the n types are still missing. The wait for the next new type is therefore geometric, and a geometric with success probability p has mean 1/p, so this stage costs n/(n-i) draws on average.

Now add the n stage means. Linearity of expectation lets you do this even though the stages are not independent of one another:

E[T] = sum_{i=0}^{n-1} n/(n-i) = n·H_n ≈ n·(ln n + 0.5772)

n = 50:   H_50 = 4.4992   ->  E[T] = 224.96      (approx: 50 × 4.489 = 224.5)
n = 365:  H_365 = 6.4785  ->  E[T] = 2,364.6

The sum reverses into n·(1/n + 1/(n-1) + ... + 1/1), which is n times H_n, the n-th harmonic number 1 + 1/2 + 1/3 + ... + 1/n. H_n grows like ln n plus the constant 0.5772.

The variance, and the two approximations hidden in the usual formula

The variance matters as much as the mean here, and the familiar version of it is two approximations deep. Both forms are worth writing down.

Exactly:

Var = n^2·H_n^(2) - n·H_n         where H_n^(2) = 1 + 1/4 + 1/9 + ... + 1/n^2

H_n^(2) is the second-order harmonic number — the same sum with the terms squared.

The asymptotic form you usually see is Var ≈ n^2·pi^2/6, giving sd ≈ 1.283·n. It makes two substitutions at once:

  1. It replaces H_n^(2) by its infinite limit pi^2/6 = 1.6449.
  2. It drops the - n·H_n term entirely.

At n = 50 those two shortcuts are worth 49.5 and 225.0 respectively, so:

exact:       Var = 3,837.9   ->  sd = 62.0
asymptotic:  Var = 4,112.3   ->  sd = 64.1

An sd of about 62 on a mean of 225. That is a long right tail, driven entirely by the wait for the last coupon, which alone costs n draws on average — the i = n-1 stage has success probability 1/n.

Three places this bites:

n·ln n is why “sample randomly until you’ve seen everything” costs ln n more than people budget.

What has to be true, and how it gets misused. The result assumes draws are independent and uniform over the n types.

Non-uniformity makes it dramatically worse, not slightly worse, because the expected time is dominated by the rarest type. Halve one coupon’s probability and you roughly double the tail.

The classic misuse is budgeting random sampling for full coverage of a long-tailed label distribution using n·ln n on the number of classes, when the rare classes have probabilities orders of magnitude below uniform. The fix is stratified sampling — drawing a fixed quota from each class rather than sampling the pool at random — not more draws.

The value of a game with re-rolls

Roll a fair die. Keep the value in dollars, or re-roll — up to three rolls. What is the game worth?

Solve it backward. The value of choosing to re-roll is exactly the value of the shorter game that follows, and the shortest game is trivial to evaluate. So start there and work up.

1 roll left:   must accept.                          V_1 = 3.5

2 rolls left:  accept x > V_1 = 3.5, i.e. {4,5,6}
               V_2 = (1/6)(4+5+6) + (3/6)(3.5) = 2.5 + 1.75      = 4.25

3 rolls left:  accept x > V_2 = 4.25, i.e. {5,6}
               V_3 = (1/6)(5+6) + (4/6)(4.25)  = 1.8333 + 2.8333 = 4.6667

Walk each line:

The threshold is the continuation value, and it falls as rolls run out: keep 5-6 first, 4-6 second, anything third. You get pickier when you have more chances left.

That structure — value in hand versus expected value of continuing — is optimal stopping. Its general form is V_k = E[max(X, V_{k-1})], and it is the same computation behind early stopping on a validation curve, a rule for when to stop exploring and start exploiting, and an agent deciding whether one more tool call is worth its cost.

What has to be true, and how it gets misused. The backward solution assumes four things: you know the distribution of the next draw, draws are independent of one another, you are maximizing expected value with no aversion to risk, and continuing is free.

Change any one and the thresholds move. A per-roll cost lowers every threshold. Risk aversion lowers them too.

The classic misuse is applying it where the distribution is being learned from the same draws you are stopping on — early stopping on a validation set you also selected the model with, where the “continuation value” you estimated is biased upward by the selection.

Monty Hall, and the variant that proves the mechanism

Three doors, one car. You pick door 1. The host, who knows where the car is, opens door 3 on a goat and offers a switch.

Switch. The reason is that the host’s action is evidence, because his choice was constrained by where the car is.

Write the likelihoods — the probability of what he did, under each possible location of the car — assuming he picks uniformly when he is free to choose. Then apply Bayes’ theorem with a prior of 1/3 on each door.

P(opens 3 | car 1) = 1/2      he could open 2 or 3
P(opens 3 | car 2) = 1        door 3 is his only legal move
P(opens 3 | car 3) = 0        he never reveals the car

P(car 1 | opens 3) = (1/3)(1/2) / [ (1/3)(1/2) + (1/3)(1) + 0 ] = (1/6)/(1/2) = 1/3
P(car 2 | opens 3) = (1/3)(1)   / (1/2)                                       = 2/3

The three likelihood lines are the whole problem. If the car is behind your door 1, the host has a free choice and picks door 3 half the time. If the car is behind door 2, he has to open door 3 — he cannot open your door and he cannot open the car. If the car is behind door 3, he never opens it.

The denominator (1/3)(1/2) + (1/3)(1) + 0 = 1/6 + 1/3 = 1/2 is P(opens 3), and each numerator divided by it gives that door’s posterior.

Monty Fall: the variant that proves the mechanism

What matters is the protocol, not the door count. Change only the protocol and the answer changes.

Monty Fall: the host trips and knocks a door open at random, and it happens to be a goat.

The observed event is identical — door 3 is open, a goat is behind it. But the likelihoods are not, because a stumbling host was equally likely to open door 3 whichever of the other two doors hid the car:

P(opens 3 | car 1) = 1/2    P(opens 3 | car 2) = 1/2    P(opens 3, goat | car 3) = 0

P(car 1 | opens 3, goat) = (1/3)(1/2) / [ (1/3)(1/2) + (1/3)(1/2) + 0 ] = 1/2

The first two likelihoods are now equal, so the evidence does not favor either door and the posterior splits evenly.

Same doors, same goat, different answer — because the likelihood of the host’s action changed. Switching is worth nothing here. The information was never in the door; it was in the rule that generated the observation.

That is exactly why missing-data mechanisms have names (Missing values three mechanisms three different correct answers). A blank field means something entirely different depending on whether it went missing at random or because of its own value, and staring at the observed data cannot distinguish the two.

You cannot infer correctly from data without a model of how the data came to be observed.

The simulation below settles both variants empirically. It plays the game trials times and reports the win rate; host_knows=False is the Monty Fall variant, where rounds in which the stumble revealed the car are discarded rather than counted, because those rounds are not the situation being asked about.

import random


def monty(switch, host_knows, trials=200_000, seed=0):
    """Empirical win rate. host_knows=False is the Monty Fall variant."""
    rng = random.Random(seed)
    wins = counted = 0
    for _ in range(trials):
        car, pick = rng.randrange(3), 0
        if host_knows:
            opened = rng.choice([d for d in range(3) if d != pick and d != car])
        else:
            opened = rng.choice([d for d in range(3) if d != pick])
            if opened == car:
                continue                       # round voided: a car was revealed
        counted += 1
        final = next(d for d in range(3) if d not in (pick, opened)) if switch else pick
        wins += (final == car)
    return wins / counted


# host_knows=True  -> switching wins ~0.667
# host_knows=False -> switching wins ~0.500

What has to be true, and how it gets misused. The 2/3 answer requires the full protocol: the host always opens a door, always opens one with a goat, never opens your door, and always offers the switch.

Drop any clause and the answer changes. A host who only offers a switch when you picked the car makes switching a guaranteed loss.

The uniform-when-free assumption matters less than people expect. It pins the answer at exactly 2/3, but any other habit the host might have still leaves switching at least as good as staying.

The classic misuse is the general one this problem exists to teach: computing a posterior from the observed outcome alone, without modelling the process that decided which outcomes you get to observe.

Cheat sheet

Every row below is derived somewhere above. This table is for recall, not for first learning. If a row is unfamiliar, go back to the section that derived it.

Three abbreviations appear in it. “LR” is likelihood ratio. “LM” is language model. “PPV” is positive predictive value — the probability that a positive prediction is correct, which is the P(D | +) computed in Bayes theorem and the question everyone gets wrong.

FactFormWhy it earns its place
Complement ruleP(at least one) = 1 - P(none)Birthday, multiple testing, hash collisions are one move
Chain ruleP(x_1..x_n) = prod P(x_t | x_<t)Exactly what an autoregressive LM computes; no approximation
Bayes, odds formposterior odds = prior odds × LRFastest path through any base-rate question
Base-rate landmark99/99 test breaks even at 1% prevalence0.1% base rate -> 9% PPV -> 10 false alarms per true one
Linearity of EE[aX + bY] = aE[X] + bE[Y], alwaysIndicator arguments work under any dependence
Uncorrelated != independentX ~ U(-1,1), Y = X^2, Cov = 0Correlation screens delete every non-monotone feature
Variance of an averagerho·s^2 + (1-rho)·s^2/BThe random-forest floor: B cannot touch the first term
Total varianceVar(Y) = E[Var(Y|X)] + Var(E[Y|X])Bias-variance is this with X = the training draw
Bias-varianceMSE = Bias^2 + Var + sigma^2Cross term vanishes by definition of the mean prediction
Bernoulli variancep(1-p), max at 0.525× smaller at p = 0.01; drives A/B sample sizes
Poisson limitBinomial(n,p) -> Poisson(np) as n->inf, p->0E = Var = lambda is a testable claim, and usually false
MemorylessnessP(T > s+t | T > s) = P(T > t)Only the exponential. Latency is not — its hazard rises
CLT, preciselysqrt(n)(Xbar-mu)/sigma -> N(0,1) in distributionNeeds finite variance and iid; says nothing about tails
Berry-Esseen<= 0.4748·rho/(sigma^3·sqrt(n)) = 0.21 at n=30, rho = E|X-mu|^3Bootstrap instead of trusting a small-n normal tail
Log spaceunderflow at 77 tokens; m + log sum exp(x - m)Accumulate log-probabilities; subtract the max in softmax
Cross-entropy = MLE-1/N sum log q(y_i) = H(p_hat, q)Three names for one objective; derive it, do not recite it
MSE = MLEGaussian noise, fixed varianceTells you what loss to use when neither assumption holds
Softmax + CE gradientdL/dz = p - yDerived in Matrix calculus derived not memorized
KL asymmetryforward = mode-covering, reverse = mode-seekingMLE minimizes forward KL, so likelihood models hedge
Perplexityexp(cross-entropy)2.0 nats = “as uncertain as a fair 7.4-way choice”
Birthday50% at n ≈ 1.177·sqrt(N); 23 for 365It is about the C(n,2) pairs, not the n people
Coupon collectorE[T] = n·H_n ≈ n·ln n, sd ≈ 1.28·nRandom coverage costs ln n more, with a long tail
Optimal stoppingV_k = E[max(X, V_{k-1})]Early stopping, bandits, agent stop/continue
Monty Hall2/3 switching; 1/2 if the host trippedThe information is in the protocol, not the observation

Next: 02 — Statistics & Inference — what changes when you only have a sample.