Ensemble Methods

The previous two lessons multiplied the same family in two different ways. Bagging built the trees in parallel and had them vote and gave 0.7897 on the test set; boosting stacked them into a sequential sum and gave 0.7738. Both were multiplying a single type of learner. The third form is to build a second layer that takes the outputs of different models as its input, and it immediately raises a question: from which set will that layer be learned.

This lesson puts the three combinations side by side. The measure is three numbers: the difference on the held-out set, the cost as the number of splits, and the candidates tried. A fourth measure is added, the similar errors of the ensembles: if three models err on the same subscribers, the combination has nothing to gain, and this is not an intuition but a ceiling that can be computed.

• NM45. The set, the split, the baseline and the features are the same as in the previous two lessons: 1,260 subscribers, 756/252/252, seed 20260218, the baseline is 0.7579 on the test set and 0.8135 on validation.
• NM46. The three base models are carried over as they are from the previous lessons and are not retuned here: a depth-6 single tree, a bagged ensemble of 41 trees (m six, depth 4) and boosting with rate 0.15 / depth 2 / 10 rounds.
• NM47. The output of each base model is in the range [0, 1]: the leaf ratio, the vote share, the boosting probability. The decision threshold is 0.5 and is fixed from the Data Preparation and Feature Engineering course.
• NM48. The second layer of stacking (the meta-learner) is a gradient-based linear model that takes the three outputs as input: three weights and a constant, 400 iterations, rate 2.0. Its settings were fixed before the measurement.
• NM49. The cost measure is the number of splits; the cost of stacking is the splits of the three base models plus the four learned numbers.
• NM50. The shared-error count is done on the test set and at the label level.
• NM51. The ceiling of a combination that looks at the label pattern is the accuracy of the rule that, in each pattern, says that pattern’s majority label; no combination that sees nothing but the three labels can pass it.
• NM52. No new tuning is searched in this lesson; the candidates tried are three. The topic’s total budget, with the carried numbers, is 1,313 candidates.

Three Combinations and Three Base Models

The setup generates the table, rebuilds two of the combinations and prints the numbers of the three base models together with their cost.

# stacking.py — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning course is rebuilt with the same seed;
# the feature set comes from the Data Preparation and Feature Engineering course and is kept fixed throughout the course.
import bisect
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, DATA = [], []
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, 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
    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)}
    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",
          "area_north", "area_south", "area_east", "area_west", "area_central"]
THRESHOLDS = {}
for a in FIELDS:
    d = sorted({x[a] for x in TRAIN})
    THRESHOLDS[a] = d[1:] if len(d) <= 24 else [d[int(i * len(d) / 25)] for i in range(1, 25)]
BINS = [len(THRESHOLDS[a]) + 1 for a in FIELDS]


def binarize(s):                            # each row: bin numbers and label
    return [(tuple(bisect.bisect_right(THRESHOLDS[a], x[a]) for a in FIELDS), x["suspicious"])
            for x in s]


B_TRAIN, B_VAL, B_TEST = binarize(TRAIN), binarize(VAL), binarize(TEST)
SETS, ALL = (B_TRAIN, B_VAL, B_TEST), list(range(len(FIELDS)))


def tree(sample, depth, columns, least=2):  # MODEL: each example is (bins, gradient, curvature)
    n, best = len(sample), None
    G, H = sum(o[1] for o in sample), sum(o[2] for o in sample)
    if depth and n >= 2 * least:
        for k in columns:
            sg, sh, sn = [0.0] * BINS[k], [0.0] * BINS[k], [0] * BINS[k]
            for b, g, h in sample:
                sg[b[k]] += g
                sh[b[k]] += h
                sn[b[k]] += 1
            cg = ch = cn = 0
            for j in range(BINS[k] - 1):
                cg, ch, cn = cg + sg[j], ch + sh[j], cn + sn[j]
                if cn < least or n - cn < least:
                    continue
                gain = cg * cg / ch + (G - cg) ** 2 / (H - ch) - G * G / H
                if best is None or gain > best[0]:
                    best = (gain, k, j)
    if best is None or best[0] <= 1e-9:
        return G / H
    _, k, j = best
    return (k, j, tree([o for o in sample if o[0][k] <= j], depth - 1, columns, least),
            tree([o for o in sample if o[0][k] > j], depth - 1, columns, least))


