M04: Supervised learning via regression, risk, and loss

CSCI 3151: Foundations of Machine Learning

Frank Rudzicz

2025-11-23

Outcomes

By the end of this module, you should be able to:

• Formulate supervised regression and its pipeline
• Write down and interpret squared loss, empirical risk, and the normal-equation solution for multivariate regression.
• Implement linear and polynomial regression in Python (NumPy + scikit-learn) and compare closed-form vs gradient descent.
• Visualize cost surfaces in weight space and interpret gradient descent trajectories.
• Explain and demonstrate why naive regression is a bad classifier and how this motivates logistic regression.
• Derive (or follow) the gradient of the logistic loss for a single example (stretch).
• Apply basic rules-of-thumb for train/val/test splits and diagnose simple forms of data leakage.

The supervised learning pipeline

1. Problem & stakeholders
   • What decision? For whom? What are the costs of errors?
2. Data & labels
   • How are \(x\) and \(y\) measured? What is the data-generating process (DGP)?
3. Split data
   • Train / validation / test (or at least train / test). No peeking.
4. Choose model family
   • Linear model? Tree? Neural net? (Capacity, interpretability, constraints.)

5. Choose loss & metric(s)
   • Loss for training, metric(s) for evaluation that match decisions.
6. Train
   • Minimize empirical loss on training data (with regularization, later).
7. Evaluate & stress-test
   • On validation / test. Look for leakage, bias, distribution shift.
8. Decide & deploy
   • Choose thresholds, guardrails, and monitoring.

Data, model, loss, metric

• Supervised dataset: We observe pairs \((x^{(i)}, y^{(i)})\) for \(i = 1,\dots,n\).
  • Collect inputs in a matrix \(X \in \mathbb{R}^{n \times d}\) and targets in a vector \(y \in \mathbb{R}^n\).
• Model family: A parametric function \(f_\theta : \mathbb{R}^d \to \mathbb{R}\) (for regression) with parameters \(\theta\).
  • For linear regression: \(f_\theta(x) = w^\top x + b = \hat{y}\), where \(\theta = \langle w, b \rangle\)
• Loss function: \(\ell(\hat{y}, y)\) measures how bad it is to predict \(\hat{y}\) when the truth is \(y\).
  • E.g. squared loss \((\hat{y}-y)^2\).
• Metric(s): Summaries used for evaluation and deciding whether a model is good enough.
  • E.g. mean squared error (MSE), mean absolute error (MAE), accuracy, recall.
• Training = choose \(\theta\) to make the average loss small on the training data.

The game against the world

• The world 🌎 samples \((X, Y)\) from some unknown distribution \(P(X, Y)\).
• On each turn,
  • We choose a model \(f_\theta\) that maps \(X \mapsto \hat{Y}\).
  • For each sample, we pay loss \(\ell(f_\theta(X), Y)\) for any error we make

We care about the average loss on new, unseen data:

• Expected risk of a model \(f_\theta\): \[ L(f_\theta) = \mathbb{E}_{(X,Y)\sim P}[\ell(f_\theta(X), Y)]. \]

But:

• We never see \(P\) directly.
• We only get a finite dataset.

So we approximate \(L(f_\theta)\) by the empirical risk.

🔬 Expected vs empirical risk

Given a training dataset \(\{(x_i, y_i)\}_{i=1}^n\), define:

• Empirical risk (i.e., ‘training objective’ or ‘cost function’):

\[ \hat{L}_n(f) = \frac{1}{n}\sum_{i=1}^n \ell(f(x_i), y_i). \]

We then (try to) choose \[ \hat{f} = \mathop{\mathrm{arg\,min}}_{f \in \mathcal{F}} \hat{L}_n(f), \] where \(\mathcal{F}\) is our model class (e.g., linear functions).

Tiny numeric example (regression):

• Suppose \(y = 10,12,9\) and our model always predicts \(\hat{y} = 11\).
• Squared error loss: \((11-10)^2 = 1\), \((11-12)^2 = 1\), \((11-9)^2 = 4\).
• Empirical risk (MSE): \((1+1+4)/3 = 2\).

So the training loss we print is really \(\hat{L}_n(\hat{f})\).

🧠 Intuition:

• If \(L(f_\theta)\) is small, then on average \(f_\theta\) makes small mistakes on new data from the same data-generating process.

• We never see \(P\) directly, so we can’t compute \(L(f_\theta)\) exactly.

• Expected \(L(f_\theta) = \mathbb{E}_{(X,Y)\sim P}[\ell(f_\theta(X), Y)]\) is theoretical.

