Classical ML From Scratch Lab

Using a library call and understanding the method are different skills. This lab has eight exercises. Each is a short function you write yourself in numpy, then check with hidden tests that run in your browser. Every exercise starts with a plain-words goal and a small hand calculation, so you know the right number before you type a line.

Recap of earlier lessons. Linear & Polynomial Regression and Logistic Regression give the first two models. K-Nearest Neighbors and Clustering give the neighbour vote and k-means. Decision Trees gives the Gini split. Dimensionality Reduction gives PCA. Naive Bayes and Model Evaluation give the last two. This lab does not teach those ideas again. It makes you build them.

Learning Objectives

After this lesson, you will be able to:

Fit a line with the normal equation and with gradient descent, and show that both reach the same weights
Write a sigmoid that cannot overflow, the log loss, and a gradient step for logistic regression
Label a new row by a vote of its k nearest rows, with a stated rule for ties
Write the two alternating steps of k-means, and handle a cluster that owns no rows
Score a split with Gini impurity and find the best single question over all features
Compute principal components from a covariance matrix, and write Gaussian naive Bayes in log space
Split rows into shuffled folds so that no row is ever in training and validation together

Prerequisites

Before You Start

Review

#How This Lab Works

Every exercise has the same shape. First comes a short goal in plain words. Then a hand trace: one small calculation done on paper, with the answer your code should give. Then a code card. The card has a starter with TODO marks, a ladder of hints that opens one step at a time, a Run button, a Check my work button, and a way to reveal the solution.

Check my work runs your code, then runs a list of hidden tests one at a time. It tells you how many passed and shows the message of the first one that failed. The starter is built to fail on purpose, so a fresh card always says some checks failed. When you finish the TODO lines, all of them pass.

Three things about the tests are worth knowing before you start.

They are tolerance based. A result close enough to the right answer passes, so a harmless rounding difference never fails you.
They are deterministic. Any random numbers come from a fixed seed, so a test that fails today fails the same way tomorrow.
They check more than the hand trace. Some compare your answer with a second route to the same number, such as numpy's own least-squares routine, an independent Newton solver, a slow loop, or a brute-force search. If your code only matches the trace, it will still fail those.

The cards use numpy only. The browser has no scikit-learn, so none of this needs an install, and each card runs in well under a second.

Words you need9 terms

The numbers a model multiplies the features by. In a line, the intercept and the slopes. Learning means finding good weights.

The direction in which the error rises fastest. Gradient descent steps the opposite way, a little at a time.

The cost of a probability. A confident wrong answer costs a lot, a confident right answer costs almost nothing.

The label that most of the chosen neighbours carry. A tie needs a stated rule.

The average position of a group of rows. K-means keeps one centre per cluster.

How mixed a group of labels is: 0 for a single label, larger when the labels are evenly mixed.

A table of how every pair of features moves together. Its diagonal holds each feature's variance.

A direction that a matrix only stretches, without turning. For a covariance matrix these are the directions of spread.

One slice of the rows. Cross-validation holds each fold out in turn and trains on the rest.

#The Conventions Used Here

Most mistakes in this area come from which of two conventions a sentence uses. Fix them before writing code. This table is the single source for the whole lab, and the tests assume it.

Quantity	Meaning in this lab
X	A numpy array of shape (n, d): n rows, d features. Row i is X[i].
y	A numpy array of n labels or targets. Class labels are whole numbers 0, 1, 2, and so on.
Intercept	Added as a column of ones in front of X. The first weight is the intercept, then one weight per feature.
Distance	Straight-line (Euclidean) distance. Squared distances give the same order, so the square root is skipped.
Ties	In a vote, the smaller label wins. In an argmin, the lower index wins. In a split, the smaller threshold wins.
Variance	Divide by n for the class variances of naive Bayes. Divide by n minus 1 for the covariance matrix in PCA.
Random numbers	`np.random.default_rng(seed)`, always with a seed, so every run is the same.
Logs	Natural logarithm, `np.log`.

Two notes. First, "divide by n or n minus 1" is a real choice that changes numbers in the fourth decimal place, and the two exercises that need it say which one. Second, each card is independent. It carries its own copy of any helper it needs, so you can do the exercises in any order.

