S090 · Random Forest & Bagging

# random_forest.py — RF from scratch + sklearn comparison, ~140 lines.
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
from sklearn.inspection import permutation_importance
 
np.random.seed(42)
 
# --- Data ---
data = load_breast_cancer()
X, y = data.data, data.target
feature_names = list(data.feature_names)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, stratify=y, random_state=42)
 
# =====================================================
# 1. Single-tree baseline
# =====================================================
tree = DecisionTreeClassifier(random_state=42).fit(X_train, y_train)
print(f"[single tree]  train={accuracy_score(y_train, tree.predict(X_train)):.4f}  "
      f"test={accuracy_score(y_test, tree.predict(X_test)):.4f}")
 
# =====================================================
# 2. Random Forest from scratch
# =====================================================
class RFScratch:
    def __init__(self, n_estimators=100, max_features="sqrt", max_depth=None, random_state=None):
        self.n_estimators = n_estimators
        self.max_features = max_features
        self.max_depth = max_depth
        self.random_state = random_state
        self.trees = []
        self.feature_subsets = []       # which features each tree was trained on
        self.oob_masks = []             # which samples were OOB for each tree
 
    def _feature_count(self, d):
        if self.max_features == "sqrt":
            return max(1, int(np.sqrt(d)))
        if isinstance(self.max_features, int):
            return self.max_features
        return d                        # "all"
 
    def fit(self, X, y):
        rng = np.random.default_rng(self.random_state)
        n, d = X.shape
        n_feat = self._feature_count(d)
        for b in range(self.n_estimators):
            # Bootstrap sample
            idx = rng.integers(0, n, size=n)
            oob = np.setdiff1d(np.arange(n), np.unique(idx))
            # Random feature subset PER TREE (simpler than per-split for scratch)
            feats = rng.choice(d, size=n_feat, replace=False)
            tree = DecisionTreeClassifier(max_depth=self.max_depth, random_state=b)
            tree.fit(X[idx][:, feats], y[idx])
            self.trees.append(tree)
            self.feature_subsets.append(feats)
            self.oob_masks.append(oob)
        return self
 
    def predict_proba(self, X):
        n_classes = 2      # binary here for simplicity
        probs = np.zeros((X.shape[0], n_classes))
        for tree, feats in zip(self.trees, self.feature_subsets):
            probs += tree.predict_proba(X[:, feats])
        return probs / len(self.trees)
 
    def predict(self, X):
        return self.predict_proba(X).argmax(axis=1)
 
    def oob_score(self, X, y):
        """OOB score: for each sample, average predictions from trees that DIDN'T see it."""
        n = X.shape[0]
        n_classes = 2
        oob_probs = np.zeros((n, n_classes))
        oob_counts = np.zeros(n)
        for tree, feats, oob in zip(self.trees, self.feature_subsets, self.oob_masks):
            if len(oob) == 0:
                continue
            preds = tree.predict_proba(X[oob][:, feats])
            oob_probs[oob] += preds
            oob_counts[oob] += 1
        mask = oob_counts > 0
        oob_probs[mask] /= oob_counts[mask, None]
        return accuracy_score(y[mask], oob_probs[mask].argmax(axis=1))
 
# --- Fit + evaluate ---
rf = RFScratch(n_estimators=200, max_features="sqrt", max_depth=None, random_state=42)
rf.fit(X_train, y_train)
print(f"\n[scratch RF]   train={accuracy_score(y_train, rf.predict(X_train)):.4f}  "
      f"test={accuracy_score(y_test, rf.predict(X_test)):.4f}  "
      f"OOB={rf.oob_score(X_train, y_train):.4f}")
 
# =====================================================
# 3. sklearn comparison
# =====================================================
sk = RandomForestClassifier(n_estimators=200, max_features="sqrt", oob_score=True, random_state=42)
sk.fit(X_train, y_train)
print(f"[sklearn RF]   train={accuracy_score(y_train, sk.predict(X_train)):.4f}  "
      f"test={accuracy_score(y_test, sk.predict(X_test)):.4f}  "
      f"OOB={sk.oob_score_:.4f}")
 
# =====================================================
# 4. Show how ensemble accuracy grows with N
# =====================================================
print("\n=== Test accuracy vs number of trees ===")
for n in [1, 5, 10, 25, 50, 100, 200]:
    m = RandomForestClassifier(n_estimators=n, max_features="sqrt", random_state=42).fit(X_train, y_train)
    print(f"  n_estimators={n:>4}: test = {accuracy_score(y_test, m.predict(X_test)):.4f}")
 
# =====================================================
# 5. Feature importance — DEFAULT (biased) vs PERMUTATION (honest)
# =====================================================
print("\n=== Feature importance ===")
mdi_imp = sk.feature_importances_
perm_imp = permutation_importance(sk, X_test, y_test, n_repeats=20, random_state=42, n_jobs=-1).importances_mean
 
# Sort by permutation importance
order = np.argsort(perm_imp)[::-1][:10]
print(f"  {'feature':<28} {'MDI':>8} {'permutation':>13}")
for i in order:
    print(f"  {feature_names[i][:27]:<28} {mdi_imp[i]:>8.4f} {perm_imp[i]:>13.4f}")