• Empirical \(\hat{L}_n(f) = \frac{1}{n}\sum_{i=1}^n \ell(f(x_i), y_i)\) is practical.

Loss functions for regression

Two common choices:

1. Squared error (L2)
   \[ \ell_{\text{SE}}(\hat{y}, y) = (\hat{y} - y)^2. \]
   • Penalizes big errors quadratically.
   • Leads to mean squared error (MSE) as the average loss.
2. Absolute error (L1)
   \[ \ell_{\text{AE}}(\hat{y}, y) = |\hat{y} - y|. \]
   • Penalizes big errors linearly.
   • Leads to mean absolute error (MAE).

Tiny numeric example:

• True \(y = 50\), models \(A\) and \(B\) predict \(\hat{y}_A = 55\), \(\hat{y}_B = 65\).
• Squared error:
  • \(A\): \((55-50)^2 = 25\), \(B\): \((65-50)^2 = 225\) (B is 9× worse).
• Absolute error:
  • \(A\): \(|55-50| = 5\), \(B\): \(|65-50| = 15\) (B is 3× worse).

🧐 Does MSE or MAE react more strongly to outliers?

1D linear regression

We start with the simplest regression problem:

• 🏹 Inputs: \(x \in \mathbb{R}\) (scalar).
• 🎯 Targets: \(t \in \mathbb{R}\) (exam score).
• 🔭 Model: a straight line \[ f_\theta(x) = w_1x + w_0, \quad \theta = \langle w_1,w_0 \rangle. \]

For each example \((x^{(i)}, t^{(i)})\):

• Prediction: \(y^{(i)} = f_\theta(x^{(i)})\).
• Residual (error): \(r^{(i)} = y^{(i)} - t^{(i)}\).

We want most residuals to be small in magnitude.

Tip

• Parameters \(\theta\) are interchangeably described as:
• ‘weight/slope’ \(w\) and ‘bias/intercept’ \(b\) and
• more generally, as a vector of weights \(\langle w_1, w_0\rangle\).
• It’s the same thing.

Squared loss and 1D empirical risk

For a single point, squared loss: \[ \ell_{\text{SE}}(y,t) = \frac{1}{2}(y-t)^2. \]

For a whole dataset \(\{(x^{(i)}, t^{(i)})\}_{i=1}^n\), empirical risk: \[ J(w_1,w_0) = \hat{L}_n(f_\theta) = \frac{1}{2n}\sum_{i=1}^n \big(w_1x^{(i)} + w_0 - t^{(i)}\big)^2. \]

Training 1D linear regression \(\equiv\) finding \((w_1,w_0)\) that minimizes \(J(w_1,w_0)\).

We’ll now see this both:

• In data space (lines through points), and
• In weight space (cost surface as a function of \(\langle w_1,w_0 \rangle\)).

Data space, weight space

import numpy as np
import matplotlib.pyplot as plt

# ----------------------------
# Synthetic data (just for context)
# ----------------------------
rng = np.random.default_rng(0)
x_data = np.linspace(0, 1, 20)
y_data = 0.5 + 1.5 * x_data + rng.normal(scale=0.1, size=x_data.shape)

# ----------------------------
# Candidate models: (w0, w1)
# ----------------------------
W = np.array([
    [0.3, 1.0],
    [0.5, 1.5],
    [0.7, 2.0],
])

colors = plt.cm.viridis(np.linspace(0, 1, len(W)))

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
ax_data, ax_weights = axes

# ----------------------------
# Left: models in data / input space
# ----------------------------
ax_data.scatter(x_data, y_data, alpha=0.6, label="data")

x_line = np.linspace(0, 1, 200)
for (w0, w1), c in zip(W, colors):
    y_line = w0 + w1 * x_line
    ax_data.plot(x_line, y_line, color=c,
                 label=fr"$w_0={w0:.1f},\, w_1={w1:.1f}$")

ax_data.set_xlabel("x")
ax_data.set_ylabel("y")
ax_data.set_title("Competing models in data space")
ax_data.legend(fontsize=8, loc="upper left")

# ----------------------------
# Right: same models in weight space
# ----------------------------
for (w0, w1), c in zip(W, colors):
    ax_weights.scatter(w0, w1, color=c)
    ax_weights.text(w0 + 0.01, w1 + 0.01,
                    fr"({w0:.1f}, {w1:.1f})",
                    fontsize=8, va="bottom")

