""" cart.py -- Classification And Regression Trees (CART) from scratch (Breiman, Friedman, Olshen & Stone 1984). Backs the notebooks in Z-ML Trees-Decision Trees (CART). A decision tree recursively splits the feature space with axis-aligned cuts, each chosen to make the two child nodes as PURE as possible. Purity is measured by an impurity criterion: Gini i(t) = 1 - sum_k p_k^2 (classification; default) Entropy i(t) = -sum_k p_k log p_k (classification) MSE i(t) = (1/n) sum (y - ybar)^2 (regression = variance) A split s sends a fraction p_L of the node's samples left and p_R right; its quality is the impurity DECREASE di = i(t) - p_L i(t_L) - p_R i(t_R), and at each node we search every feature and every threshold for the split that maximises it. The best split per feature is found in O(n log n) by sorting the feature once and sweeping the cut point with cumulative class counts (classification) or cumulative sums (regression), so the whole search is O(n_features * n log n). Grown to purity a tree OVERFITS. CART's remedy is COST-COMPLEXITY PRUNING: penalise the tree by the number of leaves, R_alpha(T) = R(T) + alpha|leaves(T)|, and for increasing alpha collapse the "weakest-link" subtree first -- the internal node t whose effective penalty alpha_eff(t) = (R(t) - R(T_t)) / (|leaves(T_t)| - 1) is smallest. This yields a nested sequence of subtrees indexed by alpha; the right alpha is picked by cross-validation. R(t) here is the node impurity weighted by the share of samples reaching it, matching scikit-learn's `ccp_alpha`. """ import numpy as np class _Node: __slots__ = ("feat", "thr", "left", "right", "value", "n", "imp") def __init__(self, value, n, imp): self.feat = None; self.thr = None; self.left = None; self.right = None self.value = value # leaf prediction: class-prob vector (clf) or mean (reg) self.n = n; self.imp = imp # samples reaching node; node impurity class DecisionTree: def __init__(self, criterion="gini", max_depth=None, min_samples_split=2, min_samples_leaf=1, ccp_alpha=0.0): self.criterion = criterion; self.max_depth = max_depth self.min_samples_split = min_samples_split; self.min_samples_leaf = min_samples_leaf self.ccp_alpha = ccp_alpha self.regression = (criterion == "mse") # ---- impurity of a set of labels ---- def _impurity(self, y): if self.regression: return np.mean((y - y.mean()) ** 2) if len(y) else 0.0 cnt = np.bincount(y, minlength=self.n_classes_).astype(float); p = cnt / cnt.sum() if self.criterion == "entropy": p = p[p > 0]; return float(-(p * np.log2(p)).sum()) return float(1.0 - (p * p).sum()) # gini def _leaf_value(self, y): if self.regression: return float(y.mean()) return np.bincount(y, minlength=self.n_classes_).astype(float) / len(y) # ---- best split over all features (cumulative sweep) ---- def _best_split(self, X, y): n, p = X.shape; parent_imp = self._impurity(y) best = (0.0, None, None) # (gain, feat, thr) for j in range(p): xs = X[:, j]; order = np.argsort(xs, kind="mergesort") xj = xs[order]; yj = y[order] cut = np.where(np.diff(xj) > 0)[0] # positions with a real threshold if len(cut) == 0: continue if self.regression: cs = np.cumsum(yj); cs2 = np.cumsum(yj * yj); tot = cs[-1]; tot2 = cs2[-1] nL = cut + 1.0; nR = n - nL sseL = cs2[cut] - cs[cut] ** 2 / nL sseR = (tot2 - cs2[cut]) - (tot - cs[cut]) ** 2 / nR child = (sseL + sseR) / n # weighted child MSE gain = parent_imp - child else: oh = np.zeros((n, self.n_classes_)); oh[np.arange(n), yj] = 1.0 cc = np.cumsum(oh, axis=0) # cumulative class counts nL = (cut + 1.0)[:, None]; nR = n - nL cL = cc[cut]; cR = cc[-1] - cL giniL = 1.0 - ((cL / nL) ** 2).sum(1); giniR = 1.0 - ((cR / nR) ** 2).sum(1) if self.criterion == "entropy": pL = cL / nL; pR = cR / nR giniL = -np.where(pL > 0, pL * np.log2(np.where(pL > 0, pL, 1)), 0).sum(1) giniR = -np.where(pR > 0, pR * np.log2(np.where(pR > 0, pR, 1)), 0).sum(1) child = (nL[:, 0] * giniL + nR[:, 0] * giniR) / n gain = parent_imp - child ok = (nL.reshape(-1) >= self.min_samples_leaf) & (nR.reshape(-1) >= self.min_samples_leaf) if not ok.any(): continue gain = np.where(ok, gain, -np.inf); k = int(np.argmax(gain)) if gain[k] > best[0]: best = (float(gain[k]), j, float((xj[cut[k]] + xj[cut[k] + 1]) / 2.0)) return best def _build(self, X, y, depth): node = _Node(self._leaf_value(y), len(y), self._impurity(y)) if (len(y) < self.min_samples_split or node.imp == 0.0 or (self.max_depth is not None and depth >= self.max_depth)): return node gain, feat, thr = self._best_split(X, y) if feat is None or gain <= 0: return node m = X[:, feat] <= thr node.feat = feat; node.thr = thr node.left = self._build(X[m], y[m], depth + 1) node.right = self._build(X[~m], y[~m], depth + 1) return node def fit(self, X, y): X = np.asarray(X, float); y = np.asarray(y) self.n_features_ = X.shape[1] if not self.regression: self.classes_, y = np.unique(y, return_inverse=True); self.n_classes_ = len(self.classes_) self.n_total_ = len(y) self.root = self._build(X, y, 0) if self.ccp_alpha > 0: self._prune(self.ccp_alpha) return self # ---- prediction ---- def _leaf(self, x, node): while node.feat is not None: node = node.left if x[node.feat] <= node.thr else node.right return node.value def predict_proba(self, X): X = np.asarray(X, float); return np.array([self._leaf(x, self.root) for x in X]) def predict(self, X): if self.regression: X = np.asarray(X, float); return np.array([self._leaf(x, self.root) for x in X]) return self.classes_[np.argmax(self.predict_proba(X), axis=1)] # ---- cost-complexity pruning ---- def _R(self, node): # weighted resubstitution risk of a node return node.imp * node.n / self.n_total_ def _collect(self, node, out): # internal nodes with (RT, nleaves) if node.feat is None: return self._R(node), 1 RTl, nl = self._collect(node.left, out); RTr, nr = self._collect(node.right, out) RT = RTl + RTr; nleaves = nl + nr out.append((node, self._R(node), RT, nleaves)) return RT, nleaves def cost_complexity_pruning_path(self): """weakest-link sequence: increasing effective alphas at which each subtree collapses.""" alphas = [0.0] import copy; root = copy.deepcopy(self.root); saved = self.root self.root = root while root.feat is not None: nodes = []; self._collect(root, nodes) nd, R, RT, nl = min(nodes, key=lambda z: (z[1] - z[2]) / (z[3] - 1)) alphas.append((R - RT) / (nl - 1)); nd.feat = None; nd.left = nd.right = None self.root = saved return np.array(alphas) def _prune(self, alpha): while self.root.feat is not None: nodes = []; self._collect(self.root, nodes) nd, R, RT, nl = min(nodes, key=lambda z: (z[1] - z[2]) / (z[3] - 1)) if (R - RT) / (nl - 1) > alpha: break nd.feat = None; nd.left = nd.right = None # ---- introspection ---- def n_leaves(self): def rec(nd): return 1 if nd.feat is None else rec(nd.left) + rec(nd.right) return rec(self.root) def get_depth(self): def rec(nd, d): return d if nd.feat is None else max(rec(nd.left, d + 1), rec(nd.right, d + 1)) return rec(self.root, 0) @property def feature_importances_(self): imp = np.zeros(getattr(self, "n_features_", self._infer_p())) def rec(nd): if nd.feat is None: return dec = nd.imp * nd.n - nd.left.imp * nd.left.n - nd.right.imp * nd.right.n imp[nd.feat] += dec; rec(nd.left); rec(nd.right) rec(self.root); s = imp.sum() return imp / s if s > 0 else imp def _infer_p(self): p = [0] def rec(nd): if nd.feat is None: return p[0] = max(p[0], nd.feat); rec(nd.left); rec(nd.right) rec(self.root); return p[0] + 1