Bias–Variance Tradeoff

The previous lesson raised capacity and showed two numbers diverging, but the entire measurement was made on a single training set. Yet those 756 rows are one sample drawn from among 1,260 subscribers. On another draw another tree is built, other thresholds are chosen and a different answer is produced for the same subscriber. Part of the error on the held-out set comes precisely from here: not from the model being wrong, but from the model swinging.

This lesson splits the error into two components. The first is how much the same procedure swings as the training sample changes: variance. The second is how far the mean prediction, taken over all the samples, deviates from the true value: bias. Both are tied to the same knob, capacity, in opposite directions, and their sum is smallest in the middle.

• GE16. In this lesson the problem is regression: estimating a numeric value. In the Statistical Analysis and Exploratory Data Analysis course of the Data Analytics curriculum, “regression” was used to interpret the coefficient of a relationship; here the measure is not a coefficient but the error of the number produced for a subscriber that has not been seen.
• GE17. The target is the subscriber’s mean consumption (m³) and is derived from the readings. The features contain no field derived from the readings: region, household size, satisfaction, number of periods.
• GE18. The measure is the root mean squared error (m³). The decomposition is made in squared units, because additivity holds only there.
• GE19. The set, the split and the seed continue unchanged from the previous lessons: 1,260 subscribers, 756 training, 252 held-out. The test set is again not opened.
• GE20. Sampling with replacement was set up in the same course of the Data Analytics curriculum and the method is not repeated here. Its only use is to measure as if there were thirty training sets while only one is at hand.
• GE21. For a fixed subscriber and $B$ models the following equality is exact: the mean squared error is the sum of the scatter of the predictions among themselves and the square of the distance of the mean prediction from the true value. The three columns in the table are this equality.
• GE22. The third column carries the bias and the irreducible noise together; with a single observation the two cannot be separated. Because the noise is independent of capacity, the differences in that column are differences in bias.

The Same Two Numbers in Regression

The setup generates the same regional measurement network (a fictional setup) with the same seed, takes the target as the subscriber’s mean consumption and builds a regression tree. The tree has the same structure as the classification tree; the only differences are that the leaf says the mean instead of the majority class and that the split is chosen not by class purity but by squared deviation.

