Probability Calibration

The whole of the previous lesson rested on the ranking. The threshold was swept, two curves were printed, the operating point was chosen with the cost calculation; at no step did the number the model says have to have a counterpart of its own. If the scores were passed through a strictly increasing transform, the ranking, the curves and the chosen set of subscribers would have stayed unchanged.

But when the model says 0.80 for a subscriber, that number claims something: about eighty percent of the subscribers given 0.80 should turn out to be suspicious. This claim is independent of the ranking and is measured separately. Calibration is the difference between the stated probability and the realized share.

• ET23. Calibration is a measure independent of the ranking; two models can give the same ranking and produce different calibration errors.
• ET24. The probabilities are split into five bins of equal width. The number of bins is a choice: as the number grows, the number of subscribers in each bin falls and the realized share becomes volatile.
• ET25. The calibration error is computed by summing, in each bin, the absolute difference between the sum of the stated probabilities and the number of suspicious subscribers realized, and dividing by the number of observations.
• ET26. The sharpened score is a strictly increasing transform; it preserves the ranking, and therefore the area under the curve, exactly.
• ET27. The calibration mapping sees no feature other than the score of the model; its only input is the score and its two parameters are the slope and the shift. As long as the slope stays positive, the ranking does not change.
• ET28. The mapping is learned from the validation set; the other two sets are printed for comparison.
• ET29. The leaf probability of a tree is the share of suspicious subscribers among the training rows that fall into that leaf; for this reason the calibration error of the tree on the training set is zero by definition.
• ET30. The baseline gives a number here too: the procedure that tells every subscriber the suspicious share of the training set is measured on calibration, and its area under the curve is 0.5000.

Stated and Realized

To measure calibration the probabilities are split into bins, and in each bin the mean probability the model states and the realized suspicious share of the subscribers in that bin are written side by side. This table is called the reliability curve, and it need not be drawn; it is read as numbers.

# evaluation.py — MODEL: the FICTIONAL subscriber table is rebuilt with the same seed.
import math

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [("north", 0.28, 21, 0.00), ("south", 0.22, 17, -0.10),
          ("east", 0.18, 26, 0.30), ("west", 0.14, 14, -0.05),
          ("central", 0.18, 23, 0.05)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


SUBSCRIBERS, READINGS = [], {}
for i in range(RAW):
    r = rng(SEED + i)
    b = REGION[pick(r(), [x[1] for x in REGION])]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() >= 0.046:
        SUBSCRIBERS.append({"id": 10001 + i, "region": b, "household": household,
                            "satisfaction": rating})
for k in SUBSCRIBERS:
    r = rng(SEED + 7000 + k["id"])
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        READINGS.setdefault(k["id"], []).append(0.0 if r() < 0.038 else math.floor(
            k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - n * 0.05) * 100 + 0.5) / 100)
        r()

DATA = []
for k in SUBSCRIBERS:
    v = READINGS.get(k["id"])
    if not v:
        continue
    mean, vol, r = sum(v) / len(v), max(v) - min(v), rng(SEED + 51000 + k["id"])
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol
         + 0.24 * (k["household"] >= 5) - 0.20 * (k["satisfaction"] >= 4) + k["region"][3]
         + 0.9 * (k["id"] - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"household": k["household"], "satisfaction": k["satisfaction"], "period": len(v),
         "mean_consumption": round(mean, 2), "volatility": round(vol, 2),
         "zero_readings": int(min(v) == 0.0), "suspicious": int(z > 1.35)}
    x["per_person"] = round(x["mean_consumption"] / x["household"], 3)
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


K = shuffle(DATA, SEED + 90000)
TRAIN, VAL, TEST = K[:756], K[756:1008], K[1008:]
FIELDS = ["mean_consumption", "volatility", "household", "satisfaction", "period",
          "zero_readings", "per_person", "area_north", "area_south", "area_east", "area_west",
          "area_central"]


def impurity(rows):
    p = sum(x["suspicious"] for x in rows) / len(rows)
    return 2 * p * (1 - p)