#Exercise 1: Linear Regression Two Ways

In plain words. A line with several features predicts y as a weighted sum of the features plus an intercept. The best weights are the ones with the smallest average squared miss. There are two ways to find them. The normal equation solves for them in one go: set the slope of the error to zero and solve the resulting equations. Gradient descent starts at all zeros and takes many small steps downhill. Both must end at the same weights.

You write two functions. With the intercept column added, called Xb, the normal equation is (Xb^T Xb) w = Xb^T y. One gradient step is w = w - lr * (2 / n) * Xb^T (Xb w - y).

Hand trace. Take the points (1, 3), (2, 5) and (3, 7), which lie on y = 1 + 2x. Then Xb^T Xb is the table [[3, 6], [6, 14]] and Xb^T y is [15, 34]. Solving gives the weights [1, 2]. For gradient descent from zeros with lr = 0.1, the first step is: the misses Xb w - y are [-3, -5, -7], so Xb^T of them is [-15, -34], times 2/3 gives [-10, -22.67], and w becomes [1.0, 2.2667]. The first step overshoots the slope, and later steps settle it.

Mini project

Linear regression: normal equation and gradient descent

~7 min

Write fit_normal_equation(X, y) and fit_gradient_descent(X, y, lr, steps). Both return the weights with the intercept first. The first solves the normal equation. The second starts from all zeros and takes steps of lr times the gradient of the mean squared error, with no change to X or y.

Expected output

normal equation: [1. 2.]
one gradient step: [1.     2.2667]

Your codeEdits save in this browser

pythonrunnable cell

import numpy as np

def add_bias(X):
  """Put a column of ones in front of X, so the first weight is the intercept."""
  X = np.asarray(X, dtype=float)
  return np.hstack([np.ones((X.shape[0], 1)), X])

def fit_normal_equation(X, y):
  """Return the weights w (intercept first) that give the smallest squared error."""
  Xb = add_bias(X)
  # TODO: solve (Xb^T Xb) w = Xb^T y  (np.linalg.solve is safer than an explicit inverse)
  return np.zeros(Xb.shape[1])

def fit_gradient_descent(X, y, lr=0.1, steps=1000):
  """The same weights, found by many small downhill steps from all zeros."""
  Xb = add_bias(X)
  n = Xb.shape[0]
  w = np.zeros(Xb.shape[1])
  for _ in range(steps):
      # TODO: gradient = (2 / n) * Xb^T (Xb w - y), then w = w - lr * gradient
      pass
  return w

X = np.array([[1.0], [2.0], [3.0]])
y = np.array([3.0, 5.0, 7.0])
print("normal equation:", np.round(fit_normal_equation(X, y), 4))
print("one gradient step:", np.round(fit_gradient_descent(X, y, lr=0.1, steps=1), 4))

Two tests deserve a note. The noisy-data test compares your normal-equation weights with numpy's own least-squares routine, a second route to the same answer. The convergence test needs gradient descent and the normal equation to agree to four decimal places, which ties the two halves of the exercise together.

#Exercise 2: Logistic Regression

In plain words. Logistic regression predicts the chance that a row belongs to class 1. It forms the same weighted sum as a line, then squeezes it through the sigmoid so the answer lies between 0 and 1. It is trained on log loss, which punishes confident wrong answers hard. The gradient of the average log loss is Xb^T (p - y) / n, where p holds the predicted chances. That is almost the linear-regression gradient, without the factor 2 and with p in place of the raw weighted sum.

The sigmoid has a trap. The formula 1 / (1 + exp(-z)) overflows for a large negative z, because exp(1000) is infinite. The fix is to split the cases: for z at or above 0 use that formula, and for z below 0 use exp(z) / (1 + exp(z)).

Hand trace. With z = 2, the sigmoid is 1 / (1 + e^-2) = 0.8808. For log loss, true labels [1, 0] and chances [0.8, 0.3] cost -(ln 0.8 + ln 0.7) / 2 = (0.2231 + 0.3567) / 2 = 0.2899. For one step from zeros, take the rows x = 0, 1, 2, 3 with labels 0, 0, 1, 1. Every p is 0.5, so p - y is [0.5, 0.5, -0.5, -0.5]. Then Xb^T (p - y) / 4 is [0, -0.5], and with lr = 0.5 the new weights are [0, 0.25].