ax_weights.set_xlabel(r"$w_0$ (intercept)")
ax_weights.set_ylabel(r"$w_1$ (slope)")
ax_weights.set_title("Same models in weight space")
ax_weights.axhline(0, linewidth=0.5)
ax_weights.axvline(0, linewidth=0.5)

fig.tight_layout()
plt.show()

Empirical loss over weight space

import numpy as np
import matplotlib.pyplot as plt  # only if you still need it for other plots
import plotly.graph_objects as go

# ----------------------------------------------------
# 1. Synthetic data (same idea as before)
# ----------------------------------------------------
rng = np.random.default_rng(0)
x_data = np.linspace(0, 1, 20)
y_data = 0.5 + 1.5 * x_data + rng.normal(scale=0.1, size=x_data.shape)

# ----------------------------------------------------
# 2. Grid in weight space (w0, w1)
# ----------------------------------------------------
w0_vals = np.linspace(0.0, 1.0, 60)   # intercept
w1_vals = np.linspace(0.0, 3.0, 60)   # slope
W0, W1 = np.meshgrid(w0_vals, w1_vals)

# ----------------------------------------------------
# 3. Compute loss surface (MSE) over the grid
# ----------------------------------------------------
loss = np.zeros_like(W0)

for xi, yi in zip(x_data, y_data):
    y_pred = W0 + W1 * xi              # broadcast over grid
    loss += (yi - y_pred) ** 2

loss /= len(x_data)  # mean over points

loss = loss + 0.02 * np.sin(5 * W0) * np.sin(5 * W1)

# ----------------------------------------------------
# 4. Find approximate minimum on the grid (just for visualization)
# ----------------------------------------------------
min_idx = np.unravel_index(np.argmin(loss), loss.shape)
best_w0 = W0[min_idx]
best_w1 = W1[min_idx]
best_loss = loss[min_idx]

# A few other candidate models to show on the surface
candidate_params = np.array([
    [0.2, 0.8],
    [0.5, 1.0],
    [0.8, 2.5],
])

def mse_for_params(w0, w1, x, y):
    y_hat = w0 + w1 * x
    return np.mean((y - y_hat) ** 2)

cand_w0 = candidate_params[:, 0]
cand_w1 = candidate_params[:, 1]
cand_loss = [mse_for_params(w0, w1, x_data, y_data)
             for w0, w1 in candidate_params]


# ----------------------------------------------------
# 5. Plotly 3D surface
# ----------------------------------------------------
fig = go.Figure()

# Loss surface
fig.add_trace(
    go.Surface(
        x=W0,
        y=W1,
        z=loss,
        colorbar=dict(title="MSE"),
        name="Loss surface"
    )
)

# Show the other candidate models
fig.add_trace(
    go.Scatter3d(
        x=cand_w0,
        y=cand_w1,
        z=cand_loss,
        mode="markers",
        marker=dict(size=4),
        name="Other models"
    )
)

# Mark the minimum
fig.add_trace(
    go.Scatter3d(
        x=[best_w0],
        y=[best_w1],
        z=[best_loss],
        mode="markers",
        marker=dict(size=6, symbol="diamond"),
        name="Approx. minimum"
    )
)

fig.update_layout(
    title="Loss surface over weight space",
    scene=dict(
        xaxis_title="w₀ (intercept)",
        yaxis_title="w₁ (slope)",
        zaxis_title="Loss (MSE)",
    ),
    margin=dict(l=0, r=0, t=40, b=0)
)

fig  

Example: 1D regression (Python)

Generate a toy exam dataset

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(3151)

n = 60
hours = rng.uniform(0, 10, size=n)
noise = rng.normal(loc=0, scale=5, size=n)

true_w = 5.0
true_b = 50.0

score = true_b + true_w * hours + noise

df_exam = pd.DataFrame({"hours_studied": hours, "exam_score": score})
df_exam.head()

	hours_studied	exam_score
0	5.346653	70.891201
1	0.698466	58.903532
2	6.454199	86.766765
3	0.607189	65.146744
4	9.652828	102.512751

Example: 1D regression (Python)

from sklearn.linear_model import LinearRegression
import plotly.express as px
import numpy as np

# Fit a simple line so we can overlay it on the scatter

X = df_exam[["hours_studied"]].values
t = df_exam["exam_score"].values

linreg = LinearRegression()
linreg.fit(X, t)

x_grid = np.linspace(df_exam["hours_studied"].min(), df_exam["hours_studied"].max(), 200)
y_grid = linreg.predict(x_grid.reshape(-1, 1))