def tree(rows, depth, fields, min_size=2):        # MODEL: decision tree
    n = len(rows)
    p = sum(x["suspicious"] for x in rows) / n
    best = None
    if depth and n >= 2 * min_size and 0.0 < p < 1.0:
        for f in fields:
            values = sorted({x[f] for x in rows})
            for v in (values[1:] if len(values) < 10 else
                      [values[int(i * len(values) / 10)] for i in range(1, 10)]):
                left = [x for x in rows if x[f] < v]
                right = [x for x in rows if x[f] >= v]
                if min(len(left), len(right)) < min_size:
                    continue
                gain = impurity(rows) - (len(left) * impurity(left)
                                         + len(right) * impurity(right)) / n
                if best is None or gain > best[0]:
                    best = (gain, f, v, left, right)
    if best is None or best[0] <= 1e-9:
        return {"p": p}
    return {"field": best[1], "threshold": best[2],
            "left": tree(best[3], depth - 1, fields, min_size),
            "right": tree(best[4], depth - 1, fields, min_size)}


def probability(node, x):
    while "p" not in node:
        node = node["left"] if x[node["field"]] < node["threshold"] else node["right"]
    return node["p"]


def linear(train, fields, w, epochs=300, step=0.5):   # MODEL: gradient-based linear model
    mu = {f: sum(x[f] for x in train) / len(train) for f in fields}
    sd = {f: (sum((x[f] - mu[f]) ** 2 for x in train) / len(train)) ** 0.5 or 1.0
          for f in fields}
    scale = lambda x: [(x[f] - mu[f]) / sd[f] for f in fields]
    coef, b = [0.0] * len(fields), 0.0
    X = [scale(x) for x in train]
    y = [x["suspicious"] for x in train]
    sw = [(w if x["suspicious"] else 1.0) for x in train]
    total = sum(sw)
    for _ in range(epochs):
        gc, gb = [0.0] * len(fields), 0.0
        for i, xi in enumerate(X):
            z = b + sum(coef[j] * xi[j] for j in range(len(coef)))
            e = (1 / (1 + math.exp(-max(-30.0, min(30.0, z)))) - y[i]) * sw[i]
            gb += e
            for j in range(len(coef)):
                gc[j] += e * xi[j]
        b -= step * gb / total
        for j in range(len(coef)):
            coef[j] -= step * gc[j] / total

    def score(x):
        v = scale(x)
        z = b + sum(coef[j] * v[j] for j in range(len(coef)))
        return 1 / (1 + math.exp(-max(-30.0, min(30.0, z))))
    return score


SCORE_FN = linear(TRAIN, FIELDS, 1.0)        # MODEL: probability-producing gradient-based linear model
BINS = 5


def reliability(score, rows, n=BINS):     # per bin: count, stated, realized
    bins = [[0, 0.0, 0] for _ in range(n)]
    for x in rows:
        p = score(x)
        k = bins[min(n - 1, int(p * n))]
        k[0], k[1], k[2] = k[0] + 1, k[1] + p, k[2] + x["suspicious"]
    return bins


def calibration_error(score, rows, n=BINS):    # bin-weighted mean absolute difference
    return sum(abs(k[1] - k[2]) for k in reliability(score, rows, n)) / len(rows)


def print_reliability(name, score, rows):
    print(f"{name}  (calibration error {calibration_error(score, rows):.4f})")
    for i, k in enumerate(reliability(score, rows)):
        if k[0]:
            print(f"  {i / BINS:.1f}-{(i + 1) / BINS:.1f}  subscribers {k[0]:>3}  "
                  f"said {k[1] / k[0]:.4f}  realized {k[2] / k[0]:.4f}  "
                  f"difference {k[1] / k[0] - k[2] / k[0]:+.4f}")


print_reliability("linear model, test set", SCORE_FN, TEST)

linear model, test set  (calibration error 0.0242)
  0.0-0.2  subscribers 152  said 0.0891  realized 0.1053  difference -0.0161
  0.2-0.4  subscribers  53  said 0.2837  realized 0.3019  difference -0.0182
  0.4-0.6  subscribers  23  said 0.4918  realized 0.4783  difference +0.0136
  0.6-0.8  subscribers  16  said 0.7079  realized 0.6250  difference +0.0829
  0.8-1.0  subscribers   8  said 0.8698  realized 1.0000  difference -0.1302

In the first three bins the difference is below the 0.02 level and the calibration error is 0.0242. The last two bins are more volatile, because there are a total of 24 subscribers there: in the 0.6-0.8 bin the model says 0.7079 while 0.6250 is realized, in the 0.8-1.0 bin it says 0.8698 while all eight of the eight subscribers turn out to be suspicious. A deviation of 0.13 in a share computed from eight subscribers is the volatility of the measure itself; as the number of bins grows this volatility grows, as it shrinks the table goes blind.

Good Ranking, Bad Calibration