Mini project

Logistic regression from scratch

~8 min

Write sigmoid(z) that never overflows, log_loss(y, p) that clips p away from 0 and 1, and fit_logistic(X, y, lr, steps) that runs gradient descent on the log loss from all zeros, intercept first.

Expected output

sigmoid(0), sigmoid(2): 0.5 0.8808
one step: [0.   0.25]

Your codeEdits save in this browser

pythonrunnable cell

import numpy as np

def add_bias(X):
  """Put a column of ones in front of X (the first weight is the intercept)."""
  X = np.asarray(X, dtype=float)
  return np.hstack([np.ones((X.shape[0], 1)), X])

def sigmoid(z):
  """Squash any number into (0, 1). Must not overflow for very large or very small z."""
  z = np.asarray(z, dtype=float)
  # TODO: 1 / (1 + exp(-z)) for z >= 0, and exp(z) / (1 + exp(z)) for z < 0
  return np.full(z.shape, 0.5)

def log_loss(y, p, eps=1e-12):
  """Average cost of the probabilities p when the true labels are y (0 or 1)."""
  y = np.asarray(y, dtype=float)
  p = np.clip(np.asarray(p, dtype=float), eps, 1.0 - eps)
  # TODO: -mean( y * log(p) + (1 - y) * log(1 - p) )
  return 0.0

def fit_logistic(X, y, lr=0.5, steps=2000):
  """Gradient descent on the log loss, starting from all zeros. Intercept first."""
  Xb = add_bias(X)
  y = np.asarray(y, dtype=float)
  n = Xb.shape[0]
  w = np.zeros(Xb.shape[1])
  for _ in range(steps):
      p = sigmoid(Xb @ w)
      # TODO: gradient = Xb^T (p - y) / n, then w = w - lr * gradient
      pass
  return w

def predict_proba(X, w):
  """Probability of class 1 for each row."""
  return sigmoid(add_bias(X) @ w)

X = np.array([[0.0], [1.0], [2.0], [3.0]])
y = np.array([0, 0, 1, 1])
print("sigmoid(0), sigmoid(2):", round(float(sigmoid(0.0)), 4), round(float(sigmoid(2.0)), 4))
print("one step:", np.round(fit_logistic(X, y, lr=0.5, steps=1), 4))

The strongest test runs a few rounds of Newton's method, a different optimiser, on the same data. Gradient descent and Newton must end at the same weights, so a sign slip or a missing intercept cannot hide. Another test flips every label and expects every weight to change sign.

#Exercise 3: K-Nearest Neighbours

In plain words. To label a new row, measure its distance to every training row, take the k closest, and let them vote. There is nothing to train: the training rows are the model. The work is in three small decisions. Which distance? Here, straight-line. Which k rows? The smallest k distances, not the largest. What if the vote is tied? Here, the smaller label wins.

For one test row, compute the squared distance to every training row, find the positions of the k smallest with np.argsort, read their labels, count them with np.bincount, and take the biggest count with np.argmax. Both np.argsort(..., kind="stable") and np.argmax keep the earlier item on a tie, which gives the rule you want.

Hand trace. The training rows are (0,0), (0,1), (1,0) with label 0 and (5,5), (5,6), (6,5) with label 1. For the test point (0.5, 0.5) the squared distances to the three 0 rows are 0.5, 0.5 and 0.5, so all three nearest are 0s, and the answer is 0. For (4, 4) the squared distances are 2 to (5,5), then 5 to (5,6) and 5 to (6,5), so the three nearest are all 1s, and the answer is 1.

Mini project

K-nearest neighbours from scratch

~6 min

Write knn_predict(X_train, y_train, X_test, k) that returns one label per test row, using straight-line distance, a majority vote of the k nearest training rows, and the smaller label on a tie.

Expected output

[0 1]

Your codeEdits save in this browser

pythonrunnable cell

One test compares your labels with a plain Python loop on 90 random rows across k = 1, 3, 5 and 7. Another builds a case where straight-line distance and the sum of the gaps (Manhattan distance) pick different nearest rows, so using the wrong distance fails it.

