{ "cells": [ { "cell_type": "markdown", "id": "25143010", "metadata": {}, "source": [ "# CS6140 Machine Learning — Fall 2026\n", "# HW1 Starter — Closed-form Regression, Decision Trees, Boosted Trees, KNN\n", "\n", "This notebook implements `HW1_26F.html`. For each problem you'll first solve it with a\n", "well-known library (a sanity-check baseline), then implement the same thing from scratch and\n", "evaluate it exactly the same way, so your two numbers should land right next to each other once\n", "everything is filled in.\n", "\n", "**Problems 1-4 (required):**\n", "1. Closed-form linear regression + Ridge (Housing, Spambase)\n", "2. Decision / regression tree from scratch (Spambase, Housing)\n", "3. Boosted regression trees, fit to residuals (Housing)\n", "4. KNN / fixed-window similarity classifiers (Spambase, MNIST)\n", "\n", "**Problems 5-6 (optional, no credit):**\n", "5. Second-order (Newton) boosting from the XGBoost paper -- a self-study mini-project (see\n", " `HW1_26F.html` for the reading + LLM-assisted implementation + no-LLM write-up structure).\n", "6. Does boosting confidence track KNN neighborhood consistency? -- a small research study that\n", " reuses models built in Problems 3-5.\n", "\n", "**How to work through this notebook.** Everything runs except the lines marked\n", "`### TODO_STUDENT ###` — those have their code replaced with a `...` placeholder and a one-line\n", "description of what to implement there. Read the surrounding (working) code and the markdown\n", "above each section for context, then replace each `...` with real code. A cell with a `...` in it\n", "will raise a `SyntaxError` or produce obviously wrong output until you fill it in — that's\n", "expected, not a bug. Work top to bottom: later cells (and later problems) depend on earlier\n", "classes (e.g. Problem 3's boosting reuses Problem 2's `DecisionTree`; Problem 6 reuses models from\n", "Problems 3-5), so finish each `TODO` before moving on to the next one that uses it. Once a\n", "\"library baseline\" cell's TODOs are filled in and run, the printed number is a target: your\n", "matching \"from scratch\" implementation just below it should land on (or very near) the same\n", "value.\n", "\n", "Datasets are read from the shared course `data/` folder (paths are relative to\n", "`1_intro_DT_RULES_REG/hw1/`, i.e. two levels up -- launch Jupyter from this notebook's own\n", "directory, or the relative path won't resolve).\n", "\n", "**Requirements (Problems 1-4 only -- Problems 5-6 are optional/no-credit and not required for\n", "either of these).**\n", "1. Fill in every `TODO_STUDENT` blank in Problems 1-4 according to the algorithm steps covered in\n", " lecture and the linked materials, and make sure the notebook runs top to bottom **through the\n", " end of Problem 4** without errors. (Problems 5-6's cells come after Problem 4 and contain their\n", " own, separate TODOs -- leaving those unfilled does not block Problems 1-4 from running; you\n", " simply stop reading/running before them if you're skipping the optional work.)\n", "2. Understand the *whole* notebook through Problem 4 — including the provided/given code, not just\n", " the blanks you filled in — well enough to explain any part of it during office hours.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ae5ec0bb", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.metrics import accuracy_score, mean_squared_error as mse_loss\n", "from sklearn.metrics.pairwise import pairwise_distances, rbf_kernel, polynomial_kernel\n", "\n", "DATA = \"../../data\"\n", "np.random.seed(42)\n" ] }, { "cell_type": "markdown", "id": "76b3d69d", "metadata": {}, "source": [ "## Data loading\n", "\n", "Housing comes pre-split (regression). Spambase is one file, so we split it ourselves. MNIST uses the 200-dim Haar-feature extraction (`mnist_haar_bingyu`) so that KNN's pairwise distances stay cheap enough to run on a laptop." ] }, { "cell_type": "code", "execution_count": null, "id": "c185b5e4", "metadata": {}, "outputs": [], "source": [ "# ---- Housing (regression): already split into train/test by the course ----\n", "train_housing = np.loadtxt(f\"{DATA}/housing_train.txt\")\n", "test_housing = np.loadtxt(f\"{DATA}/housing_test.txt\")\n", "\n", "X_train_housing, y_train_housing = train_housing[:, :-1], train_housing[:, -1]\n", "X_test_housing, y_test_housing = test_housing[:, :-1], test_housing[:, -1]\n", "\n", "# standardize using train statistics only\n", "housing_scaler = StandardScaler().fit(X_train_housing)\n", "X_train_housing = housing_scaler.transform(X_train_housing)\n", "X_test_housing = housing_scaler.transform(X_test_housing)\n", "\n", "print(\"Housing train:\", X_train_housing.shape, \" test:\", X_test_housing.shape)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "62d09c59", "metadata": {}, "outputs": [], "source": [ "# ---- Spambase (binary classification): one file, we split it ----\n", "spam_raw = np.loadtxt(f\"{DATA}/spambase/spambase.data\", delimiter=\",\")\n", "X_spam, y_spam = spam_raw[:, :-1], spam_raw[:, -1]\n", "\n", "X_train_spam, X_test_spam, y_train_spam, y_test_spam = train_test_split(\n", " X_spam, y_spam, test_size=0.2, random_state=42)\n", "\n", "spam_scaler = StandardScaler().fit(X_train_spam)\n", "X_train_spam = spam_scaler.transform(X_train_spam)\n", "X_test_spam = spam_scaler.transform(X_test_spam)\n", "\n", "print(\"Spambase train:\", X_train_spam.shape, \" test:\", X_test_spam.shape)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "040f6c42", "metadata": {}, "outputs": [], "source": [ "# ---- MNIST (Haar features, for KNN): pre-extracted 200-dim features ----\n", "# Subsample the training set (KNN is O(n_train) per query, and this is a from-scratch\n", "# implementation, not a KD-tree) so a laptop can run every (k, distance) combination below.\n", "X_train_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/training_image.txt\", delimiter=\",\")\n", "y_train_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/training_label.txt\", delimiter=\",\")\n", "X_test_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/testing_image.txt\", delimiter=\",\")\n", "y_test_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/testing_label.txt\", delimiter=\",\")\n", "\n", "n_train_digits, n_test_digits = 3000, 1000\n", "rng = np.random.RandomState(0)\n", "train_idx = rng.choice(len(X_train_digits_full), size=n_train_digits, replace=False)\n", "test_idx = rng.choice(len(X_test_digits_full), size=n_test_digits, replace=False)\n", "\n", "X_train_digits = X_train_digits_full[train_idx]\n", "y_train_digits = y_train_digits_full[train_idx]\n", "X_test_digits = X_test_digits_full[test_idx]\n", "y_test_digits = y_test_digits_full[test_idx]\n", "\n", "print(\"MNIST (Problem 4 split) train:\", X_train_digits.shape, \" test:\", X_test_digits.shape)\n" ] }, { "cell_type": "markdown", "id": "bb8e3ddc", "metadata": {}, "source": [ "---\n", "## Problem 1 — Closed-form Linear Regression + Ridge\n", "\n", "### (a) Library baseline: `sklearn.linear_model.LinearRegression` / `Ridge`\n", "\n", "Both are already closed-form solvers under the hood (`LinearRegression` calls a least-squares\n", "solver, `Ridge` solves the penalized normal equations), so they are the natural baseline for\n", "\"did I implement the normal equations correctly?\" For Spambase we reuse the real-valued\n", "regression output and threshold it at 0.5 to get a class label, exactly as the from-scratch\n", "version will.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "23003134", "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import LinearRegression, Ridge as SkRidge\n", "\n", "# Housing: library linear regression\n", "sk_lr = LinearRegression... ### TODO_STUDENT ### fit sklearn's LinearRegression as the library baseline for closed-form OLS on Housing\n", "mse_train = mse_loss(y_train_housing, sk_lr.predict(X_train_housing))\n", "mse_test = mse_loss(y_test_housing, sk_lr.predict(X_test_housing))\n", "print(f\"[library] Housing LinearRegression: train MSE = {mse_train:.3f}, test MSE = {mse_test:.3f}\")\n", "\n", "# Spambase: library linear regression, thresholded at 0.5 for classification\n", "sk_lr_spam = LinearRegression... ### TODO_STUDENT ### fit sklearn's LinearRegression as the library baseline on Spambase (thresholded below into a classifier)\n", "acc_train = accuracy_score(y_train_spam, sk_lr_spam.predict(X_train_spam) >= 0.5)\n", "acc_test = accuracy_score(y_test_spam, sk_lr_spam.predict(X_test_spam) >= 0.5)\n", "print(f\"[library] Spambase LinearRegression (threshold=0.5): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "bbb8bcfc", "metadata": {}, "outputs": [], "source": [ "# Housing: library ridge regression (same lambda=1 as our from-scratch version below)\n", "sk_ridge = SkRidge... ### TODO_STUDENT ### fit sklearn's Ridge(alpha=1.0) as the library baseline for closed-form ridge regression on Housing\n", "mse_train = mse_loss(y_train_housing, sk_ridge.predict(X_train_housing))\n", "mse_test = mse_loss(y_test_housing, sk_ridge.predict(X_test_housing))\n", "print(f\"[library] Housing Ridge(alpha=1): train MSE = {mse_train:.3f}, test MSE = {mse_test:.3f}\")\n", "\n", "# Spambase: library ridge, thresholded at 0.5\n", "sk_ridge_spam = SkRidge... ### TODO_STUDENT ### fit sklearn's Ridge(alpha=1.0) as the library baseline on Spambase\n", "acc_train = accuracy_score(y_train_spam, sk_ridge_spam.predict(X_train_spam) >= 0.5)\n", "acc_test = accuracy_score(y_test_spam, sk_ridge_spam.predict(X_test_spam) >= 0.5)\n", "print(f\"[library] Spambase Ridge(alpha=1, threshold=0.5): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "9caf3098", "metadata": {}, "source": [ "### (b) From scratch: Normal Equations / Ridge implemented by hand\n", "\n", "Normal equations: $w = (X^TX)^{-1}X^Ty$, and the Ridge version $w = (X^TX + \\lambda I)^{-1}X^Ty$.\n", "Both are the exact linear-algebra solution, no iteration needed. Same datasets, same threshold,\n", "so the printed numbers should land right next to the library ones above.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e157a277", "metadata": {}, "outputs": [], "source": [ "class NormalEquationRegressor:\n", " \"\"\"Linear regression via the normal equations (closed form).\"\"\"\n", "\n", " def __init__(self, bias: bool = True):\n", " self.bias = bias\n", "\n", " def fit(self, X, y):\n", " n = X.shape[0]\n", " X_aug = np.concatenate([X, np.ones((n, 1))], axis=1) if self.bias else X\n", " self.w = np.linalg.... ### TODO_STUDENT ### implement the normal-equation closed form w = (XtX)^-1 Xty\n", "\n", " def predict(self, X):\n", " n = X.shape[0]\n", " X_aug = np.concatenate([X, np.ones((n, 1))], axis=1) if self.bias else X\n", " return ... ### TODO_STUDENT ### compute predictions as X @ w\n", "\n", "\n", "class NormalEquationClassifier(NormalEquationRegressor):\n", " \"\"\"Same closed-form fit, but predictions are thresholded into {0, 1}.\"\"\"\n", "\n", " def __init__(self, bias: bool = True, threshold: float = 0.5):\n", " super().__init__(bias)\n", " self.threshold = threshold\n", "\n", " def predict(self, X):\n", " scores = super().predict(X)\n", " return ... ### TODO_STUDENT ### threshold the regression score into a 0/1 class label\n" ] }, { "cell_type": "code", "execution_count": null, "id": "be48b434", "metadata": {}, "outputs": [], "source": [ "housing_ols = NormalEquationRegressor()\n", "housing_ols.fit(X_train_housing, y_train_housing)\n", "\n", "mse_train = mse_loss(y_train_housing, housing_ols.predict(X_train_housing))\n", "mse_test = mse_loss(y_test_housing, housing_ols.predict(X_test_housing))\n", "print(f\"[scratch] Housing (Normal Equations): train MSE = {mse_train:.3f}, test MSE = {mse_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c2c09c28", "metadata": {}, "outputs": [], "source": [ "spam_ols = NormalEquationClassifier(threshold=0.5)\n", "spam_ols.fit(X_train_spam, y_train_spam)\n", "\n", "acc_train = accuracy_score(y_train_spam, spam_ols.predict(X_train_spam))\n", "acc_test = accuracy_score(y_test_spam, spam_ols.predict(X_test_spam))\n", "print(f\"[scratch] Spambase (Normal Equations, threshold=0.5): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b5e84481", "metadata": {}, "outputs": [], "source": [ "class RidgeRegressor:\n", " \"\"\"Ridge regression via the closed-form penalized normal equations.\"\"\"\n", "\n", " def __init__(self, lmbda: float = 1.0, bias: bool = True):\n", " self.lmbda = lmbda\n", " self.bias = bias\n", "\n", " def fit(self, X, y):\n", " n, m = X.shape\n", " X_aug = np.concatenate([X, np.ones((n, 1))], axis=1) if self.bias else X\n", " m_aug = X_aug.shape[1]\n", " # don't penalize the bias term\n", " penalty = np.eye(m_aug) * self.lmbda\n", " if self.bias:\n", " penalty[-1, -1] = ... ### TODO_STUDENT ### don't penalize the bias/intercept term when building the ridge penalty matrix\n", " self.w = np.linalg.... ### TODO_STUDENT ### implement the closed-form ridge solution w = (XtX + lambda*I)^-1 Xty\n", "\n", " def predict(self, X):\n", " n = X.shape[0]\n", " X_aug = np.concatenate([X, np.ones((n, 1))], axis=1) if self.bias else X\n", " return ... ### TODO_STUDENT ### compute predictions as X @ w\n", "\n", "\n", "class RidgeClassifier(RidgeRegressor):\n", " def __init__(self, lmbda: float = 1.0, bias: bool = True, threshold: float = 0.5):\n", " super().__init__(lmbda, bias)\n", " self.threshold = threshold\n", "\n", " def predict(self, X):\n", " scores = super().predict(X)\n", " return ... ### TODO_STUDENT ### threshold the ridge regression score into a 0/1 class label\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9d27f56e", "metadata": {}, "outputs": [], "source": [ "housing_ridge = RidgeRegressor(lmbda=1.0)\n", "housing_ridge.fit(X_train_housing, y_train_housing)\n", "\n", "mse_train = mse_loss(y_train_housing, housing_ridge.predict(X_train_housing))\n", "mse_test = mse_loss(y_test_housing, housing_ridge.predict(X_test_housing))\n", "print(f\"[scratch] Housing (Ridge, lambda=1): train MSE = {mse_train:.3f}, test MSE = {mse_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e0163b6b", "metadata": {}, "outputs": [], "source": [ "spam_ridge = RidgeClassifier(lmbda=1.0, threshold=0.5)\n", "spam_ridge.fit(X_train_spam, y_train_spam)\n", "\n", "acc_train = accuracy_score(y_train_spam, spam_ridge.predict(X_train_spam))\n", "acc_test = accuracy_score(y_test_spam, spam_ridge.predict(X_test_spam))\n", "print(f\"[scratch] Spambase (Ridge, lambda=1, threshold=0.5): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "58b5c8da", "metadata": {}, "source": [ "### Bonus plot: Ridge regularization path (Housing)\n", "\n", "Sweep `lambda` across several orders of magnitude instead of the single fixed value above, and\n", "watch the classic bias-variance trade-off: training MSE rises monotonically as the penalty grows,\n", "while test MSE dips before rising again once the penalty overwhelms the signal.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dcadbbb8", "metadata": {}, "outputs": [], "source": [ "lambdas = np.logspace(-2, 3, 12)\n", "\n", "lib_ridge_train, lib_ridge_test = [], []\n", "scratch_ridge_train, scratch_ridge_test = [], []\n", "for lmbda in lambdas:\n", " sk_r = SkRidge(alpha=lmbda).fit(X_train_housing, y_train_housing)\n", " lib_ridge_train.append(mse_loss(y_train_housing, sk_r.predict(X_train_housing)))\n", " lib_ridge_test.append(mse_loss(y_test_housing, sk_r.predict(X_test_housing)))\n", "\n", " r = RidgeRegressor(lmbda=lmbda)\n", " r.fit(X_train_housing, y_train_housing)\n", " scratch_ridge_train.append(mse_loss(y_train_housing, r.predict(X_train_housing)))\n", " scratch_ridge_test.append(mse_loss(y_test_housing, r.predict(X_test_housing)))\n", "\n", "plt.figure(figsize=(6, 4))\n", "plt.plot(lambdas, lib_ridge_train, color=\"#2a78d6\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library train\")\n", "plt.plot(lambdas, lib_ridge_test, color=\"#e34948\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library test\")\n", "plt.plot(lambdas, scratch_ridge_train, color=\"#2a78d6\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch train\")\n", "plt.plot(lambdas, scratch_ridge_test, color=\"#e34948\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch test\")\n", "plt.xscale(\"log\")\n", "plt.xlabel(\"lambda (log scale)\")\n", "plt.ylabel(\"MSE\")\n", "plt.title(\"Ridge regularization path on Housing: library vs. from scratch\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "cc66e2b5", "metadata": {}, "source": [ "---\n", "## Problem 2 — Decision / Regression Tree from scratch\n", "\n", "### (a) Library baseline: `sklearn.tree.DecisionTreeClassifier` / `DecisionTreeRegressor`\n", "\n", "Same split criteria (`entropy` for Spambase, `squared_error`/variance for Housing) and the same\n", "depth sweep, so the train/test numbers can be read side by side with the from-scratch tree.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "15b9b752", "metadata": {}, "outputs": [], "source": [ "from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor\n", "\n", "depths = [2, 4, 6, 8]\n", "lib_tree_spam_train, lib_tree_spam_test = [], []\n", "for depth in depths:\n", " clf = DecisionTreeClassifier(criterion=\"entropy\", max_depth=depth, random_state=0)\n", " clf... ### TODO_STUDENT ### fit sklearn's DecisionTreeClassifier(criterion='entropy') as the library baseline\n", " acc_train = accuracy_score(y_train_spam, clf.predict(X_train_spam))\n", " acc_test = accuracy_score(y_test_spam, clf.predict(X_test_spam))\n", " lib_tree_spam_train.append(acc_train)\n", " lib_tree_spam_test.append(acc_test)\n", " print(f\"[library] Spambase tree depth={depth}: train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "317b62a3", "metadata": {}, "outputs": [], "source": [ "lib_tree_house_train, lib_tree_house_test = [], []\n", "for depth in depths:\n", " reg = DecisionTreeRegressor(criterion=\"squared_error\", max_depth=depth, random_state=0)\n", " reg... ### TODO_STUDENT ### fit sklearn's DecisionTreeRegressor(criterion='squared_error') as the library baseline\n", " mse_train = mse_loss(y_train_housing, reg.predict(X_train_housing))\n", " mse_test = mse_loss(y_test_housing, reg.predict(X_test_housing))\n", " lib_tree_house_train.append(mse_train)\n", " lib_tree_house_test.append(mse_test)\n", " print(f\"[library] Housing tree depth={depth}: train MSE = {mse_train:.3f}, test MSE = {mse_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "7333e7b2", "metadata": {}, "source": [ "### (b) From scratch: entropy / variance-reduction tree\n", "\n", "A single tree class handles both cases: **entropy / information gain** for classification\n", "(Spambase), **variance / MSE reduction** for regression (Housing). Since features are\n", "continuous, candidate splits are the midpoints between consecutive sorted unique values of a\n", "feature (this is the standard trick — it's exactly the threshold that could possibly separate\n", "two adjacent points, so we never need to check more than that).\n", "\n", "This tree is written as a reusable module (`fit`/`predict`): HW2's gradient-boosting problem and\n", "HW1 Problem 3 below both reuse it as the weak learner.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dd3e7f7d", "metadata": {}, "outputs": [], "source": [ "def entropy(y):\n", " _, counts = np.unique(y, return_counts=True)\n", " p = counts / len(y)\n", " return -np.... ### TODO_STUDENT ### implement the Shannon entropy formula, -sum(p * log2(p))\n", "\n", "\n", "def variance(y):\n", " if len(y) == 0:\n", " return 0.0\n", " return np.... ### TODO_STUDENT ### implement the variance (mean squared deviation from the mean) formula\n", "\n", "\n", "def candidate_thresholds(column):\n", " \"\"\"Midpoints between consecutive sorted unique values -- the only splits that can matter.\"\"\"\n", " values = np.unique(column)\n", " return ... ### TODO_STUDENT ### compute the midpoints between consecutive sorted unique values as candidate split thresholds\n", "\n", "\n", "class TreeNode:\n", " def __init__(self, feature=None, threshold=None, left=None, right=None, value=None):\n", " self.feature = feature\n", " self.threshold = threshold\n", " self.left = left\n", " self.right = right\n", " self.value = value # only set on leaves\n", "\n", " def is_leaf(self):\n", " return self.value is not None\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b1ce8ab7", "metadata": {}, "outputs": [], "source": [ "class DecisionTree:\n", " \"\"\"Binary decision/regression tree with threshold splits on numeric features.\n", "\n", " mode='classification' splits on information gain (entropy) and leaves store the\n", " majority class; mode='regression' splits on variance reduction and leaves store the mean.\n", " \"\"\"\n", "\n", " def __init__(self, mode=\"classification\", max_depth=5, min_samples_split=2):\n", " assert mode in (\"classification\", \"regression\")\n", " self.mode = mode\n", " self.max_depth = max_depth\n", " self.min_samples_split = min_samples_split\n", " self.impurity = ... ### TODO_STUDENT ### pick the impurity function for this tree: entropy for classification splits, variance for regression splits\n", "\n", " def fit(self, X, y):\n", " self.root = self._grow(X, y, depth=0)\n", " return self\n", "\n", " def _leaf_value(self, y):\n", " if self.mode == \"classification\":\n", " values, counts = np.unique(y, return_counts=True)\n", " return ... ### TODO_STUDENT ### return the majority class label among this leaf's points\n", " return np.... ### TODO_STUDENT ### return the mean target value among this leaf's points\n", "\n", " def _grow(self, X, y, depth):\n", " n, m = X.shape\n", "\n", " # stopping conditions: max depth reached, pure node, or too few points to split\n", " if ...: ### TODO_STUDENT ### decide when to stop splitting and turn this node into a leaf\n", " return TreeNode(value=self._leaf_value(y))\n", "\n", " parent_impurity = self.impurity(y)\n", " best_gain, best_feature, best_threshold = 0.0, None, None\n", "\n", " for feature in range(m):\n", " for threshold in candidate_thresholds(X[:, feature]):\n", " left_mask = ... ### TODO_STUDENT ### define which points fall left vs. right of this candidate threshold\n", " right_mask = ~left_mask\n", " if ...: ### TODO_STUDENT ### skip candidate thresholds that would leave one side empty\n", " continue\n", "\n", " y_left, y_right = y[left_mask], y[right_mask]\n", " weighted_child_impurity = ... ### TODO_STUDENT ### weight each child's impurity by the fraction of the parent's points it received\n", " gain = ... ### TODO_STUDENT ### compute the information/variance gain of this candidate split (parent impurity minus weighted child impurity)\n", "\n", " if gain > best_gain:\n", " best_gain, best_feature, best_threshold = ... ### TODO_STUDENT ### keep track of the best (feature, threshold) split seen so far\n", "\n", " # no split improves things -> make a leaf\n", " if best_feature is None:\n", " return TreeNode(value=self._leaf_value(y))\n", "\n", " left_mask = ... ### TODO_STUDENT ### recompute the partition using the winning (best_feature, best_threshold) before recursing\n", " left = self._grow... ### TODO_STUDENT ### recursively grow the left subtree on the left split\n", " right = self._grow... ### TODO_STUDENT ### recursively grow the right subtree on the right split\n", " return TreeNode(feature=best_feature, threshold=best_threshold, left=left, right=right)\n", "\n", " def _predict_one(self, x, node):\n", " if node.is_leaf():\n", " return node.value\n", " branch = ... ### TODO_STUDENT ### walk left or right at this node based on the split threshold\n", " return self._predict_one(x, branch)\n", "\n", " def predict(self, X):\n", " return np.array([self._predict_one(x, self.root) for x in X])\n" ] }, { "cell_type": "code", "execution_count": null, "id": "97887d97", "metadata": {}, "outputs": [], "source": [ "# Spambase: classification tree (information gain), a few depths to see over/under-fitting\n", "scratch_tree_spam_train, scratch_tree_spam_test = [], []\n", "for depth in depths:\n", " tree = DecisionTree(mode=\"classification\", max_depth=depth).fit(X_train_spam, y_train_spam)\n", " acc_train = accuracy_score(y_train_spam, tree.predict(X_train_spam))\n", " acc_test = accuracy_score(y_test_spam, tree.predict(X_test_spam))\n", " scratch_tree_spam_train.append(acc_train)\n", " scratch_tree_spam_test.append(acc_test)\n", " print(f\"[scratch] Spambase tree depth={depth}: train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f390c57d", "metadata": {}, "outputs": [], "source": [ "# Housing: regression tree (variance reduction), same depth sweep\n", "scratch_tree_house_train, scratch_tree_house_test = [], []\n", "for depth in depths:\n", " tree = DecisionTree(mode=\"regression\", max_depth=depth).fit(X_train_housing, y_train_housing)\n", " mse_train = mse_loss(y_train_housing, tree.predict(X_train_housing))\n", " mse_test = mse_loss(y_test_housing, tree.predict(X_test_housing))\n", " scratch_tree_house_train.append(mse_train)\n", " scratch_tree_house_test.append(mse_test)\n", " print(f\"[scratch] Housing tree depth={depth}: train MSE = {mse_train:.3f}, test MSE = {mse_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1ff9afc0", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "\n", "axes[0].plot(depths, lib_tree_spam_train, color=\"#2a78d6\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library train\")\n", "axes[0].plot(depths, lib_tree_spam_test, color=\"#e34948\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library test\")\n", "axes[0].plot(depths, scratch_tree_spam_train, color=\"#2a78d6\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch train\")\n", "axes[0].plot(depths, scratch_tree_spam_test, color=\"#e34948\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch test\")\n", "axes[0].set_xlabel(\"Tree depth\")\n", "axes[0].set_ylabel(\"Accuracy\")\n", "axes[0].set_title(\"Spambase classification tree\")\n", "axes[0].legend()\n", "\n", "axes[1].plot(depths, lib_tree_house_train, color=\"#2a78d6\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library train\")\n", "axes[1].plot(depths, lib_tree_house_test, color=\"#e34948\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library test\")\n", "axes[1].plot(depths, scratch_tree_house_train, color=\"#2a78d6\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch train\")\n", "axes[1].plot(depths, scratch_tree_house_test, color=\"#e34948\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch test\")\n", "axes[1].set_xlabel(\"Tree depth\")\n", "axes[1].set_ylabel(\"MSE\")\n", "axes[1].set_title(\"Housing regression tree\")\n", "axes[1].legend()\n", "\n", "fig.suptitle(\"Train/test error vs. tree depth: library vs. from scratch\")\n", "fig.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "c64808f0", "metadata": {}, "source": [ "Both sweeps (library and from-scratch) should show the classic pattern: deeper trees keep improving train error while test error flattens out or gets worse. The right depth is the one where *test* performance stops improving (the bias-variance tradeoff), not the one where train and test happen to be numerically closest -- a very shallow, underfit tree can have a tiny train/test gap while both numbers are simply bad, which is not what we want. In practice you'd pick a final depth via a validation split or cross-validation, not by reading the test curve directly (that's tuning on the test set); the sweep here is a diagnostic showing you the tradeoff, not a tuning procedure. Check that your from-scratch numbers track the library ones closely at every depth on Spambase -- that's the sanity check to look for. On Housing, don't be surprised if the two start to diverge on *test* MSE at deeper depths (4+) even though both still look reasonable by *training* loss: with only 433 training rows and 13 correlated features, many splits tie or nearly tie on variance reduction, so `sklearn`'s tie-breaking and yours can legitimately pick different (but similarly good) trees deeper down. That would be expected small-sample variance, not a bug -- and it's itself a nice illustration of the instability of deep single trees, which Problem 3's boosting exists to tame." ] }, { "cell_type": "markdown", "id": "9a532e88", "metadata": {}, "source": [ "---\n", "## Problem 3 — Boosted Trees for Regression (Housing)\n", "\n", "### (a) Library baseline: `sklearn.ensemble.GradientBoostingRegressor`\n", "\n", "With `loss=\"squared_error\"`, this *is* residual-fitting boosting (each new tree fits the\n", "current residual). We set `learning_rate=1.0` to match our from-scratch version, which takes a\n", "full step at each round rather than shrinking it.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "104ab109", "metadata": {}, "outputs": [], "source": [ "from sklearn.ensemble import GradientBoostingRegressor\n", "\n", "n_rounds = 60\n", "sk_gbr = GradientBoostingRegressor(\n", " loss=\"squared_error\", n_estimators=n_rounds, max_depth=2, learning_rate=1.0, random_state=0\n", ")... ### TODO_STUDENT ### fit sklearn's GradientBoostingRegressor(learning_rate=1.0) as the residual-boosting library baseline\n", "\n", "lib_train_mses = [mse_loss(y_train_housing, pred) for pred in sk_gbr.staged_predict(X_train_housing)]\n", "lib_test_mses = [mse_loss(y_test_housing, pred) for pred in sk_gbr.staged_predict(X_test_housing)]\n", "\n", "print(f\"[library] Final (T={n_rounds}): train MSE = {lib_train_mses[-1]:.3f}, test MSE = {lib_test_mses[-1]:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "d55d2be9", "metadata": {}, "source": [ "### (b) From scratch: residual boosting with our own tree\n", "\n", "Fit a sequence of shallow regression trees to the *residuals* of the running prediction (no\n", "re-weighting of points, unlike AdaBoost). Each new tree nudges the ensemble's prediction a\n", "little closer to the true labels. We reuse the `DecisionTree` class from Problem 2 (depth 2) as\n", "the weak learner, exactly as the assignment asks.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4217c7d3", "metadata": {}, "outputs": [], "source": [ "class BoostedRegressionTrees:\n", " def __init__(self, n_rounds=100, tree_depth=2):\n", " self.n_rounds = n_rounds\n", " self.tree_depth = tree_depth\n", "\n", " def fit(self, X, y):\n", " self.trees = []\n", " residual = y.copy()\n", " for _ in range(self.n_rounds):\n", " tree = DecisionTree... ### TODO_STUDENT ### fit a shallow regression tree to the current residuals\n", " self.trees.append(tree)\n", " residual = ... ### TODO_STUDENT ### update the residuals by subtracting this round's tree predictions\n", " return self\n", "\n", " def predict(self, X):\n", " return ... ### TODO_STUDENT ### sum every tree's prediction to form the boosted ensemble output\n", "\n", " def staged_predict(self, X):\n", " \"\"\"Yield the ensemble prediction after each additional tree (for the T-sweep plot).\"\"\"\n", " total = np.zeros(X.shape[0])\n", " for tree in self.trees:\n", " total = ... ### TODO_STUDENT ### accumulate this round's tree prediction into the running ensemble total\n", " yield total\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a5336d9e", "metadata": {}, "outputs": [], "source": [ "boosted = BoostedRegressionTrees(n_rounds=n_rounds, tree_depth=2).fit(X_train_housing, y_train_housing)\n", "\n", "train_mses, test_mses = [], []\n", "for train_pred, test_pred in zip(boosted.staged_predict(X_train_housing),\n", " boosted.staged_predict(X_test_housing)):\n", " train_mses.append(mse_loss(y_train_housing, train_pred))\n", " test_mses.append(mse_loss(y_test_housing, test_pred))\n", "\n", "print(f\"[scratch] Final (T={n_rounds}): train MSE = {train_mses[-1]:.3f}, test MSE = {test_mses[-1]:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2c2c1dc9", "metadata": {}, "outputs": [], "source": [ "# Library vs. scratch, overlaid: both should tell the same overfitting story as T grows.\n", "plt.figure(figsize=(6, 4))\n", "rounds = range(1, n_rounds + 1)\n", "plt.plot(rounds, lib_train_mses, \"b--\", label=\"library train MSE\")\n", "plt.plot(rounds, lib_test_mses, \"r--\", label=\"library test MSE\")\n", "plt.plot(rounds, train_mses, \"b-\", label=\"scratch train MSE\")\n", "plt.plot(rounds, test_mses, \"r-\", label=\"scratch test MSE\")\n", "plt.xlabel(\"Number of boosting rounds T\")\n", "plt.ylabel(\"MSE\")\n", "plt.title(\"Boosted regression trees on Housing: library vs. from scratch\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "2c3f301c", "metadata": {}, "source": [ "---\n", "## Problem 4 — Similarity & KNN from scratch\n", "\n", "### (a) Library baseline: `sklearn.neighbors.KNeighborsClassifier` / `RadiusNeighborsClassifier`\n", "\n", "`KNeighborsClassifier` supports `'euclidean'` and `'cosine'` directly. For the RBF and\n", "degree-2-polynomial *kernel* similarities there's no ready-made \"kernel-KNN\" in sklearn, so we\n", "precompute the kernel matrix with `sklearn.metrics.pairwise` (still 100% library code) and hand\n", "it to `KNeighborsClassifier(metric=\"precomputed\")` — the library still owns the neighbor search\n", "and the majority vote, which is the part we're benchmarking.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "eba004f0", "metadata": {}, "outputs": [], "source": [ "from sklearn.neighbors import KNeighborsClassifier, RadiusNeighborsClassifier\n", "\n", "# (A) Fixed k -- Spambase with Euclidean distance\n", "for k in [1, 3, 7]:\n", " knn = KNeighborsClassifier... ### TODO_STUDENT ### fit sklearn's KNeighborsClassifier with Euclidean distance as the library baseline\n", " acc_train = accuracy_score(y_train_spam, knn.predict(X_train_spam))\n", " acc_test = accuracy_score(y_test_spam, knn.predict(X_test_spam))\n", " print(f\"[library] Spambase kNN (euclidean, k={k}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4815c7a6", "metadata": {}, "outputs": [], "source": [ "# (A) Fixed k -- MNIST with cosine distance (native metric) and RBF / poly-2 kernels (precomputed)\n", "for k in [1, 3, 7]:\n", " knn = KNeighborsClassifier... ### TODO_STUDENT ### fit sklearn's KNeighborsClassifier with cosine distance as the library baseline\n", " acc_train = accuracy_score(y_train_digits, knn.predict(X_train_digits))\n", " acc_test = accuracy_score(y_test_digits, knn.predict(X_test_digits))\n", " print(f\"[library] MNIST kNN (cosine, k={k}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n", "\n", "kernel_specs = {\n", " \"rbf\": lambda A, B: ..., ### TODO_STUDENT ### compute the RBF (Gaussian) kernel similarity matrix between two point sets\n", " \"poly\": lambda A, B: ..., ### TODO_STUDENT ### compute the degree-2 polynomial kernel similarity matrix between two point sets\n", "}\n", "for name, kernel_fn in kernel_specs.items():\n", " K_train = kernel_fn(X_train_digits, X_train_digits)\n", " K_test = kernel_fn(X_test_digits, X_train_digits)\n", " # sklearn's metric=\"precomputed\" requires a non-negative *distance* matrix; \"larger kernel\n", " # value = closer\" flips to \"smaller distance = closer\" via this order-reversing shift\n", " shift = ... ### TODO_STUDENT ### pick one global shift so both kernel blocks convert to comparable non-negative distances\n", " D_train = shift - K_train\n", " D_test = shift - K_test\n", " for k in [1, 3, 7]:\n", " knn = KNeighborsClassifier... ### TODO_STUDENT ### fit sklearn's KNeighborsClassifier on the precomputed kernel-distance matrix\n", " acc_train = accuracy_score(y_train_digits, knn.predict(D_train))\n", " acc_test = accuracy_score(y_test_digits, knn.predict(D_test))\n", " print(f\"[library] MNIST kNN ({name}, k={k}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c903ce2b", "metadata": {}, "outputs": [], "source": [ "# (B) Fixed window -- Spambase with Euclidean distance\n", "# most_frequent fallback for points with an empty radius, matching the spirit of our\n", "# from-scratch fallback below (predict *something* rather than erroring out)\n", "spam_radii = [4.0, 6.0, 8.0, 10.0]\n", "lib_window_spam_train, lib_window_spam_test = [], []\n", "for r in spam_radii:\n", " rnn = RadiusNeighborsClassifier(radius=r, metric=\"euclidean\", outlier_label=\"most_frequent\")\n", " rnn... ### TODO_STUDENT ### fit sklearn's RadiusNeighborsClassifier as the library baseline for the fixed-window classifier on Spambase\n", " acc_train = accuracy_score(y_train_spam, rnn.predict(X_train_spam))\n", " acc_test = accuracy_score(y_test_spam, rnn.predict(X_test_spam))\n", " lib_window_spam_train.append(acc_train)\n", " lib_window_spam_test.append(acc_test)\n", " print(f\"[library] Spambase window (euclidean, r={r}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "265d8af9", "metadata": {}, "outputs": [], "source": [ "# (B) Fixed window -- MNIST with cosine distance\n", "digits_radii = [0.05, 0.1, 0.15, 0.2]\n", "lib_window_digits_train, lib_window_digits_test = [], []\n", "for r in digits_radii:\n", " rnn = RadiusNeighborsClassifier(radius=r, metric=\"cosine\", algorithm=\"brute\", outlier_label=\"most_frequent\")\n", " rnn... ### TODO_STUDENT ### fit sklearn's RadiusNeighborsClassifier as the library baseline for the fixed-window classifier on MNIST\n", " acc_train = accuracy_score(y_train_digits, rnn.predict(X_train_digits))\n", " acc_test = accuracy_score(y_test_digits, rnn.predict(X_test_digits))\n", " lib_window_digits_train.append(acc_train)\n", " lib_window_digits_test.append(acc_test)\n", " print(f\"[library] MNIST window (cosine, r={r}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "d7ec16b6", "metadata": {}, "source": [ "At the tightest radius (`r=0.05`), expect a meaningful fraction of the 1000 MNIST test points to have *zero* training points within range (cosine distance is a tight metric on 200-d Haar features). The library's `outlier_label=\"most_frequent\"` dumps every one of those into the single most common training digit -- watch what that does to test accuracy at this radius. Our from-scratch fallback below instead backs off to the single nearest point for those cases, which should be a much better rule for a well-separated multiclass problem: look for a noticeably higher test accuracy than the library's fallback at this same tight radius, with the gap closing up as the radius grows and empty neighborhoods become rare." ] }, { "cell_type": "markdown", "id": "ed59f379", "metadata": {}, "source": [ "### (b) From scratch: KNN and fixed-window classifiers\n", "\n", "No training phase: at query time we compute the distance/similarity from each test point to\n", "every training point, then either take a **majority vote over the k nearest** (part A) or a\n", "**majority vote over everyone inside a fixed radius/window** (part B).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4f0ee98a", "metadata": {}, "outputs": [], "source": [ "class KNNClassifier:\n", " \"\"\"k-NN (or fixed-window) classifier with a choice of distance/similarity function.\n", "\n", " If `radius` is set, `k` is ignored: every training point within `radius` votes.\n", " `metric` is one of {'euclidean', 'cosine', 'rbf', 'poly'} (the last two are kernel\n", " similarities, so larger = closer -- we negate them to reuse the same \"smallest distance\n", " wins\" logic).\n", " \"\"\"\n", "\n", " def __init__(self, k=1, radius=None, metric=\"euclidean\", gamma=1.0, degree=2, coef0=0.0):\n", " self.k = k\n", " self.radius = radius\n", " self.metric = metric\n", " self.gamma = gamma\n", " self.degree = degree\n", " self.coef0 = coef0\n", "\n", " def fit(self, X, y):\n", " self.X_train = X\n", " self.y_train = y\n", " return self\n", "\n", " def _distances(self, X_test):\n", " if self.metric == \"euclidean\":\n", " return ... ### TODO_STUDENT ### compute the Euclidean distance from every test point to every training point\n", " if self.metric == \"cosine\":\n", " return ... ### TODO_STUDENT ### compute the cosine distance from every test point to every training point\n", " if self.metric == \"rbf\":\n", " return -... ### TODO_STUDENT ### compute the (negated) RBF kernel similarity as a distance surrogate\n", " if self.metric == \"poly\":\n", " return -... ### TODO_STUDENT ### compute the (negated) degree-d polynomial kernel similarity as a distance surrogate\n", " raise ValueError(f\"unknown metric {self.metric}\")\n", "\n", " def _vote(self, neighbor_labels):\n", " values, counts = np.unique(neighbor_labels, return_counts=True)\n", " return ... ### TODO_STUDENT ### return the majority label among the neighbor votes\n", "\n", " def predict(self, X_test):\n", " D = self._distances(X_test)\n", " predictions = np.empty(X_test.shape[0], dtype=self.y_train.dtype)\n", "\n", " for i in range(X_test.shape[0]):\n", " row = D[i]\n", " if self.radius is None:\n", " nearest = np.... ### TODO_STUDENT ### select the indices of the k nearest training points\n", " else:\n", " nearest = np.... ### TODO_STUDENT ### select the indices of every training point within the fixed radius\n", " if len(nearest) == 0:\n", " # nothing in the window: fall back to the single nearest point\n", " nearest = np.... ### TODO_STUDENT ### when nothing falls inside the radius, back off to the single nearest point rather than refusing to predict\n", " predictions[i] = self._vote(self.y_train[nearest])\n", "\n", " return predictions\n" ] }, { "cell_type": "code", "execution_count": null, "id": "81b4f693", "metadata": {}, "outputs": [], "source": [ "# (A) Fixed k -- Spambase with Euclidean distance\n", "for k in [1, 3, 7]:\n", " knn = KNNClassifier(k=k, metric=\"euclidean\").fit(X_train_spam, y_train_spam)\n", " acc_train = accuracy_score(y_train_spam, knn.predict(X_train_spam))\n", " acc_test = accuracy_score(y_test_spam, knn.predict(X_test_spam))\n", " print(f\"[scratch] Spambase kNN (euclidean, k={k}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "56b49848", "metadata": {}, "outputs": [], "source": [ "# (A) Fixed k -- MNIST with cosine distance, Gaussian (RBF) kernel, degree-2 polynomial kernel\n", "k_values = [1, 3, 7]\n", "digits_knn_by_metric = {}\n", "for metric, kwargs in [(\"cosine\", {}), (\"rbf\", {\"gamma\": 1e-3}), (\"poly\", {\"degree\": 2, \"gamma\": 1e-3})]:\n", " train_accs, test_accs = [], []\n", " for k in k_values:\n", " knn = KNNClassifier(k=k, metric=metric, **kwargs).fit(X_train_digits, y_train_digits)\n", " acc_train = accuracy_score(y_train_digits, knn.predict(X_train_digits))\n", " acc_test = accuracy_score(y_test_digits, knn.predict(X_test_digits))\n", " train_accs.append(acc_train)\n", " test_accs.append(acc_test)\n", " print(f\"[scratch] MNIST kNN ({metric}, k={k}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n", " digits_knn_by_metric[metric] = (train_accs, test_accs)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "574e76de", "metadata": {}, "outputs": [], "source": [ "metric_colors = {\"cosine\": \"#2a78d6\", \"rbf\": \"#eb6834\", \"poly\": \"#1baf7a\"}\n", "\n", "plt.figure(figsize=(6, 4))\n", "for metric, (train_accs, test_accs) in digits_knn_by_metric.items():\n", " color = metric_colors[metric]\n", " plt.plot(k_values, train_accs, color=color, linestyle=\"--\", marker=\"o\", markersize=5, label=f\"{metric} train\")\n", " plt.plot(k_values, test_accs, color=color, linestyle=\"-\", marker=\"o\", markersize=5, label=f\"{metric} test\")\n", "plt.xlabel(\"k\")\n", "plt.ylabel(\"Accuracy\")\n", "plt.title(\"MNIST kNN accuracy vs. k, by distance/kernel (from scratch)\")\n", "plt.xticks(k_values)\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "1ef152c9", "metadata": {}, "source": [ "Expect cosine and RBF to land well above the degree-2 polynomial kernel here (matching the library numbers above) -- and if the polynomial kernel looks surprisingly weak in both the library and scratch versions, that's not necessarily a bug, it may be the kernel: both use `coef0=0` (a *homogeneous* poly kernel), and for `coef0=0` with any positive `gamma`, rescaling `gamma` only rescales every entry by the same monotone function, so it can never change the neighbor *ranking*. (This is specific to `coef0=0` -- an *inhomogeneous* kernel, `coef0>0`, does not have this invariance, since adding a constant before raising to a power is not a monotone rescaling of the dot product.) If your numbers come out this way, it's worth reporting as a real finding about this kernel/feature combination, not something to tune away." ] }, { "cell_type": "code", "execution_count": null, "id": "8a360f71", "metadata": {}, "outputs": [], "source": [ "# (B) Fixed window -- Spambase with Euclidean distance\n", "scratch_window_spam_train, scratch_window_spam_test = [], []\n", "for r in spam_radii:\n", " knn = KNNClassifier(radius=r, metric=\"euclidean\").fit(X_train_spam, y_train_spam)\n", " acc_train = accuracy_score(y_train_spam, knn.predict(X_train_spam))\n", " acc_test = accuracy_score(y_test_spam, knn.predict(X_test_spam))\n", " scratch_window_spam_train.append(acc_train)\n", " scratch_window_spam_test.append(acc_test)\n", " print(f\"[scratch] Spambase window (euclidean, r={r}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5b3340b7", "metadata": {}, "outputs": [], "source": [ "# (B) Fixed window -- MNIST with cosine distance\n", "scratch_window_digits_train, scratch_window_digits_test = [], []\n", "for r in digits_radii:\n", " knn = KNNClassifier(radius=r, metric=\"cosine\").fit(X_train_digits, y_train_digits)\n", " acc_train = accuracy_score(y_train_digits, knn.predict(X_train_digits))\n", " acc_test = accuracy_score(y_test_digits, knn.predict(X_test_digits))\n", " scratch_window_digits_train.append(acc_train)\n", " scratch_window_digits_test.append(acc_test)\n", " print(f\"[scratch] MNIST window (cosine, r={r}): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "20aa59e4", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "\n", "axes[0].plot(spam_radii, lib_window_spam_train, color=\"#2a78d6\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library train\")\n", "axes[0].plot(spam_radii, lib_window_spam_test, color=\"#e34948\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library test\")\n", "axes[0].plot(spam_radii, scratch_window_spam_train, color=\"#2a78d6\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch train\")\n", "axes[0].plot(spam_radii, scratch_window_spam_test, color=\"#e34948\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch test\")\n", "axes[0].set_xlabel(\"Radius (Euclidean)\")\n", "axes[0].set_ylabel(\"Accuracy\")\n", "axes[0].set_title(\"Spambase fixed-window\")\n", "axes[0].legend()\n", "\n", "axes[1].plot(digits_radii, lib_window_digits_train, color=\"#2a78d6\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library train\")\n", "axes[1].plot(digits_radii, lib_window_digits_test, color=\"#e34948\", linestyle=\"--\", marker=\"o\", markersize=5, label=\"library test\")\n", "axes[1].plot(digits_radii, scratch_window_digits_train, color=\"#2a78d6\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch train\")\n", "axes[1].plot(digits_radii, scratch_window_digits_test, color=\"#e34948\", linestyle=\"-\", marker=\"o\", markersize=5, label=\"scratch test\")\n", "axes[1].set_xlabel(\"Radius (cosine)\")\n", "axes[1].set_ylabel(\"Accuracy\")\n", "axes[1].set_title(\"MNIST fixed-window\")\n", "axes[1].legend()\n", "\n", "fig.suptitle(\"Accuracy vs. radius: library vs. from scratch\")\n", "fig.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "55969d0d", "metadata": {}, "source": [ "### (C) Kernel density estimation [optional, no credit]\n", "\n", "Fit a per-class kernel density $P(z \\mid C)$ using only that class's training points, then classify by $P(C \\mid z) \\propto P(C) \\cdot P(z \\mid C)$." ] }, { "cell_type": "code", "execution_count": null, "id": "a7ee3a0a", "metadata": {}, "outputs": [], "source": [ "from sklearn.naive_bayes import GaussianNB\n", "from sklearn.neighbors import KernelDensity\n", "\n", "# library baseline: per-class sklearn KernelDensity, combined with class priors by hand\n", "# (this *is* what GaussianNB does, except with a parametric Gaussian instead of a KDE --\n", "# included here too since it's the library's closest ready-made \"generative density\" classifier)\n", "gnb = GaussianNB... ### TODO_STUDENT ### fit sklearn's GaussianNB as a reference parametric-density classifier\n", "acc_train = accuracy_score(y_train_spam, gnb.predict(X_train_spam))\n", "acc_test = accuracy_score(y_test_spam, gnb.predict(X_test_spam))\n", "print(f\"[library] Spambase GaussianNB (parametric density, for reference): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n", "\n", "bandwidth = np.... ### TODO_STUDENT ### convert the RBF kernel's gamma into the equivalent Gaussian bandwidth, so this library KDE matches our scratch kernel\n", "classes = np.unique(y_train_spam)\n", "class_kdes = {c: KernelDensity... ### TODO_STUDENT ### fit one sklearn KernelDensity estimator per class on that class's training points\n", " for c in classes}\n", "priors = {c: np.mean(y_train_spam == c) for c in classes}\n", "\n", "def kde_predict(X):\n", " log_scores = np.zeros((X.shape[0], len(classes)))\n", " for j, c in enumerate(classes):\n", " log_scores[:, j] = np.... ### TODO_STUDENT ### combine the class prior and the per-class KDE log-density into a log-posterior score\n", " return ... ### TODO_STUDENT ### predict the class with the highest posterior score\n", "\n", "acc_train = accuracy_score(y_train_spam, kde_predict(X_train_spam))\n", "acc_test = accuracy_score(y_test_spam, kde_predict(X_test_spam))\n", "print(f\"[library] Spambase KernelDensity (Gaussian kernel): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "42be8866", "metadata": {}, "outputs": [], "source": [ "class KDEClassifier:\n", " \"\"\"Naive-Bayes-style classifier: per-class kernel density estimate, Gaussian kernel.\"\"\"\n", "\n", " def __init__(self, gamma=1.0):\n", " self.gamma = gamma\n", "\n", " def fit(self, X, y):\n", " self.classes = np.unique(y)\n", " self.class_data = ... ### TODO_STUDENT ### partition the training data by class, so each class's density estimate only sees its own points\n", " self.priors = ... ### TODO_STUDENT ### estimate each class's prior P(C) as its frequency in the training data\n", " return self\n", "\n", " def predict(self, X_test):\n", " log_scores = np.zeros((X_test.shape[0], len(self.classes)))\n", " for j, c in enumerate(self.classes):\n", " # average Gaussian-kernel similarity to this class's training points == density estimate\n", " density = ... ### TODO_STUDENT ### estimate the class-conditional density as the average Gaussian-kernel similarity to that class's training points\n", " log_scores[:, j] = np.... ### TODO_STUDENT ### combine the class prior and the estimated density into a log-posterior score\n", " return ... ### TODO_STUDENT ### predict the class with the highest posterior score\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d324ac26", "metadata": {}, "outputs": [], "source": [ "kde = KDEClassifier(gamma=1e-3).fit(X_train_spam, y_train_spam)\n", "acc_train = accuracy_score(y_train_spam, kde.predict(X_train_spam))\n", "acc_test = accuracy_score(y_test_spam, kde.predict(X_test_spam))\n", "print(f\"[scratch] Spambase KDE (Gaussian kernel): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "cb202790", "metadata": {}, "source": [ "---\n", "## Problem 5 (optional, substantial) — Second-order (Newton) Boosting, from the XGBoost paper\n", "\n", "Problem 3's boosting fits trees to raw residuals -- a *first-order* method (it only uses the\n", "gradient of the loss). XGBoost (Chen & Guestrin, 2016) generalizes this to a *second-order*\n", "method: it also uses the loss's Hessian, which lets the regularized **structure score** decide\n", "splits and gives a closed-form optimal weight for each leaf. We redo Problem 3's boosting as\n", "**classification** on Spambase with log-loss, where the gradient/Hessian are no longer trivial\n", "(unlike squared loss, where the Hessian is just 1).\n", "\n", "### (a) Library baseline: `sklearn.ensemble.GradientBoostingClassifier`\n", "\n", "`xgboost` itself isn't in the course environment, so we use sklearn's own gradient-boosted\n", "classifier (`loss=\"log_loss\"`) as the reference -- it isn't *exactly* the XGBoost algorithm (no\n", "per-leaf L2 term, no explicit `gamma`), but it optimizes the same loss with the same \"trees on\n", "pseudo-residuals\" family, so its train/test log-loss curve is the right shape to check ours\n", "against.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3ab1c5fc", "metadata": {}, "outputs": [], "source": [ "from sklearn.ensemble import GradientBoostingClassifier\n", "from sklearn.metrics import log_loss\n", "\n", "n_rounds5 = 60\n", "sk_gbc = GradientBoostingClassifier(\n", " loss=\"log_loss\", n_estimators=n_rounds5, max_depth=3, learning_rate=0.3, random_state=0\n", ")... ### TODO_STUDENT ### fit sklearn's GradientBoostingClassifier(loss='log_loss') as the library baseline\n", "\n", "lib_train_logloss = [log_loss(y_train_spam, p[:, 1])\n", " for p in sk_gbc.staged_predict_proba(X_train_spam)]\n", "lib_test_logloss = [log_loss(y_test_spam, p[:, 1])\n", " for p in sk_gbc.staged_predict_proba(X_test_spam)]\n", "\n", "print(f\"[library] Final (T={n_rounds5}): train log-loss = {lib_train_logloss[-1]:.3f}, \"\n", " f\"test log-loss = {lib_test_logloss[-1]:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "a5ec86e5", "metadata": {}, "source": [ "### (b) From scratch: gradients, Hessians, and the regularized structure score\n", "\n", "For log-loss with prediction `p = sigmoid(F(x))` and label `y`, the per-point gradient and Hessian\n", "(w.r.t. the raw score `F`, not `p`) are:\n", "\n", "$$g_i = p_i - y_i \\qquad h_i = p_i(1-p_i)$$\n", "\n", "A candidate split on node data with gradient/Hessian sums `G, H` (and `G_L,H_L` / `G_R,H_R` for\n", "the two children) is scored by the **structure score gain**:\n", "\n", "$$\\text{gain} = \\frac{G_L^2}{H_L+\\lambda} + \\frac{G_R^2}{H_R+\\lambda} - \\frac{G^2}{H+\\lambda} - \\gamma$$\n", "\n", "and a leaf's optimal weight, for fixed tree structure, is the closed form `w* = -G/(H+lambda)`.\n", "`lambda` shrinks every leaf weight toward 0 (L2 regularization); `gamma` is a flat cost charged\n", "per split, so a candidate split is only taken if it improves the score by more than `gamma`. We\n", "reuse `candidate_thresholds` from Problem 2 for the split search.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5ce1c3fc", "metadata": {}, "outputs": [], "source": [ "def sigmoid(z):\n", " return ... ### TODO_STUDENT ### implement the sigmoid function that maps a raw score to a probability\n", "\n", "\n", "class XGBoostStyleTree:\n", " \"\"\"A regression tree whose splits are scored by the XGBoost structure score (gain), using\n", " per-point gradients/Hessians instead of variance or entropy.\"\"\"\n", "\n", " def __init__(self, max_depth=3, min_samples_split=2, lmbda=1.0, gamma=0.0):\n", " self.max_depth = max_depth\n", " self.min_samples_split = min_samples_split\n", " self.lmbda = lmbda\n", " self.gamma = gamma\n", "\n", " def fit(self, X, g, h):\n", " self.root = self._grow(X, g, h, depth=0)\n", " return self\n", "\n", " def _leaf_weight(self, g, h):\n", " return ... ### TODO_STUDENT ### compute the closed-form optimal leaf weight w* = -G/(H+lambda)\n", "\n", " def _grow(self, X, g, h, depth):\n", " n, m = X.shape\n", " G, H = g.sum(), h.sum()\n", "\n", " if ...: ### TODO_STUDENT ### decide when to stop splitting and turn this node into a leaf\n", " return TreeNode(value=self._leaf_weight(g, h))\n", "\n", " parent_score = ... ### TODO_STUDENT ### compute the parent node's structure score G^2/(H+lambda)\n", " best_gain, best_feature, best_threshold = 0.0, None, None\n", "\n", " for feature in range(m):\n", " for threshold in candidate_thresholds(X[:, feature]):\n", " left_mask = X[:, feature] <= threshold\n", " right_mask = ~left_mask\n", " if left_mask.sum() == 0 or right_mask.sum() == 0:\n", " continue\n", "\n", " G_L, H_L = ... ### TODO_STUDENT ### sum the gradients/Hessians of the points that would go left\n", " G_R, H_R = ... ### TODO_STUDENT ### sum the gradients/Hessians of the points that would go right\n", "\n", " gain = ... ### TODO_STUDENT ### compute the structure-score gain of this split, charging the flat cost gamma per split\n", "\n", " if gain > best_gain:\n", " best_gain, best_feature, best_threshold = ... ### TODO_STUDENT ### keep track of the best (feature, threshold) split seen so far\n", "\n", " if best_feature is None:\n", " return TreeNode(value=self._leaf_weight(g, h))\n", "\n", " left_mask = X[:, best_feature] <= best_threshold\n", " left = self._grow... ### TODO_STUDENT ### recursively grow the left subtree on the left split's gradients/Hessians\n", " right = self._grow... ### TODO_STUDENT ### recursively grow the right subtree on the right split's gradients/Hessians\n", " return TreeNode(feature=best_feature, threshold=best_threshold, left=left, right=right)\n", "\n", " def _predict_one(self, x, node):\n", " if node.is_leaf():\n", " return node.value\n", " branch = node.left if x[node.feature] <= node.threshold else node.right\n", " return self._predict_one(x, branch)\n", "\n", " def predict(self, X):\n", " return np.array([self._predict_one(x, self.root) for x in X])\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d9305b09", "metadata": {}, "outputs": [], "source": [ "class SecondOrderBoostedClassifier:\n", " \"\"\"Binary classification boosting on log-loss using XGBoost-style Newton steps: each round\n", " fits an XGBoostStyleTree to the current gradients/Hessians instead of raw residuals.\"\"\"\n", "\n", " def __init__(self, n_rounds=60, tree_depth=3, learning_rate=0.3, lmbda=1.0, gamma=0.0):\n", " self.n_rounds = n_rounds\n", " self.tree_depth = tree_depth\n", " self.learning_rate = learning_rate\n", " self.lmbda = lmbda\n", " self.gamma = gamma\n", "\n", " def fit(self, X, y):\n", " self.trees = []\n", " F = np.zeros(X.shape[0])\n", "\n", " for _ in range(self.n_rounds):\n", " p = sigmoid(F)\n", " g = ... ### TODO_STUDENT ### compute the log-loss gradient g_i = p_i - y_i\n", " h = ... ### TODO_STUDENT ### compute the log-loss Hessian h_i = p_i(1-p_i)\n", " tree = XGBoostStyleTree(max_depth=self.tree_depth, lmbda=self.lmbda, gamma=self.gamma)\n", " tree.fit(X, g, h)\n", " self.trees.append(tree)\n", " F = ... ### TODO_STUDENT ### take a Newton step: add this round's (learning-rate-scaled) tree output to the running score\n", " return self\n", "\n", " def staged_decision_function(self, X):\n", " F = np.zeros(X.shape[0])\n", " for tree in self.trees:\n", " F = ... ### TODO_STUDENT ### accumulate this round's Newton step into the running score, for the staged log-loss curve\n", " yield F\n", "\n", " def predict_proba(self, X):\n", " F = np.zeros(X.shape[0])\n", " for tree in self.trees:\n", " F = F + self.learning_rate * tree.predict(X)\n", " return sigmoid(F)\n", "\n", " def predict(self, X, threshold=0.5):\n", " return (self.predict_proba(X) >= threshold).astype(int)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7c5cfec9", "metadata": {}, "outputs": [], "source": [ "second_order = SecondOrderBoostedClassifier(n_rounds=n_rounds5, tree_depth=3, learning_rate=0.3,\n", " lmbda=1.0, gamma=0.0).fit(X_train_spam, y_train_spam)\n", "\n", "scratch_train_logloss = [log_loss(y_train_spam, sigmoid(F))\n", " for F in second_order.staged_decision_function(X_train_spam)]\n", "scratch_test_logloss = [log_loss(y_test_spam, sigmoid(F))\n", " for F in second_order.staged_decision_function(X_test_spam)]\n", "\n", "print(f\"[scratch] Final (T={n_rounds5}): train log-loss = {scratch_train_logloss[-1]:.3f}, \"\n", " f\"test log-loss = {scratch_test_logloss[-1]:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0fd6560b", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(6, 4))\n", "rounds5 = range(1, n_rounds5 + 1)\n", "plt.plot(rounds5, lib_train_logloss, color=\"#2a78d6\", linestyle=\"--\", label=\"library train\")\n", "plt.plot(rounds5, lib_test_logloss, color=\"#e34948\", linestyle=\"--\", label=\"library test\")\n", "plt.plot(rounds5, scratch_train_logloss, color=\"#2a78d6\", linestyle=\"-\", label=\"scratch train\")\n", "plt.plot(rounds5, scratch_test_logloss, color=\"#e34948\", linestyle=\"-\", label=\"scratch test\")\n", "plt.xlabel(\"Number of boosting rounds T\")\n", "plt.ylabel(\"Log-loss\")\n", "plt.title(\"Second-order boosting on Spambase: library (first-order) vs. scratch (second-order)\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "c628b66a", "metadata": {}, "source": [ "Both curves should show the same overfitting shape (test log-loss bottoming out then creeping\n", "back up) -- that's the sanity check to look for: your Newton-step trees are optimizing the same\n", "loss the library's first-order boosting is, just with an extra (Hessian-aware) way of scoring\n", "splits and weighting leaves, so the qualitative curve shape should match even though the two\n", "methods differ.\n", "\n", "### Reference notes for the no-LLM write-up (grading aid, not a student answer)\n", "\n", "- **What second-order gets you over Problem 3:** Problem 3's boosting only knows the residual's\n", " *sign and size*; it can't tell \"confidently wrong\" from \"barely wrong\" the way a Hessian can.\n", " For log-loss, `h_i = p_i(1-p_i)` is largest when the model is most unsure (`p≈0.5`) and shrinks\n", " toward 0 as the model becomes confident -- so Newton steps are naturally smaller/more cautious\n", " exactly where the loss surface is curving fastest, which first-order gradient boosting can't see.\n", "- **`g_i`, `h_i` for log-loss:** `g_i = p_i - y_i` (residual on the probability scale), `h_i =\n", " p_i(1-p_i)` (bounded in `[0, 0.25]`, i.e. a per-point \"how flat is the loss here\" curvature term).\n", "- **Why `/(H+lambda)` regularizes:** it's a Newton step `-G/H` (like Newton's method for\n", " optimization) shrunk by adding `lambda` to the denominator -- exactly like ridge regression adds\n", " `lambda` to `XtX` before inverting. Small `H` (an uncertain/sparse leaf) gets shrunk hardest.\n", " - **Why `w*=-G/(H+lambda)` is optimal:** it's the minimizer of the node's *quadratic* (2nd-order\n", " Taylor) approximation to the loss plus an L2 penalty `0.5*lambda*w^2` on the leaf weight --\n", " literally solving `d/dw [G*w + 0.5*(H+lambda)*w^2] = 0`.\n", "- **`lambda` vs `gamma`:** `lambda` regularizes *how big* a leaf's weight can be (shrinks toward\n", " 0); `gamma` regularizes *how many* leaves you're allowed to have (a split must clear a minimum\n", " gain bar). Set both to 0 and this degenerates toward unregularized Newton boosting -- every split\n", " that helps even a hair gets taken, and leaf weights are unshrunk, so it overfits faster than\n", " Problem 3's boosting (which had no split criterion tunable at all beyond depth).\n", "- **Where a real implementation typically first diverges from the library:** forgetting that `g`\n", " and `h` are derivatives *of the loss with respect to the raw score `F`*, not with respect to\n", " `p` -- plugging in `y - p` instead of `p - y` (sign flip) or using `1` instead of `p(1-p)` for the\n", " Hessian (i.e. accidentally reimplementing Problem 3's boosting) are the two most common bugs, and\n", " both show up as a *plausible-looking but wrong* log-loss curve rather than a crash.\n" ] }, { "cell_type": "markdown", "id": "367213a4", "metadata": {}, "source": [ "### (c) Bonus: more than two classes (library only, MNIST)\n", "\n", "Our from-scratch `XGBoostStyleTree`/`SecondOrderBoostedClassifier` assume binary log-loss.\n", "Generalizing them to multiclass (softmax loss, a gradient/Hessian *per class*, one tree per class\n", "per round) is a substantially bigger build and isn't required here -- we just check how the\n", "*library* behaves with more than two classes, on a fresh MNIST split. (This split and model are\n", "reused as-is in Problem 6's MNIST analysis.)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5096a0ba", "metadata": {}, "outputs": [], "source": [ "# Fresh random 90/10 split, stratified per class (no folds), from a subsample of the pooled\n", "# train+test Haar-feature arrays already loaded for Problem 4.\n", "mnist_pool_X = np.vstack([X_train_digits_full, X_test_digits_full])\n", "mnist_pool_y = np.concatenate([y_train_digits_full, y_test_digits_full])\n", "\n", "rng_mnist = np.random.RandomState(2)\n", "pool_idx = rng_mnist.choice(len(mnist_pool_X), size=4000, replace=False)\n", "X_mnist, y_mnist = mnist_pool_X[pool_idx], mnist_pool_y[pool_idx]\n", "\n", "X_train_mnist, X_test_mnist, y_train_mnist, y_test_mnist = train_test_split(\n", " X_mnist, y_mnist, test_size=0.10, stratify=y_mnist, random_state=2)\n", "print(\"MNIST (fresh 90/10 split) train:\", X_train_mnist.shape, \" test:\", X_test_mnist.shape)\n", "\n", "n_rounds5_mnist = 60\n", "sk_gbc_mnist = GradientBoostingClassifier(\n", " loss=\"log_loss\", n_estimators=n_rounds5_mnist, max_depth=3, learning_rate=0.3, random_state=0\n", ")... ### TODO_STUDENT ### fit sklearn's multiclass GradientBoostingClassifier(loss='log_loss') on the MNIST split\n", "\n", "mnist_train_logloss = [log_loss(y_train_mnist, p)\n", " for p in sk_gbc_mnist.staged_predict_proba(X_train_mnist)]\n", "mnist_test_logloss = [log_loss(y_test_mnist, p)\n", " for p in sk_gbc_mnist.staged_predict_proba(X_test_mnist)]\n", "acc_mnist = sk_gbc_mnist.score(X_test_mnist, y_test_mnist)\n", "\n", "print(f\"MNIST (10-class, library only): test accuracy = {acc_mnist:.3f}, \"\n", " f\"final train log-loss = {mnist_train_logloss[-1]:.3f}, \"\n", " f\"test log-loss = {mnist_test_logloss[-1]:.3f}\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f6af59d7", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(6, 4))\n", "plt.plot(range(1, n_rounds5_mnist + 1), mnist_train_logloss, color=\"#2a78d6\", label=\"train\")\n", "plt.plot(range(1, n_rounds5_mnist + 1), mnist_test_logloss, color=\"#e34948\", label=\"test\")\n", "plt.xlabel(\"Number of boosting rounds T\")\n", "plt.ylabel(\"Log-loss\")\n", "plt.title(\"Multiclass gradient boosting on MNIST (library only, 10 classes)\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "d31b0655", "metadata": {}, "source": [ "With 10 classes instead of 2, expect log-loss to start much higher (guessing uniformly among 10\n", "classes costs `ln(10) ≈ 2.30` nats vs. `ln(2) ≈ 0.69` for binary) -- but the overall *shape* of\n", "the curve (a fast initial drop, then a slow climb back up as it overfits) should tell the same\n", "qualitative story as Spambase. The multiclass case just fits *K* trees per round internally (one\n", "per class) instead of one, which is exactly the generalization our from-scratch second-order tree\n", "would need to support this properly.\n" ] }, { "cell_type": "markdown", "id": "d03691dc", "metadata": {}, "source": [ "---\n", "## Problem 6 (optional) — Does boosting confidence track KNN neighborhood consistency?\n", "\n", "**Hypothesis.** Shallow trees split the same feature space k-NN measures proximity in, so a test\n", "point sitting in a locally *pure/consistent* neighborhood should tend to land on the \"easy\" side\n", "of most trees in a boosted ensemble -- their contributions reinforce rather than cancel, giving a\n", "large-magnitude, confident prediction. In a mixed/sparse neighborhood, different trees disagree,\n", "partially cancel, and the ensemble's output stays close to \"unsure.\"\n", "\n", "We test this on **three** datasets, reusing already-fitted models and the same\n", "`pairwise_distances` machinery from Problems 3/4/5 throughout -- no new data, and (for Spambase and\n", "Housing) essentially no new training:\n", "\n", "- **Spambase (classification):** confidence = boosting margin `|p(z) - 0.5| * 2` from the already\n", " -fit `second_order` classifier (Problem 5) -- zero retraining. Neighborhood consistency = k-NN\n", " **purity** (majority-class fraction among the k nearest Spambase training points), reusing the\n", " exact `pairwise_distances` call from Problem 4.\n", "- **Housing (regression):** there's no \"decision boundary\" to be far from, so a large prediction\n", " doesn't mean \"confident\" -- it just means \"predicted to be big.\" The regression-appropriate\n", " notion of confidence is *predictive uncertainty*: we bootstrap-resample the training set, refit\n", " Problem 3's `BoostedRegressionTrees` a handful of times, and use the **std of predictions across\n", " resamples** as an uncertainty score. Neighborhood consistency = k-NN **target spread** (std of y\n", " among the k nearest Housing training points) -- the direct regression analog of \"purity.\"\n", "- **MNIST (10-class classification):** reuses the fresh, class-stratified 90/10 split (no folds)\n", " and the library `sk_gbc_mnist` classifier from Problem 5's bonus section -- no new training here\n", " either. Confidence generalizes the Spambase margin to multiclass: `p_top1(z) - p_top2(z)`, the\n", " gap between the top two predicted-class probabilities (this collapses to exactly `|p-0.5|*2` in\n", " the binary case, so it's the same quantity, not an analogy). `knn_purity` is already\n", " class-count-generic, so it needs no change.\n", "\n", "All three pairs are framed the same direction (confidence-vs-confidence for Spambase/MNIST,\n", "uncertainty-vs-uncertainty for Housing) so a *positive* correlation is the predicted sign in every\n", "case.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "45b6d47f", "metadata": {}, "outputs": [], "source": [ "from scipy.stats import pearsonr, spearmanr\n", "\n", "def knn_purity(D, y_train, k):\n", " \"\"\"Majority-class fraction among the k nearest neighbors, for each query row of D.\"\"\"\n", " purities = np.empty(D.shape[0])\n", " for i in range(D.shape[0]):\n", " nearest = np.argsort(D[i])[:k]\n", " _, counts = np.unique(y_train[nearest], return_counts=True)\n", " purities[i] = ... ### TODO_STUDENT ### compute the majority-class fraction among these k neighbors (the purity score)\n", " return purities\n", "\n", "\n", "def knn_target_spread(D, y_train, k):\n", " \"\"\"Std of the target value among the k nearest neighbors, for each query row of D.\"\"\"\n", " spreads = np.empty(D.shape[0])\n", " for i in range(D.shape[0]):\n", " nearest = np.argsort(D[i])[:k]\n", " spreads[i] = ... ### TODO_STUDENT ### compute the std of the target among these k neighbors (the spread score)\n", " return spreads\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a1660ff3", "metadata": {}, "outputs": [], "source": [ "# ---- Spambase: boosting margin vs. k-NN purity ----\n", "D_spam_test_train = pairwise_distances(X_test_spam, X_train_spam, metric=\"euclidean\")\n", "margin_spam = np.... ### TODO_STUDENT ### define the boosting confidence score as the margin |p-0.5| rescaled to [0,1]\n", "\n", "print(\"Spambase: margin vs. k-NN purity, at a few k\")\n", "for k in [5, 15, 31]:\n", " purity_spam_k = knn_purity(D_spam_test_train, y_train_spam, k)\n", " r_p, _ = pearsonr(margin_spam, purity_spam_k)\n", " r_s, _ = spearmanr(margin_spam, purity_spam_k)\n", " print(f\" k={k:>2}: Pearson r = {r_p:.3f}, Spearman rho = {r_s:.3f}\")\n", "\n", "k_confidence = 15\n", "purity_spam = knn_purity(D_spam_test_train, y_train_spam, k_confidence)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3ff56874", "metadata": {}, "outputs": [], "source": [ "# ---- Housing: bootstrap-ensemble disagreement vs. k-NN target spread ----\n", "n_bootstrap = 15\n", "rng6 = np.random.RandomState(1)\n", "bootstrap_preds = np.zeros((n_bootstrap, X_test_housing.shape[0]))\n", "\n", "for b in range(n_bootstrap):\n", " idx = rng6.... ### TODO_STUDENT ### draw a bootstrap resample (sampling with replacement) of the training set\n", " boosted_b = BoostedRegressionTrees(n_rounds=n_rounds, tree_depth=2)\n", " boosted_b.fit(X_train_housing[idx], y_train_housing[idx])\n", " bootstrap_preds[b] = boosted_b.predict(X_test_housing)\n", "\n", "boost_std_housing = bootstrap_preds.... ### TODO_STUDENT ### define the uncertainty score as the std of predictions across bootstrap refits\n", "\n", "D_house_test_train = pairwise_distances(X_test_housing, X_train_housing, metric=\"euclidean\")\n", "print(\"Housing: bootstrap-ensemble std vs. k-NN target spread, at a few k\")\n", "for k in [5, 15, 31]:\n", " spread_housing_k = knn_target_spread(D_house_test_train, y_train_housing, k)\n", " r_p, _ = pearsonr(boost_std_housing, spread_housing_k)\n", " r_s, _ = spearmanr(boost_std_housing, spread_housing_k)\n", " print(f\" k={k:>2}: Pearson r = {r_p:.3f}, Spearman rho = {r_s:.3f}\")\n", "\n", "spread_housing = knn_target_spread(D_house_test_train, y_train_housing, k_confidence)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "70fa7192", "metadata": {}, "outputs": [], "source": [ "# ---- MNIST (10-class): margin vs. k-NN purity ----\n", "# Reuse the fresh 90/10 split and the library GradientBoostingClassifier from Problem 5's bonus\n", "# section -- no new data, no new training.\n", "proba_mnist = sk_gbc_mnist.predict_proba(X_test_mnist)\n", "proba_sorted = np.sort(proba_mnist, axis=1)\n", "margin_mnist = ... ### TODO_STUDENT ### generalize the margin to multiclass as the gap between the top two predicted-class probabilities\n", "\n", "D_mnist_test_train = pairwise_distances(X_test_mnist, X_train_mnist, metric=\"euclidean\")\n", "print(\"MNIST: margin vs. k-NN purity, at a few k\")\n", "for k in [5, 15, 31]:\n", " purity_mnist_k = knn_purity(D_mnist_test_train, y_train_mnist, k)\n", " r_p, _ = pearsonr(margin_mnist, purity_mnist_k)\n", " r_s, _ = spearmanr(margin_mnist, purity_mnist_k)\n", " print(f\" k={k:>2}: Pearson r = {r_p:.3f}, Spearman rho = {r_s:.3f}\")\n", "\n", "purity_mnist = knn_purity(D_mnist_test_train, y_train_mnist, k_confidence)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "eef3fdb7", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n", "\n", "axes[0].scatter(purity_spam, margin_spam, s=12, alpha=0.4, color=\"#2a78d6\")\n", "trend = np.polyfit(purity_spam, margin_spam, 1)\n", "xs = np.linspace(purity_spam.min(), purity_spam.max(), 50)\n", "axes[0].plot(xs, np.polyval(trend, xs), color=\"#e34948\", linewidth=2)\n", "axes[0].set_xlabel(f\"{k_confidence}-NN purity\")\n", "axes[0].set_ylabel(\"Boosting margin |p - 0.5| * 2\")\n", "axes[0].set_title(\"Spambase: margin vs. purity\")\n", "\n", "axes[1].scatter(spread_housing, boost_std_housing, s=12, alpha=0.4, color=\"#2a78d6\")\n", "trend = np.polyfit(spread_housing, boost_std_housing, 1)\n", "xs = np.linspace(spread_housing.min(), spread_housing.max(), 50)\n", "axes[1].plot(xs, np.polyval(trend, xs), color=\"#e34948\", linewidth=2)\n", "axes[1].set_xlabel(f\"{k_confidence}-NN target spread (std of y)\")\n", "axes[1].set_ylabel(\"Bootstrap-ensemble prediction std\")\n", "axes[1].set_title(\"Housing: uncertainty vs. spread\")\n", "\n", "axes[2].scatter(purity_mnist, margin_mnist, s=12, alpha=0.4, color=\"#2a78d6\")\n", "trend = np.polyfit(purity_mnist, margin_mnist, 1)\n", "xs = np.linspace(purity_mnist.min(), purity_mnist.max(), 50)\n", "axes[2].plot(xs, np.polyval(trend, xs), color=\"#e34948\", linewidth=2)\n", "axes[2].set_xlabel(f\"{k_confidence}-NN purity\")\n", "axes[2].set_ylabel(\"Margin p_top1 - p_top2\")\n", "axes[2].set_title(\"MNIST: margin vs. purity\")\n", "\n", "fig.suptitle(\"Boosting confidence/uncertainty vs. KNN neighborhood consistency\")\n", "fig.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "9a107ef8", "metadata": {}, "source": [ "### Does neighborhood-based confidence also predict correctness?\n", "\n", "The more useful claim isn't just \"margin correlates with purity\" -- it's \"purity (a\n", "model-free, geometric quantity) predicts where the boosted model is likely to be *wrong*,\" which\n", "is the practically interesting payoff (e.g. for flagging low-trust predictions without needing the\n", "model at all).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f1afbbf1", "metadata": {}, "outputs": [], "source": [ "pred_spam = second_order.predict(X_test_spam)\n", "correct_spam = (pred_spam == y_test_spam)\n", "\n", "# equal-*frequency* bins by rank, not equal-width value bins: purity_spam is highly discrete\n", "# (multiples of 1/k_confidence), so value-based quantile edges can collide and leave an empty bin.\n", "n_bins = 5\n", "order = np.argsort(purity_spam)\n", "bin_of_rank = np.... ### TODO_STUDENT ### split test points into equal-frequency (not equal-width) bins by purity rank, so ties/duplicates can't create an empty bin\n", "\n", "print(\"Spambase test accuracy and mean margin, by k-NN-purity quintile (low purity -> high purity):\")\n", "acc_by_bin, margin_by_bin, range_by_bin = [], [], []\n", "for ranks in bin_of_rank:\n", " idx = order[ranks]\n", " acc_by_bin.append(correct_spam[idx].mean())\n", " margin_by_bin.append(margin_spam[idx].mean())\n", " range_by_bin.append(f\"[{purity_spam[idx].min():.2f}, {purity_spam[idx].max():.2f}]\")\n", " print(f\" purity {range_by_bin[-1]}: accuracy = {acc_by_bin[-1]:.3f}, mean margin = {margin_by_bin[-1]:.3f}\")\n", "\n", "acc_by_bin = np.array(acc_by_bin)\n", "margin_by_bin = np.array(margin_by_bin)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "80ce865f", "metadata": {}, "outputs": [], "source": [ "pred_mnist = sk_gbc_mnist.predict(X_test_mnist)\n", "correct_mnist = (pred_mnist == y_test_mnist)\n", "\n", "order_m = np.argsort(purity_mnist)\n", "bin_of_rank_m = np.array_split(np.arange(len(purity_mnist)), n_bins)\n", "\n", "print(\"MNIST test accuracy and mean margin, by k-NN-purity quintile (low purity -> high purity):\")\n", "acc_by_bin_m, margin_by_bin_m, range_by_bin_m = [], [], []\n", "for ranks in bin_of_rank_m:\n", " idx = order_m[ranks]\n", " acc_by_bin_m.append(correct_mnist[idx].mean())\n", " margin_by_bin_m.append(margin_mnist[idx].mean())\n", " range_by_bin_m.append(f\"[{purity_mnist[idx].min():.2f}, {purity_mnist[idx].max():.2f}]\")\n", " print(f\" purity {range_by_bin_m[-1]}: accuracy = {acc_by_bin_m[-1]:.3f}, mean margin = {margin_by_bin_m[-1]:.3f}\")\n", "\n", "acc_by_bin_m = np.array(acc_by_bin_m)\n", "margin_by_bin_m = np.array(margin_by_bin_m)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "af7c32e7", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "\n", "axes[0].bar(range_by_bin, acc_by_bin, color=\"#2a78d6\")\n", "axes[0].axhline(correct_spam.mean(), color=\"#e34948\", linestyle=\"--\", label=\"overall test accuracy\")\n", "axes[0].set_xlabel(f\"{k_confidence}-NN purity bin (low -> high)\")\n", "axes[0].set_ylabel(\"Test accuracy\")\n", "axes[0].set_title(\"Spambase: accuracy by neighborhood-purity bin\")\n", "axes[0].tick_params(axis=\"x\", rotation=20)\n", "axes[0].legend()\n", "\n", "axes[1].bar(range_by_bin_m, acc_by_bin_m, color=\"#2a78d6\")\n", "axes[1].axhline(correct_mnist.mean(), color=\"#e34948\", linestyle=\"--\", label=\"overall test accuracy\")\n", "axes[1].set_xlabel(f\"{k_confidence}-NN purity bin (low -> high)\")\n", "axes[1].set_ylabel(\"Test accuracy\")\n", "axes[1].set_title(\"MNIST: accuracy by neighborhood-purity bin\")\n", "axes[1].tick_params(axis=\"x\", rotation=20)\n", "axes[1].legend()\n", "\n", "fig.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "7b3c9763", "metadata": {}, "source": [ "### Questions to answer from your three correlation/scatter plots\n", "\n", "This is the actual point of Problem 6 -- work these out from *your* numbers and plots, not from\n", "a pre-written answer key:\n", "\n", "- **Does the hypothesis hold equally on all three datasets, or does it hold more strongly on\n", " some than others?** Compare Pearson r and Spearman rho for Spambase, MNIST, and Housing side by\n", " side. Do the two classification datasets (binary vs. 10-class) behave similarly to each other,\n", " or does class count itself seem to matter? (Think about how many *ways* a neighborhood can be\n", " \"impure\" as the number of classes grows.)\n", "- **Do Pearson and Spearman ever disagree noticeably for the same dataset?** If so, look at the\n", " scatter plot for that case -- is there a single outlier point dragging the linear correlation\n", " around while the rank correlation barely moves? This is exactly why the assignment asks for\n", " both, not just one.\n", "- **Read the payoff plot (test accuracy/error by k-NN-purity quintile).** Is there a real,\n", " monotonic trend from the least- to most-consistent neighborhoods? How big is the swing between\n", " the extremes, and does that swing size track the correlation strength you computed above?\n", "- **Housing is a different kind of confidence signal (bootstrap std, not a margin) -- does it\n", " behave differently from the two classification datasets?** If the relationship looks weaker\n", " there, is that more likely a property of *this specific dataset* (small test set, noisy\n", " estimates) or of *bootstrap-resampling variance as a confidence proxy* in general? What would\n", " you check to tell those two explanations apart?\n", "- **Bottom line:** state, in your own words and citing your own numbers, whether the hypothesis\n", " from the top of this problem held up -- and if it held up more on some datasets than others,\n", " what you think is actually driving that difference.\n" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 5 }