07 — Imbalanced Data, Calibration & Drift

Three problems that look separate are the same problem underneath:

• what to do when the class you care about is rare,
• what to do when the probabilities your model prints are wrong,
• what to do when the live data stops resembling the data you trained on.

This chapter shows how to tell the three apart, gives the one-line correction that fixes two of them exactly, and states the assumption that correction rests on so you can recognize when it stops holding.

Every number below is derived and worked through.

Everything in this chapter concerns a classifier: a model that takes a row of features x — one transaction, one loan application, one image — and returns p(y=1|x).

Read p(y=1|x) aloud as “the probability that y equals one given x.” It is a number between 0 and 1, and it answers “how likely is it that this row belongs to the class I care about?”

The label y is 1 for the class of interest (fraud, click, disease) and 0 otherwise. By convention the y = 1 class is the positive class and y = 0 is the negative class.

That single output feeds two entirely different kinds of consumer:

• The ranking is the order the scores put the rows in. A top-k queue, a search results page, or a fixed-capacity review list uses nothing but the order.
• The probability is the number taken at face value. An expected-loss calculation, an insurance price, or a rule like “auto-block above 0.99” reads the number itself.

Which of the two your system consumes is the organizing question of the whole chapter. The two break for different reasons and are repaired by different tools.

All three topics of the title are the same topic underneath: a prior that moved. The prior here is p(y=1), the fraction of all rows that are positive before you look at any features at all — the base rate.

• Imbalance is a prior that is extreme.
• Resampling — deliberately changing the mix of positives and negatives in the training set — is a prior you moved on purpose.
• Label shift is a prior that moved on its own, with no change you made.

In all three cases the repair is the same correction in log-odds space. The log-odds of a probability p is ln(p/(1-p)), read aloud as “the natural log of p divided by one minus p.” It is the scale on which a change of prior stops being a messy rescaling and becomes plain addition.

And in all three cases the failure presents as “the model got worse” when nothing about p(x|y) has changed at all. Read p(x|y) as “p of x given y”: the distribution of feature values within a class.

In short: a classifier’s ranking and its probabilities are separate objects, damaged by separate things, and repaired by separate tools — so the first question in every section of this chapter is whether anything downstream reads the number, or only the order.

1. Imbalance is usually a threshold problem wearing a data problem’s clothes

Diagnose a class-imbalance complaint before treating it. Class imbalance means one class is far more common than the other in your data: 1% fraud and 99% legitimate is a 1:99 imbalance.

The reflex on hearing “1% positives” is to rebalance the data. First ask what actually broke. In most cases the model is fine and something downstream of it is not.

The triage

The diagram below asks two questions in order.

The first question is how many positive rows you have in absolute terms. That separates a genuine shortage of data from a mere ratio.

The second question, asked only if you have plenty of positives, is which specific downstream component is misbehaving. Each of its four branches names its own fix.

flowchart TD
    S["1% positive class"] --> Q1{"How many positives<br/>in absolute terms?"}
    Q1 -->|"under ~1,000"| DATA["A real DATA problem.<br/>The model cannot estimate<br/>the minority conditional at all.<br/>Get labels, add features,<br/>use a simpler model."]
    Q1 -->|"thousands or more"| Q2{"What is broken?"}
    Q2 --> M1["Metric is accuracy<br/>-> switch metric<br/>ml/06 section 5"]
    Q2 --> M2["Decisions use argmax<br/>at p = 0.5<br/>-> derive threshold from cost<br/>ml/06 section 11"]
    Q2 --> M3["Optimizer is drowned by<br/>easy negatives<br/>-> class weights or focal loss"]
    Q2 --> M4["Nothing. Log loss already<br/>fits a low-prior model fine.<br/>-> do nothing"]

    style DATA fill:#9d0208,color:#fff
    style M2 fill:#2d6a4f,color:#fff
    style M4 fill:#1d3557,color:#fff

The second question, “what is broken?”, has four answers:

• The metric is accuracy. Switch to a metric that is not dominated by the majority class, as worked through in The confusion matrix is the root everything else is a ratio of its cells.
• Decisions use argmax at p = 0.5. Argmax means “take whichever class has the larger score,” which for two classes is exactly the same thing as thresholding the probability at 0.5. The fix is to derive the threshold from the cost of each kind of mistake, as in Choosing a threshold from the cost matrix.
• The optimizer is drowned by easy negatives. The training signal is dominated by the enormous number of obviously-negative rows. Fix it with class weights or focal loss, both derived later in this chapter.
• Nothing is broken. Log loss already fits a low-prior model perfectly well. The correct action is to do nothing.

Relative rarity and absolute rarity are different problems

Only one of them is about imbalance.

1% of 10,000,000 rows is 100,000 positives. That is more minority examples than most datasets have rows, and nothing is wrong.

1% of 2,000 rows is 20 positives, and no resampling technique creates information that 20 examples do not contain. SMOTE is the Synthetic Minority Over-sampling Technique: it manufactures new minority rows by interpolating between existing ones, and Resampling and the calibration it breaks dissects it in detail. Run on 20 points, SMOTE produces convex combinations of 20 points. A convex combination is any point you can reach by sitting somewhere on the straight line between two of them.

Where the diagram’s ~1,000 boundary comes from

It is not a constant anyone measured. It is roughly where the minority conditional — the distribution of feature values among the positive rows, p(x|y=1) — acquires enough support to estimate anything at the feature count in play. Support means “enough observed examples, spread over enough of the feature space, that the estimate is not mostly guesswork.”

Two independent readings land on the same order of magnitude.

Reading one: events per variable. Logistic regression is the standard linear model for two-class problems, and it fits one coefficient per feature. The events-per-variable rule of thumb asks for 10–20 positive events per fitted coefficient. A 50-feature model therefore wants 500–1,000 positives.

Reading two: how noisy your metric is. Sampling noise on any minority-class metric scales as 1/sqrt(n_pos), where n_pos is the number of positives you have.

The block below reads that scaling for recall, the fraction of true positives your model actually catches. Each line gives recall’s standard error — the typical amount a measurement varies purely because you measured it on a finite sample rather than on the whole population — and then doubles it, because plus-or-minus two standard errors is the conventional 95% interval. The right-hand column falls from 22 points to 0.3 as n_pos grows:

recall of 0.60, standard error sqrt(0.6 x 0.4 / n_pos), reported at 2 SE:

  n_pos =      20    SE = 0.1095   ->  +- 21.9 points
  n_pos =     100    SE = 0.0490   ->  +-  9.8 points
  n_pos =     200    SE = 0.0346   ->  +-  6.9 points
  n_pos =   1,000    SE = 0.0155   ->  +-  3.1 points
  n_pos = 100,000    SE = 0.0015   ->  +-  0.3 points

Below roughly a thousand positives you cannot measure a change you make, never mind learn one. At 20 positives the error bar on recall is wider than the entire range of recalls you would argue about. At the 100,000 positives of the 1%-of-10M case it is three tenths of a point.

That is the scale the boundary marks. It moves with your feature count and with how much metric resolution you need, so read ~1,000 as an order of magnitude, not a test.

The model is fine; what breaks is downstream

A 1% positive rate is not by itself a problem. The next question is whether the model is even damaged by it. It is not.

Take a logistic regression, or a GBDT — a gradient-boosted decision tree, an ensemble that adds up many small decision trees, each fitted to the errors of the ones before it. Train either with plain log loss and it converges to a good estimate of p(y|x): a correct estimate that happens to be small nearly everywhere.

Log loss, also called cross-entropy, charges you -ln(p) when the true label is positive and -ln(1-p) when it is negative. So it punishes confident mistakes savagely and is minimized by telling the truth.

That last property has a name. Log loss is a proper scoring rule: the prediction that minimizes it in expectation is the true probability itself and nothing else (Proper scoring rules log loss and brier). Its minimizer — the prediction that makes the loss as small as it can go — is the true conditional p(y|x), whatever the prior is.

So the model is fine. What breaks is downstream. The table below takes the five complaints people file under “the imbalance problem” and names what each one actually is. Only the last row is a problem with the data:

Five terms in that table deserve a definition before you meet them again.

• Prevalence is the fraction of rows that are positive — 1% in the running example. A metric being prevalence-weighted means its value depends on that fraction rather than only on the model.
• PR-AUC is the area under the precision-recall curve. It ignores the true negatives entirely, so it does not flatter a model for correctly dismissing 99% of the data.
• MCC, the Matthews correlation coefficient, is a single number between -1 and +1 that uses all four cells of the confusion matrix. It is only high when the model does well on both classes.
• Expected cost is what you get by assigning a dollar figure to each kind of mistake — C_fp for a false positive, C_fn for a false negative — and adding up what the model’s errors would actually bill you.
• Anomaly detection, in the last row, means giving up on learning the positive class from examples. Instead you model only what “normal” looks like and flag whatever falls far outside it. That is the honest move when 20 positives is all you will ever have.

Expected cost is where the optimal threshold C_fp/(C_fp+C_fn) comes from. It is the cutoff at which one more false positive costs exactly what one more false negative costs.

What interviewers probe: “You have a 1:99 imbalance — what do you do?” The answer that separates candidates is “first I check whether I have an imbalance problem at all.” Count absolute positives, name the metric, name the decision threshold and where it came from. Reaching straight for SMOTE is the junior answer.

2. Resampling, and the calibration it breaks

What happens to a model’s output when you change the class mix it trains on?

Resampling means altering the composition of the training set on purpose — duplicating minority rows, discarding majority rows, or inventing synthetic minority rows — so that the two classes appear in a ratio you chose rather than the one nature supplied.

The consequence, derived in full below: resampling leaves the model’s ranking untouched and destroys its probabilities by a fixed, computable amount.

Calibration is the property that gets destroyed. A model is calibrated when its stated probabilities match observed reality, so that among the rows it labels “70% likely,” about 70% really are positive. Calibration what it means and when it matters defines calibration properly. Here you only need to know that resampling breaks it, and by exactly how much.

Walking the argument before the algebra

The diagram below is the whole argument in one picture. Walk it now; the derivation later fills in the algebra.

Box 1 to box 2. You train on a rebalanced set whose prior is pi' (read “pi prime”) instead of the true pi. All you changed was how many rows of each class you kept, so p(x|y) is unchanged and only the prior moved.

