Decision Trees, Forests, and Boosting

One access-change request breaks the obvious rule. R3 has an ambiguity score of only 2.0, yet its label says it needed human review. Several higher-ambiguity requests also needed review, but R8 at 3.2 didn't. How can a model catch exceptions without sending every request to a human?

The previous lesson fitted a smooth review score over ambiguity and policy support, with a straight decision boundary. A decision tree asks yes/no questions instead. Each answer leads to another question or a leaf, which stores a prediction. Fit one split, then investigate whether extra branches or multiple trees handle R3 more reliably.

Keep the same fictional error costs: $120 for a missed required review and $18 for an unnecessary review. These aren't measured incident costs. Authorization rules still apply independently of the classifier. A tree's feature threshold chooses a branch; a separate operating threshold turns its final score into a routing action.

View larger

The same eight requests, now as rectangles

Keep the same two evidence scores from the logistic example. ambiguity rises when request text and retrieved evidence conflict. auto_policy_support rises when policy P-7 backs automatic handling. needs_review is 1 when a human had to approve the permission change.

In this small synthetic dataset, each row is one request. Which single feature cut would leave the fewest mixed labels?

Four requests needed review and four didn't. Ambiguity looks promising because R4, R5, and R6 are high, but R3 and R8 are counterexamples. Keep eight separate constructed requests frozen for validation, using the same train / validation split as logistic regression. They represent a later batch in the story; their results aren't evidence about real traffic.

A decision tree splits on one feature at a time using an orthogonal, axis-aligned cut ($x_j \le \theta$). In our two-feature plane, each split draws a strictly horizontal or vertical line, partitioning the feature space into hyper-rectangles.

This geometry has a huge advantage and a subtle trap:

• Ordering rather than distance: An exact numeric split searches gaps between distinct values. Strictly increasing transformations preserve those candidate training partitions, so ordinary feature standardization isn't required. They don't necessarily preserve predictions inside unseen gaps. Ties, quantization, finite precision, and approximate split-search implementations also need checking.
• Diagonal boundary failure: If the true policy boundary is diagonal (say, $\text{ambiguity} + \text{support} > 6$), a single linear model draws it with one equation. An axis-aligned tree must approximate that diagonal using a jagged staircase of many orthogonal cuts, requiring deeper trees and fragmenting training data into tiny leaves.

Impurity scores how mixed a leaf is

Four of the eight requests needed review. Before choosing a cut, we need a mathematical score that measures how mixed or pure a leaf is: all automatic, all review, or at least less mixed than the parent.

For a leaf with $K$ classes where class $k$ has proportion $p_k$, Gini impurity measures the probability of misclassifying a randomly chosen sample if it were labeled randomly according to the leaf's class distribution:

For binary review routing with review proportion $p$, this simplifies to:

A pure leaf ($p = 0$ or $p = 1$) has $G = 0$. A leaf split evenly between classes ($p = 0.5$) has $G = 0.50$. Gini is a fast, bounded measure of label disagreement.