#Exercise 4: K-Means

In plain words. K-means groups unlabelled rows into k clusters. Start with k guessed centres. Then repeat two steps. First, give every row to its nearest centre. Second, move every centre to the average of the rows it was given. Stop when no row changes cluster. Each round can only lower or keep the inertia, the total squared distance from rows to their centres, so the loop always settles.

There is one awkward case. A centre can end up owning no rows, and the average of nothing is not a number. The rule here is that such a centre stays where it was.

Hand trace. The rows are 1, 2, 9, 10 and the starting centres are 0 and 5. Rows 1 and 2 are nearer 0 than 5, and rows 9 and 10 are nearer 5, so the labels are [0, 0, 1, 1]. The new centres are the averages, 1.5 and 9.5. Assigning again changes nothing, so the loop stops. Every row is 0.5 from its centre, so the inertia is 4 x 0.25 = 1.0.

Mini project

K-means from scratch

~8 min

Write assign_clusters(X, centers), update_centers(X, labels, centers) where an empty cluster keeps its old centre, and kmeans(X, init_centers, max_iter) that returns the final centres, the labels and the inertia. Do not change the starting centres you are given.

Expected output

[1.5 9.5] [0 0 1 1] 1.0

Your codeEdits save in this browser

pythonrunnable cell

import numpy as np

def assign_clusters(X, centers):
  """For each row of X, the index of the nearest centre (a tie goes to the lower index)."""
  X = np.asarray(X, dtype=float)
  centers = np.asarray(centers, dtype=float)
  # TODO: squared distance from every row to every centre -> an (n, k) table, then argmin along axis 1
  return np.zeros(len(X), dtype=int)

def update_centers(X, labels, centers):
  """Move each centre to the mean of the rows assigned to it. An empty cluster keeps its old centre."""
  X = np.asarray(X, dtype=float)
  new_centers = np.array(centers, dtype=float)
  for j in range(len(new_centers)):
      # TODO: if any rows have label j, set new_centers[j] to their mean
      pass
  return new_centers

def kmeans(X, init_centers, max_iter=100):
  """Alternate assign and update until the labels stop changing. Return (centers, labels, inertia)."""
  X = np.asarray(X, dtype=float)
  centers = np.array(init_centers, dtype=float)
  labels = assign_clusters(X, centers)
  for _ in range(max_iter):
      # TODO: new_centers = update_centers(...), new_labels = assign_clusters(...);
      #       store them, and stop early when new_labels equals the old labels
      break
  inertia = 0.0  # TODO: sum of squared distances from each row to its own centre
  return centers, labels, inertia

X = np.array([[1.0], [2.0], [9.0], [10.0]])
centers, labels, inertia = kmeans(X, [[0.0], [5.0]])
print(centers.ravel(), labels, inertia)

The tests check several properties that must hold for any correct k-means. Running more rounds never makes the inertia worse. At the end every row is labelled with its nearest final centre. And at convergence each centre equals the mean of its own rows.

#Exercise 5: A Decision Stump by Gini

In plain words. A decision stump is a tree with one question, such as "is feature 2 at most 0.45?". To choose the question, try every sensible cut and keep the one that leaves the two sides least mixed. Mixing is measured by Gini impurity: 1 minus the sum of the squared class shares. A pure group scores 0. A fifty-fifty group of two classes scores 0.5.

A cut is scored by the weighted Gini: (n_left x Gini_left + n_right x Gini_right) / n. The candidate cuts are the midpoints between neighbouring different values of a feature. Rows with a value at or below the cut go left. When two cuts tie, the smaller threshold wins.

Hand trace. The labels [0, 0, 0, 1] have shares 0.75 and 0.25, so Gini = 1 - (0.5625 + 0.0625) = 0.375. Now take x = [1, 2, 3, 4, 5, 6] with labels [0, 0, 0, 1, 1, 1]. The cut at 3.5 sends [0, 0, 0] left and [1, 1, 1] right. Both are pure, so the score is 0. For x = [1, 2, 3, 4] with labels [0, 1, 1, 0], a cut at 1.5 gives left [0] and right [1, 1, 0], and (1 x 0 + 3 x 0.4444) / 4 = 0.3333. The cut at 3.5 gives the same score, so the smaller threshold, 1.5, wins the tie.