fig = px.scatter(
df_exam,
x="hours_studied",
y="exam_score",
opacity=0.7,
labels={"hours_studied": "hours_studied", "exam_score": "exam_score"},
title="Exam score vs hours studied (interactive)",
)

fig.add_scatter(x=x_grid, y=y_grid, mode="lines", name="linear fit")

fig

Data space: scatter and fitted line

In weight space, for 1D regression, the cost \[ J(w_1,w_0) = \frac{1}{2n}\sum_{i=1}^n (w_qx^{(i)} + w_0 - t^{(i)})^2 \] is a quadratic and convex function of \((w_1,w_0)\).

We can visualize \(J(w_1,w_0)\) over a grid of \((w_0,w_1)\) values:

# Build a grid over w and b
w_vals = np.linspace(linreg.coef_[0] - 4, linreg.coef_[0] + 4, 80)
b_vals = np.linspace(linreg.intercept_ - 20, linreg.intercept_ + 20, 80)

W, B = np.meshgrid(w_vals, b_vals)
J_grid = np.zeros_like(W)

for i in range(W.shape[0]):
    for j in range(W.shape[1]):
        w_ij = W[i, j]
        b_ij = B[i, j]
        y_ij = w_ij * hours + b_ij
        J_grid[i, j] = 0.5 * np.mean((y_ij - score) ** 2)

fig, ax = plt.subplots()
cs = ax.contour(W, B, J_grid, levels=20)
ax.clabel(cs, inline=True, fontsize=8)
ax.scatter([linreg.coef_[0]], [linreg.intercept_], color="red", label="optimum (sklearn)")
ax.set_xlabel("w")
ax.set_ylabel("b")
ax.set_title("Cost surface in weight space (contours of J(w_1,w_0))")
ax.legend()
plt.tight_layout()
plt.show()

Why is \(J(w_1,w_0)\) quadratic and convex?

Why is \(J(w_1,w_0)\) quadratic?

• Each term \[ \big(w_1x^{(i)} + w_0 - t^{(i)}\big)^2 \] is a square of a linear function in \((w_1,w_0)\).
• If you expand it, you only get terms like \(w_1^2\), \(w_1\cdot w_0\), \(w_0^2\), \(w_1\), \(w_0\), and a constant.
• Averaging over \(i\) keeps it a polynomial of degree at most 2 in \((w_1,w_0)\) → quadratic.

Why is \(J(w_1,w_0)\) convex?

• Function \(z \mapsto \tfrac{1}{2} z^2\) is 🔗 convex on \(\mathbb{R}\).
• Each term \(w_1x^{(i)} + w_0 - t^{(i)}\) is an affine (linear + constant) function of \((w_1,w_0)\).
• A convex function composed with an affine map is convex, and an average of convex functions is convex.
• Therefore the average \[ J(w,b) = \frac{1}{n}\sum_{i=1}^n \tfrac{1}{2}\big(w_1x^{(i)} + w_0 - t^{(i)}\big)^2 \] is convex in \((w_1,w_0)\).

This “single-bowl” shape is why gradient descent works so reliably for linear regression: any local minimum is a global minimum.

Gradient descent visualized

👉 The path moves roughly orthogonal to level sets, towards the optimum.

Once Loop Reflect

Gradient descent mathemetized

We can minimize \(J(w_1,w_0)\) via gradient descent:

• Gradients: \[ \frac{\partial J}{\partial w_1} = \frac{1}{n}\sum_{i=1}^n (w_1x^{(i)} + w_0 - t^{(i)})x^{(i)}, \quad \frac{\partial J}{\partial w_0} = \frac{1}{n}\sum_{i=1}^n (w_1x^{(i)} + w_0 - t^{(i)}). \]
• Update rule: \[ w_1 \leftarrow w_1 - \eta \frac{\partial J}{\partial w_1}, \quad w_0 \leftarrow w_0 - \eta \frac{\partial J}{\partial b_0}. \]

Let’s see a few steps from a bad initial guess.

def cost_and_grads(w, b, x, t):
    y = w * x + b
    r = y - t
    J = 0.5 * np.mean(r**2)
    dJ_dw = np.mean(r * x)
    dJ_db = np.mean(r)
    return J, dJ_dw, dJ_db

w, b = -2.0, 20.0  # deliberately bad starting point
eta = 1e-3
history = []

for step in range(120):
    J, dJ_dw, dJ_db = cost_and_grads(w, b, hours, score)
    history.append((w, b, J))
    w = w - eta * dJ_dw
    b = b - eta * dJ_db

print("w_1, w_0, J")
history[-5:]

w_1, w_0, J