An alternative criterion from information theory is Shannon entropy^{[1]Reference 1A Mathematical Theory of Communicationhttps://people.math.harvard.edu/~ctm/home/text/others/shannon/entropy/entropy.pdf}:

For our binary case with base-two logarithms, entropy measures uncertainty in bits:

Treat a zero-probability term as zero using $\lim_{p\to0^+}p\log_2p=0$. An evenly mixed binary leaf has $H=1$ bit; a pure leaf has $H=0$. Both criteria prefer purer children, but their weighted rankings can differ. Entropy includes logarithms; actual fitting speed also depends on the implementation and data.

What if the target is continuous (such as handling duration in minutes) rather than binary? In regression trees, CART minimizes squared error through variance reduction:

where each leaf predicts the regional mean $\bar{y}$. Here $\text{Var}(S)=\frac{1}{|S|}\sum_{i\in S}(y_i-\bar y_S)^2$: the denominator is the number of rows, not one less. We're measuring mean training loss, not estimating an unknown population variance. A singleton leaf therefore has zero impurity. Minimizing weighted child variance is mathematically equivalent to maximizing the reduction in Mean Squared Error (MSE).^{[2]Reference 2Decision Trees: Missing Values Supporthttps://scikit-learn.org/stable/modules/tree.html}

For targets [10, 20, 60], the parent mean is 30 and its MSE is 466.67. Splitting into [10, 20] and [60] gives leaf means 15 and 60, with child impurities 25 and 0. Their weighted impurity is (2/3) × 25 = 16.67, so the split reduces MSE by 450. Using the unbiased sample-variance denominator would change this objective and leave the singleton undefined.

For classification, a split's impurity decrease is parent impurity minus weighted child impurity. When entropy is the metric, that decrease is information gain. When Gini is the metric, it's Gini gain. We'll use Gini gain because the arithmetic stays clean and fast.

At our training root, four of eight rows need review:

Now test the tempting cut ambiguity <= 3.35. Predict the outcome before checking the arithmetic: the three high-ambiguity rows (R4, R5, R6) are all reviews, so only the low side remains mixed.

• Left labels [0, 0, 1, 0, 0], so $p = 1/5 = 0.20$ and $G_{\text{left}} = 1 - 0.20^2 - 0.80^2 = 0.320$.
• Right labels [1, 1, 1], so $p = 3/3 = 1.00$ and $G_{\text{right}} = 0$.
• Weighted child impurity: $(5/8)(0.320) + (3/8)(0.000) = 0.200$.
• Gini gain: $0.500 - 0.200 = 0.300$.

The cut removes most label disagreement, but not all of it. R3 remains a required-review request inside the mostly automatic leaf.

from math import log2

TRAIN_X = [
    (1.0, 1.0),
    (1.5, 2.0),
    (2.0, 3.0),
    (3.5, 6.0),
    (4.0, 7.0),
    (4.5, 5.0),
    (2.8, 4.0),
    (3.2, 8.0),
]
TRAIN_Y = [0, 0, 1, 1, 1, 1, 0, 0]
VALID_X = [
    (1.2, 2.0),
    (2.1, 4.0),
    (2.4, 3.0),
    (3.0, 5.0),
    (3.6, 4.0),
    (3.2, 7.0),
    (4.2, 5.0),
    (2.8, 6.0),
]
VALID_Y = [0, 0, 1, 0, 1, 0, 1, 1]
NAMES = ["ambiguity", "auto_policy_support"]

def gini(labels):
    n = len(labels)
    if n == 0:
        return 0.0
    p = sum(labels) / n
    return 1.0 - p * p - (1.0 - p) * (1.0 - p)

def entropy(labels):
    n = len(labels)
    if n == 0:
        return 0.0
    p = sum(labels) / n
    terms = []
    for q in (p, 1.0 - p):
        if q > 0:
            terms.append(-q * log2(q))
    return sum(terms)

def midpoints(values):
    unique = sorted(set(values))
    return [(a + b) / 2 for a, b in zip(unique, unique[1:])]

root = TRAIN_Y
left = [y for x, y in zip(TRAIN_X, TRAIN_Y) if x[0] <= 3.35]
right = [y for x, y in zip(TRAIN_X, TRAIN_Y) if x[0] > 3.35]
weighted = (len(left) * gini(left) + len(right) * gini(right)) / len(root)

print(f"root gini={gini(root):.3f} entropy={entropy(root):.3f}")
print(f"left gini={gini(left):.3f} right gini={gini(right):.3f}")
print(f"weighted={weighted:.3f} gain={gini(root) - weighted:.3f}")

root gini=0.500 entropy=1.000
left gini=0.320 right gini=0.000
weighted=0.200 gain=0.300

Two criteria can prefer different cuts

Use a separate 24-row counterexample with twelve labels of each class. Candidate A isolates four negatives; candidate B produces two equally sized, moderately pure leaves. Which kind of purity does each criterion reward more?

parent = [0] * 12 + [1] * 12
partitions = {
    "A": ([0] * 4, [0] * 8 + [1] * 12),
    "B": ([0] * 3 + [1] * 9, [0] * 9 + [1] * 3),
}

for name, children in partitions.items():
    gains = []
    for criterion in (gini, entropy):
        weighted = sum(len(child) * criterion(child) for child in children) / len(parent)
        gains.append(criterion(parent) - weighted)
    print(f"{name}: Gini gain={gains[0]:.4f}; entropy gain={gains[1]:.4f}")

A: Gini gain=0.1000; entropy gain=0.1909
B: Gini gain=0.1250; entropy gain=0.1887

Gini prefers B; entropy prefers A. Neither result declares a better model on unseen requests. Return to our eight-row table and search its actual feature gaps.

A stump tests midpoints

A decision stump is a tree with one split. To find its rule, sort one feature, place a candidate between each pair of adjacent values, and score each partition by impurity decrease. The winner is the feature-threshold pair that leaves the children most uniform.

For ambiguity, 3.35 sits between 3.2 and 3.5, separating R4-R6 from the lower-ambiguity rows. With a <= rule, any threshold from 3.2 inclusive to 3.5 exclusive makes the same training partition. The midpoint is convenient; any value in this gap gives the same partition. Only gaps between distinct values matter because this feature cannot separate equal values.

def split_gain(values, labels, threshold):
    left = [y for x, y in zip(values, labels) if x <= threshold]
    right = [y for x, y in zip(values, labels) if x > threshold]
    weighted = (len(left) * gini(left) + len(right) * gini(right)) / len(labels)
    return gini(labels) - weighted

for column, name in enumerate(NAMES):
    values = [row[column] for row in TRAIN_X]
    scored = [(threshold, split_gain(values, TRAIN_Y, threshold)) for threshold in midpoints(values)]
    best_threshold, best_gain = max(scored, key=lambda item: item[1])
    print(f"{name:20} best threshold={best_threshold:.2f} gain={best_gain:.3f}")

ambiguity            best threshold=3.35 gain=0.300
auto_policy_support  best threshold=2.50 gain=0.167

The search selects ambiguity <= 3.35. Its left leaf reports the observed review rate, $1/5 = 0.20$; its right leaf reports $3/3 = 1.00$. The number in each leaf comes from labels that reached that leaf, not from a smooth curve.

View larger

A score of 0.20 isn't automatically "safe." The routing threshold still has to price a missed review against an extra queue item, as it did for logistic regression.

Fit the stump

Read the next cell as the stump's complete loop: scan every feature, try every midpoint, keep the largest gain, and store the two leaf rates. Then use that stored rule to score the training rows.

def fit_stump(rows, labels):
    root_impurity = gini(labels)
    best = None
    n_features = len(rows[0])
    for feature in range(n_features):
        values = [row[feature] for row in rows]
        for threshold in midpoints(values):
            left = [y for row, y in zip(rows, labels) if row[feature] <= threshold]
            right = [y for row, y in zip(rows, labels) if row[feature] > threshold]
            weighted = (len(left) * gini(left) + len(right) * gini(right)) / len(labels)
            gain = root_impurity - weighted
            if best is None or gain > best["gain"]:
                best = {
                    "feature": feature,
                    "threshold": threshold,
                    "gain": gain,
                    "left_probability": sum(left) / len(left),
                    "right_probability": sum(right) / len(right),
                }
    return best

def predict_proba(stump, rows):
    feature = stump["feature"]
    threshold = stump["threshold"]
    scores = []
    for row in rows:
        if row[feature] <= threshold:
            scores.append(stump["left_probability"])
        else:
            scores.append(stump["right_probability"])
    return scores

stump = fit_stump(TRAIN_X, TRAIN_Y)
scores = predict_proba(stump, TRAIN_X)
predictions = [int(score >= 0.50) for score in scores]
correct = sum(pred == label for pred, label in zip(predictions, TRAIN_Y))

print(
    f"rule: {NAMES[stump['feature']]} <= {stump['threshold']:.2f}; "
    f"gain={stump['gain']:.3f}"
)
print(
    f"leaf probabilities: left={stump['left_probability']:.2f}, "
    f"right={stump['right_probability']:.2f}"
)
print("training predictions:", predictions)
print(f"training accuracy={correct / len(TRAIN_Y):.3f}")

rule: ambiguity <= 3.35; gain=0.300
leaf probabilities: left=0.20, right=1.00
training predictions: [0, 0, 0, 1, 1, 1, 0, 0]
training accuracy=0.875

The returned scores make the tradeoff visible. At routing threshold 0.50, the stump misses R3. A deeper tree could repair this training row; evaluate later requests to measure whether the branch transfers. This minimal implementation assumes nonempty numeric data with at least one varying feature; a general estimator must also handle constant features, missing values, and invalid inputs.

The simple search rescans labels for every threshold. With $D$ features and up to $N-1$ cuts per feature, that costs roughly $O(DN^2)$ plus sorting. Efficient exact implementations sort and update cumulative class counts; histogram methods search a smaller set of bins. Don't extrapolate this eight-row trainer's performance to large datasets.

Multiplying ambiguity by ten preserves this stump's partitions and rescales its midpoint. A nonlinear transform gives a subtler result. Take two training values, 1 and 4. The original midpoint is 2.5; after squaring, it's 8.5. Which side receives an unseen value of 2.8?

training = [1.0, 4.0]
original_cut = sum(training) / 2
transformed_cut = sum(value**2 for value in training) / 2
unseen = 2.8

print("original training sides:", [value <= original_cut for value in training])
print("squared training sides:", [value**2 <= transformed_cut for value in training])
print("unseen 2.8 goes left originally:", unseen <= original_cut)
print("unseen 2.8 goes left after squaring:", unseen**2 <= transformed_cut)

original training sides: [True, False]
squared training sides: [True, False]
unseen 2.8 goes left originally: False
unseen 2.8 goes left after squaring: True

The training partition is identical; the unseen route changes. This counterexample concerns the chosen midpoint, not a new learned relationship. Scikit-learn also converts tree inputs to float32 internally, so extreme scaling can collapse distinct representable values.^{[2]Reference 2Decision Trees: Missing Values Supporthttps://scikit-learn.org/stable/modules/tree.html}

Held-out cost still picks the threshold

Freeze the fitted stump before looking at the later labels. These eight requests didn't choose its split and therefore give the first check on transfer:

Use the same operational costs as logistic regression:

A 0.50 threshold doesn't minimize cost by definition. Before running the cell, predict which threshold wins when a missed review costs almost seven times an unnecessary review. Compare candidates on validation, then confirm the selected policy on a fresh test period before shipping.

def metrics(y_true, y_pred):
    tp = sum(yt == 1 and yp == 1 for yt, yp in zip(y_true, y_pred))
    fp = sum(yt == 0 and yp == 1 for yt, yp in zip(y_true, y_pred))
    fn = sum(yt == 1 and yp == 0 for yt, yp in zip(y_true, y_pred))
    precision = tp / (tp + fp) if tp + fp else 0.0
    recall = tp / (tp + fn) if tp + fn else 0.0
    f1 = 2 * precision * recall / (precision + recall) if precision + recall else 0.0
    cost = 18 * fp + 120 * fn
    return precision, recall, f1, cost, fp, fn

probabilities = predict_proba(stump, VALID_X)
for threshold in (0.20, 0.50):
    predicted = [int(score >= threshold) for score in probabilities]
    precision, recall, f1, cost, fp, fn = metrics(VALID_Y, predicted)
    print(
        f"threshold={threshold:.2f} precision={precision:.3f} "
        f"recall={recall:.3f} f1={f1:.3f} fp={fp} fn={fn} cost=${cost}"
    )

threshold=0.20 precision=0.500 recall=1.000 f1=0.667 fp=4 fn=0 cost=$72
threshold=0.50 precision=1.000 recall=0.500 f1=0.667 fp=0 fn=2 cost=$240

Both thresholds tie on F1 here, but they make different mistakes. At 0.20, four false positives cost $72; at 0.50, two false negatives cost $240. The lower threshold buys recall with queue work, which is cheaper on this batch.

Extra depth can memorize R3

A one-split tree is easy to audit and can underfit. An unconstrained tree can isolate each of these distinct training inputs; conflicting labels at identical inputs couldn't be separated. Extra depth becomes overfitting when branches capture quirks that don't transfer. Before growing it, predict the two scores: training accuracy should rise, but V2 may inherit R3's exception leaf.

Grow the same Gini tree until every leaf is pure:

def fit_tree(rows, labels, max_depth, min_leaf=1, rng=None, mtry=None, depth=0):
    n = len(labels)
    p = sum(labels) / n if n else 0.0
    n_features = len(rows[0]) if rows else 0
    if depth >= max_depth or n < 2 * min_leaf or len(set(labels)) == 1 or n_features == 0:
        return {"leaf": True, "p": p, "n": n}

    features = list(range(n_features))
    if rng is not None and mtry is not None:
        rng.shuffle(features)
        features = features[:mtry]

    best = None
    for feature in features:
        values = [row[feature] for row in rows]
        for threshold in midpoints(values):
            left_idx = [i for i, row in enumerate(rows) if row[feature] <= threshold]
            left_set = set(left_idx)
            right_idx = [i for i in range(n) if i not in left_set]
            if len(left_idx) < min_leaf or len(right_idx) < min_leaf:
                continue
            left_y = [labels[i] for i in left_idx]
            right_y = [labels[i] for i in right_idx]
            gain = gini(labels) - (
                len(left_y) * gini(left_y) + len(right_y) * gini(right_y)
            ) / n
            if best is None or gain > best["gain"]:
                best = {
                    "feat": feature,
                    "threshold": threshold,
                    "gain": gain,
                    "left_idx": left_idx,
                    "right_idx": right_idx,
                }
    if best is None:
        return {"leaf": True, "p": p, "n": n}

    left_rows = [rows[i] for i in best["left_idx"]]
    right_rows = [rows[i] for i in best["right_idx"]]
    left_y = [labels[i] for i in best["left_idx"]]
    right_y = [labels[i] for i in best["right_idx"]]
    return {
        "leaf": False,
        "feat": best["feat"],
        "threshold": best["threshold"],
        "gain": best["gain"],
        "left": fit_tree(left_rows, left_y, max_depth, min_leaf, rng, mtry, depth + 1),
        "right": fit_tree(right_rows, right_y, max_depth, min_leaf, rng, mtry, depth + 1),
    }

def predict_one(tree, row):
    node = tree
    while not node["leaf"]:
        if row[node["feat"]] <= node["threshold"]:
            node = node["left"]
        else:
            node = node["right"]
    return node["p"]

def report(name, tree):
    train_scores = [predict_one(tree, row) for row in TRAIN_X]
    valid_scores = [predict_one(tree, row) for row in VALID_X]
    train_pred = [int(score >= 0.50) for score in train_scores]
    valid_pred = [int(score >= 0.50) for score in valid_scores]
    train_acc = sum(pred == y for pred, y in zip(train_pred, TRAIN_Y)) / len(TRAIN_Y)
    _, _, f1, cost, fp, fn = metrics(VALID_Y, valid_pred)
    print(
        f"{name}: train_accuracy={train_acc:.3f} val_f1={f1:.3f} "
        f"fp={fp} fn={fn} cost=${cost}"
    )
    print(f"{name} val scores: {[round(score, 2) for score in valid_scores]}")

stump_tree = fit_tree(TRAIN_X, TRAIN_Y, max_depth=1)
deep_tree = fit_tree(TRAIN_X, TRAIN_Y, max_depth=8)
report("stump", stump_tree)
report("deep tree", deep_tree)

stump: train_accuracy=0.875 val_f1=0.667 fp=0 fn=2 cost=$240
stump val scores: [0.2, 0.2, 0.2, 0.2, 1.0, 0.2, 1.0, 0.2]
deep tree: train_accuracy=1.000 val_f1=0.750 fp=1 fn=1 cost=$138
deep tree val scores: [0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0]

View larger

The extra splits isolate R3 between ambiguity 1.75 and 2.40. That fits one training exception; whether the interval captures a reusable rule requires more evidence.

V2 has ambiguity 2.1 and label 0, so it inherits the review leaf and incurs an $18 false-positive cost. V8 has ambiguity 2.8 and label 1, so it falls into the R7/R8 automatic leaf and incurs a $120 false-negative cost.

Training accuracy reaches 1.000. At the same operating threshold 0.50, validation cost actually falls from the stump's $240 to $138; perfect training fit doesn't automatically imply worse held-out performance. The deep tree still costs more than the stump's $72 at threshold 0.20. Greedy, hierarchical splits can be unstable under changes to training data. Tune depth, minimum leaf size, and pruning with representative validation or cross-validation, then verify the selected policy independently.^{[3]Reference 3The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/[4]Reference 4Scikit-learn: Machine Learning in Python.https://jmlr.org/papers/v12/pedregosa11a.html}

Forests average perturbed trees

Bagging (bootstrap aggregating) attacks that high variance by fitting many trees independently in parallel. Each tree receives a bootstrap sample: draw $N$ rows from our $N$-row table with replacement. Some requests repeat; others are left out entirely.

At a fixed input $x$, consider $B$ identically distributed tree scores across training-data and tree-building randomness. Assume each has variance $\sigma^2$ and every distinct pair has correlation $\rho$. Summing their variances and covariances gives:

As $B$ grows, the second term shrinks; the shared variation remains for this assumed building procedure. If $\rho=1$, averaging doesn't reduce variance. This floor isn't irreducible outcome noise or a formula for classification error. Bootstrap randomizations can be independent conditional on one fixed training set while the fitted scores remain correlated across shared random training sets.

Random forests perturb both rows and features to encourage useful disagreement between members.^{[5]Reference 5Random Forestshttps://doi.org/10.1023/A:1010933404324}

1. Bootstrap sampling: Each tree sees a different perturbed slice of data.
2. Feature subsampling: At each split, restrict the candidate features. This can reduce dependence on a dominant predictor, but may also weaken individual trees; tune the tradeoff. In scikit-learn 1.9.1, RandomForestClassifier defaults to max_features="sqrt", while RandomForestRegressor defaults to 1.0, considering all features. The older $D/3$ regression rule is a heuristic, not that library's current default.^{[6]Reference 6Ensembles: Gradient Boosting, Random Forests, Bagging, Voting, Stackinghttps://scikit-learn.org/stable/modules/ensemble.html}

Bootstrap sampling also provides a built-in validation mechanism: Out-of-Bag (OOB) validation. For a dataset of size $N$, the probability that a specific row isn't selected in a bootstrap draw of size $N$ is:

The 36.8% value is a large-$N$ limit; for our eight rows the expected fraction is about 34.4%. Score each row using only trees whose bootstrap omitted it. OOB evaluation can estimate same-distribution performance without a separate fitting holdout, but isn't automatically unbiased: each row uses fewer trees, and a small forest may leave some rows without any OOB prediction. Breiman explicitly discusses the finite-forest effect.^{[5]Reference 5Random Forestshttps://doi.org/10.1023/A:1010933404324} Ordinary row bootstrapping also doesn't isolate a future period or related requests. Tuning against OOB scores still requires independent verification of the selected policy.

Here we average leaf class-frequency scores, as scikit-learn's random-forest classifier does.^{[7]Reference 7RandomForestClassifierhttps://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html} Breiman's original classifier used hard majority voting. Averages preserve more score detail than hard votes; neither makes the probabilities calibrated automatically.

With two features, use one random feature per split (mtry = 1) to make that disagreement obvious. Predict R3's forest score before running the cell: if two trees never draw R3 and one tree over-samples it, the mean should land around 0.33, rather than becoming a hard majority vote of 0.

View larger

import random
from statistics import mean

rng = random.Random(0)
forest = []
bootstrap_members = []
for tree_index in range(3):
    sample_idx = [rng.randrange(len(TRAIN_X)) for _ in range(len(TRAIN_X))]
    sampled_rows = [TRAIN_X[i] for i in sample_idx]
    sampled_labels = [TRAIN_Y[i] for i in sample_idx]
    tree = fit_tree(sampled_rows, sampled_labels, max_depth=2, rng=rng, mtry=1)
    forest.append(tree)
    bootstrap_members.append(set(sample_idx))
    r3_score = predict_one(tree, TRAIN_X[2])
    feat = NAMES[tree["feat"]] if not tree["leaf"] else "leaf"
    threshold = tree.get("threshold", float("nan"))
    print(
        f"tree {tree_index + 1}: {feat} <= {threshold:.2f}; "
        f"bootstrap={sample_idx}; R3 p={r3_score:.2f}"
    )

r3_mean = mean(predict_one(tree, TRAIN_X[2]) for tree in forest)
valid_scores = [mean(predict_one(tree, row) for tree in forest) for row in VALID_X]
print(f"R3 forest mean={r3_mean:.2f}")
print("validation probabilities:", [round(score, 3) for score in valid_scores])
for threshold in (0.20, 0.50):
    predicted = [int(score >= threshold) for score in valid_scores]
    _, _, f1, cost, fp, fn = metrics(VALID_Y, predicted)
    print(f"threshold={threshold:.2f} f1={f1:.3f} fp={fp} fn={fn} cost=${cost}")

r3_oob = [
    predict_one(tree, TRAIN_X[2])
    for tree, members in zip(forest, bootstrap_members)
    if 2 not in members
]
print("R3 OOB trees:", len(r3_oob), "mean score:", round(mean(r3_oob), 2))
print("expected OOB fraction at N=8:", round((1 - 1 / 8)**8, 4))

tree 1: ambiguity <= 3.60; bootstrap=[6, 6, 0, 4, 7, 6, 4, 7]; R3 p=0.00
tree 2: auto_policy_support <= 2.50; bootstrap=[3, 2, 4, 2, 1, 4, 2, 4]; R3 p=1.00
tree 3: ambiguity <= 3.35; bootstrap=[1, 5, 7, 1, 5, 6, 5, 3]; R3 p=0.00
R3 forest mean=0.33
validation probabilities: [0.0, 0.333, 0.333, 0.333, 0.667, 0.333, 1.0, 0.333]
threshold=0.20 f1=0.727 fp=3 fn=0 cost=$54
threshold=0.50 f1=0.667 fp=0 fn=2 cost=$240
R3 OOB trees: 2 mean score: 0.0
expected OOB fraction at N=8: 0.3436

Tree 1 never draws R3, so its low-ambiguity leaf is pure automatic and scores R3 as 0.00. Tree 2 over-samples R3 and splits on policy support, scoring it as 1.00. Tree 3 also drops R3.

The mean 0.33 is the forest's review score for that exception: two trees never saw it, and one tree treated it as a review case. It's an average of tree scores, not a new label.

Our toy considers exactly the sampled features at each node and stops if none yields a valid partition. Library split search can inspect additional features to find a valid partition. This demonstrates row/feature perturbation and averaging; it isn't a full clone of RandomForestClassifier.

Validation check: On validation, threshold 0.20 costs $54. V1's score is 0.00, a true negative, while the other low-ambiguity rows go to review. The threshold improves cost on this eight-row slice. Eight rows and three trees demonstrate the mechanism; a production comparison needs representative data and broader evaluation.

R3's OOB score is 0.00, using only trees 1 and 3; tree 2 must be excluded because its bootstrap contains R3. The full-forest mean 0.33 isn't an OOB score. With eight rows and three trees, neither supports a release estimate.

More trees stabilize the average of this tree-building procedure; they don't repair a bad feature set or a biased dataset. Individual forest trees often fit training data closely, and their average can still generalize well. Select depth and leaf size using validation rather than assuming that deeper members necessarily worsen the ensemble.

The forest emits more varied scores than the stump. Its threshold still belongs to the routing policy and must be evaluated with the same costs. Averaging trees that don't see one another's errors leaves a different question: what if the next tree was trained to fix the mistakes the current score still makes?

Boosting fits the leftover error

A forest averages trees that were fit separately. The next tree can instead be trained on the mistakes the current score still makes. Gradient boosting does that in order: each new tree targets the leftover error of the ensemble so far.

In standard gradient descent, we update model parameter weights $w$ in the direction of negative parameter gradients: $w \leftarrow w - \eta \nabla_w L$. In gradient boosting, we perform gradient descent in function space^{[8]Reference 8Greedy Function Approximation: A Gradient Boosting Machinehttps://projecteuclid.org/journals/annals-of-statistics/volume-29/issue-5/Greedy-Function-Approximation--A-Gradient-Boosting-Machine/10.1214/aos/1013203451.full}. Our "parameters" are the predicted outputs $F(x_i)$ at each training example.

To minimize an arbitrary differentiable loss function $L(y, F(x))$, the steepest descent step is the negative gradient of the loss with respect to the current predictions:

We call $r_{im}$ a pseudo-residual.

For squared-error regression $L(y, F) = \frac{1}{2}(y - F)^2$, the negative gradient is the ordinary difference between the label and the current prediction:

The target for each new tree is that residual, $y - \hat{y}$.

Before returning to binary review routing, predict one correction step on continuous handling time. For ambiguity values [1, 2, 5] and handling times [10, 20, 60], the initial prediction is their mean, 30. The first two rows need downward corrections, while the last needs a large upward correction.

from statistics import mean

ambiguity = [1.0, 2.0, 5.0]
minutes = [10.0, 20.0, 60.0]
prediction = [mean(minutes)] * len(minutes)
residual = [target - current for target, current in zip(minutes, prediction)]

threshold = 3.5
left = [value for value, flag in zip(residual, ambiguity) if flag <= threshold]
right = [value for value, flag in zip(residual, ambiguity) if flag > threshold]
correction = [mean(left) if value <= threshold else mean(right) for value in ambiguity]
learning_rate = 0.10
updated = [current + learning_rate * delta for current, delta in zip(prediction, correction)]

print("base prediction:", [round(value, 1) for value in prediction])
print("residual targets:", [round(value, 1) for value in residual])
print("stump correction:", [round(value, 1) for value in correction])
print("updated prediction:", [round(value, 1) for value in updated])

base prediction: [30.0, 30.0, 30.0]
residual targets: [-20.0, -10.0, 30.0]
stump correction: [-15.0, -15.0, 30.0]
updated prediction: [28.5, 28.5, 33.0]

View larger

Why scale the update by learning rate (also known as shrinkage) $\eta = 0.10$?

At $\eta=1$, this step gives [15, 15, 60]: it fits the last row exactly but still misses the first two. Shrinkage scales each added correction and changes the learning path; it isn't the explicit $L_2$ leaf penalty introduced below. Smaller rates often need more rounds, and don't guarantee diverse trees or prevent noise fitting. Tune the rate together with tree complexity and the stopping round on appropriate validation data.^{[6]Reference 6Ensembles: Gradient Boosting, Random Forests, Bagging, Voting, Stackinghttps://scikit-learn.org/stable/modules/ensemble.html}

The next cell repeats that update four times. Watch whether each residual stump keeps the same threshold and whether training MSE falls:

from statistics import mean

def residual_stump(values, residual):
    unique = sorted(set(values))
    candidates = [(a + b) / 2 for a, b in zip(unique, unique[1:])]
    best = None
    for threshold in candidates:
        left = [res for value, res in zip(values, residual) if value <= threshold]
        right = [res for value, res in zip(values, residual) if value > threshold]
        prediction = [mean(left) if value <= threshold else mean(right) for value in values]
        mse = mean((res - pred) ** 2 for res, pred in zip(residual, prediction))
        if best is None or mse < best["error"]:
            best = {"threshold": threshold, "prediction": prediction, "error": mse}
    return best

x = [1.0, 2.0, 5.0]
y = [10.0, 20.0, 60.0]
prediction = [mean(y)] * len(y)
learning_rate = 0.10

print(
    f"round=0 predictions={[round(value, 2) for value in prediction]} "
    f"mse={mean((target - pred) ** 2 for target, pred in zip(y, prediction)):.2f}"
)
for round_number in range(1, 5):
    residual = [target - pred for target, pred in zip(y, prediction)]
    stump = residual_stump(x, residual)
    prediction = [
        current + learning_rate * delta
        for current, delta in zip(prediction, stump["prediction"])
    ]
    mse = mean((target - pred) ** 2 for target, pred in zip(y, prediction))
    print(
        f"round={round_number} threshold={stump['threshold']:.2f} "
        f"predictions={[round(value, 2) for value in prediction]} mse={mse:.2f}"
    )

round=0 predictions=[30.0, 30.0, 30.0] mse=466.67
round=1 threshold=3.50 predictions=[28.5, 28.5, 33.0] mse=381.17
round=2 threshold=3.50 predictions=[27.15, 27.15, 35.7] mse=311.91
round=3 threshold=3.50 predictions=[25.93, 25.93, 38.13] mse=255.82
round=4 threshold=3.50 predictions=[24.84, 24.84, 40.32] mse=210.38

Classification boosting uses y minus p

For access routing, we need probabilities rather than handling minutes. Let $F$ be the current additive logit score, and let $p = \sigma(F) = \frac{1}{1 + e^{-F}}$ be the review probability.

For a single request with true label $y \in \{0, 1\}$, the binary cross-entropy loss is:

Differentiating with respect to the logit $F$ using the chain rule ($\frac{d\sigma}{dF} = p(1-p)$):

The pseudo-residual (the negative gradient) is therefore:

When $y = 1$ and the current probability is low (like $p = 0.20$), the residual is positive ($+0.80$), asking for a higher logit. When $y = 0$ and $p = 0.80$, the residual is negative ($-0.80$), asking for a lower logit.

For an approximate Newton leaf update, also compute the second derivative:

Our gradient-boosting fixture fits split regions by squared error on $y-p$, then takes one Newton step in each fixed leaf:^{[8]Reference 8Greedy Function Approximation: A Gradient Boosting Machinehttps://projecteuclid.org/journals/annals-of-statistics/volume-29/issue-5/Greedy-Function-Approximation--A-Gradient-Boosting-Machine/10.1214/aos/1013203451.full}

Start at $F=0$, the best constant unpenalized logit for these balanced training labels. In general, a constant baseline uses the training positive fraction $\pi$ through $\log(\pi/(1-\pi))$ when $0<\pi<1$. Here $p=0.5$ and every curvature term is $0.25$; the residual stump splits at ambiguity <= 3.35.

• Left leaf ($5$ rows: R1, R2, R3, R7, R8): One review ($y=1$) and four automatic ($y=0$). The residual sum is $1 \times (1 - 0.5) + 4 \times (0 - 0.5) = 0.5 - 2.0 = -1.5$. $$\gamma_{\text{left}} = \frac{-1.5}{5 \times 0.25} = -1.20$$
• Right leaf ($3$ rows: R4, R5, R6): Three reviews ($y=1$). The residual sum is $3 \times (1 - 0.5) = +1.5$. $$\gamma_{\text{right}} = \frac{1.5}{3 \times 0.25} = +2.00$$

With learning rate $\eta = 0.10$, the updated logits become $0 + 0.10(-1.20) = -0.12$ on the left and $0 + 0.10(2.00) = +0.20$ on the right. The corresponding probabilities shift to $\sigma(-0.12) \approx 0.470$ and $\sigma(0.20) \approx 0.550$. Notice the contrast with an isolated decision stump: boosting doesn't overwrite predictions with extreme 0.20 and 1.00 probabilities. It nudges the ensemble in the direction of the gradient.

The Newton values aren't exact minimizers of the original log loss. Starting at zero, the left leaf's exact constant optimum would be $\operatorname{logit}(0.20)\approx-1.386$, not $-1.20$. The all-positive right leaf has no finite unpenalized optimum, despite its finite first Newton step of $+2$. A quadratic approximation, a step-size choice, and an exact optimum are different objects.

XGBoost integrates gradient and curvature statistics directly into its split score, rather than selecting the tree solely by residual squared error.^{[9]Reference 9XGBoost: A Scalable Tree Boosting System.https://www.kdd.org/kdd2016/papers/files/rfp0697-chenAemb.pdf[10]Reference 10Introduction to Boosted Treeshttps://xgboost.readthedocs.io/en/stable/tutorials/model.html} For this binary log-loss objective, its second-order approximation is:

where $g_i = p_i - y_i$ is the gradient, $h_i = p_i(1 - p_i)$ is the Hessian, $T$ is the number of leaves, and $\lambda$ is an $L_2$ leaf regularizer.

Minimizing that quadratic approximation for a fixed tree gives:

With $\lambda=0$ and positive summed curvature, this is the same Newton step. Positive $\lambda$ shrinks its magnitude for a fixed numerator and curvature. Check calibration and extreme predictions separately. Define $G_L,G_R$ as summed gradients and $H_L,H_R$ as summed Hessians. The proposed split's improvement, including the extra-leaf charge, is:

For our first stump with $\lambda=0$, the improvement before the charge is 2.4; a charge greater than 2.4 makes this proposal unfavorable. Once the charge has been subtracted, compare gain with zero, not with the charge again. Actual growth also obeys estimator constraints such as depth and minimum child weight.^{[10]Reference 10Introduction to Boosted Treeshttps://xgboost.readthedocs.io/en/stable/tutorials/model.html}

Predict the validation consequence before running the four-round boosting code. The probability loss should fall, but if V3 and V8 stay on the low-ambiguity side, threshold 0.50 will still miss them.

from math import exp, log1p
from statistics import mean

def sigmoid(logit):
    if logit >= 0:
        return 1.0 / (1.0 + exp(-logit))
    exponential = exp(logit)
    return exponential / (1.0 + exponential)

def log_loss_from_logits(labels, logits):
    terms = []
    for label, logit in zip(labels, logits):
        signed_error = (1 - 2 * label) * logit
        terms.append(max(signed_error, 0.0) + log1p(exp(-abs(signed_error))))
    return mean(terms)

def fit_residual_stump(rows, residual):
    best = None
    n_features = len(rows[0])
    for feature in range(n_features):
        values = [row[feature] for row in rows]
        for threshold in midpoints(values):
            left_idx = [i for i, row in enumerate(rows) if row[feature] <= threshold]
            left_set = set(left_idx)
            right_idx = [i for i in range(len(rows)) if i not in left_set]
            left_mean = mean(residual[i] for i in left_idx)
            right_mean = mean(residual[i] for i in right_idx)
            prediction = [
                left_mean if row[feature] <= threshold else right_mean for row in rows
            ]
            mse = mean((res - pred) ** 2 for res, pred in zip(residual, prediction))
            if best is None or mse < best["mse"]:
                best = {
                    "feature": feature,
                    "threshold": threshold,
                    "left_idx": left_idx,
                    "right_idx": right_idx,
                    "mse": mse,
                }
    return best

def newton_value(idxs, residual, hessian, leaf_l2=0.0):
    if leaf_l2 < 0:
        raise ValueError("leaf_l2 must be nonnegative")
    numerator = sum(residual[i] for i in idxs)
    denominator = sum(hessian[i] for i in idxs) + leaf_l2
    if denominator <= 0:
        raise ValueError("nonpositive leaf curvature")
    return numerator / denominator

logits = [0.0] * len(TRAIN_Y)
stumps = []
learning_rate = 0.10

def boosted_logit(row):
    return sum(
        learning_rate * (
            fitted["left"] if row[fitted["feature"]] <= fitted["threshold"] else fitted["right"]
        )
        for fitted in stumps
    )

def boosted_probability(row):
    return sigmoid(boosted_logit(row))

print(f"round=0 val_logloss={log_loss_from_logits(VALID_Y, [0.0] * len(VALID_Y)):.3f}")
for round_number in range(1, 5):
    probabilities = [sigmoid(logit) for logit in logits]
    residual = [label - probability for label, probability in zip(TRAIN_Y, probabilities)]
    hessian = [probability * (1.0 - probability) for probability in probabilities]
    candidate = fit_residual_stump(TRAIN_X, residual)
    stump = {
        "feature": candidate["feature"],
        "threshold": candidate["threshold"],
        "left": newton_value(candidate["left_idx"], residual, hessian),
        "right": newton_value(candidate["right_idx"], residual, hessian),
    }
    stumps.append(stump)
    logits = [
        logit
        + learning_rate
        * (stump["left"] if row[stump["feature"]] <= stump["threshold"] else stump["right"])
        for logit, row in zip(logits, TRAIN_X)
    ]

    valid_logits = [boosted_logit(row) for row in VALID_X]
    print(
        f"round={round_number} {NAMES[stump['feature']]} <= {stump['threshold']:.2f} "
        f"left={stump['left']:.3f} right={stump['right']:.3f} "
        f"val_logloss={log_loss_from_logits(VALID_Y, valid_logits):.3f}"
    )

valid_scores = [boosted_probability(row) for row in VALID_X]
print("validation probabilities:", [round(score, 3) for score in valid_scores])
for threshold in (0.20, 0.50):
    predicted = [int(score >= threshold) for score in valid_scores]
    _, _, f1, cost, fp, fn = metrics(VALID_Y, predicted)
    print(f"threshold={threshold:.2f} f1={f1:.3f} fp={fp} fn={fn} cost=${cost}")

print("stable wrong-logit loss:", round(log_loss_from_logits([0], [1000.0]), 1))
try:
    newton_value([0], [-1.0], [0.0])
except ValueError as error:
    print("unpenalized saturated leaf:", error)
print("same leaf with L2=1:", newton_value([0], [-1.0], [0.0], leaf_l2=1.0))

round=0 val_logloss=0.693
round=1 ambiguity <= 3.35 left=-1.200 right=2.000 val_logloss=0.656
round=2 ambiguity <= 3.35 left=-1.084 right=1.819 val_logloss=0.626
round=3 ambiguity <= 3.35 left=-0.985 right=1.683 val_logloss=0.603
round=4 ambiguity <= 3.35 left=-0.900 right=1.577 val_logloss=0.584
validation probabilities: [0.397, 0.397, 0.397, 0.397, 0.67, 0.397, 0.67, 0.397]
threshold=0.20 f1=0.667 fp=4 fn=0 cost=$72
threshold=0.50 f1=0.667 fp=0 fn=2 cost=$240
stable wrong-logit loss: 1000.0
unpenalized saturated leaf: nonpositive leaf curvature
same leaf with L2=1: -1.0

Validation log loss falls from 0.693 to 0.584; all four rounds retain the original split. Cost at 0.50 stays $240 because V3 and V8 remain on the low-ambiguity side. A better probability loss needn't improve a selected policy. Choose rounds and thresholds on validation, freeze both, then evaluate a locked test period.

The stress check computes loss from finite logits, preserving a confidently wrong loss of 1000 even when sigmoid rounds to 1. Rounded probabilities can also make curvature zero; the unpenalized helper reports that condition instead of silently returning a zero update. An L2 penalty can make this leaf's quadratic problem well defined, but changes the fitted objective. This small trainer handles valid 0/1 labels, finite features, and a varying feature. A production estimator also needs controls for invalid inputs and broader data conditions.

The default unweighted Gini and log-loss criteria don't price a false negative at $120 while fitting. When false negatives cost far more than false positives, per-row sample weights can make high-harm rows influence splits. Sweep the decision threshold on validation anyway: weights change the fit, while the threshold turns scores into actions.

Leaf frequencies and forest averages are ranking scores, not automatically honest probabilities. The logistic chapter measured expected calibration error on these same eight validation rows. Treat ensemble predict_proba as a score until a reliability audit on a reserved split says otherwise.

A tree prediction is constant outside the training support of each split. If production ambiguity climbs past every training value, the leaf still returns the training rate for that region. Monitor feature ranges and don't expect the smooth extrapolation a linear model might give.

What bias and variance can tell us

The familiar decomposition is for squared prediction error, not every classification metric. At a fixed input $x$, let $Y=f_*(x)+\epsilon$, with zero-mean noise independent of the training data and variance $\sigma_\epsilon^2$. Across random training sets and fitting seeds $D$:

Perfect training accuracy doesn't measure statistical bias. A singleton tree can fit all eight labels while producing systematically wrong scores elsewhere. A stump can also be unstable if small data changes switch its best feature or threshold.

Averaging identically distributed member predictions preserves their common expectation and can reduce variance through imperfect correlation. Changing bootstrap, feature, depth, or leaf-size rules changes the member distribution too, so a forest isn't guaranteed to have the same bias as an arbitrary deep tree.

Boosting expands an additive function to improve a chosen loss. It can reduce approximation error, but its fitted variance can also change. More rounds can overfit; eventual memorization isn't guaranteed, especially with restricted learners or contradictory labels at identical inputs. Depth, learning rate, sampling, regularization, and stopping interact. For binary review routing, compare the resulting validation costs and calibration directly.^{[6]Reference 6Ensembles: Gradient Boosting, Random Forests, Bagging, Voting, Stackinghttps://scikit-learn.org/stable/modules/ensemble.html}

Those controls change the fit. They don't say which feature the forest will claim mattered, or whether that claim survives on later rows.

MDI can crown a useless ID

Tree libraries often report mean decrease in impurity (MDI) as default feature importance. Each split contributes its impurity decrease weighted by the fraction of training samples reaching that node. Contributions are grouped by feature, normalized, and averaged across trees. It's fast, but high-cardinality features offer many candidate thresholds and more chances to fit accidental patterns. A numeric request ID can receive large importance without predicting later labels.^{[11]Reference 11Permutation Feature Importancehttps://scikit-learn.org/stable/modules/permutation_importance.html}

To isolate that failure, use a separate synthetic experiment. A binary feature agrees with the label 70% of the time; a continuous feature is independent uniform noise. Both training and validation come from the same process. Scikit-learn supplies the forest and importance calculations here, so we don't need a second implementation of the tree builder. Forty training rows leave plenty of opportunity to overfit; 2,000 held-out rows give a less noisy comparison.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.metrics import f1_score

rng = np.random.default_rng(1)

def make_rows(n):
    labels = rng.integers(0, 2, size=n)
    binary = np.where(rng.random(n) < 0.70, labels, 1 - labels)
    noise = rng.random(n)
    return np.column_stack([binary, noise]), labels

train_x, train_y = make_rows(40)
valid_x, valid_y = make_rows(2000)
model = RandomForestClassifier(
    n_estimators=20, max_depth=5, max_features=1, random_state=0
).fit(train_x, train_y)
importance = permutation_importance(
    model, valid_x, valid_y, scoring="f1", n_repeats=15, random_state=0
)

print(f"held-out F1={f1_score(valid_y, model.predict(valid_x)):.3f}")
for index, name in enumerate(["weak binary", "continuous noise"]):
    print(
        f"{name}: MDI={model.feature_importances_[index]:.3f}; "
        f"permutation F1 drop={importance.importances_mean[index]:.3f} "
        f"(shuffle SD={importance.importances_std[index]:.3f})"
    )

held-out F1=0.575
weak binary: MDI=0.080; permutation F1 drop=0.097 (shuffle SD=0.007)
continuous noise: MDI=0.920; permutation F1 drop=0.003 (shuffle SD=0.007)

Permutation importance asks how this fitted model's selected metric changes when one feature is shuffled across held-out rows. Marginal shuffling breaks that feature's relationship with labels and other features, potentially creating unusual combinations. It doesn't measure intrinsic feature quality or establish a causal effect.^{[11]Reference 11Permutation Feature Importancehttps://scikit-learn.org/stable/modules/permutation_importance.html}

In this seeded experiment, MDI assigns 0.920 to continuous noise, while its held-out permutation F1 drop is only 0.003. The binary predictor's drop is 0.097. This exposes a misleading training-importance ranking. The reported 0.007 shuffle SD describes repeated permutations on the fixed model and data; it isn't a confidence interval over new datasets, refits, or production traffic.

With strongly correlated features, another feature may retain redundant information, making individual permutation effects small.^{[12]Reference 12Permutation Importance with Multicollinear or Correlated Featureshttps://scikit-learn.org/stable/auto_examples/inspection/plot_permutation_importance_multicollinear.html} A fitted tree that relies entirely on one copy can still suffer when that copy is shuffled. Consider a joint permutation of a feature group or a justified conditional method, and state which question it answers.

An importance score doesn't certify model quality. For access routing, include permutation effects on the declared cost at the frozen threshold, alongside recall and queue capacity. F1 alone can't express the two error prices.

A local explanation isn't a verdict

Global importances tell you which features matter across an entire dataset. They don't explain why the model made a specific decision on request R3.

SHAP (SHapley Additive exPlanations) solves this by adapting Shapley values from cooperative game theory to machine learning models^{[13]Reference 13A Unified Approach to Interpreting Model Predictionshttps://arxiv.org/abs/1705.07874}. It decomposes an individual prediction $f(x)$ into an additive sum of a baseline expectation plus individual feature contributions:

For our chosen background game, $\phi_0=\mathbb E[f(X)]$. The original SHAP paper identifies three properties: local accuracy, missingness in the simplified feature-presence representation, and consistency. Consistency requires the feature's marginal contribution not to decrease for any coalition when the model changes; then its attribution must not decrease. The reconstruction above expresses local accuracy. These properties concern an attribution game; they don't validate labels or establish causal effects.^{[13]Reference 13A Unified Approach to Interpreting Model Predictionshttps://arxiv.org/abs/1705.07874}

The baseline and contributions depend on the background data and how absent features are replaced. Here we explain the stump's review probability using an interventional background: hold selected features at the request's values and fill the others from each training row, without conditioning on the held features.^{[14]Reference 14shap.TreeExplainerhttps://shap.readthedocs.io/en/latest/generated/shap.TreeExplainer.html}

With two features, a Shapley contribution averages the change from adding that feature first and adding it second. Our stump only uses ambiguity, so revealing policy support never changes its output under this replacement rule. That gives a clean, exact check:

• Baseline review score across the training background: $0.50$.
• Left-leaf score: $0.20$, so ambiguity contributes $-0.30$.
• Right-leaf score: $1.00$, so ambiguity contributes $+0.50$.
• auto_policy_support contributes 0.00 under the stated interventional calculation.

from statistics import mean

attribution_stump = fit_stump(TRAIN_X, TRAIN_Y)

def background_score(row, known_features):
    completed = [
        [row[j] if j in known_features else background[j] for j in range(2)]
        for background in TRAIN_X
    ]
    return mean(predict_proba(attribution_stump, completed))

for name, row in [("R5 high ambiguity", TRAIN_X[4]), ("R3 missed exception", TRAIN_X[2])]:
    base_score = background_score(row, set())
    only_ambiguity = background_score(row, {0})
    only_policy = background_score(row, {1})
    prediction = background_score(row, {0, 1})
    ambiguity_contribution = 0.5 * (
        only_ambiguity - base_score + prediction - only_policy
    )
    policy_contribution = 0.5 * (
        only_policy - base_score + prediction - only_ambiguity
    )
    rebuilt = base_score + ambiguity_contribution + policy_contribution
    assert abs(rebuilt - prediction) < 1e-12
    print(
        f"{name}: base={base_score:.2f} ambiguity={ambiguity_contribution:+.2f} "
        f"policy={policy_contribution:+.2f} prediction={rebuilt:.2f}"
    )

R5 high ambiguity: base=0.50 ambiguity=+0.50 policy=+0.00 prediction=1.00
R3 missed exception: base=0.50 ambiguity=-0.30 policy=+0.00 prediction=0.20

R3 exposes the divide between explanation and decision correctness. The attribution faithfully reports why the stump scored the request low: its ambiguity of 2.0 placed it in the low-ambiguity leaf, contributing $-0.30$. But that doesn't prove automatic processing was safe! R3's true label is $y = 1$ (it required human review). Local explanations audit what a model computed; they don't validate real-world ground truth.

If instead the explainer conditions on correlated features, revealing policy support can convey information about ambiguity even though the tree never splits on support. Attribution then answers a different question. Always document whether your SHAP values use interventional or conditional expectations rather than assuming one attribution is uniquely defined.

With TreeExplainer, also record model_output and feature_perturbation. Its default raw output depends on the estimator; for binary XGBoost, raw attributions reconstruct log odds, not a percentage. Current auto chooses interventional replacement when a background is supplied and tree-path-dependent replacement otherwise. Tree-path-dependent replacement isn't a general conditional-expectation calculation.^{[14]Reference 14shap.TreeExplainerhttps://shap.readthedocs.io/en/latest/generated/shap.TreeExplainer.html}

Missing values are a library contract

"Trees handle missing data automatically" is a half-truth that breaks pipelines. Before moving an estimator, check which component owns the missing-value route. The behavior depends on library, estimator, and version:

The following contracts were checked in scikit-learn 1.9.1, XGBoost 3.4.2, and LightGBM 4.7.0 documentation on September 21, 2026. Our hand-written trainers don't implement native missing-value routing.

The trap is assuming the behavior is uniform. Moving a pipeline between estimators or versions can turn native routing into an input error.

The quieter trap changes data without raising one. zero_as_missing and sparse-matrix conventions mean a literal 0 and a missing value can be treated identically or differently depending on the tool, which changes splits. Decide the missing-value contract explicitly: impute, encode a missingness flag, or rely on documented native behavior. A tree won't infer that contract for you.

Pick a model the same way you picked logistic

Use the same later-request contract for every candidate. Compare a regularized tree or ensemble with logistic regression under identical splits and costs, then inspect whether added flexibility repairs meaningful failures. Read the table below as a set of decisions, not a leaderboard: the cheapest row still rests on eight validation requests.

The forest improves on our stump, but it doesn't beat the logistic baseline on this batch. Greater flexibility isn't a reason to replace a working model. We used these rows to inspect and select policies, so even the logistic result needs a fresh verification set before making a performance claim.

Trees expose conditional splits. Forests reduce single-tree instability by averaging perturbed members. Boosting accumulates targeted corrections. Each mechanism changes a different failure mode, and none replaces held-out evaluation.

Compute the gain of a cut, and watch R3 either stay in a mixed leaf or disappear into a one-row leaf. Compare those scores with the logistic baseline on the same later requests and the same dollar costs.

Change the example, then check your reasoning