def value(d, b):
    while type(d) is tuple:
        d = d[2] if b[d[0]] <= d[1] else d[3]
    return d


def squash(z):
    return 1.0 / (1.0 + math.exp(-z)) if z > -30 else 0.0


def split_count(d):
    return 0 if type(d) is not tuple else 1 + split_count(d[2]) + split_count(d[3])


def forest(seed, count, m, depth):          # COMBINATION 1: bagging, parallel vote
    r, trees = rng(seed), []
    for _ in range(count):
        chosen = [int(r() * len(B_TRAIN)) for _ in range(len(B_TRAIN))]
        pool, columns = list(ALL), []
        for _ in range(min(m, len(ALL))):
            columns.append(pool.pop(int(r() * len(pool))))
        trees.append(tree([(B_TRAIN[i][0], B_TRAIN[i][1], 1.0) for i in chosen],
                          depth, sorted(columns)))
    return trees


def boost(rate, rounds, depth):             # COMBINATION 2: boosting, sequential sum
    p = sum(y for _, y in B_TRAIN) / len(B_TRAIN)
    F, trees = [[math.log(p / (1 - p))] * len(k) for k in SETS], []
    for _ in range(rounds):
        q = [squash(z) for z in F[0]]
        d = tree([(b, y - q[i], max(q[i] * (1 - q[i]), 1e-6))
                  for i, (b, y) in enumerate(B_TRAIN)], depth, ALL)
        trees.append(d)
        for u, subset in enumerate(SETS):
            for i, (b, _) in enumerate(subset):
                F[u][i] += rate * value(d, b)
    return trees, [[squash(z) for z in f] for f in F]


SINGLE = tree([(b, y, 1.0) for b, y in B_TRAIN], 6, ALL)
FOREST = forest(SEED, 41, 6, 4)
BOOSTED, BOOSTED_P = boost(0.15, 10, 2)
BASE = [("single tree", [[value(SINGLE, b) for b, _ in k] for k in SETS], split_count(SINGLE)),
        ("bagging", [[sum(value(x, b) > 0.5 for x in FOREST) / 41 for b, _ in k]
                     for k in SETS], sum(map(split_count, FOREST))),
        ("boosting", BOOSTED_P, sum(map(split_count, BOOSTED)))]


def accuracy(output, subset):
    return sum((c > 0.5) == y for c, (_, y) in zip(output, subset)) / len(subset)


print(f"{'model':<12} {'validation':>10} {'test':>8} {'splits':>7}")
print(f"{'baseline':<12} {sum(y == 0 for _, y in B_VAL) / len(B_VAL):>10.4f} "
      f"{sum(y == 0 for _, y in B_TEST) / len(B_TEST):>8.4f} {0:>7}")
for name, output, cost in BASE:
    print(f"{name:<12} {accuracy(output[1], B_VAL):>10.4f} "
          f"{accuracy(output[2], B_TEST):>8.4f} {cost:>7}")

model        validation     test  splits
baseline         0.8135   0.7579       0
single tree      0.8056   0.7778      43
bagging          0.8413   0.7897     568
boosting         0.8373   0.7738      30

The three rows show three different trade-offs. Bagging gives the highest number, 0.7897, and is the most expensive: 568 splits. Boosting gives 0.7738 with 30 splits, that is, nineteen times cheaper. The single tree stands at 0.7778 with 43 splits and is above boosting. The real question is this: how differently do these three models err.

def labels(u):
    return [tuple(c[u][i] > 0.5 for _, c, _ in BASE) for i in range(len(SETS[u]))]


E = labels(2)
ERRORS = [sum(e[j] != y for j in range(3)) for e, (_, y) in zip(E, B_TEST)]
PATTERNS = {}
for e, (_, y) in zip(E, B_TEST):
    PATTERNS.setdefault(e, [0, 0])[y] += 1