Mini project

A Gini decision stump

~8 min

Write gini(y), best_split(x, y) that returns (threshold, weighted Gini) for one feature, and best_stump(X, y) that returns (feature index, threshold, weighted Gini) over all features. Cut only between different values. A feature with one value cannot split.

Expected output

0.5 0.375 0.0
(3.5, 0.0)

Your codeEdits save in this browser

pythonrunnable cell

import numpy as np

def gini(y):
  """Gini impurity of a list of labels: 1 minus the sum of squared class shares. Empty gives 0.0."""
  y = np.asarray(y)
  if len(y) == 0:
      return 0.0
  # TODO: shares = count of each class / len(y), then 1 - sum(shares ** 2)
  return 1.0

def best_split(x, y):
  """Best threshold for one feature. Rows with x <= threshold go left, the rest go right.

Try the midpoint between each pair of neighbouring DISTINCT values. Score a split by the
  weighted Gini: (n_left * gini_left + n_right * gini_right) / n.
  Return (threshold, score). Ties go to the smaller threshold.
  If the feature has one distinct value, there is nothing to split: return (None, gini(y)).
  """
  x = np.asarray(x, dtype=float)
  y = np.asarray(y)
  best_t, best_score = None, gini(y)
  # TODO: loop over the midpoints of np.unique(x); keep the lowest weighted Gini
  return best_t, best_score

def best_stump(X, y):
  """Best single-question rule over all features. Return (feature_index, threshold, score).

Ties go to the lower feature index. If no feature can split, return (None, None, gini(y)).
  """
  X = np.asarray(X, dtype=float)
  # TODO: call best_split on each column and keep the lowest score
  return None, None, gini(y)

print(gini([0, 0, 1, 1]), gini([0, 0, 0, 1]), gini([1, 1, 1]))
print(best_split([1, 2, 3, 4, 5, 6], [0, 0, 0, 1, 1, 1]))

The most useful test recomputes the weighted Gini of your returned threshold from scratch and compares it with the score you returned. A function that splits by position in the sorted list, rather than by value, gives a score that does not match its own threshold when values repeat, and this test catches it. A brute-force search on random data with repeated values catches the rest.

#Exercise 6: PCA by Eigendecomposition

In plain words. PCA finds the directions along which the data spreads most. First centre the data by subtracting each column's mean. Then build the covariance matrix, Xc^T Xc / (n - 1). Its eigenvectors are the directions, and its eigenvalues are the variances along them. The direction with the biggest eigenvalue is the first principal component. Projecting the centred data onto the top directions gives the scores.

np.linalg.eigh handles symmetric matrices like a covariance matrix. It returns the eigenvalues in increasing order and the eigenvectors as columns, so you must reverse the order and turn the chosen columns into rows. A direction and its negative are the same direction, so the sign of a component is arbitrary, and the tests ignore it.

Hand trace. The rows (1,1), (2,2), (3,3), (4,4) have column means 2.5, so the centred rows are -1.5, -0.5, 0.5, 1.5 in both columns. The covariance is [[5/3, 5/3], [5/3, 5/3]]. Its eigenvalues are 10/3 = 3.3333 and 0, and the first eigenvector is (1, 1) / sqrt(2) = (0.7071, 0.7071). The scores are the centred rows times that direction: -2.1213, -0.7071, 0.7071 and 2.1213, up to one overall sign.

Mini project

PCA from scratch

~7 min

Write pca(X, n_components) that returns the components (one unit-length direction per row, biggest spread first), their variances (dividing by n minus 1), and the scores of the centred data. Do not change X.

Expected output

[[0.7071 0.7071]] [3.3333] [-2.1213 -0.7071  0.7071  2.1213]

Your codeEdits save in this browser

pythonrunnable cell

import numpy as np