If the scores are passed through a strictly increasing transform, no subscriber’s position changes. Every measure that rests on the ranking stays constant; calibration does not.

def area_roc(score, rows):                  # rank sum averaging ties
    p = sorted(((score(x), x["suspicious"]) for x in rows))
    P = sum(y for _, y in p)
    t, i = 0.0, 0
    while i < len(p):
        j = i
        while j < len(p) and p[j][0] == p[i][0]:
            j += 1
        t += (i + j + 1) / 2.0 * sum(y for _, y in p[i:j])
        i = j
    return (t - P * (P + 1) / 2) / (P * (len(p) - P))


def sharpen(p, g=3.0):                      # strictly increasing, rank-preserving transform
    return p ** g / (p ** g + (1 - p) ** g)


SHARP = lambda x: sharpen(SCORE_FN(x))      # MODEL: sharpened score
TREE_D8 = tree(TRAIN, 8, FIELDS)            # MODEL: depth 8 decision tree
BASELINE = sum(x["suspicious"] for x in TRAIN) / len(TRAIN)
print(f"{'model':<25}{'ROC':>8}{'cal.error':>10}{'mean score':>11}{'suspicious share':>17}")
for name, s in (("linear", SCORE_FN), ("sharpened", SHARP),
                ("tree d8", lambda x: probability(TREE_D8, x)),
                ("baseline (training rate)", lambda x: BASELINE)):
    print(f"{name:<25}{area_roc(s, TEST):>8.4f}{calibration_error(s, TEST):>10.4f}"
          f"{sum(s(x) for x in TEST) / len(TEST):>11.4f}"
          f"{sum(x['suspicious'] for x in TEST) / len(TEST):>17.4f}")
print()
print_reliability("sharpened score, test set", SHARP, TEST)

model                         ROC cal.error mean score suspicious share
linear                     0.8058    0.0242     0.2309           0.2421
sharpened                  0.8058    0.1512     0.1504           0.2421
tree d8                    0.7391    0.1435     0.2586           0.2421
baseline (training rate)   0.5000    0.0198     0.2222           0.2421

sharpened score, test set  (calibration error 0.1512)
  0.0-0.2  subscribers 204  said 0.0199  realized 0.1569  difference -0.1369
  0.2-0.4  subscribers  10  said 0.3033  realized 0.3000  difference +0.0033
  0.4-0.6  subscribers   7  said 0.4756  realized 0.8571  difference -0.3816
  0.6-0.8  subscribers   8  said 0.7001  realized 0.3750  difference +0.3251
  0.8-1.0  subscribers  23  said 0.9509  realized 0.7391  difference +0.2117

The first two rows of the first table are the proof of the claim: the area under the curve of the two models is the same, 0.8058, because both order the 252 subscribers in exactly the same sequence. The calibration error, on the other hand, rises from 0.0242 to 0.1512, to more than six times.

The third row carries a separate warning. The mean score of the depth 8 tree is 0.2586 and the real suspicious share on the test set is 0.2421 — the difference between them is only 0.0165. A model that looks right in total gives an error of 0.1435 at the bin level. A mean that holds is not calibration; it hides overconfidence and underconfidence that cancel each other.

The fourth row says what calibration alone does not measure. The baseline that tells every subscriber the suspicious share of the training set, 0.2222, has a calibration error of 0.0198, the smallest number in the table; its area under the curve is 0.5000, that is, its ranking is the same as random. A procedure that has learned nothing can come first in calibration. The two measures answer separate questions and both need to be written.

The second table shows where the sharpened score breaks. In the most crowded bin the model says 0.0199 and 0.1569 is realized — more than seven times. At the top end it says 0.9509 and 0.7391 is realized. The model is overly timid on one side and overly confident on the other; the two errors cancel each other somewhat on average, at the bin level both stand.

The Correction: A Single-Variable Mapping

The breakage is not in the ranking but in the numeric value of the score. The correction should not touch the ranking either: a mapping that takes the score as input, is strictly increasing and has two parameters is enough.

def logit(p, e=1e-6):
    p = min(1 - e, max(e, p))
    return math.log(p / (1 - p))


def mapping(score, rows, epochs=600, step=0.5):  # MODEL: single-variable calibration mapping
    z = [logit(score(x)) for x in rows]
    y = [x["suspicious"] for x in rows]
    a, b = 1.0, 0.0
    for _ in range(epochs):
        ga = gb = 0.0
        for i, zi in enumerate(z):
            e = 1 / (1 + math.exp(-max(-30.0, min(30.0, a * zi + b)))) - y[i]
            ga, gb = ga + e * zi, gb + e
        a, b = a - step * ga / len(z), b - step * gb / len(z)
    return lambda x: 1 / (1 + math.exp(-max(-30.0, min(30.0, a * logit(score(x)) + b)))), a, b