print("models wrong on test  " +
      "  ".join(f"{k}: {ERRORS.count(k):>3}" for k in range(4)))
print(f"ceiling of a combination that sees the label pattern "
      f"{sum(max(v) for v in PATTERNS.values()) / len(B_TEST):.4f}, distinct patterns {len(PATTERNS)}")
for i in range(3):
    for j in range(i + 1, 3):
        differ = sum(e[i] != e[j] for e in E)
        both = sum(e[i] != y and e[j] != y for e, (_, y) in zip(E, B_TEST))
        print(f"{BASE[i][0]:>11} - {BASE[j][0]:<11} classified differently {differ:>3}, "
              f"both wrong {both:>3}")

models wrong on test  0: 185  1:  13  2:   9  3:  45
ceiling of a combination that sees the label pattern 0.7897, distinct patterns 7
single tree - bagging     classified differently  17, both wrong  46
single tree - boosting    classified differently  21, both wrong  46
    bagging - boosting    classified differently   6, both wrong  52

The first line fixes the ceiling. On 185 of the 252 subscribers all three are right, on 45 all three are wrong; the remaining 22 subscribers are the only region where a combination can work. Those 45 subscribers cannot be rescued by any vote, any sum or any second layer, because all three inputs point in the same wrong direction.

The second line turns this into an exact number. Three labels could form eight patterns, only seven are seen on the test set, and the rule that says the majority label in each pattern gives 0.7897. No combination that sees nothing but the three labels can pass this, and this number is equal to bagging’s own number.

The last three lines say why. Bagging and boosting classify only 6 of the 252 subscribers differently. One is a parallel vote and the other a sequential sum; their structures are opposites, their decisions nearly the same. Bagging errs on 53 subscribers, boosting on 57, and the two err together on 52: nearly every subscriber boosting misses, bagging misses too. The single tree is a little further away (17 and 21 subscribers), but it too is fed from the same family.

From Which Set the Second Layer Is Learned

Stacking builds a model that takes the three outputs as input. The data this model will learn from is the outputs of the base models; on which set those outputs were produced is the lesson’s question.

def second_layer(u, rate=2.0, rounds=400):  # COMBINATION 3: stacking
    subset, w, bias = SETS[u], [0.0] * len(BASE), 0.0
    for _ in range(rounds):
        gw, gs = [0.0] * len(BASE), 0.0
        for i, (_, y) in enumerate(subset):
            h = squash(sum(w[j] * BASE[j][1][u][i]
                           for j in range(len(BASE))) + bias) - y
            for j in range(len(BASE)):
                gw[j] += h * BASE[j][1][u][i]
            gs += h
        for j in range(len(BASE)):
            w[j] -= rate * gw[j] / len(subset)
        bias -= rate * gs / len(subset)
    return w, bias


def layer_accuracy(model, u):
    w, bias = model
    return sum((sum(w[j] * BASE[j][1][u][i] for j in range(len(BASE))) + bias > 0)
               == y for i, (_, y) in enumerate(SETS[u])) / len(SETS[u])


RIGHT, WRONG = second_layer(1), second_layer(0)
print(f"{'second layer':<22} {'train':>8} {'validation':>10} {'test':>8}   weights")
for name, model, u in (("fit on validation", RIGHT, 1), ("fit on training", WRONG, 0)):
    print(f"{name:<22}" + "".join(f"{layer_accuracy(model, v):>{(8, 10, 8)[v]}.4f}"
                                  for v in range(3))
          + "   " + " ".join(f"{q:+.2f}" for q in model[0]))
FINAL = [("bagging", accuracy(BASE[1][1][1], B_VAL), accuracy(BASE[1][1][2], B_TEST),
          BASE[1][2], 31),
         ("boosting", accuracy(BASE[2][1][1], B_VAL), accuracy(BASE[2][1][2], B_TEST),
          BASE[2][2], 1279),
         ("stacking", layer_accuracy(RIGHT, 1), layer_accuracy(RIGHT, 2),
          sum(t[2] for t in BASE) + 4, 3)]