[(8.900952575260046, 22.568603162487364, 123.89854231186739),
 (8.907953940299901, 22.577338417235847, 123.77433362139037),
 (8.914690043533586, 22.586030971432944, 123.65444212230794),
 (8.92116953659019, 22.594682154627684, 123.53859294767682),
 (8.927400789228505, 22.603293253068422, 123.426528848044)]

Multidimensions & vectorization

Now let \(x \in \mathbb{R}^d\) with \(d>1\). \(x\) is a multidimensional vector

• Stack data into a matrix \(X \in \mathbb{R}^{n \times d}\).
• Parameters: \(w_1 \in \mathbb{R}^d\), bias \(w_0 \in \mathbb{R}\).
• Predictions on all \(n\) examples: \[ y = Xw_1 + w_0, \] or with the bias trick (augment \(x\) with a 1): \[ \tilde{x} = \begin{bmatrix} x \\\\ 1 \end{bmatrix}, \quad \tilde{w} = \begin{bmatrix} w_1 \\\\ w_0 \end{bmatrix}, \quad y = X_{\text{aug}}\tilde{w}. \]

Empirical risk: \[ J(w_1,w_0) = \frac{1}{2n} \|Xw_1 + w_0\mathbf{1} - t\|_2^2. \]

This is still a convex quadratic in \((w_1,w_0)\).

🔬 Closed-form solution (normal equations)

Assume we use the bias trick and write \(y = Xw_1\) with \(X \in \mathbb{R}^{n \times d}\) (including a column of 1s).

Then: \[ J(w_1) = \frac{1}{2n}\|Xw_1 - t\|_2^2. \]

Setting the gradient to zero: \[ \begin{align} \nabla_{w_1} J(w_1) &= \frac{1}{n} X^\top (Xw_1 - t) \\ &= 0\\ \quad \Rightarrow \quad X^\top X w_1 &= X^\top t. \end{align} \]

If \(X^\top X\) is invertible: \[ \boxed{w^* = (X^\top X)^{-1} X^\top t} \]

This is the normal equation solution for linear regression.

In practice:

• For moderate \(d\), this is very fast.
• For large \(d\), we often prefer iterative methods (GD, SGD).

Multivariate example in Python

1. Create a synthetic multivariate dataset and write the closed-form vector solution

# Create a synthetic multivariate regression problem
rng = np.random.default_rng(3151)
n = 200
d = 3

X_mv = rng.normal(size=(n, d))
true_w_mv = np.array([2.0, -1.0, 0.5])
true_b_mv = 3.0

noise = rng.normal(scale=0.5, size=n)
t_mv = X_mv @ true_w_mv + true_b_mv + noise

# Closed-form solution via normal equations (with bias trick)
X_aug = np.hstack([X_mv, np.ones((n, 1))])
w_closed = np.linalg.inv(X_aug.T @ X_aug) @ (X_aug.T @ t_mv)

w_closed

array([ 2.04729233, -0.96085044,  0.46363859,  3.03477073])

2. Solve it with sklearn

from sklearn.linear_model import LinearRegression

linreg_mv = LinearRegression()
linreg_mv.fit(X_mv, t_mv)

print("True w,b      :", true_w_mv, true_b_mv)
print("Closed-form   :", w_closed[:-1], w_closed[-1])
print("sklearn       :", linreg_mv.coef_, linreg_mv.intercept_)

True w,b      : [ 2.  -1.   0.5] 3.0
Closed-form   : [ 2.04729233 -0.96085044  0.46363859] 3.034770726256364
sklearn       : [ 2.04729233 -0.96085044  0.46363859] 3.034770726256364

👉 Our closed-form and sklearn solutions match (up to numerical noise)!

Nonlinear regression

Real relationships are often not linear in \(x\).

Trick: use a feature map \(\phi(x)\) to build nonlinear models while keeping linear parameters.

Example: polynomial features up to degree 3 in 1D: \[ \phi(x) = [x, x^2, x^3], \quad f_\theta(x) = w_1 x + w_2 x^2 + w_3 x^3 + w_0. \]

More generally:

• Choose \(\phi : \mathbb{R}^d \to \mathbb{R}^k\).
• Model: \(f_\theta(x) = w^\top \phi(x)\).
• All our linear-regression math still works, but the function of the original \(x\) can be nonlinear.

Polynomial regression (parody)

• Tiny dataset, smooth underlying function, and increasingly complex polynomial models:
  • Red dashed = true underlying function.
  • Dots = noisy samples (training data).
  • Coloured lines = polynomial fits.