def pca(X, n_components):
  """Principal components by eigendecomposition of the covariance matrix.

Return (components, variances, scores):
    components  shape (n_components, d), one unit-length direction per row, biggest spread first
    variances   shape (n_components,), the variance along each direction (divide by n - 1)
    scores      shape (n, n_components), the centred data expressed along those directions
  The sign of a component is arbitrary (v and -v are the same direction).
  """
  X = np.asarray(X, dtype=float)
  n, d = X.shape
  # TODO 1: centre the columns (subtract the column means)
  Xc = X
  # TODO 2: covariance = Xc^T Xc / (n - 1)
  cov = np.eye(d)
  # TODO 3: np.linalg.eigh(cov) gives eigenvalues in INCREASING order. Reorder from largest to smallest.
  values, vectors = np.linalg.eigh(cov)
  # TODO 4: keep the first n_components; components are the COLUMNS of vectors, so transpose them into rows
  components = vectors[:, :n_components].T
  variances = values[:n_components]
  scores = Xc @ components.T
  return components, variances, scores

X = np.array([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [4.0, 4.0]])
components, variances, scores = pca(X, 1)
print(np.round(components, 4), np.round(variances, 4), np.round(scores.ravel(), 4))

Several tests approach the answer by a different route. One compares your variances with np.cov. Another compares your directions with the singular value decomposition of the centred data. A third rebuilds the original data from the scores with all components kept, and a fourth checks that the variances add up to the trace of the covariance matrix.

#Exercise 7: Naive Bayes

In plain words. Naive Bayes picks the class with the highest posterior chance: the class prior times the chance of seeing these feature values in that class. It is called naive because it treats the features as independent given the class, so the chance of all the features is the product of the chance of each one. For numeric features, each chance comes from a bell curve fitted to that class's values: a mean and a variance per feature per class.

Multiplying many small numbers underflows to exactly 0 in floating point, and then every class ties. The fix is to add logs instead of multiplying. For one feature the log of the bell-curve density is -0.5 * (log(2 * pi * v) + (x - m) ** 2 / v). A score for a class is its log prior plus the sum of those over the features.

Hand trace. The class 0 rows are 1, 2, 3, so its mean is 2 and its variance (dividing by 3) is 2/3. The class 1 rows are 7, 8, 9: mean 8, variance 2/3. The priors are both 0.5. For x = 4 the squared gap over the variance is 4 / (2/3) = 6 for class 0 and 16 / (2/3) = 24 for class 1. The terms with 2 pi v are equal, so class 0 scores higher by 0.5 x (24 - 6) = 9, and the prediction is 0. For x = 6 the gaps swap, and the prediction is 1.

Mini project

Gaussian naive Bayes in log space

~8 min

Write fit_gnb(X, y) that returns a dict with classes, priors, means and variances (divide by the count, then add 1e-9). Write log_scores(model, X) that returns an (n, classes) table of log prior plus summed log densities, and predict_gnb(model, X) that returns the original labels.

Expected output

[2. 8.] [0.6667 0.6667]
[0 1]

Your codeEdits save in this browser

pythonrunnable cell

import numpy as np

def fit_gnb(X, y):
  """Learn a Gaussian naive Bayes model. y holds whole-number class labels.

Return a dict with:
    classes    the sorted distinct labels
    priors     share of rows in each class
    means      shape (n_classes, d), the mean of each feature in each class
    variances  shape (n_classes, d), the variance of each feature in each class
               (divide by the count, not count - 1) plus 1e-9 so it is never zero
  """
  X = np.asarray(X, dtype=float)
  y = np.asarray(y)
  classes = np.unique(y)
  # TODO: for each class c, use only the rows X[y == c] to get prior, mean and variance
  priors = np.zeros(len(classes))
  means = np.zeros((len(classes), X.shape[1]))
  variances = np.ones((len(classes), X.shape[1]))
  return {"classes": classes, "priors": priors, "means": means, "variances": variances}

def log_scores(model, X):
  """Shape (n, n_classes): log(prior) + the sum over features of log N(x; mean, variance).

log N(x; m, v) = -0.5 * (log(2 * pi * v) + (x - m) ** 2 / v). Work in logs, never multiply densities.
  """
  X = np.asarray(X, dtype=float)
  # TODO: build the table, one column per class
  return np.zeros((len(X), len(model["classes"])))

def predict_gnb(model, X):
  """The class with the highest log score for each row."""
  scores = log_scores(model, X)
  return model["classes"][np.argmax(scores, axis=1)]