print(f"\n{'combination':<12} {'reported':>13} {'test':>8} {'inflation':>9} "
      f"{'splits':>7} {'candidates':>10}")
for name, d, s, cost, candidates in FINAL:
    print(f"{name:<12} {d:>13.4f} {s:>8.4f} {d - s:>9.4f} {cost:>7} {candidates:>10}")
print(f"total tuning budget of the topic {sum(t[4] for t in FINAL)} candidates")

second layer              train validation     test   weights
fit on validation       0.8333    0.8254  0.7778   +1.99 +1.54 +1.45
fit on training         0.8532    0.7937  0.7778   +5.04 +2.38 +0.54

combination       reported     test inflation  splits candidates
bagging             0.8413   0.7897    0.0516     568         31
boosting            0.8373   0.7738    0.0635      30       1279
stacking            0.8254   0.7778    0.0476     645          3
total tuning budget of the topic 1313 candidates

The last column of the two rows is the same: 0.7778. Learning from the wrong set did not change the delivered number in this table. There are two things it did change, and both are visible.

The first is the reported number. When the layer is learned from the training set it reads 0.8532 on that set; the real number is 0.7778, so the inflation is 0.0754. The layer learned from validation reads 0.8254 on its own set and its inflation is 0.0476. Both numbers are inflated, because both were read from the set the layer learned from; the wrong set only enlarged the inflation by 0.0278, that is, seven subscribers.

The second is the weights. The layer learned from validation gives balanced weights: +1.99, +1.54, +1.45. The layer learned from training writes +5.04 to the single tree and lowers boosting to +0.54. The cause is in the training column: the single tree reads 0.8426 on the training set, the bagged ensemble 0.8175, boosting 0.8082. The second layer takes the model that memorized most to be the best model. Because the memorization here is mild, the harm stayed mild; had one of the base models learned the training set exactly, the layer would have given all the weight to it and the delivered number would have dropped too.

Learning from validation has a cost as well: that set is now spent. The way that spends nothing is to learn the second layer from the base models’ outputs on rows that did not train them; the out-of-bag prediction is exactly this form.

The last table is the topic’s bill. The inflation of the three combinations is 0.0516, 0.0635 and 0.0476; the smallest budget gives the highest number. Bagging reaches 0.7897 with 31 candidates, boosting stays at 0.7738 with 1,279. The topic’s total budget is 1,313 candidates and its return is a number 0.0318 above the baseline.

Summary

• The three base models give 0.7778, 0.7897 and 0.7738 on the test set; their costs are 43, 568 and 30 splits, that is, the cheapest combination is one nineteenth of the most expensive.
• On 185 of the 252 subscribers all three are right, on 45 all three are wrong; the region where a combination can work is 22 subscribers.
• The ceiling of a combination that looks at the three labels is 0.7897 and equals bagging’s number; bagging and boosting classify only 6 subscribers differently.
• When the second layer is learned from training it reports 0.8532 and gives 0.7778 (inflation 0.0754); when it is learned from validation it reports 0.8254 and again 0.7778 (inflation 0.0476).
• The wrong set shifts the weight to the model that memorizes most: +5.04 to the single tree against +1.99 in the correct setup; the smallest tuning budget of the three combinations gives the highest number.

Next Step

This topic tried seven structures in seven lessons: an example-based procedure, a probabilistic classifier built on a conditional independence assumption, a margin-based classifier, a decision tree, a bagged tree ensemble, sequential boosting and stacking. All were compared with the same table, the same split and the same baseline, and next to every number was written how many candidates were tried. But the comparison itself was made from start to finish with a single number: accuracy. What that number hides was never asked. Of the 45 subscribers on whom all three were wrong, how many were really suspicious and how many meant a team sent to a door for nothing was never counted. It is possible for two models that give the same accuracy to make entirely different errors, and the costs of the two kinds of error are not the same. The next topic, model evaluation, opens with this question.