• Good generalization usually comes from compromise. We’ll return to this concept in a later module.

import numpy as np
import plotly.graph_objects as go

rng = np.random.default_rng(3151)

# Small dataset: x in [0, 1], y = sin(2πx) + noise
n = 10
x_train = np.linspace(0.0, 1.0, n)
y_true = np.sin(2 * np.pi * x_train)
y_train = y_true + rng.normal(scale=0.2, size=n)

degrees = [0, 1, 3, 9]

def fit_poly(x, y, degree):
    """
    Fit a polynomial of given degree by least squares.
    Uses a simple Vandermonde design matrix and closed-form solution.
    """
    X = np.vander(x, N=degree + 1, increasing=True)  # [1, x, x^2, ...]
    # Pseudo-inverse for numerical stability
    w = np.linalg.pinv(X) @ y
    return w

# Dense grid for plotting curves
x_plot = np.linspace(0.0, 1.0, 400)
y_plot_true = np.sin(2 * np.pi * x_plot)

fig = go.Figure()

# Training data
fig.add_scatter(
    x=x_train,
    y=y_train,
    mode="markers",
    name="data",
    marker=dict(size=8),
)

# True underlying function
fig.add_scatter(
    x=x_plot,
    y=y_plot_true,
    mode="lines",
    name="true function",
    line=dict(dash="dash"),
)

# Polynomial fits of various degrees
for degree in degrees:
    w = fit_poly(x_train, y_train, degree)
    Xp = np.vander(x_plot, N=degree + 1, increasing=True)
    y_fit = Xp @ w

    # Only show degree 0 initially; others start as legend-only
    visible = True if degree == 0 else "legendonly"

    fig.add_scatter(
        x=x_plot,
        y=y_fit,
        mode="lines",
        name=f"degree {degree}",
        visible=visible,
    )

fig.update_layout(
    title="Under- and over-fitting with polynomial regression",
    xaxis_title="x",
    yaxis_title="y",
    legend_title="Click to toggle fits",
)

fig

• This is based on Bishop (2006), who shows this phenomenon in a slightly exaggerated fashion:

Polynomial regression example

Let’s solve it with sklearn:

from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

rng = np.random.default_rng(3151)
n = 200
x = rng.uniform(-3, 3, size=n)

def true_fun(x):
    return 2 + 0.5 * x - 1.5 * x**2 + 0.5 * x**3

y = true_fun(x) + rng.normal(scale=2.0, size=n)

X1 = x.reshape(-1, 1)
X_train1, X_val1, y_train1, y_val1 = train_test_split(
    X1, y, test_size=0.3, random_state=3151
)

degrees = [1, 3, 9]
models = {}

for deg in degrees:
    poly = PolynomialFeatures(degree=deg, include_bias=False)
    Xtr = poly.fit_transform(X_train1)
    Xva = poly.transform(X_val1)
    reg = LinearRegression()
    reg.fit(Xtr, y_train1)
    ytr_hat = reg.predict(Xtr)
    yva_hat = reg.predict(Xva)
    models[deg] = {
        "poly": poly,
        "reg": reg,
        "train_mse": mean_squared_error(y_train1, ytr_hat),
        "val_mse": mean_squared_error(y_val1, yva_hat),
    }

models

{1: {'poly': PolynomialFeatures(degree=1, include_bias=False),
  'reg': LinearRegression(),
  'train_mse': 25.86667098071478,
  'val_mse': 27.213356959052955},
 3: {'poly': PolynomialFeatures(degree=3, include_bias=False),
  'reg': LinearRegression(),
  'train_mse': 3.8888390926469065,
  'val_mse': 4.0729021289950165},
 9: {'poly': PolynomialFeatures(degree=9, include_bias=False),
  'reg': LinearRegression(),
  'train_mse': 3.7344714418676097,
  'val_mse': 4.410209605455015}}

Visualizing under/overfitting

x_plot = np.linspace(-3, 3, 400).reshape(-1, 1)
y_true_plot = true_fun(x_plot.flatten())

fig, ax = plt.subplots()
ax.scatter(x, y, alpha=0.3, label="data")
ax.plot(x_plot, y_true_plot, "k--", label="true function")

for deg, info in models.items():
    y_plot = info["reg"].predict(info["poly"].transform(x_plot))
    ax.plot(x_plot, y_plot, label=f"deg={deg}")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Polynomial regression: model complexity")
ax.legend()
plt.tight_layout()
plt.show()

👉 Deg 1 underfits, deg 3 is pretty good, deg 9 might overfit noise.

Train vs validation error vs degree