X = np.array([[1.0], [2.0], [3.0], [7.0], [8.0], [9.0]])
y = np.array([0, 0, 0, 1, 1, 1])
model = fit_gnb(X, y)
print(np.round(model["means"].ravel(), 4), np.round(model["variances"].ravel(), 4))
print(predict_gnb(model, np.array([[4.0], [6.0]])))

Two tests are aimed at the log trick. One uses 1000 features, where a product of densities underflows to 0 for every class, so a multiplying solution labels everything as class 0. The other uses a feature that is constant inside a class, which has variance 0 unless you add the small constant.

#Exercise 8: A Cross-Validation Split

In plain words. To judge a model fairly, test it on rows it did not learn from. K-fold cross-validation does that k times. Cut the rows into k slices, called folds. For each fold in turn, train on the other k minus 1 folds and score on the held-out fold. The final score is the average of the k fold scores.

Two details matter. The rows must be shuffled first, because a file sorted by label would give folds that hold only one class. And no row may be in both the training part and the validation part of the same split. This leak is easy to cause and hard to see, because the score just looks excellent.

Hand trace. Ten rows in 3 folds cannot be cut evenly. np.array_split makes the sizes 4, 3, 3, with the larger folds first. The validation sets then have sizes 4, 3, 3, and the matching training sets have sizes 6, 7, 7. Together the three validation folds hold each of the rows 0 to 9 exactly once.

Mini project

K-fold cross-validation

~7 min

Write kfold_indices(n, k, seed) that returns k pairs (train_idx, val_idx) from a seeded shuffle, and cross_val_accuracy(predict_fn, X, y, k, seed) that returns the plain mean of the k fold accuracies. predict_fn(X_train, y_train, X_val) is any function that returns predicted labels.

Expected output

number of folds: 3
validation sizes: [4, 3, 3]

Your codeEdits save in this browser

pythonrunnable cell

import numpy as np

def kfold_indices(n, k, seed=0):
  """Split row numbers 0 .. n-1 into k folds. Return a list of k pairs (train_idx, val_idx).

Shuffle the row numbers first with np.random.default_rng(seed). Fold sizes differ by at most 1,
  the larger folds first (np.array_split does this). Every row is in exactly one validation fold.
  """
  # TODO 1: shuffled = np.random.default_rng(seed).permutation(n)
  # TODO 2: folds = np.array_split(shuffled, k)
  # TODO 3: for each i, val = folds[i] and train = every other fold joined together (np.concatenate)
  return []

def cross_val_accuracy(predict_fn, X, y, k=5, seed=0):
  """Mean of the k fold accuracies.

predict_fn(X_train, y_train, X_val) returns predicted labels for X_val. It may only learn from the
  training part of each split.
  """
  X = np.asarray(X)
  y = np.asarray(y)
  scores = []
  # TODO: for each (train, val) pair, predict the validation rows from the training rows and
  #       append the share of correct predictions to scores
  return float(np.mean(scores)) if scores else 0.0

folds = kfold_indices(10, 3, seed=0)
print("number of folds:", len(folds))
print("validation sizes:", [len(v) for _, v in folds])

The leak test passes a spy function as predict_fn. The spy raises an error if any validation row also appears in the training rows it receives, so a split that reuses rows is caught at the moment it happens. Another test builds an uneven split, with 23 rows in 4 folds, and checks that you report the plain mean of the fold scores and not a mean weighted by fold size.

#Where the Traps Are

Every exercise has one common wrong answer that looks right. They fall into four families.

Scale and centring. Distances and covariances depend on the units. PCA without centring finds the average, not the spread. Nearest neighbours on raw features lets the largest unit win. This lab fixes the data, so the traps show up in the tests, but on real data you rescale first, and you fit the scaler on the training rows only.

Divide by what. The covariance for PCA divides by n minus 1. The class variances in naive Bayes divide by n. The loss gradient divides by n so the step size does not grow with the data. Each choice is defensible, and each changes the numbers, so the one you use needs to be stated.

Numbers that quietly break. An overflowing sigmoid, a log of exactly 0, a product of 1000 small densities, and a variance of exactly 0 all produce infinities or zeros that then spread through the answer. Clip, work in logs, and add a tiny constant.