# regression.py — MODEL. The regional measurement network (a fictional setup) of the Data Analytics curriculum is generated with the same seed;
# the target is the subscriber's mean consumption and is derived from the readings.
import math

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


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 = []
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})
DATA = []
for k in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + k["id"]), []
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.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()
    if not v:
        continue
    x = {"id": k["id"], "household": k["household"], "satisfaction": k["satisfaction"],
         "period": len(v), "target": round(sum(v) / len(v), 2)}
    for b in REGION:
        x["reg_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)

FIELDS = ["household", "satisfaction", "period",
          "reg_north", "reg_south", "reg_east", "reg_west", "reg_central"]


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, HELD = K[:756], K[756:1008]


def mean_sq_dev(s):
    m = sum(x["target"] for x in s) / len(s)
    return sum((x["target"] - m) ** 2 for x in s) / len(s)


def reg_tree(s, depth, min_size=5):
    best = None
    if depth and len(s) >= 2 * min_size:
        for f in FIELDS:
            values = sorted({x[f] for x in s})
            for v in (values[1:] if len(values) < 9 else
                      [values[int(i * len(values) / 9)] for i in range(1, 9)]):
                left = [x for x in s if x[f] < v]
                right = [x for x in s if x[f] >= v]
                if min(len(left), len(right)) < min_size:
                    continue
                k = (mean_sq_dev(s)
                     - (len(left) * mean_sq_dev(left) + len(right) * mean_sq_dev(right)) / len(s))
                if best is None or k > best[0]:
                    best = (k, f, v, left, right)
    if best is None or best[0] <= 1e-9:
        return {"prediction": sum(x["target"] for x in s) / len(s)}
    return {"field": best[1], "threshold": best[2],
            "left": reg_tree(best[3], depth - 1, min_size),
            "right": reg_tree(best[4], depth - 1, min_size)}


def predict(d, x):
    while "prediction" not in d:
        d = d["left"] if x[d["field"]] < d["threshold"] else d["right"]
    return d["prediction"]


def rmse(f, s):
    return math.sqrt(sum((f(x) - x["target"]) ** 2 for x in s) / len(s))


MEAN = sum(x["target"] for x in TRAIN) / len(TRAIN)
print(f"train {len(TRAIN)}, held-out {len(HELD)}, train mean {MEAN:.3f} m3")
print(f"\n{'procedure':<14} {'train':>9} {'held-out':>9}")
print(f"{'baseline':<14} {rmse(lambda x: MEAN, TRAIN):>9.3f} {rmse(lambda x: MEAN, HELD):>9.3f}")
for dp in (1, 3, 12):
    m = reg_tree(TRAIN, dp)
    print(f"{'depth ' + str(dp):<14} {rmse(lambda x: predict(m, x), TRAIN):>9.3f} "
          f"{rmse(lambda x: predict(m, x), HELD):>9.3f}")

train 756, held-out 252, train mean 19.767 m3

procedure          train  held-out
baseline           6.200     6.419
depth 1            5.696     5.979
depth 3            5.105     5.380
depth 12           4.639     5.594

The baseline here is “saying the mean”: every subscriber is assigned 19.767 m³, the mean of the training set. It gives an error of 6.419 m³ on the held-out set. The depth 3 tree lowers this to 5.380 m³. The depth 12 tree is best on the training set with 4.639 m³ but rises to 5.594 m³ on the held-out set — the shape of the previous lesson, this time in cubic meters.

The table does not say what has gone wrong. Is the depth 12 tree aiming at a wrong place, or does it aim at the right place and land somewhere different every time? The corrections for the two cases are different, and the two must be separated.

Splitting the Error in Two

The way to separate them is to multiply the single training set. Sampling with replacement builds 30 different sets of 756 rows from the same training set; on each of them the same procedure is run and 30 separate answers are produced for the subscribers in the same held-out set.

B = 30
SAMPLES = []
for b in range(B):
    r = rng(SEED + 61000 + b * 13)
    SAMPLES.append([TRAIN[int(r() * len(TRAIN))] for _ in range(len(TRAIN))])


def decompose(name, fit, predictor):
    MODELS = [fit(s) for s in SAMPLES]
    PRED = [[predictor(m, x) for x in HELD] for m in MODELS]
    total = variance = bias = 0.0
    for i, x in enumerate(HELD):
        p = [PRED[b][i] for b in range(B)]
        m = sum(p) / B
        total += sum((q - x["target"]) ** 2 for q in p) / B
        variance += sum((q - m) ** 2 for q in p) / B
        bias += (m - x["target"]) ** 2
    n = len(HELD)
    print(f"{name:<14} {total / n:>9.3f} {variance / n:>9.3f} {bias / n:>9.3f} "
          f"{math.sqrt(total / n):>7.3f}")


print(f"{B} samples drawn with replacement, {len(TRAIN)} rows each")
print(f"\n{'procedure':<14} {'total':>9} {'variance':>9} {'bias':>9} {'rmse':>7}")
decompose("baseline", lambda s: sum(y["target"] for y in s) / len(s), lambda m, x: m)
for dp in (1, 2, 3, 5, 8, 12):
    decompose(f"depth {dp}", lambda s, dp=dp: reg_tree(s, dp), predict)

30 samples drawn with replacement, 756 rows each

procedure          total  variance      bias    rmse
baseline          41.242     0.040    41.202   6.422
depth 1           36.158     2.138    34.020   6.013
depth 2           31.840     2.218    29.622   5.643
depth 3           29.768     1.321    28.448   5.456
depth 5           30.745     2.274    28.471   5.545
depth 8           32.922     3.983    28.939   5.738
depth 12          33.810     4.409    29.400   5.815

The second and third columns give the first exactly in every row; the equality is not an approximation but the definition itself. The columns move in opposite directions. The variance rises from 0.040 to 4.409, that is, by more than a hundredfold. The bias and noise fall from 41.202 to 28.448 and stop there: after depth 3 this column does not fall, it rises slightly.

The baseline is one end of the tradeoff and its place in the table shows it. Its variance is almost zero at 0.040 — whichever sample it is trained on, it says almost the same number. The price of this is a bias of 41.202: it always stands in the same place and is always wrong in the same place. Stability is not a virtue on its own.

The floor of the third column is the most important number of the lesson. Whatever is done with capacity, this column does not go below 28.4 m³ squared; this is the remaining noise that no feature can explain. The bias that capacity actually removes is 12.754 m³ squared between 41.202 and 28.448, that is, 30.9% of the error of the baseline. The remaining 69.1% is removed neither by a deeper tree nor by more data; it is the measurement network’s own scatter.

The third column has a second reading, and it gives the rationale of a family of methods. That column is the error on the held-out set of a procedure that takes the mean of the thirty models: when a mean is taken, by definition the variance component vanishes and only the third column remains. For depth 12, a single tree makes a squared error of 33.810 m³ while the mean of thirty trees makes 29.400, that is, 13.0% of the error is erased merely by averaging; on the root scale, 5.422 m³ instead of 5.815 m³. This number is smaller even than the 5.456 m³ of a single depth 3 tree. Ensemble methods do exactly this: instead of trying to lower the bias, they multiply high-capacity, swinging models and leave the variance to the average. The way to lower the error without touching the capacity knob at all passes through here.

The variance column is not a one-directional function of capacity, and the reason can be read. At depth 3 the variance, 1.321, is below the 2.218 of depth 2. At depth 1 and 2 the tree has to separate five regions with one or two splits, and which region it will separate changes from sample to sample; at depth 3 the five regions get their own leaf in almost every sample and the thirty models resemble each other. When capacity is insufficient, the swinging comes not from what the model has learned but from what it has been forced to choose.

The Same Subscriber, Thirty Models

The components are averages; what they look like in a single subscriber is asked separately.

print(f"{'capacity':<12} {'mean range':>17} {'widest subscriber':>17}")
MODELS = {}
for dp in (1, 3, 12):
    MODELS[dp] = [reg_tree(s, dp) for s in SAMPLES]
    a = [max(t) - min(t) for t in
         ([predict(m, x) for m in MODELS[dp]] for x in HELD)]
    print(f"{'depth ' + str(dp):<12} {sum(a) / len(a):>17.3f} {max(a):>17.3f}")

print(f"\n{'subscriber':>10} {'region':<8} {'actual':>7}", end="")
for dp in (1, 3, 12):
    print(f" {'d' + str(dp) + ' mean':>12} {'range':>7}", end="")
print()
for x in (HELD[3], HELD[11], HELD[40]):
    b = [k for k in ("north", "south", "east", "west", "central") if x["reg_" + k]][0]
    print(f"{x['id']:>10} {b:<8} {x['target']:>7.2f}", end="")
    for dp in (1, 3, 12):
        p = [predict(m, x) for m in MODELS[dp]]
        print(f" {sum(p) / len(p):>12.2f} {max(p) - min(p):>7.2f}", end="")
    print()

capacity            mean range widest subscriber
depth 1                  3.587             5.246
depth 3                  4.595            16.362
depth 12                 8.172            24.236

subscriber region    actual      d1 mean   range      d3 mean   range     d12 mean   range
     10227 north      18.92        20.25    2.87        20.52    2.92        21.36    6.54
     10683 west       11.99        15.66    5.25        13.80    5.16        13.72    5.48
     10424 south      26.22        20.25    2.87        16.94    4.05        16.15    9.17

In the first row the bias stands bare. At depth 1, subscriber 10227 from the north and subscriber 10424 from the south get the same number: 20.25 m³. Although the true values of the two subscribers are 18.92 and 26.22 m³, the model cannot tell them apart, because its single split leaves these two regions in the same leaf. Its range is narrow at 2.87 — it is stably wrong.

At depth 3 the prediction for 10424 drops to 16.94 m³, and this is not a worse but a more accurate number: it is close to the true mean of the south region. That the subscriber’s own value is 26.22 m³ is not a sign of bias but the very thing that noise is — that subscriber’s readings came out above those of its region, and no feature in the set announces this. To expect a model to know this is to expect it to remove the floor of the third column.

At depth 12 the mean prediction stays almost the same at 16.15, but the range rises from 4.05 to 9.17 m³. The same thing happens over the whole held-out set: the mean range is 3.587 at depth 1 and 8.172 m³ at depth 12. In the widest subscriber the answers of the thirty models spread over a range of 24.236 m³, and this is larger than the target’s own mean of 19.767 m³. The answer given by high capacity is not wrong, it is random; it depends on which training sample was drawn.

Summary

• The mean squared error on the held-out set is the exact sum of two components: the scatter of the predictions from sample to sample and the distance of the mean prediction from the true value.
• As capacity rises the variance goes from 0.040 up to 4.409 and the bias and noise from 41.202 down to 28.448; the total is smallest at depth 3 with 29.768 and its root measure is 5.456 m³.
• The floor of the third column is 28.4 m³ squared and is irreducible noise; the bias removed by capacity is 12.754 m³ squared, that is, 30.9% of the baseline’s error. The same column is the error of an ensemble method that takes the mean of thirty models and brings 33.810 down to 29.400 at depth 12.
• The baseline is one end of the tradeoff: its variance is almost zero at 0.040 and its bias is the largest at 41.202; stability alone is not a virtue.
• The variance is not a one-directional function of capacity; at depth 1 and 2 the swinging is greater than at depth 3, because which region the tree will separate changes from sample to sample.
• At the level of a single subscriber, depth 1 gives the subscribers of two different regions the same number, 20.25 m³, while depth 12 produces a range of 9.17 m³ among thirty models.

Next Step

The tradeoff has been found but the way it was found is crude: depth is an integer and there is no intermediate value between 3 and 5. Instead of cutting capacity off, penalizing it gives a finer knob — the complexity of the model is added to the error measure as a cost and the penalty coefficient is tuned continuously. The next lesson raises the penalty coefficient step by step, prints that the training number gets worse without interruption while the held-out number first improves and then gets worse, and shows from which set the best point is read.