max_deg = 12
train_mses = []
val_mses = []
degrees_full = list(range(1, max_deg+1))

for deg in degrees_full:
    poly = PolynomialFeatures(degree=deg, include_bias=False)
    Xtr = poly.fit_transform(X_train1)
    Xva = poly.transform(X_val1)
    reg = LinearRegression()
    reg.fit(Xtr, y_train1)
    ytr_hat = reg.predict(Xtr)
    yva_hat = reg.predict(Xva)
    train_mses.append(mean_squared_error(y_train1, ytr_hat))
    val_mses.append(mean_squared_error(y_val1, yva_hat))

fig, ax = plt.subplots()
ax.plot(degrees_full, train_mses, marker="o", label="train MSE")
ax.plot(degrees_full, val_mses, marker="o", label="val MSE")
ax.set_xlabel("polynomial degree")
ax.set_ylabel("MSE")
ax.set_title("Model complexity vs train/validation error")
ax.legend()
plt.tight_layout()
plt.show()

• Train MSE often keeps going down with degree.
• Validation MSE typically has a sweet spot: bias–variance tradeoff in action.
• We’ll return to this in the evaluation cluster (M12–M14).

Loss functions for classification

Now switch to a binary label \(y \in \{0,1\}\) and a predicted probability \(\hat{p} \in [0,1]\) for class 1.

• 0–1 loss \[ \ell_{0/1}(\hat{y}, y) = \mathbf{1}[\hat{y} \neq y]. \]
  • 🤷🏿‍♂️ Intuitive, but non-differentiable and hard to optimize directly.
• Logistic (cross-entropy) loss \[ \ell_{\text{log}}(\hat{p}, y) = -\big(y \log \hat{p} + (1-y)\log(1-\hat{p})\big). \]
  • 👍 Smooth and convex in the logit → good for gradient-based optimization.
  • 👍 Punishes overconfident wrong predictions very strongly.

Tiny numeric example (binary):

• True label \(y=1\).
  • Model A: \(\hat{p}_A = 0.9\) → loss \(\approx -\log 0.9 \approx 0.105\).
  • Model B: \(\hat{p}_B = 0.6\) → loss \(\approx -\log 0.6 \approx 0.511\).
  • Model C: \(\hat{p}_C = 0.1\) → loss \(\approx -\log 0.1 \approx 2.302\) 😬
• Models A and B both get zero 0–1 loss (correct label with threshold 0.5), but we still prefer A (0.9) over B (0.6).
  • An approximation of ‘confidence’ in a prediction can be helpful.

Naïve regression is a bad classifier

• Scenario (ought be re-used later):
  • Features:
    • \(x_1\) = videos_watched in the first few days.
    • \(x_2\) = avg_quiz_score (%).
  • Label:
    • \(y \in \{0,1\}\) where 1 = A+, 0 = otherwise.
• Naïve idea:
  1. Treat \(y \in \{0,1\}\) as a real-valued target.
  2. Fit a linear regression model with squared loss.
  3. Threshold predictions at 0.5.
• This gives some classifier, but:
  • Predictions can be \(< 0\) or \(>1\) → not valid probabilities.
  • Squared loss is misaligned with the underlying Bernoulli structure.
  • Overconfident wrong predictions aren’t punished enough.

Let’s see this concretely.

Naïve regression classifier in Python

from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

rng = np.random.default_rng(3151)

n = 400
videos = rng.normal(loc=6, scale=2.0, size=n)
quizzes = rng.normal(loc=70, scale=10.0, size=n)

# Ground-truth logistic model for activity
w_true = np.array([0.8, 0.12])
b_true = -13.0

logits = w_true[0]*videos + w_true[1]*quizzes + b_true
probs = 1 / (1 + np.exp(-logits))

active = rng.binomial(1, probs, size=n)

import pandas as pd
df = pd.DataFrame({
    "videos_watched": videos,
    "avg_quiz_score": quizzes,
    "active_14_days": active
})
df["videos_watched"] = df["videos_watched"].clip(0, 20)
df["avg_quiz_score"] = df["avg_quiz_score"].clip(0, 100)

X_cls = df[["videos_watched", "avg_quiz_score"]].values
y_cls = df["active_14_days"].values

X_train_c, X_test_c, y_train_c, y_test_c = train_test_split(
    X_cls, y_cls, test_size=0.3, random_state=3151, stratify=y_cls
)

linreg_cls = LinearRegression()
linreg_cls.fit(X_train_c, y_train_c)