Box 2 to box 4. That makes the posterior odds multiplied by a constant c. Posterior odds means p/(1-p) computed after seeing the features. Taking logarithms turns that multiplication into logit(p) = logit(p') + ln(c), a pure intercept shift. The logit of a probability is that same log-odds ln(p/(1-p)) from the introduction, and the intercept is the constant term added to every prediction.

Box 4 splits two ways, and the two consequences point in opposite directions.

The first is that AUC is unchanged. A constant shift is a monotone map — one that stretches the numbers without ever reordering them — and AUC reads only the ordering.

AUC is the area under the ROC curve. ROC stands for receiver operating characteristic, a name inherited from wartime radar and safely ignorable. The quantity itself has an exact plain-English meaning: the probability that a randomly chosen positive row scores higher than a randomly chosen negative one (Roc auc is a probability and here is the derivation). It knows nothing about the face value of the scores.

The second is that every probability is wrong. A printed 0.90 really means 0.083, which the worked table below reproduces.

The repair is to add ln(c) to every logit. It is exact and it costs one constant.

The dashed branch is the exception that matters. SMOTE breaks the assumption in box 2, because it changes p(x|y=1) itself. No single constant can undo it, so you must fit a calibrator instead. A calibrator is a small extra model, fitted on held-out data, whose only job is to map the raw score to an honest probability.

flowchart LR
    A["Train on a rebalanced set<br/>prior pi' instead of pi"] --> B["p of x given y unchanged<br/>only the prior moved"]
    B --> C["posterior ODDS multiplied<br/>by a constant c"]
    C --> D["logit p = logit p' + ln c<br/>a pure INTERCEPT shift"]
    D --> E["AUC unchanged<br/>monotone map, ranks preserved"]
    D --> F["Every probability wrong<br/>0.90 really means 0.083"]
    F --> G["Add ln c to every logit<br/>-- exact, one constant"]
    B -.->|"SMOTE breaks<br/>this assumption"| H["p of x given y=1 CHANGED<br/>-> no constant undoes it<br/>-> must fit a calibrator"]

    style D fill:#1d3557,color:#fff
    style E fill:#2d6a4f,color:#fff
    style F fill:#bc6c25,color:#fff
    style H fill:#9d0208,color:#fff

The three moves, the one non-move, and what each distorts

People actually do four things about imbalance, and each introduces its own specific distortion.

The heading says three moves and one non-move because the first three rows change the data and the fourth changes only the loss function. Read the Distorts column as the price of each. The fourth row’s price is the punchline two headings below: class weights are the exact objective that the other three only approximate.

The vocabulary in that table, unpacked:

• The minority class is whichever class is rarer — the positives, in every example here. The majority class is the common one.
• A high-variance learner is a model flexible enough to fit the training rows almost exactly. That is what makes duplicate rows dangerous: a decision tree can isolate a row it saw 40 times into its own leaf (a terminal node of the tree, holding one final prediction for everything that reaches it) and score it perfectly. Training loss falls while nothing is learned that transfers to new data.
• SMOTE’s k nearest minority neighbours are the k other positive rows closest to a given positive row in feature space. The synthetic point is placed at a random position on the line segment joining the row to one of them.
• w_c in the last row is a number multiplying the loss contribution of every row belonging to class c, so the optimizer treats one rare row as if it were w_c rows.

Every one of those four techniques rests on an assumption, and each fails differently when it is false.

Random oversampling assumes duplicated rows carry the same information as fresh ones. That is false for any model that can memorize. It survives with strongly regularized linear models and fails with deep trees and neural networks.

Random undersampling assumes the discarded negatives were redundant. That is true when you have millions of them and false the moment the negative class has structure — rare-but-legitimate patterns you threw away come back as false positives.

Class weights assume the only thing wrong is the relative emphasis of the two classes. That is exactly right, and it is why the technique distorts nothing structural.

SMOTE makes the strongest assumption of the four, which the rest of this subsection takes apart.

SMOTE’s assumption is that the segment between two minority points is minority territory, and there are three standard ways that is false:

1. A non-convex or multi-modal minority class. Multi-modal means the positives form two or more separate clusters rather than one blob, and non-convex means the straight line between two positives can leave positive territory. Card-testing fraud lives at small amounts; account-takeover fraud lives at large ones. Interpolating between the two clusters manufactures “medium-amount fraud” that does not exist, and the model learns a decision region straddling ordinary purchases. False positives rise in the region SMOTE invented.

2. Categorical features. A categorical feature takes one of a fixed set of unordered values, such as a country code. It is usually one-hot encoded: one column per possible value, holding 1 for the value the row has and 0 for all the others.

   Interpolating those columns produces country_US = 0.5, which corresponds to no country at all. It sits off the data manifold — the thin region of the feature space where real rows actually live.

   SMOTE-NC, the variant for nominal-and-continuous data, substitutes the most frequent category among the neighbours instead. That erases exactly the rare-category structure you were trying to model.

3. High dimensions. As the number of features p grows, all nearest neighbours become nearly equidistant from any given point (Knn and the curse of dimensionality), so “nearest minority neighbour” degenerates toward “random minority point,” and SMOTE collapses into random convex combinations of the whole class — which contracts the minority cloud toward its centroid, the average position of all the positives, and reduces recall on the tails you cared about.

The fourth failure is procedural rather than statistical, and it shows up in every code review.

Cross-validation (CV) is the standard way to estimate how a model will do on new data: split the rows into folds, train on all but one, score on the held-out one, and repeat.

The block below compares the CV estimate against a truly untouched holdout set for three setups. The 0.274 gap in the first row is the size of the error:

                             CV AUC    holdout AUC
SMOTE applied, then split      0.987       0.713
split, then SMOTE in-fold      0.844       0.826
class weights, no resampling   0.851       0.839

The first row’s 0.987 is spurious, and the mechanism is simple.

A synthetic point is a convex combination of two real minority points, so it carries information about both of them. Oversample before splitting and a validation row ends up partially present in training.

That is a leak: evaluation data influencing the fitted model, so the reported score measures memory rather than generalization. Its severity grows with the oversampling ratio.

Any resampling, weighting, or calibration step belongs strictly inside the cross-validation fold.

The calibration break, derived

Exactly how wrong are a resampled model’s probabilities? The derivation from first principles is short and sets the exact size of the error.

The one assumption. Resampling changes the class prior from pi to pi' while leaving the class-conditional densities f_1(x) and f_0(x) untouched.

A class-conditional density is the distribution of feature values within one class. f_1 describes what fraudulent transactions look like; f_0 describes what legitimate ones look like. Keeping or dropping whole rows at random changes how many of each class you have without changing what either class looks like.

That holds for random oversampling, random undersampling and class weights. It is emphatically not true for SMOTE.

Under that assumption, write the posterior odds under each prior. Three lines, then a division:

p(y=1|x)  =  pi·f_1(x) / (pi·f_1(x) + (1-pi)·f_0(x))

odds:      o  = p/(1-p)   = [ pi  / (1-pi ) ] · f_1(x)/f_0(x)      true world
           o' = p'/(1-p') = [ pi' / (1-pi') ] · f_1(x)/f_0(x)      resampled world

divide:    o  = o' · [ pi/(1-pi) ] · [ (1-pi')/pi' ]
                     \___________________________/
                              call this c

Read that block one line at a time.

Line 1 is Bayes’ rule for two classes: the probability that a row is positive equals the prior times the positive class’s density, divided by the same quantity summed over both classes.

Lines 2 and 3 convert probability to odds. The point of doing so is that the denominator cancels, and what survives is the ratio f_1(x)/f_0(x). That ratio is called the likelihood ratio and means “how much more typical this row is of positives than of negatives.” It appears identically in both worlds, because resampling did not touch f_1 or f_0.

Line 4 divides the true-world odds by the resampled-world odds. The likelihood ratio cancels entirely and nothing is left but prior terms.

That leftover is c. The critical property is that c does not depend on x — it is the same number for every row in the dataset, which is why one constant can undo the damage.

Take logs and the whole correction collapses to one line:

logit(p) = logit(p') + ln(c)

Resampling shifts the log-odds by a constant. It is an intercept bug, not a model bug.

That is why applying the correction leaves AUC exactly unchanged: a constant shift in logit space is a strictly increasing map, and AUC reads only ranks (Roc auc is a probability and here is the derivation). Meanwhile it destroys every probability the model emits.

One thing not to confuse: applying the correction does not move AUC, but resampling itself is a different fit on less data, and its AUC does move — which is what the CV trace above shows.

Worked: a model trained on a 50/50 balanced set

This is the commonest case. The true prevalence is pi = 0.01. Rebalancing to equal class sizes makes pi' = 0.5. Substituting into c = [pi/(1-pi)]·[(1-pi')/pi']:

c    = (0.01 / 0.99) · (0.5 / 0.5) = 0.010101
ln c = -4.5951

Now run five model outputs through it. Convert each p' to odds, multiply the odds by c = 0.010101, and convert back. The right-hand column is what the model should have said:

Take the third row line by line. The model says p' = 0.900. Odds are 0.900/0.100 = 9.00. Corrected odds are 0.010101 × 9.00 = 0.09091. Back to a probability: 0.09091/1.09091 = 0.0833.

A model trained on balanced data that says “90% confident” means 8.3%. Feed that 0.90 into an expected-loss calculation, an insurance price, or a threshold of 0.5 derived from a cost matrix — the small table stating what a false positive and a false negative each cost you — and every downstream number is wrong by an order of magnitude.

Note the break-even row at p' = 0.99 -> p = 0.50. The balanced model must be 99% sure before the real-world posterior even reaches a coin flip.

Worked: negative downsampling, the case you will actually meet

In negative downsampling you keep every positive row and keep each negative row with probability w. A w of 0.05 means you retain one negative in twenty.

That is a prior move like any other, and it is the move every large-scale system makes. The reason is compute, not statistics: keeping all the negatives would mean training on tens of billions of rows for no statistical gain.

Substituting pi' = pi / (pi + w(1-pi)) into c makes the whole expression collapse. Follow the third line, where the 0.01 in the numerator meets the 0.01 in the denominator:

pi = 0.01,  w = 0.05

pi'  = 0.01 / (0.01 + 0.05 x 0.99) = 0.01 / 0.0595 = 0.16807
       rows kept = 0.0595 N   ->  a 16.8x smaller training set

(1 - pi')/pi'  =  0.05 x 0.99 / 0.01  =  4.95
c              =  (0.01/0.99) x 4.95  =  0.05000        <-  exactly w

logit(p) = logit(p') + ln(w) = logit(p') - 2.9957

c collapses to w itself. The pi cancels, so the offset is ln(w) at every prevalence.

That is the practically important result of the whole section: you do not need to know the true prevalence to fix a downsampled model, only the sampling rate you chose — which you always know, because you chose it.

Worked on one score: a downsampled model emitting p' = 0.95 carries odds 0.95/0.05 = 19.0. Corrected odds are 19.0 × 0.05 = 0.95, so p = 0.95/1.95 = 0.487. Half what it said.

Extreme imbalance negative downsampling and the correction it forces runs exactly this through an ad auction — same w = 0.05, same -2.9957 — and prices the missing offset at a 1.95x inflation of the ad’s eCPM, the effective cost per thousand impressions, which is the currency an ad auction ranks on.

That chapter also gives the Fisher-information argument for why 16.8x less data costs only 9% in standard errors. Fisher information is the formal measure of how much a sample tells you about a parameter, and the point is that the discarded rows were the least informative ones.

Extending it from two classes to K classes

The derivation never used the fact that there were only two classes. It used only the statement that the class-conditional densities did not move, and that statement makes sense for any number of classes.

Applying the same odds argument to every class against a common reference class gives a per-class offset on the logits. A logit here means the raw pre-softmax score z_k the network emits for class k. Softmax is the function that turns those raw scores into probabilities that sum to one, by exponentiating each and dividing by the total.

z_k_corrected  =  z_k + ln( pi_new_k / pi_old_k )        then re-softmax

The correction reads: add the log prior ratio to each class’s logit, then re-softmax.

Here it is worked on three classes. The model was trained where the class priors were (0.5, 0.3, 0.2) and deployed into a market where they are (0.2, 0.3, 0.5). On some particular row it emits the probabilities (0.50, 0.20, 0.30).

Multiply each class’s probability by the ratio of its new prior to its old one, then divide by the sum to renormalize. Read the block downward, one row per step:

prior ratio      0.2/0.5 = 0.4      0.3/0.3 = 1.0      0.5/0.2 = 2.5
log ratio          -0.91629            0.00000            +0.91629
reweighted     0.50x0.4 = 0.200   0.20x1.0 = 0.200   0.30x2.5 = 0.750
                                                       sum      = 1.150
corrected           0.1739             0.1739             0.6522

The argmax moved from class 1 to class 3. Before the correction the model’s best guess was class 1 at 0.50. Afterwards it is class 3 at 0.65.

The two-class formula from earlier is just the K = 2 case of this one, after the constant shared by both classes cancels.

That argmax move is the fine print on the “argmax over classes needs no calibration” row in the table at the end of When calibration matters and when it is irrelevant. Softmax is invariant to adding the same constant to every logit, and temperature scaling (defined in Platt scaling vs isotonic regression) is monotone, so neither of those moves the argmax. But a prior correction adds a different constant per class, and that does move it.

Prior correction is not calibration, and you owe it whether or not anything downstream reads the number.

Both forms of the correction, as code. prior_shift_correction takes a probability and returns a probability; logit_offset returns the constant you add if you are working in logit space:

from math import log, exp

def prior_shift_correction(p_resampled, pi_true, pi_train):
    """Undo a change of class prior. Exact when p(x|y) is unchanged.

    Valid for random over/undersampling and class weights.
    NOT valid for SMOTE, which alters p(x|y=1) itself.
    """
    c = (pi_true / (1 - pi_true)) * ((1 - pi_train) / pi_train)
    odds = p_resampled / (1 - p_resampled)
    return (c * odds) / (1 + c * odds)


def logit_offset(pi_true, pi_train):
    """The same correction as a constant to add to every logit."""
    return log((pi_true / (1 - pi_true)) * ((1 - pi_train) / pi_train))

The SMOTE caveat is the practical payoff of the derivation.

The formula assumed f_1(x) was unchanged. That assumption is what makes c a single constant rather than a function of x.

SMOTE violates it by construction, because the synthetic points are new probability mass in places where no real positive ever sat.

So SMOTE’s miscalibration varies from row to row. It is not a constant in logit space and cannot be undone by an intercept. It has to be repaired by fitting a calibrator on untouched held-out data, which is the subject of Calibration what it means and when it matters.

Class weights are the exact objective; resampling is a noisy estimate of it

Class weights and resampling are the same statistical operation — weights compute it exactly, resampling only estimates it, and both therefore need the identical calibration fix.

Why weights are exact and resampling is not. Weighting each positive row by w in the loss is, in expectation, identical to replicating that row w times.

But replication samples which rows get duplicated, and undersampling discards rows outright. Both introduce randomness that the weighted version does not have.

The weighted objective is computed exactly, with every distinct row contributing and no Monte Carlo noise in between. Monte Carlo is the general name for estimating a quantity by random sampling instead of computing it directly — fine on average, noisy on any particular run.

What scikit-learn’s “balanced” setting actually does. The block below substitutes a 1% dataset into its formula. The ratio of the two weights is 99, which is exactly (1-pi)/pi.

sklearn "balanced": w_c = n / (k · n_c)
  n = 10,000, n_1 = 100, n_0 = 9,900, k = 2
  w_1 = 10,000 / 200   = 50.00
  w_0 = 10,000 / 19,800 =  0.5051
  ratio = 99  =  (1 - pi) / pi          <- exactly the ratio that makes pi' = 0.5

In that block n is the total number of rows, n_c the number in class c, and k the number of classes. So class_weight="balanced" gives each class a weight inversely proportional to how common it is.

That makes class_weight="balanced" and “resample to 50/50” the same statistical operation. They break calibration in exactly the same way and need exactly the same correction, with pi' = w·pi / (w·pi + (1 - pi)) in general. XGBoost’s scale_pos_weight parameter is the same knob under a different name.

Two traps with the letter w

Trap 1: this is a different w from the one in the previous subsection. They collide because each is, in its own setting, the single number the whole correction depends on.

• In negative downsampling, w was a keep rate: a probability below 1 that shrinks the negative class. It shifts every logit by +ln(w), and since w < 1 that is a negative offset.
• Here, w is an up-weight ratio: a number above 1 that inflates the positive class. It shifts every logit by -ln(w), and since w > 1 that is also a negative offset.

The two are reciprocals of each other as prior moves — c = w in the first case, c = 1/w in the second — which is exactly why one takes +ln and the other -ln. Carry the wrong one and you get the offset backwards, which doubles the error rather than removing it.

Whenever you see w, ask first whether it multiplies rows or multiplies loss.

Trap 2: w is the positive-to-negative weight ratio, not either weight on its own. The block above prints two numbers and only their quotient matters, because scaling both weights by any constant scales the whole loss and moves nothing about where its minimum sits.

Feed the salient-looking number, w_1 = 50.00, into pi' = w·pi/(w·pi + 1 - pi) and you get pi' = 0.336 — not the 0.5 that the balanced setting is supposed to produce. Feed the ratio w = w_1/w_0 = 99 and it lands exactly on 0.500.

The block below works the balanced case, then a deliberately unbalanced one at w = 12, to show that the collapse to c = 1/w is not special to 50/50:

w = 99, pi = 0.01:   pi' = 0.99 / (0.99 + 0.99)  = 0.500       <- the balanced case
w = 12, pi = 0.01:   pi' = 0.12 / (0.12 + 0.99)  = 0.12/1.11 = 0.10811

  c    = (0.01/0.99) x (0.89189/0.10811) = 0.083333  =  1/12
  ln c = -2.4849                                     = -ln(12)

So c = 1/w exactly, at any prevalence. It is the same cancellation as the downsampling case, running the other way. Downsampling negatives at rate w shifts every logit by +ln(w); up-weighting positives by w shifts them by -ln(w).

A model trained with scale_pos_weight = 12 needs -2.485 added to every logit, and you never have to know pi to apply it.

Which of the four to pick

The table below puts the four techniques in order of preference. It is the decision you should actually make.

Its third row mentions bagging, short for bootstrap aggregating: train several copies of the same model on different random subsets of the data and average their outputs, which cancels much of the noise any single copy has. BalancedRandomForest and EasyEnsemble are two library implementations of that idea applied to undersampling specifically.

Focal loss — rebalancing by difficulty instead of by class

One rebalancing tool remains, and it works on a different axis entirely: difficulty rather than class membership.

Focal loss is a modified log loss that automatically pays less attention to examples the model already gets right. The training signal ends up dominated by the hard cases rather than by the numerous easy ones.

Its formula introduces p_t, read “p sub t.” That is a bookkeeping trick meaning “the probability the model assigned to the correct answer for this row.” So p_t is high when the model is right and low when it is wrong, whichever class the row belongs to:

FL(p_t) = -alpha_t · (1 - p_t)^gamma · log(p_t)        p_t = p if y=1 else 1-p

Reading that formula, right to left:

• -log(p_t) is ordinary log loss.
• alpha_t is an optional per-class weight of the kind the previous subsection covered.
• (1 - p_t)^gamma is the new part. This factor is near zero when p_t is near 1 and near 1 when p_t is near 0, so it shrinks the loss of easy examples and leaves hard ones almost alone.
• gamma (the Greek letter gamma) controls how aggressively it does that. gamma = 0 recovers plain log loss; gamma = 2 is the standard setting.

The table below tabulates that factor at gamma = 2. The third column is 1/(1-p_t)^2 — how much each row’s loss is shrunk. The fourth column divides each row’s factor by the bottom row’s 0.81:

The fourth column is the one that answers the question focal loss exists to answer. Not “how much is this example down-weighted” but “how much less does an easy example matter than a hard one.”

At gamma = 2 an example at p_t = 0.9 contributes 81x less than one at p_t = 0.1. One at p_t = 0.99 contributes 8,100x less.

That column is a direct multiplier on the loss split. Easy examples only keep control of the gradient while they outnumber the hard ones by more than that ratio.

Worked: dense object detection, the problem focal loss was invented for

A dense detector does not propose a handful of candidate regions and then classify them. It lays a fixed grid of roughly 100,000 candidate boxes — anchor boxes — over every image and asks of each one “is there an object here?”

About 10 of those 100,000 contain an object. The rest — 99,990 exactly, which the arithmetic below rounds to 100,000 since the ten make no difference to a four-digit total — are obviously empty background.

Suppose the easy background boxes sit at p_t = 0.99, meaning the model already assigns 99% probability to “empty,” and the positives sit at p_t = 0.10, meaning the model is still mostly wrong about them. Compute the total loss contributed by each group under both losses:

cross-entropy:
  100,000 easy negatives x -ln(0.99) = 100,000 x 0.01005 = 1,005.0
       10 hard positives  x -ln(0.10) =      10 x 2.30259 =    23.0
  ratio  ->  background owns 43.6x more of the loss than the objects do

focal, gamma = 2:
  100,000 x 0.0001 x 0.01005 =   0.10
       10 x 0.8100 x 2.30259 =  18.65
  ratio  ->  objects now own 186x more of the loss than the background

The gradient signal flips from 44:1 background to 186:1 foreground without touching the data. The gradient is the direction the optimizer moves the weights on each step, and it is a sum over rows, so whichever rows contribute most of the loss also steer most of the learning.

The two ratios are connected by exactly the modulation factor from the table: 43.648 / 8,100 = 1/185.6. Focal loss divided the background’s share by the 8,100x from the fourth column, and that one division is the whole flip.

Focal loss’s own assumption is that your majority-class rows are overwhelmingly easy, not merely numerous. That is why it belongs to dense detection and rarely helps on a tabular 1:99 problem. Tabular negatives are genuinely hard to tell apart from positives; they are not sitting at p_t = 0.99. Focal loss rebalances difficulty, not counts, so when difficulty is not the problem it does nothing useful.

The cost, stated precisely: focal loss is not a proper scoring rule. Its minimizer is not p(y|x), so its output is a score to be ranked or thresholded, not a probability to be multiplied by money.

It happens to pull deep networks toward higher entropy — less peaked, more hedged output distributions — which usually reduces the overconfidence described in Why modern deep nets are overconfident the mechanism. But that is a coincidence of direction and not a guarantee. If you need probabilities, calibrate afterwards.

When to do nothing

More often than the literature suggests, the correct action on an imbalanced dataset is no action.

Do nothing when all four of these hold:

1. The metric is proper (log loss or Brier) or rank-based (AUC, PR-AUC).
2. The threshold came from a cost matrix rather than defaulting to 0.5.
3. Positives number in the thousands.
4. The downstream system consumes either a rank or a probability you have already corrected.

That is a large fraction of real problems. Resampling adds a distortion and a correction step; if you were not going to be wrong without it, you are adding two ways to be wrong.

The one thing “do nothing” does not excuse you from

Adding ln(c) to every logit and subtracting ln(c) from the decision threshold in logit space are the same decision expressed two ways. The rule logit(p) > t and the rule logit(p') > t - ln(c) select exactly the same rows.

So a rank-only consumer does not escape this section’s correction. It applies the correction to its cutoff instead of to its scores.

Work it. A cost matrix gives an optimal threshold of 0.20 on true probabilities. On a balanced-trained model at pi = 0.01, that cutoff becomes something very different:

ln c            = -4.59512
logit(0.20)     = -1.38629
shifted         = -1.38629 - (-4.59512) = 3.20883
threshold on p' = sigmoid(3.20883) = 0.9612

The sigmoid in that last line is the function that undoes a logit, turning a log-odds number back into a probability: sigmoid(z) = 1/(1 + e^-z).

Threshold the balanced-trained model at 0.961, not at 0.20. Same rows, no probabilities touched.

The two cutoffs differ by a factor of exactly 99 in the odds — 0.20 has odds 0.25, 0.9612 has odds 24.75, and 24.75/0.25 = 99, which is 1/c. Take the 0.20 at face value and you flag almost everything.

3. Calibration: what it means, and when it matters

Everything in §2 traded on calibration — the property resampling destroys, the repair SMOTE forces — without ever pinning it down. Start with the definition.

A model is calibrated if, among all inputs it assigns probability q, exactly a fraction q are positive.

Concretely: gather every row the model scored at 0.70 and count how many were actually positive. Calibration says that count should be 70% of them.

for all q in [0,1]:   P(Y = 1 | p_hat(X) = q) = q

Read that line aloud as: for every probability q between 0 and 1, the probability that the true label Y is 1, given that the model’s estimate p_hat(X) equals q, is itself q. The hat on p_hat is standard notation for “an estimate of,” as opposed to the unknown true quantity.

Calibration is a weak property on its own, and that is the first thing to internalize. A model that outputs the base rate 0.03 for every input is perfectly calibrated and completely useless.

The reliability - resolution + uncertainty decomposition of the Brier score is what makes that precise (Proper scoring rules log loss and brier). The Brier score is the mean squared difference between the predicted probability and the 0-or-1 outcome, and like log loss it is a proper scoring rule. Its three terms:

• reliability is miscalibration. You want it small.
• resolution is how far the model’s predictions spread away from the base rate. You want it large.
• uncertainty is a property of the data alone. You cannot change it.

A constant predictor has reliability = 0 and resolution = 0, so it looks flawless on the calibration term and contributes nothing. Calibration is one of three terms; never report it without a discrimination number beside it.

Discrimination is the other axis, and it is what that word names: the model’s ability to rank positives above negatives, measured by ROC-AUC, PR-AUC, or the resolution term.

Calibration is about whether the number is right. Discrimination is about whether the ordering is. A monotone recalibration can fix the first and cannot touch the second, and no amount of the first substitutes for the second.

Reliability diagrams and ECE, with numbers

In practice calibration is measured with a picture and a number, and both leave out more of the story than they tell.

A reliability diagram, also called a calibration curve, is the standard picture. Sort the predictions into bins by their predicted value. In each bin, plot the mean predicted confidence on one axis against the fraction actually positive on the other.

A perfectly calibrated model traces the 45-degree line. A bin sitting below the line is overconfident; one above it is underconfident.

ECE, the expected calibration error, is that picture collapsed into one number: average the vertical gaps, weighted by how many rows each bin holds.

The table below holds 10,000 predictions from an overconfident network. n_b is the number of rows in bin b, “mean conf” is the average probability the model assigned in that bin, and “observed acc” is the fraction that were actually positive. Every gap is positive, meaning every bin is overconfident — that fact does real work two subsections from now:

ECE = sum over bins of (n_b / n) · |acc_b - conf_b| = 0.1011
MCE = max over bins of |acc_b - conf_b|            = 0.125   (the 0.8-0.9 bin)
mean confidence - accuracy = 0.898 - 0.797         = 0.101

Check one row by hand: the 0.9-1.0 bin holds 6,300 of the 10,000 rows, its gap is 0.978 - 0.874 = 0.104, and its contribution is (6300/10000) × 0.104 = 0.06552 — nearly two thirds of the total ECE.

MCE in that block is the maximum calibration error: the single worst bin rather than the weighted average.

When miscalibration runs in one direction, ECE equals |mean confidence - accuracy|. That scalar takes one line of code and requires no binning decision at all, which makes it the fastest overconfidence check you can run.

MCE is a weaker signal. Being the worst bin, it is usually dominated by whichever bin is sparsest, so it reports noise as often as it reports a problem.

The same 10,000 rows, decomposed into all three Brier terms

ECE summarises only one of the three, and running all three over this table shows how little of the story that one term carries.

Take bin b to hold n_b rows at mean confidence p_b against observed rate o_b. Write pbar (read “p bar”) for the overall positive rate across all 10,000 rows:

pbar        = 0.797211                              <- the accuracy total
uncertainty = pbar(1 - pbar)                         = 0.161666
reliability = sum_b (n_b/n)(p_b - o_b)^2             = 0.010585
resolution  = sum_b (n_b/n)(o_b - pbar)^2            = 0.012185

Brier = reliability - resolution + uncertainty
      = 0.010585 - 0.012185 + 0.161666               = 0.160065

reliability is the weighted mean squared gap; ECE is the weighted mean absolute gap. The two say the same thing in different units: sqrt(0.010585) = 0.1029 against ECE = 0.1011. They nearly coincide only because every bin here errs in the same direction.

Now make the comparison ECE cannot make: against the dumbest possible model.

A constant predictor emitting 0.797 on every row scores Brier 0.161666. That is the uncertainty term on its own, since a constant has zero reliability and zero resolution.

This model, with its five distinct confidence levels, scores 0.160065. The gap is 0.161666 - 0.160065 = 0.001601 — it beats a constant by 1%.

Recalibrate it perfectly (set each bin’s output to its own observed rate) and reliability drops to zero, landing Brier at 0.161666 - 0.012185 = 0.149480. That is a 6.6% gain, and the resolution = 0.012185 that remains is the sum total of what the model knows.

Recalibrating this model buys 0.010585 of Brier. Owning the model at all rather than a constant buys 0.001601 — a factor of 6.6 smaller.

The confidence numbers are badly wrong and the ranking is barely there. No ECE value at any bin count could have told you the second half of that sentence.

ECE’s binning sensitivity, and why it is a lower bound

ECE is a property of the binning choice you made, not a fixed property of the model, and the error that choice introduces always runs in the same direction. Consider first when the binning choice can move the number at all.

Two binning schemes to know. Equal-width binning cuts the interval [0,1] into slices of the same size, so bins can hold wildly different numbers of rows. Equal-mass binning cuts it so that every bin holds the same number of rows, so the slices have different widths.

Run the same 10,000 predictions through four different schemes. All four return the identical number:

 5 equal-width bins   ECE = 0.1011
10 equal-width bins   ECE = 0.1011
15 equal-width bins   ECE = 0.1011
10 equal-mass bins    ECE = 0.1011

That invariance is the point. Here is why it happens.

Merge two bins A and B. The merged bin reports |w_A·g_A + w_B·g_B|, where g is a bin’s signed gap and w its share of the rows. The two separate bins had reported w_A·|g_A| + w_B·|g_B|.

Those two expressions are equal whenever the two gaps carry the same sign — because then the absolute value on the outside can be pushed inside without changing anything.

Every bin of that model is overconfident, as the reliability table showed, so every gap is positive. No regrouping of them can move the total off

|mean confidence - accuracy|  =  0.898 - 0.797  =  0.1011

That is the same one-line shortcut the previous subsection offered as the fastest overconfidence check, now doing double duty as a proof. When the error runs one way everywhere, ECE is that scalar — and a scalar has no bin count in it.

No rebinning of an everywhere-overconfident model can change its ECE. Coarsen these five bins all the way down to one and you still get 0.1011. Cut them finer and every piece is still overconfident, so you still get 0.1011.

Run the chapter’s own ece() below at n_bins = 5 and n_bins = 10 on these predictions and it returns the same 0.1011 twice.

When bin count does matter: gaps of opposite sign

Bin count starts mattering the moment a single bin holds errors of opposite sign, and then it matters enormously.

Consider one wide bin covering [0.0, 0.5) that happens to contain two very different groups of rows. One group is badly underconfident, the other badly overconfident, and they are the same size:

100 predictions at 0.10, of which 40 are positive   (true rate 0.40)
100 predictions at 0.40, of which 10 are positive   (true rate 0.10)

bin mean confidence = 0.25       bin observed rate = 0.25       gap = 0.000

The bin’s mean confidence is (0.10 + 0.40)/2 = 0.25 and its observed rate is (0.40 + 0.10)/2 = 0.25, so the gap is exactly zero.

ECE contribution: zero. Actual miscalibration: 0.30 on every one of those 200 rows.

Within-bin errors of opposite sign cancel, and coarser bins give them more room to cancel. Add one bin edge at 0.25 and the cancellation has nowhere left to go:

two bins:  [0.00, 0.25)   conf 0.10   rate 0.40   gap 0.30   weight 0.5
           [0.25, 0.50)   conf 0.40   rate 0.10   gap 0.30   weight 0.5

ECE contribution  =  0.5(0.30) + 0.5(0.30)  =  0.300        (it was 0.000)

One extra bin edge moved the reported error from 0.000 to 0.300 without changing a single prediction.

That is the bin-count dependence in its purest form, and it is exactly the sign-cancellation mechanism the overconfident model above never triggers. There, every gap pointed the same way, so there was nothing for a coarse bin to hide and every scheme reported 0.1011. Here, two gaps point in opposite directions, one bin edge decides whether they cancel, and the answer swings the entire width of the error.

Bin count is harmless on a model that errs one way and decisive on a model that errs both ways — and you cannot tell which you have from the ECE number alone. Two consequences:

• ECE is a lower bound on miscalibration, and it is monotone in bin coarseness. “Our ECE is 0.02” is not a claim until the binning scheme, the bin count, and the bin-mass policy are stated alongside it.
• ECE is not a proper scoring rule and can be gamed — a model can lower its ECE while getting worse. Log loss and Brier cannot be gamed this way. Report ECE for interpretability and a proper score for the decision.

The implementation below is deliberately short so you can see where the binning decision enters. Both n_bins and equal_mass change the answer on some models and not on others, which is precisely why they belong in the number you report:

def ece(probs, labels, n_bins=10, equal_mass=False):
    """Expected calibration error. Report n_bins and the binning mode with it."""
    pairs = sorted(zip(probs, labels))
    n = len(pairs)
    if equal_mass:
        edges = [i * n // n_bins for i in range(n_bins + 1)]
        groups = [pairs[edges[i]:edges[i + 1]] for i in range(n_bins)]
    else:
        groups = [[] for _ in range(n_bins)]
        for p, y in pairs:
            groups[min(int(p * n_bins), n_bins - 1)].append((p, y))
    total = 0.0
    for g in groups:
        if not g:
            continue
        conf = sum(p for p, _ in g) / len(g)
        acc = sum(y for _, y in g) / len(g)
        total += (len(g) / n) * abs(acc - conf)
    return total

Platt scaling vs isotonic regression

Two calibrators are standard, each resting on its own assumption — and when the popular one’s assumption fails, it makes things worse.

Platt scaling fits a sigmoid curve from the model’s raw score to a probability, using just two fitted numbers: a slope and an offset.

Isotonic regression instead fits any staircase that never goes down, with as many steps as the data supports.

Both fit a monotone map from raw score to probability on held-out data — a map that never reorders two rows — so neither can change AUC by more than tie-breaking effects. They differ in how much shape they can express, and therefore in what they assume about the distortion they are correcting.

The row of that table to read first is Fails when. Everything else follows from it:

Four pieces of shorthand in that table need expanding.

• PAVA is the pool-adjacent-violators algorithm: the procedure that fits an isotonic regression by repeatedly merging any neighbouring pair that runs the wrong way. It is worked by hand below.
• NLL is the negative log-likelihood, which for a classifier is the same quantity as log loss under another name.
• One-vs-rest means fitting one binary model per class — “is it class 1 or not,” “is it class 2 or not” — and combining the answers afterwards.
• s in Platt’s sigmoid(a·s + b) is the model’s raw score, with a and b the two fitted numbers. In practice s is usually the logit rather than the probability, which is what the worked example below uses.

Worked: when Platt scaling is worse than doing nothing

That is the most decision-relevant entry in the table, so here it is with numbers.

Take a model that is exactly right at the bottom of its range and badly overconfident only at the top. That is an asymmetric distortion, not a sigmoidal one, which is exactly the case Platt’s assumption excludes.

There are 1,000 held-out rows sitting at three distinct score levels. Only the top level is wrong:

n = 400   raw 0.10   true rate 0.10      gap 0.000
n = 400   raw 0.50   true rate 0.50      gap 0.000
n = 200   raw 0.90   true rate 0.55      gap 0.350

raw ECE  =  0.2 x 0.350  =  0.0700

Only the 200 rows at 0.90 are miscalibrated, and they are off by 0.350. Weighted by their 20% share, that gives raw ECE = 0.0700.

Fitting Platt scaling here means choosing a and b in sigmoid(a·logit(s) + b) to minimize log loss on these held-out rows. A two-parameter sigmoid has no way to bend only at the top, so it tilts everywhere instead:

fitted  a = 0.5382,  b = -0.4822

  0.10  ->  0.1591     gap 0.0591      (was 0.000)
  0.50  ->  0.3817     gap 0.1183      (was 0.000)
  0.90  ->  0.6683     gap 0.1183      (was 0.350)

Platt ECE     =  0.0946      <- worse than the 0.0700 it started from
isotonic ECE  =  0.0000      <- the three levels are already monotone; PAVA reproduces them exactly

Platt cut log loss from 0.6261 to 0.5683 — it is fitted to do precisely that — while raising calibration error from 0.0700 to 0.0946. It did that by smearing one bin’s error across the 800 rows that had been exactly right.

Two parameters buy low variance, meaning the fit barely moves if you resample the calibration set. They pay for it in bias, meaning the fitted shape is systematically wrong no matter how much data you give it.

When the true distortion is not sigmoid-shaped, that bias has to land somewhere. Here it lands on the rows that were already correct.

Isotonic assumes only monotonicity, so it has room to absorb the shape instead. The three levels here are already increasing, so PAVA leaves them alone and its ECE is 0. That is what the Fixes: any monotone distortion cell means in practice.

State the assumptions side by side. Platt assumes the distortion is a sigmoid and needs only a few hundred rows. Isotonic assumes only that the distortion never reorders anything and needs a few thousand. Choosing between them is choosing which of those two assumptions your situation can afford.

The two multiclass forms are not symmetric either

Temperature scaling fits one scalar across all K logits at once, so the output is still a distribution by construction.

One-vs-rest isotonic fits K independent step functions, and nothing makes K separately calibrated numbers sum to 1. You renormalize afterwards, and that renormalization is itself an uncalibrated operation that can undo part of what each fit achieved.

That asymmetry, not accuracy, is why temperature is the deep-learning default and one-vs-rest isotonic is a fallback for small K with plenty of held-out data.

PAVA worked by hand

It is short enough to follow completely.

Take eight held-out rows, sorted by model score from lowest to highest, whose true labels turn out to be [0, 1, 0, 0, 1, 1, 0, 1].

Isotonic regression wants a fitted value for each position that never decreases as you move right. PAVA gets there by repeatedly finding a neighbouring pair that decreases — a “violator” — and replacing both with their average. Repeat until nothing decreases anywhere:

start                       0 | 1 | 0 | 0 | 1 | 1 | 0 | 1
1 > 0 at positions 2,3   -> 0 | (1,0)=0.500 | 0 | 1 | 1 | 0 | 1
0.5 > 0 at position 4    -> 0 | (1,0,0)=0.333 | 1 | 1 | 0 | 1
1 > 0 at positions 5,6,7 -> 0 | 0.333 x3 | (1,1,0)=0.667 x3 | 1
check: 0 <= 0.333 <= 0.667 <= 1                        monotone, done

fitted map: [0.000, 0.333, 0.333, 0.333, 0.667, 0.667, 0.667, 1.000]

Notice what isotonic did: positions 2, 3 and 4 held three distinct raw scores and all three collapsed to the single value 0.333.

That is the ties problem. Inside a pooled block, ranking information is destroyed, which is why isotonic can shave a fraction off AUC even though it is monotone.

Temperature scaling

This is the deep-learning default and deserves its own heading.

You fit a single positive number T, the temperature, by minimizing log loss on a validation set. Then you divide every logit by it before the softmax. A T above 1 flattens the distribution and reduces confidence; a T below 1 sharpens it.

Because z/T is strictly increasing in z, argmax is unchanged, accuracy is unchanged, and AUC is unchanged. It changes only the numbers. With one parameter fitted on thousands of rows, it essentially cannot overfit.

Its assumption is that the whole logit vector is miscalibrated by one common factor. That is a strong assumption. It happens to describe the zero-training-error mechanism of the next subsection almost exactly, and it fails when different classes are miscalibrated by different amounts.

The rule that catches everyone

The calibrator must be fitted on data the model did not train on.

A calibrator fitted on training predictions learns the model’s training-set confidence, which is far higher than the confidence it will have on new inputs. So you will “calibrate” a model into being worse.

Use a dedicated calibration split held out from the start. If data is scarce, use cross-fitted calibration instead: split into folds, fit the model on all but one fold, collect that fold’s out-of-sample predictions, repeat, and fit one calibrator on the pooled out-of-sample predictions.

Why modern deep nets are overconfident — the mechanism

Why does a large neural network reliably claim more confidence than it has earned? The explanation is a chain of six steps, laid out in the diagram below. The step that matters most is the second-to-last one.

flowchart TD
    A["Training error reaches zero"] --> B["Every training point is<br/>on the correct side"]
    B --> C["NLL is still positive,<br/>and scaling all logits by k > 1<br/>strictly lowers it"]
    C --> D["No finite minimizer:<br/>the optimizer keeps growing<br/>the logit magnitude"]
    D --> E["Softmax saturates<br/>-> confidence near 1.0<br/>on training inputs"]
    E --> F["That confidence transfers<br/>to test inputs where the model<br/>is right only 80% of the time"]
    F --> G["mean confidence 0.90<br/>accuracy 0.80<br/>ECE 0.10"]

    style C fill:#bc6c25,color:#fff
    style D fill:#9d0208,color:#fff
    style G fill:#1d3557,color:#fff

Read that chain forward, because each box is forced by the one before it.

Training error reaching zero means every training point already sits on the correct side of the boundary. The classification is finished. And yet NLL is still strictly positive, because -log p vanishes only in the limit.

That leaves exactly one direction in which the loss can still fall: scale all the logits by some k > 1. That moves no point across the boundary and strictly lowers NLL on every correctly-classified row.

So there is no finite minimizer, and the optimizer keeps growing the logit magnitude until the softmax saturates and confidence on training inputs sits at 1.0.

Then comes the step that matters. That confidence is a property of the weights, not of the training set, so it transfers intact to test inputs where the model is right only 80% of the time. It arrives as mean confidence 0.90 against accuracy 0.80.

The load-bearing step is the third one, and it is one line of arithmetic. Same two logit vectors, one scaled up by 3:

logits (2, 0)  ->  p = sigmoid(2) = 0.8808   ->  loss = 0.1269
logits (6, 0)  ->  p = sigmoid(6) = 0.9975   ->  loss = 0.0025

same decision, same ranking, 51x lower loss

Once training error is zero, increasing confidence is the only remaining direction in which the loss can decrease, so gradient descent keeps walking that way for as long as you let it.

Three modern practices make this worse rather than better:

• Capacity large enough to reach zero training error is now routine.
• Weight decay — a penalty on the size of the weights, which is what would otherwise stop the logits from growing — is typically dialled down because it costs accuracy.
• Architectures built around residual connections and normalization layers make reaching zero training error easy.

Put those together and overconfidence is the predicted outcome rather than a surprise.

The trace below shows the resulting gap on train and test, then the same model rescored after fitting a temperature of 1.9. Compare the accuracy column before and after:

                accuracy   mean conf   ECE     after T = 1.9
train              1.000       0.999     -      -
test               0.797       0.898   0.1011   ECE 0.0121, accuracy 0.797

Accuracy did not move — it cannot, because temperature is monotone. The calibration error fell by a factor of eight for the cost of fitting one number.

Why boosted trees are under-confident at the extremes

Boosted trees fail in the opposite direction, by a completely different mechanism — and the difference matters because the fix is different too.

Where a deep network’s confidence runs away, a gradient-boosted tree’s confidence is structurally prevented from ever getting there.

How a boosted model builds its score. It is additive: F = F_0 + eta · sum_m h_m. Start from a base score F_0, then add the output of tree m, each multiplied by a learning rate eta that is deliberately small (0.05 is typical) so that no single tree can dominate.

Every leaf of every tree carries the value w = -G/(H + lambda), where:

• G is the summed gradient of the loss over the rows in that leaf,
• H is the summed second derivative,
• lambda is an L2 regularization term, a constant added to the denominator specifically to keep leaf values small (Xgboost second order and a regularizer inside the split criterion).

So each leaf value is shrunk twice: once by lambda in the denominator and again by eta outside.

Why that caps confidence. To emit p = 0.99 the model needs F = ln(99) = 4.595. With eta = 0.05 and typical leaf values, that takes a hundred-plus rounds all pushing the same row in the same direction.

Two things cut the run short before it gets there. Early stopping halts training when validation loss stops improving, long before the confident rows saturate. And min_child_weight independently blocks the tiny leaves that would carry extreme values in the first place.

Where “a hundred-plus rounds” comes from

That should not be an assertion either, so here is the arithmetic.

Take a leaf holding n = 20 rows, all of them positive, in a problem with 1% prevalence, so the base score is F_0 = logit(0.01) = -4.595.

For logistic loss the per-row gradient is g_i = p_i - y_i and the per-row second derivative is h_i = p_i(1 - p_i). Summing over a leaf where every y_i = 1 gives G = -n(1 - p) and H = n·p(1 - p), so the leaf value works out to the expression below. Watch eta·w in the last column shrink from 0.83 to 0.008 across the run — nearly 100x smaller steps by the end:

w  =  n(1 - p) / ( n·p(1-p) + lambda )

round   1   p = 0.0100   G = -19.800   H = 0.198   w = 16.528   eta·w = +0.8264
round  19   p = 0.4988                             w =  1.671   eta·w = +0.0835
round  61   p = 0.8999                             w =  0.714   eta·w = +0.0357
round 198   p = 0.9899                             w =  0.168   eta·w = +0.0084

Now iterate F <- F + eta·w at eta = 0.05, lambda = 1, and count how many rounds each level of confidence costs. The last two lines rerun the same loop at two other values of lambda:

rounds to reach p = 0.50      19
rounds to reach p = 0.90      61
rounds to reach p = 0.99     198     <- the last 0.09 of probability costs 137 rounds
rounds to reach p = 0.999  1,146
rounds to reach p = 0.9999 10,193

same run to p = 0.99 at lambda = 0     75
same run to p = 0.99 at lambda = 5    631

The step size collapses exactly where you need it not to. As p -> 1 both G and H go to zero, so w -> n(1-p)/lambda. The increments shrink in proportion to 1-p.

That shows up as a round count that multiplies by a growing factor for each additional nine of confidence: 3.2x going from 0.9 to 0.99 (61 to 198), then 5.8x (198 to 1,146), then 8.9x (1,146 to 10,193).

It also changes which term controls the leaf value. lambda is negligible against H in the middle of the range and becomes the entire denominator at the top. That is why lambda alone stretches the run from 75 rounds to 631 while changing nothing about the fit’s direction.

And early stopping watches validation loss, which is dominated by the many rows in the middle. So it cuts the run off in the flat part of that curve.

The parameter min_child_weight closes the other door. It is a floor on H = sum_i p_i(1 - p_i) within a leaf: a minimum amount of second-derivative mass a leaf must contain before the algorithm is allowed to create it.

That sum is precisely the quantity that collapses as predictions saturate. The block below gives per-row h at four confidence levels, then divides 1 by each to get how many rows a leaf needs:

h per row:      p = 0.50 -> 0.2500     p = 0.90 -> 0.0900
                p = 0.95 -> 0.0475     p = 0.99 -> 0.0099

rows needed for H >= 1 (the XGBoost default):
                p = 0.50 ->     4      p = 0.99 ->   102

A leaf that would push already-confident rows further needs about 25x more rows to be allowed to exist than one operating in the middle — 102 against 4. The regularizer that keeps boosting from overfitting is the same regularizer that keeps it from ever saying 0.99.

What that looks like on a real reliability table is below. Read the mean pred and observed rate columns against each other: too high at the bottom, honest in the middle, too low at the top:

predicted bin   mean pred   observed rate
0.00-0.10         0.041         0.012      <- predicts too high at the bottom
0.10-0.30         0.196         0.131
0.30-0.70         0.492         0.501      <- fine in the middle
0.70-0.90         0.804         0.881
0.90-1.00         0.938         0.983      <- predicts too low at the top

observed score range: [0.021, 0.961]       the model never emits 0.99

That is the classic sigmoid-shaped reliability curve. It is exactly the shape Platt scaling was designed to invert, and it is why Platt is the traditional recommendation for boosted trees and for SVMs — support vector machines, a classical model that outputs an uncalibrated distance from a decision boundary rather than a probability at all.

The production failure this causes is specific and silent. A business rule written as “auto-block when p > 0.99” never fires, because the model’s maximum output is 0.961. Nothing errors. The rule is simply dead, and it stays dead until someone plots the score histogram.

Random forests are under-confident too, for an unrelated reason

The reason is averaging.

A random forest is a collection of B decision trees. Each is grown on a bootstrap resample of the rows — a sample drawn with replacement, so each tree sees a slightly different dataset — and each is allowed to consider only m randomly chosen features at every split.

The forest’s output for a row is the fraction of its trees that voted positive. So emitting exactly 1.0 requires unanimity, and how hard unanimity is depends entirely on how correlated the trees are.

Trees in a forest are not independent, because they are all grown on resamples of the same data.

Bagging a variance fix and why it needs unstable base learners derives the variance floor rho·sigma_t^2. Here rho (the Greek letter rho) is the average correlation between any two trees’ votes and sigma_t^2 is the variance of a single tree’s vote. The floor says averaging can never drive the variance below rho times the single-tree variance, no matter how many trees you add.

Random forests decorrelation which matters more than the bagging tabulates the correlations that different feature-sampling settings produce: rho = 0.35 when every feature is available to every split, 0.12 when each split sees m = sqrt(p) of the p features, and 0.03 at the extreme of m = 1.

Turning that correlation into a statement about the vote. Model the B tree votes as exchangeable: any two of them have the same correlation rho with each other, and none is special.

The standard distribution with that property is the beta-binomial. Instead of each tree flipping an independent coin with fixed bias, all the trees in one forest share a single bias drawn at random from a beta distribution, which is what makes their votes agree more than chance.

Set that shared bias to have mean equal to the per-tree accuracy and intra-class correlation rho. The block below gives the probability that all 100 trees vote positive. Compare each rho line against the rho = 0 line above it — that is the independence figure people quote:

P(all 100 vote positive), per-tree accuracy 0.95
   rho = 0      (independent)   0.95^100  =  0.0059
   rho = 0.03                             =  0.097
   rho = 0.12                             =  0.362
   rho = 0.35                             =  0.669

same, per-tree accuracy 0.75  (a borderline row)
   rho = 0                                =  3e-13
   rho = 0.12                             =  0.005

So the independence figure is not just imprecise, it points at the wrong rows.

On an easy row (per-tree accuracy 0.95) a real forest reaches 1.0 between a third and two thirds of the time at realistic correlations. The extremes are not bounded away there at all — the 0.95^100 = 0.0059 figure is off by a factor of 60 at rho = 0.12.

Where averaging genuinely traps the output is on borderline rows. Bootstrap resampling keeps a persistent dissenting minority there no matter how many trees you add: at per-tree accuracy 0.75 the forest lands near 0.75 and essentially never at 1.0.

The full vote distribution makes the point more sharply than the unanimity figure does. Take the same beta-binomial model with B = 100 trees, rho = 0.12, and per-tree accuracy 0.75. Read off the standard deviation of the vote fraction, then how much of its probability mass sits away from the boundaries:

sd of the vote fraction    0.1554      (0.0433 if the trees were independent)
P(vote in (0.05, 0.95))    0.9240      (0.9428 if the endpoints are included)
P(vote = 1.00)             0.0050

Correlation multiplies the spread of the vote by 3.6 — 0.1554 against the independent 0.0433 — and still leaves 92% of the mass in the interior.

And adding trees does not help. The variance of the vote fraction is p(1-p)·(rho + (1 - rho)/B). Drive B to infinity and the (1 - rho)/B term vanishes, leaving the vote scattered around 0.75 with a standard deviation of sqrt(0.75 x 0.25 x 0.12) = 0.15.

At B = 1,000 it is already 0.1505, so you are at the floor long before you notice. That is the variance floor from Bagging a variance fix and why it needs unstable base learners showing up as a calibration symptom instead of an accuracy one.

That is the honest version of the symptom, and it is still the one that costs you money, because borderline rows are exactly where a decision threshold sits.

The mechanism is different from boosting’s shrinkage. The fix is the same: fit a monotone calibrator, or read the output as a rank and stop pretending it is a probability.

When calibration matters, and when it is irrelevant

The practical payoff of everything above fits in one table: what your model’s output feeds into, and whether calibration is worth any effort in each case.

The rule underneath the table is one sentence. Calibration matters exactly when something reads the number’s face value rather than its position in an ordering.

NDCG in the first row is normalized discounted cumulative gain, the standard quality measure for a ranked list. Like AUC, it depends only on the order.

If nothing downstream reads the number, calibration is a metric you can afford to lose. If anything downstream multiplies it, compares it to a constant, or shows it to a person, it is the whole game.

4. Drift: three different failures with three different signals

There are three kinds of drift. They need three different fixes, and only two of them can be detected before the labels arrive — so telling them apart is worth real effort.

Drift is the general name for the live data drifting away from the data the model was trained on. It is the reason a model’s accuracy decays over time even though nobody touched it.

The taxonomy is not arbitrary. Write the joint distribution p(x, y) — the probability of seeing a particular input together with a particular label — and factor it the two ways it can be factored.

Each factoring isolates a different pair of things that can move independently. Every drift type is one of those pieces moving while its partner stays put:

p(x, y) = p(y | x) · p(x)          <- covariate shift lives here
        = p(x | y) · p(y)          <- label shift lives here

Covariate shift is p(x) moving while p(y|x) holds. The inputs changed, but the relationship between input and answer did not, so a customer with a given profile is still just as risky as before. A covariate is an input feature, so covariate shift means “the feature distribution moved.”

Label shift is p(y) moving while p(x|y) holds. The mix of classes changed, but what each class looks like did not: twice as much fraud, and fraud still looks like fraud. It also travels under the name prior shift, since p(y) is the prior of Resampling and the calibration it breaks. The two terms mean the same thing; this chapter says label shift throughout.

Concept drift is p(y|x) itself moving. The same input now deserves a different answer. This is the case no reweighting can repair.

flowchart TD
    D["Model degraded in production"] --> Q1{"Did p of x move?"}
    Q1 -->|yes| CS["COVARIATE SHIFT<br/>inputs moved, p of y given x fixed"]
    Q1 -->|no| Q2{"Did the predicted<br/>positive rate move?"}
    Q2 -->|yes| LS["LABEL SHIFT<br/>the prior moved, p of x given y fixed"]
    Q2 -->|no| CD["CONCEPT DRIFT<br/>the relationship itself changed"]
    CS --> F1["Detect: feature PSI,<br/>domain classifier - NO LABELS<br/>Fix: importance weighting,<br/>retrain, drop unstable features"]
    LS --> F2["Detect: predicted-rate monitor,<br/>BBSE - NO LABELS<br/>Fix: constant logit shift<br/>-- the section 2 formula again"]
    CD --> F3["Detect: rolling loss on<br/>delayed or proxy labels - LABELS REQUIRED<br/>Fix: retrain. Nothing else works."]

    style CS fill:#2d6a4f,color:#fff
    style LS fill:#40916c,color:#fff
    style CD fill:#9d0208,color:#fff

That diagram is the triage you run from the entry box, model degraded in production, when you do not yet know why.

The first question — did p of x move?, meaning did the distribution of the inputs themselves change — is answerable immediately and without labels. It splits covariate shift off from the other two.

The second question, whether the predicted-positive rate moved, splits label shift from concept drift.

Concept drift is the residual diagnosis: every distribution you can see is stable and the model is wrong anyway.

Each of the three boxes names its detector and its fix. Both are worth reading now, because the rest of the section elaborates them.

• Covariate shift is detected by feature PSI — the population stability index, a per-feature measure of how far a distribution has moved, computed in Psi and kl computed — or by a domain classifier. Neither needs labels. It is fixed by importance weighting, by retraining, or by dropping the unstable features.
• Label shift is detected by watching the predicted-positive rate, or by BBSE, black-box shift estimation, derived at the end of this section. Neither needs labels. It is fixed by the same constant logit shift derived in Resampling and the calibration it breaks.
• Concept drift can only be detected by a rolling loss computed on real or proxy labels, which is why it requires labels. The only fix is to retrain.

The table below is the same three columns side by side, which is the form to memorize for an interview:

Why covariate shift hurts at all

In theory it should not hurt. If your model were correctly specified — meaning the true relationship is inside the family of functions the model can express — and well-estimated everywhere, covariate shift would cost you nothing. p(y|x) is unchanged, so the same function is still right on the new inputs.

It hurts for two concrete reasons.

1. The model is misspecified, so the fit was a compromise across regions — and that compromise was weighted by the old p(x), which has now changed.
2. The new inputs land in regions where there was little training data, so the model is extrapolating: making predictions in a region it never saw.

That is why “PSI is high” is not by itself a reason to retrain. It is a reason to check whether the moved mass sits in a region the model was ever good at.

The three covariate-shift fixes, in detail

“Drop unstable features” needs its criterion stated or it is just a slogan. The instability that earns a drop is not high PSI on its own. It is high PSI and meaningful permutation importance and no reason to expect the movement to stop.

Permutation importance measures how much a model’s accuracy falls when you randomly shuffle one feature’s values. It is the standard way to ask “does the model actually lean on this?”

A feature that drifts and the model leans on is a scheduled outage. A feature that drifts and the model ignores is noise in your alerting, and the answer there is to fix the alert.

The trade is explicit and permanent: dropping a moved feature costs accuracy on the whole population that did not move, in order to buy stability against the part that did. Take it when the drift is structural — an upstream vendor you do not control, a field whose definition changes with every product release. Do not take it for a one-off you could have handled with a retrain.

Importance weighting means retraining or re-evaluating on the old data with each row weighted by how much more common it is now than it used to be, so that the old data is made to stand in for the new distribution.

The obvious way to compute that weight is to estimate the two densities p_new(x) and p_old(x) and divide. That is hopeless in more than a few dimensions.

There is an estimator that avoids density estimation entirely. Train a domain classifier d(x) = P(x came from the new period): an ordinary binary classifier whose label is not the outcome but simply which time period the row came from, fitted on old rows against new rows.

Bayes’ rule then converts its odds directly into the density ratio you wanted:

d(x)/(1 - d(x)) = [ p_new(x)·P(new) ] / [ p_old(x)·P(old) ]

  ->  p_new(x)/p_old(x) = [ d(x)/(1 - d(x)) ] · (n_old / n_new)

Here P(new) and P(old) are the shares of the pooled dataset coming from each period, which is why they turn into the row counts n_new and n_old in the second line.

Work the simplest case first. Pool 100,000 rows from each period and label them old or new. Then n_old/n_new = 1, the correction factor disappears, and the weight is just the classifier’s own odds:

d(x) = 0.50   ->   p_new/p_old = 1.00      this region did not move
d(x) = 0.60   ->   p_new/p_old = 1.50      half again as common now
d(x) = 0.75   ->   p_new/p_old = 3.00      three times as common
d(x) = 0.90   ->   p_new/p_old = 9.00      nine times; a retrain weights these rows 9x

The n_old/n_new factor starts mattering the moment the two periods are different sizes.

Score four weeks of history against one week of production and the factor is 4. Now d(x) = 0.50 no longer means “unchanged.” It means “as many new rows as old rows land here,” and with 4x fewer new rows overall that is a 4x enrichment:

n_old = 200,000, n_new = 50,000:   d(x) = 0.50  ->  1.00 x 4 = 4.00
                                   d(x) = 0.75  ->  3.00 x 4 = 12.0

For a reweighted retrain the constant is harmless. Multiplying every weight by the same number scales the whole loss and does not move its argmin, the parameter setting at which the loss is smallest.

But the moment you read a weight as “how much more common is this region now,” dropping the ratio is the difference between reporting 1x and 4x.

Importance weighting’s own assumption is that the new distribution’s support is contained in the old one. Every region that has new rows must have had at least some old rows.

Where that fails, the density ratio is infinite, the estimated weight is whatever your classifier’s saturated output happens to be, and a handful of old rows end up carrying the entire retrain. That failure mode is the reason for the AUC thresholds below.

The same classifier is a drift detector for free. Its held-out AUC is the drift magnitude: 0.5 means the two periods are indistinguishable, 0.85 means severe. Its feature importances name which features moved. One model gives you the alarm, the diagnosis, and the correction weights.

This chapter quotes two AUC figures for that classifier — 0.85 here, and an alert threshold of 0.70 in the monitoring table of Detecting drift without labels the real production problem. They are not in conflict. 0.70 is where you look; 0.85 is where you act.

At 0.70 the periods are distinguishable but heavily overlapping, which is the ordinary condition of any system with weekly seasonality. At 0.85 the classifier can nearly separate them, which means most new rows sit in regions the training data barely covered — so importance weighting is running on ratios estimated from almost no old data.

Label shift, and why it is the §2 formula again

Label shift needs no new mathematics at all — only one genuinely new problem comes with it.

Label shift is a change of prior with p(x|y) fixed. That is the identical assumption behind resampling in Resampling and the calibration it breaks: the classes changed their proportions but not their appearance.

So the identical correction applies, with the old and new prevalences in place of the training and true ones:

logit(p_new) = logit(p_model) + ln( [pi_new/(1-pi_new)] · [(1-pi_old)/pi_old] )

The only new problem is estimating the new prevalence pi_new when you have no labels from the new period.

BBSE, black-box shift estimation, does it from the model’s predictions alone. It rests on two numbers you already measured on the old validation set, at your operating threshold:

• TPR, the true positive rate: the fraction of actual positives the model flags.
• FPR, the false positive rate: the fraction of actual negatives it flags by mistake.

Both are properties of p(x|y), which label shift leaves alone by definition, so they carry over to the new period unchanged.

The predicted-positive rate you can observe today is then a simple linear function of the unknown prevalence. Every flagged row is either a true positive (a pi fraction of rows, caught at rate TPR) or a false positive (a 1-pi fraction, flagged at rate FPR):

q = P(predict positive) = TPR·pi + FPR·(1 - pi)

You know q, TPR and FPR. Solve for pi.

Here is BBSE worked end to end. On the old validation set the model had TPR = 0.80 and FPR = 0.05 at the threshold in use, and the prevalence then was pi_old = 0.10. In the new period the model has scored 10,000 rows and flagged 2,300 of them, so the observed predicted-positive rate is q = 2300/10000 = 0.23:

q = 0.23 = 0.80·pi + 0.05·(1 - pi) = 0.05 + 0.75·pi
pi_new = 0.18 / 0.75 = 0.240                       <- prevalence 2.4x, zero labels used

logit shift = ln( (0.240/0.760) · (0.900/0.100) ) = ln(2.842) = 1.0445

a row scoring 0.30:  odds 0.4286 -> 0.4286 × 2.842 = 1.2180 -> p = 0.549

You corrected the model today, from predictions only, while the true labels are still 45 days out.

Contrast the alternative: wait for chargebacks — the disputed-transaction reversals that are a card issuer’s ground-truth fraud label, and which take six weeks to settle — notice accuracy fell, then retrain. That path takes two months.

The assumption you are betting on is that p(x|y) held. If the attackers’ modus operandi changed rather than merely their volume, then what you have is concept drift: p(x|y) moved, the TPR and FPR you carried over are stale, and the correction is wrong in an unknown direction.

Sanity-check it before you apply it. Confirm that feature PSI stayed low while the predicted-positive rate moved. Same-looking inputs plus more positive predictions is the signature of genuine label shift; moved inputs is the signature that BBSE’s assumption has failed.

5. PSI and KL, computed

The industry-standard drift number is worth computing by hand at least once, because underneath it is a symmetrized KL divergence (defined below) — and three facts about it determine whether your alert threshold is meaningful or theatre.

PSI, the population stability index, answers one question about one feature: how far has this feature’s distribution moved since training?

You compute it by cutting the feature’s range into bins whose edges are frozen at training time — usually the training deciles, the ten cut points that split the training data into equal tenths. Then you compare what share of rows falls in each bin now against what share fell there during training:

PSI = sum over bins of  (a_i - e_i) · ln(a_i / e_i)

Here e_i is the expected share — bin i’s share of the training rows — and a_i is the actual share, bin i’s share of today’s production rows.

Because the bins were cut at training deciles, every e_i is 0.10 by construction. That is why the e_i column below is a constant.

Freezing the edges is not a detail but the entire method. Re-binning on production data would compare a distribution against itself and always return approximately zero.

In the worked table below, the feature has drifted upward: the low bins have emptied out (bin 1 fell from 0.10 to 0.03) and the high bins have filled. Watch the rightmost column — bin 1 alone contributes more than a third of the total:

bin    e_i    a_i    a/e    ln(a/e)     (a-e)·ln(a/e)
 1     0.10   0.03   0.30   -1.20397      0.084278
 2     0.10   0.05   0.50   -0.69315      0.034657
 3     0.10   0.06   0.60   -0.51083      0.020433
 4     0.10   0.08   0.80   -0.22314      0.004463
 5     0.10   0.09   0.90   -0.10536      0.001054
 6     0.10   0.11   1.10    0.09531      0.000953
 7     0.10   0.12   1.20    0.18232      0.003646
 8     0.10   0.14   1.40    0.33647      0.013459
 9     0.10   0.16   1.60    0.47000      0.028200
10     0.10   0.16   1.60    0.47000      0.028200
                                        -----------
                                    PSI = 0.219344

The conventional thresholds people read that number against:

So this feature, at PSI 0.219, lands in “investigate.”

Those thresholds are folklore rather than statistics. The three facts below explain what the number actually is, which part of it is doing the work, and why the thresholds cannot be right as stated.

Fact 1: PSI is symmetric KL divergence

KL divergence, named for Kullback and Leibler, measures how much information you lose by using one distribution in place of another. It is zero when they match and positive otherwise, and it is asymmetric: KL(a||e) and KL(e||a) are different numbers.

Splitting the PSI sum shows it is exactly the two directions added together. The middle step uses -ln(a/e) = ln(e/a):

sum (a - e)·ln(a/e) = sum a·ln(a/e) - sum e·ln(a/e)
                    = sum a·ln(a/e) + sum e·ln(e/a)
                    = KL(a || e) + KL(e || a)          <- Jeffreys divergence

here:  0.219344 = 0.1011 + 0.1182

That sum of the two directions has a name of its own: the Jeffreys divergence. Every term in it is non-negative, so PSI is zero exactly when the two distributions match and positive otherwise.

PSI is not an arbitrary index. It is a symmetrized information distance with a folklore threshold table bolted on.

Fact 2: it is dominated by the bins that emptied, not the bins that filled

Bin 1 alone contributes 0.084278 of the 0.219344 total. That is 38% of the number coming from one bin that holds 3% of the mass.

The reason is that ln(a/e) runs off to negative infinity as a approaches zero, while it grows only logarithmically as a grows. Emptying is punished far harder than filling.

A bin that empties completely would make PSI infinite, so every implementation quietly floors a_i at some small epsilon, a placeholder minimum value, with 0.0001 the common choice.

That arbitrary epsilon silently sets your alert threshold. A feature with one empty bin reports whatever number the epsilon dictates, and nothing else about the data matters.

Fact 3: the thresholds are independent of sample size, which cannot be right

Even with no drift whatsoever, a finite sample will not reproduce the training shares exactly. So PSI is positive by luck alone.

How positive is calculable. Under pure multinomial sampling noise — the randomness of throwing n independent rows into K bins — the quantity 2n·KL follows a chi-square distribution with K-1 degrees of freedom. A chi-square is the standard distribution of a sum of squared random deviations, and its degrees of freedom is the number of independent ways the counts can vary.

The mean of that distribution is K-1, so one KL direction alone has mean (K-1)/(2n).

The factor of 2 comes back because PSI is not one KL but two. By fact 1 above, PSI = KL(a||e) + KL(e||a). To leading order in the deviation a - e both directions equal the same chi2/(2n), so the two halves add, the 2s cancel, and the noise floor is twice the single-direction figure:

E[KL(a||e)]  ~  (K-1)/(2n)      one direction
E[KL(e||a)]  ~  (K-1)/(2n)      the other, equal to leading order
             ----------------------------------------------------
E[PSI]       ~  2 · (K-1)/(2n)  =  (K - 1)/n

That doubling is the step people drop, and dropping it halves every floor below.

Simulation agrees with the doubled version and not the halved one. 200,000 no-drift draws at K = 10, n = 10,000 give a mean PSI of 0.00090, which matches (K-1)/n = 9/10,000 and is twice (K-1)/(2n).

Now read the floor off at three batch sizes. The third line is the one that should change how you set alerts:

K = 10 bins, n = 10,000 rows  ->  E[PSI] ≈ 0.0009    (noise floor, negligible)
K = 10 bins, n =    500 rows  ->  E[PSI] ≈ 0.018
K = 10 bins, n =    100 rows  ->  E[PSI] ≈ 0.090     <- already at "investigate"

On a small daily slice, the 0.1 threshold fires on pure noise. On 10 million rows a PSI of 0.05 is overwhelmingly significant and operationally meaningless.

The fix is to calibrate the threshold to your own batch size. Bootstrap the PSI: resample your stable historical windows with replacement, at the batch size you actually score, and recompute PSI many times. Then set the alert at the p99 of that distribution — the value exceeded by only 1% of the no-drift runs. That is the mature answer, and it takes twenty lines of code.

Two limits PSI cannot check about itself

PSI is univariate. It looks at one feature at a time, so it cannot see a change in how features relate to each other. income and debt can each hold a perfectly calm PSI of 0.02 while their joint relationship inverts — and a model that splits on both is destroyed while every monitor stays green.

PSI is blind to importance. A drifting feature the model barely uses will page you at 3 a.m. for nothing.

Weight PSI by permutation importance to fix the second. Use the domain classifier from Drift three different failures with three different signals to fix the first: that classifier is inherently multivariate and would notice the income-versus-debt inversion immediately.

The two functions below are the section in code. psi() computes the index, and psi_noise_floor() gives the number you must compare it against before alerting:

from math import log

def psi(expected_counts, actual_counts, eps=1e-4):
    """Population Stability Index = KL(a||e) + KL(e||a) over fixed bins.

    Bin edges must be frozen at training time; re-binning on production
    data compares a distribution against itself and always returns ~0.
    """
    e_tot, a_tot = sum(expected_counts), sum(actual_counts)
    total = 0.0
    for e_c, a_c in zip(expected_counts, actual_counts):
        e = max(e_c / e_tot, eps)          # eps here sets your alert threshold
        a = max(a_c / a_tot, eps)
        total += (a - e) * log(a / e)
    return total


def psi_noise_floor(n_rows, n_bins):
    """E[PSI] under no drift. Compare before alerting.

    2n*KL -> chi2(K-1), so ONE direction has mean (K-1)/(2n). PSI is both
    directions and they are equal to leading order, so the floor is twice
    that: (K-1)/n. Dropping the doubling halves every threshold you set.
    """
    return (n_bins - 1) / n_rows

6. Detecting drift without labels — the real production problem

In a running system, the thing you would most like to measure — whether the model is still right — is the thing you cannot have.

Labels arrive late (chargebacks take 45 to 60 days to settle), or partially (only the loans you approved ever produce a repayment outcome), or never (nobody writes in to say the recommendation was bad).

Every drift signal that arrives in time is a signal that needs no labels, and none of them can see concept drift.

The table below is the full menu. Read it as a ladder: each row is slower and more informative than the one above it, and the Misses column is the reason you cannot rely on any single one.

Two terms in it. Latency means how long after the problem starts the signal can tell you about it. p95 in the second row means the 95th percentile — the score value that 95% of predictions fall below — and it is watched alongside the mean because a distribution can shift in its tail while its average holds still.

The row that needs unpacking is the embedding one, because it is the only signal in this chapter that works on unstructured data.

PSI needs bins on a scalar feature, and a sentence or an image does not have one.

The workaround is an embedding: a fixed-length vector of numbers that the model’s own encoder produces as its internal summary of an input, and in which similar inputs land near each other.

You freeze a reference sample of training embeddings at release time — 10,000 to 100,000 vectors is plenty — and monitor how far each production batch sits from that frozen reference. Three standard distances do the job:

• Mean cosine distance to the k nearest reference vectors per row. Cosine distance measures the angle between vectors, so it ignores their length.
• Mahalanobis distance under a Gaussian fitted to the reference. That is a Euclidean distance rescaled by how much the reference actually varies in each direction.
• Maximum mean discrepancy between the two sets, which compares the two clouds as wholes rather than row by row.

What this catches that no per-feature monitor can is the input that is unremarkable in every individual coordinate and jointly sits somewhere the model has never been: a new language, a new document template, a camera with a different colour response.

What it misses is exactly what everything else in the table misses. An in-domain input whose correct answer changed sits in the same place in embedding space as it always did.

The strategic conclusion: you cannot buy your way out of this with a better detector. No detector can see concept drift without labels. The information is not there, because p(y|x) is what changed and y is what you do not have.

So the investment that pays is engineering a fast proxy label: a cheap, quick signal that stands in for the slow true one. A click at 60 seconds that correlates 0.7 with a conversion at 30 days converts a two-month detection latency into a one-day one, and no amount of input monitoring does that.

Feedback loops: the model changes the distribution it is measured on

One failure evades every monitor above, because the input distribution p(x) never moves at all.

What moves is which labels you are allowed to observe. The mechanism: the model’s own decisions determine which outcomes ever get recorded.

The trace below follows one fraud pattern, called P, through a single cycle. Note that the only monitor that would have fired — feature PSI — reads 0.03 in week 9, which is nothing:

week 0   model blocks 95% of fraud pattern P.  P is 2.0% of labeled fraud.
week 4   labeled data now contains only the 5% that slipped through;
         P's apparent share of labeled fraud has fallen to 0.4%.
week 6   scheduled retrain. P's learned risk drops from 0.71 to 0.28.
week 7   deploy. Threshold is 0.30, so P now passes cleanly.
week 9   losses on P up 6x.  Feature PSI across every input: 0.03. No alarm.
week 14  chargebacks settle; aggregate accuracy finally moves 0.4 points.

The model suppressed the evidence for its own decision, and the retrain believed the evidence.

Credit scoring has exactly the same structure: only the applicants you approved ever produce a repayment outcome, and the problem of guessing what the rejected ones would have done is called reject inference.

So does every recommender system, where an item that is never shown accumulates no clicks and therefore looks permanently irrelevant (Why offline ranking metrics disagree with online ctr).

The fix is structural rather than statistical, because no amount of cleverness recovers information that was never recorded:

• Randomized hold-back. Deliberately let 1-2% of would-be blocks through, logged and flagged as an experiment. It costs a known, budgeted amount of fraud and buys an unbiased stream of labels covering the rows the model would otherwise have hidden from itself. This is the only fix that actually works.
• Log the propensity of every action taken — the propensity being the probability that the system would have taken that action for that row, which you know because you chose the policy. Recording it lets you reweight the historical data later by inverse propensity, giving rarely-taken actions proportionally more weight so the biased log can stand in for an unbiased one.
• Monitor the action distribution, not just the input distribution — block rate, approval rate, and score distribution per segment, where a segment is any slice of traffic you can name in advance, such as a country, a device type, or a merchant category. In the trace above, P’s block rate was the one series that moved in week 7, and it moved seven weeks before the accuracy number did.

Monitoring design, and why aggregate accuracy is too slow

The metric everyone reaches for first — overall accuracy on real labels — is the wrong thing to alert on.

There are three independent reasons it is too slow, and each one has arithmetic behind it.

1. Label latency. If chargebacks settle at 45 days, then aggregate accuracy computed today is a report on the model as it existed a month and a half ago. Whatever it tells you, you cannot act on it in time.

2. Dilution. An aggregate is a weighted average, so a catastrophe confined to a small slice barely registers.

Suppose a segment holding 5% of traffic collapses from 0.92 accuracy to 0.50 while everything else holds steady at 0.91. The two lines below are the same weighted average with only that segment’s number changed:

before:  0.95(0.91) + 0.05(0.92) = 0.9105
after:   0.95(0.91) + 0.05(0.50) = 0.8895
headline moves 2.1 points for a total failure in a real segment

3. Statistical power. Statistical power is the ability of a test to notice a real effect. It is set by how much data you have relative to how small the effect is.

To detect a 1-point drop from 0.91 at roughly two standard errors, use se_diff = sqrt(2·p(1-p)/n) for the standard error of a difference between two accuracy measurements each on n rows, set 2·se_diff = 0.01, and solve for n:

2 · sqrt(2(0.91)(0.09)/n) = 0.01
       sqrt(0.1638/n)     = 0.005
                       n  = 6,552 labeled rows per period

At 500 labeled rows a day, 6,552 rows is 13 days per period — so 26 days to compare two of them.

A one-point regression takes about a month to become statistically visible and is diluted 20:1 if it lives in one segment. Aggregate accuracy is an audit metric, not an alert metric.

The latency ladder

Alert on a ladder instead, ordered cheapest and fastest first. The diagram below names six layers, each sitting above the next because it is faster and more certain.

• Pipeline checks the data itself: schema, nulls and ranges. Latency of minutes, and its certainty is absolute rather than statistical, since a violated range check is a fact and not an inference.
• Input runs feature PSI and the domain classifier over hours. Statistical rather than certain.
• Output watches the predicted-positive rate, the mean score and the mean confidence.
• Action watches the block and approval rates per segment.
• Proxy outcome uses a fast label standing in for the slow one, at a latency of days.
• True outcome computes log loss and AUC on real labels at a latency of weeks to months, which makes it an audit and never a page.

flowchart TD
    P["PIPELINE<br/>schema · nulls · ranges<br/>latency: minutes<br/>certainty: absolute"] --> I["INPUT<br/>feature PSI · domain classifier<br/>latency: hours<br/>certainty: statistical"]
    I --> O["OUTPUT<br/>predicted-positive rate · mean score<br/>mean confidence<br/>latency: hours"]
    O --> AC["ACTION<br/>block and approval rate<br/>PER SEGMENT<br/>latency: hours"]
    AC --> PX["PROXY OUTCOME<br/>fast label standing in<br/>for the slow one<br/>latency: days"]
    PX --> T["TRUE OUTCOME<br/>log loss · AUC on real labels<br/>latency: weeks to months<br/>AUDIT, never a page"]

    style P fill:#2d6a4f,color:#fff
    style PX fill:#40916c,color:#fff
    style T fill:#9d0208,color:#fff

The diagram gives the order; the table gives the cadence and the rule. Same six layers — what only the table carries is how often each one runs and what specifically fires it, and those two columns are the entire difference between a monitoring design and a monitoring diagram.

Two conventions in that Alert rule column are worth naming.

Sigma is the Greek letter for a standard deviation, so “3 sigma versus the trailing 4 weeks” means “further from the recent norm than three standard deviations of that norm.” Under ordinary noise that happens about once in 370 observations.

Hour-of-week comparison means you compare Tuesday 3 p.m. against previous Tuesdays at 3 p.m. rather than against yesterday, so that ordinary weekly seasonality does not read as drift.

Alert on the pipeline before the model. The most common cause of a “drift” alert is an upstream change, and those are catchable in minutes with certainty rather than in weeks with statistics.

Here is one, minute by minute, with the same incident replayed underneath as it would have gone without the range check:

09:14  upstream job changes `days_since_signup` from days to seconds
09:15  null rate 0.0% (no alarm), range check FAILS: max 315,360,000 vs training max 3,650
09:15  range alert -> rollback.  Total exposure: 1 minute.

without the range check:
  next day    feature PSI = 4.1 (flagged, 15 hours later)
  +45 days    accuracy moves

The unit change made days_since_signup 86,400x larger, which the range check catches instantly because 315,360,000 is far outside the training range of 3,650. PSI needs a day to notice; accuracy needs 45.

Finally, prefer scheduled retraining on a cadence shorter than your drift timescale over drift-triggered retraining. Triggered retraining chases noise, and it creates a feedback loop of its own, since each retrain is fitted on data that the previous model shaped.

Whatever the cadence, always compare the candidate model against the incumbent on two holdout sets rather than one: a frozen historical holdout, which catches regressions, and a fresh recent one, which catches staleness.

A candidate that wins on the fresh set and loses on the frozen one has learned the drift. That may be exactly what you wanted, or it may be overfitting to a two-week anomaly, and the two look identical until you ask which it is.

7. Cheat sheet

The whole chapter, compressed into a lookup table: a symptom you might observe on the left, the mechanism that produces it in the middle, and the fix on the right.

Every row is derived somewhere above. If a mechanism reads as an assertion, the section that proves it is linked in the corresponding fix.

The chapter in three sentences:

1. Imbalance, resampling, class weights and label shift are all the same prior moving, and all are repaired by logit(p) = logit(p') + ln(c).
2. Calibration is a separate axis from ranking — broken by zero-training-error logit growth in deep nets and by shrinkage in boosted trees, repaired by a monotone map that cannot change AUC, and worth nothing unless something downstream reads the number.
3. Drift splits into three types whose only practical distinction is that two are detectable without labels and the third — the dangerous one — is not, which is why the thing to build is a fast proxy label rather than a better detector.

Next: 08 — Reinforcement Learning — decisions that change the distribution they are evaluated on, taken deliberately this time.