Leakage. A validation row inside the training set, a scaler fitted on all the rows, or a feature built from the label all give a score that is too good. The split function is where leakage is easiest to prevent and hardest to notice.

Symptom	Likely cause
Gradient descent blows up to huge numbers	The gradient used a sum where it needed a mean, or the step size is too large
Logistic weights have the opposite sign	The gradient used y minus p instead of p minus y
Nearest neighbours gives one label for everything	The farthest rows were taken, or one feature has much larger units
K-means gives nan for a centre	The mean of an empty cluster was taken
A stump returns a threshold equal to a data value	The cut was placed on a point, not between two different values
PCA says the first direction is the data's average	The data was not centred
Naive Bayes always predicts the first class	Densities were multiplied and underflowed to 0, or the prior was left out
Cross-validation accuracy is near perfect	A validation row also sat in the training part

#Capstone Checklist

Use this list to check that you can read and write the basic models without a library. Tick each item only if you can do it without looking.

I can write the normal equation and a gradient descent loop for a line, and say why both reach the same weights
I can write a sigmoid that cannot overflow, and the log loss with clipping
I can write the gradient of the log loss and say which way the sign goes
I can write a nearest-neighbour vote and state how ties are broken
I can write the assign and update steps of k-means, and say what happens to an empty cluster
I can compute the Gini impurity and the weighted Gini of a cut, and place cuts between values, not on them
I can compute principal components from a covariance matrix, and say why the data is centred first
I can write Gaussian naive Bayes in log space and say why
I can split rows into shuffled folds and name one way a validation row could leak into training
Given a library call, I can name the choices it hides: tie rules, n versus n minus 1, centring, and clipping

If you can tick every line, you can read the source of a classical ML library and explain its defaults. If you cannot, go back to the exercise for that line, open the solution, and write the function again without it.

#Key Takeaways

Linear regression has two routes to one answer. The normal equation solves (Xb^T Xb) w = Xb^T y, and gradient descent steps w by lr times (2 / n) Xb^T (Xb w - y). They must agree
Logistic regression is the same loop with a sigmoid and the gradient Xb^T (p - y) / n. Split the sigmoid on the sign of z so it cannot overflow, and clip p so the log stays finite
Nearest neighbours has no training step. Take the k smallest distances, vote, and state the tie rule
K-means alternates assign and update. The inertia never rises, an empty cluster keeps its old centre, and at the end each centre is the mean of its rows
A stump scores each cut between different values by the weighted Gini, with ties going to the smaller threshold and the lower feature index
PCA centres the data, builds the covariance matrix, and takes its eigenvectors from the largest eigenvalue down. The variances add up to the trace
Naive Bayes adds a log prior to a sum of per-feature log densities, which survives a thousand features. Cross-validation shuffles row numbers, cuts them into folds, and never lets a row be in both parts

#Quick Check

Quick Check1 / 4

A logistic regression has been trained for one step from all-zero weights. Every predicted chance is 0.5. Which expression gives the gradient of the average log loss?

Where to next. Revisit Model Evaluation and read its cross-validation code with the leakage trap in mind. When you are ready for networks that learn their own features, continue to Track 4, Deep Learning.

Classical ML From Scratch Lab

Learning Objectives

After this lesson, you will be able to:

Fit a line with the normal equation and with gradient descent, and show that both reach the same weights
Write a sigmoid that cannot overflow, the log loss, and a gradient step for logistic regression
Label a new row by a vote of its k nearest rows, with a stated rule for ties
Write the two alternating steps of k-means, and handle a cluster that owns no rows
Score a split with Gini impurity and find the best single question over all features
Compute principal components from a covariance matrix, and write Gaussian naive Bayes in log space
Split rows into shuffled folds so that no row is ever in training and validation together

Prerequisites

Before You Start

Review

#How This Lab Works

Three things about the tests are worth knowing before you start.

They are tolerance based. A result close enough to the right answer passes, so a harmless rounding difference never fails you.
They are deterministic. Any random numbers come from a fixed seed, so a test that fails today fails the same way tomorrow.
They check more than the hand trace. Some compare your answer with a second route to the same number, such as numpy's own least-squares routine, an independent Newton solver, a slow loop, or a brute-force search. If your code only matches the trace, it will still fail those.