yhat_test_reg = linreg_cls.predict(X_test_c)
yhat_test_reg_cls = (yhat_test_reg >= 0.5).astype(int)

acc_reg = accuracy_score(y_test_c, yhat_test_reg_cls)
yhat_test_reg[:5], acc_reg

(array([0.51946392, 0.66849211, 0.25890321, 0.11860755, 0.96139261]), 0.825)

Regression predictions on binary labels

• Many predictions lie outside [0,1].
• Threshold at 0.5 is arbitrary and not tied to any probability semantics.

Regression predictions on binary labels

import numpy as np
import pandas as pd
import plotly.graph_objects as go
from sklearn.linear_model import LinearRegression

np.random.seed(123)

# Simulate x
n = 200
x = np.linspace(0, 10, n)

# Probability of y = 1 increases with x (logistic-ish)
p = 1 / (1 + np.exp(-(x - 5)))
y = np.random.binomial(1, p)

df = pd.DataFrame({"x": x, "y": y})

# Fit linear regression y ~ x (deliberately wrong)
model = LinearRegression()
model.fit(df[["x"]], df["y"])

x_grid = np.linspace(df["x"].min(), df["x"].max(), 200)
y_hat  = model.predict(x_grid.reshape(-1, 1))

# Jitter y so the points aren’t perfectly on 0 and 1
rng = np.random.default_rng(123)
df["y_jitter"] = df["y"] + rng.uniform(-0.05, 0.05, size=len(df))

# Build interactive figure
fig = go.Figure()

# Scatter of observed data
fig.add_trace(
    go.Scatter(
        x=df["x"],
        y=df["y_jitter"],
        mode="markers",
        marker=dict(size=6),
        name="Observed 0/1 outcomes",
        text=[f"x = {xi:.2f}<br>y = {yi}" for xi, yi in zip(df["x"], df["y"])],
        hoverinfo="text"
    )
)

# Linear regression line
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=y_hat,
        mode="lines",
        name="Linear regression fit (y ~ x)"
    )
)

fig.update_layout(
    title="Why linear regression is a bad model for binary data",
    xaxis_title="Predictor x",
    yaxis=dict(
        title="Outcome (0/1)",
        range=[-0.3, 1.3],
        tickmode="array",
        tickvals=[0, 1],
        ticktext=["0", "1"]
    )
)

fig

🧐: regression & loss at scale

The ideas in this module scale up:

• Large language models (LLMs) are trained with token-level cross-entropy loss (a big logistic-like loss over huge vocabularies).
• Image classifiers optimize cross-entropy over thousands of classes.
• Deep nets are compositions of linear maps + nonlinearities, trained by gradient-based optimization on an average loss over a dataset.

If you understand:

• Data → model → loss → empirical risk.
• Linear regression & logistic regression.
• Gradient descent as “follow negative gradient in weight space”.

…you’ve understood a core template behind much of modern ML.

Quick checks

1. For 1D linear regression, the cost \[ J(w,b) = \frac{1}{2n}\sum_{i=1}^n (wx^{(i)} + b - t^{(i)})^2 \] is, as a function of \((w,b)\):

   A. Linear.
   B. Convex quadratic.
   C. Non-convex.
   D. Discontinuous.

2. In the polynomial regression example, increasing the degree tends to:

   A. Always improve validation error.
   B. Always worsen validation error.
   C. Decrease train error, but validation error may eventually increase.
   D. Have no effect on either train or validation error.

3. In the early-warning system, which metric is most aligned with the decision to catch at-risk students?

   A. Overall accuracy.
   B. Recall for \(y=1\) (at-risk).
   C. Train MSE.
   D. Number of parameters.

4. Using test data to choose hyperparameters is best described as:

   A. Regularization.
   B. Cross-validation.
   C. Data leakage / test-set overfitting.
   D. Semi-supervised learning.

Looking ahead

Next steps in the course:

• M05 will continue the supervised story and set us up for:
  • Deeper probabilistic views (MLE, likelihood).
  • More sophisticated classification setups.
• M06 will introduce unsupervised objectives (clustering, density).

Later clusters ought revisit today’s ideas in more sophisticated forms:

• C03 MLE & probabilistic modelling: formalize LR as a likelihood model.
• C04 Evaluation & bias–variance: build on the polynomial regression curves.
• C05 Regularization & optimization: revisit linear/logistic regression with penalties and more advanced optimizers.

For now: you’ve seen how regression, loss, risk, and optimization fit together in a scientifically grounded ML pipeline.

References

Bishop, Christopher M. 2006. Pattern Recognition and Machine Learning. Springer.