print(f"{'mapping learned on':<28}{'a':>8}{'b':>9}{'ROC':>8}{'cal.error':>10}")
print(f"{'no correction':<28}{'-':>8}{'-':>9}{area_roc(SHARP, TEST):>8.4f}"
      f"{calibration_error(SHARP, TEST):>10.4f}")
for name, rows in (("train", TRAIN), ("validation", VAL), ("test", TEST)):
    d, a, b = mapping(SHARP, rows)
    print(f"{name:<28}{a:>8.4f}{b:>9.4f}{area_roc(d, TEST):>8.4f}"
          f"{calibration_error(d, TEST):>10.4f}")
TREE_FN = lambda x: probability(TREE_D8, x)
d8, a8, b8 = mapping(TREE_FN, TRAIN)
print(f"\ntree d8 calibration error: train {calibration_error(TREE_FN, TRAIN):.4f}, "
      f"validation {calibration_error(TREE_FN, VAL):.4f}, test {calibration_error(TREE_FN, TEST):.4f}")
print(f"mapping learned on train: a {a8:.4f}, b {b8:.4f}, calibration error on test "
      f"{calibration_error(d8, TEST):.4f}")

mapping learned on                 a        b     ROC cal.error
no correction                      -        -  0.8058    0.1512
train                         0.3340   0.0015  0.8058    0.0244
validation                    0.3013  -0.3695  0.8058    0.0463
test                          0.3038  -0.0043  0.8058    0.0208

tree d8 calibration error: train 0.0000, validation 0.1624, test 0.1435
mapping learned on train: a 1.0002, b 0.0001, calibration error on test 0.1435

The slope of the mapping is around 0.30 on all three sets: as far as the sharpening pushed the scores toward the extremes, the mapping takes that much back. The area under the curve is 0.8058 on all four rows — as long as the slope is positive the ranking does not change and the correction does not move any subscriber from its place. The calibration error falls from 0.1512 to 0.0463 with the validation set, that is, to a third.

The mapping learned from the test set gives 0.0208; this number is not a performance but the result of the correction having been learned from the set it is measured on, and the 0.0255 difference between it and 0.0463 is not an estimate but an optimism.

The training row comes out well on this set, 0.0244; because the linear model does not fit the training set tightly, its probabilities there are not overconfident either. The reason for the rule is in the last two lines. The calibration error of the depth 8 tree is 0.0000 on the training set and 0.1435 on the test set — the leaf probability of the tree is the training share of that leaf, so on the training set the stated and the realized are equal by definition. The mapping learned from training finds a slope of 1.0002 and a shift of 0.0001 and carries the error on the test set from 0.1435 to 0.1435: it has seen nothing to correct, because there is nothing to correct in the set it looked at.

Summary

• Calibration is the difference between the stated probability and the realized share and is measured by splitting the probabilities into bins; the calibration error of the linear model on the test set is 0.0242.
• A strictly increasing transform preserves the ranking: the area under the curve of the sharpened score is also 0.8058, while its calibration error is 0.1512, more than six times.
• It is not enough for the mean score to be close to the real share: the mean of the depth 8 tree is 0.2586, the real share is 0.2421, and its error at the bin level is 0.1435. The baseline that says a constant 0.2222, on the other hand, is the best in calibration with 0.0198 and the worst in the area under the curve with 0.5000.
• The single-variable mapping learned from validation lowers the calibration error from 0.1512 to 0.0463 and leaves the area under the curve at 0.8058; because the slope is positive, no position in the ranking changes.
• The calibration error of the tree on the training set is 0.0000 by definition; the mapping learned from there converges to the identity mapping and leaves the error of 0.1435 on the test set as it is.

Next Step

In four lessons each of the measures was for classification and all of them came out of a single operation: counting. The confusion matrix was four counts, precision and recall two ratios, the curves a pair of counts per threshold, calibration a comparison of counts per bin. Because the target was a class, right and wrong were two discrete boxes. In the same fictional set, when the target is the mean consumption of period 3, there is no box left to count: the model says 22.40, the truth is 26.15, and this estimate is neither right nor wrong. The next lesson turns the measure into distance and rebuilds the baseline as the training mean.