40 - Linear regression: hypothesis, mean squared error and the cost function

You want to estimate the price of a house from its square metres. You have twenty listings, you put them on a sheet with the square metres on the horizontal axis and the price on the vertical one, and it comes naturally to take a ruler and draw the line that "fits best" among the points. You have just done, by eye, a linear regression. What is missing to make a computer do it is a precise definition of "fits best": a cost function, a number that measures how wrong a line is. The whole rest of the module, from logistic regression to neural networks, is built on the same idea: choose a model with parameters, define a cost, search for the parameters that minimise it.

This lesson is the first building block. You derive the cost function from an assumption about noise, find the optimal line with mathematics, check it with code, and discover why summary numbers (like R²) can mislead.

The words we will use (and will no longer take for granted)

An example (or sample, or row) is a single observation, for instance a house listing: we have m of them. The feature is an input value that describes the example (the square metres) and is written x; the target (or label) is the value we want to predict (the price), written y, and since it is a number it takes us into the territory of regression. The parameters are the numbers the model learns (θ₀ and θ₁): you do not choose them, the algorithm searches for them. The model (or hypothesis) is the rule that turns features into a prediction, made specific by those parameters.

The last term holds up all the rest: the cost function is a number, computed over all the data, that says how much the model is wrong; the lower it is, the better. Training means finding the parameters that make it as low as possible. Keep the distinction between the three levels in mind, because it will return in every lesson: the data are given, the parameters are what changes, the cost is the yardstick by which we decide what to change.

The hypothesis: a line

With a single feature x (square metres) and a target y (price), the hypothesis is a linear function with two parameters:

h(x) = θ₀ + θ₁ · x

θ₁ is the slope (how much the price grows for each additional square metre), θ₀ the intercept (the value at x = 0). Every pair (θ₀, θ₁) gives a different line. Learning means choosing the right pair. For each example i, the error (or residual) is the difference between prediction and true value:

e_i = h(x_i) − y_i

How to read the drawing: each blue point is a listing (square metres horizontally, price vertically). The black line is the hypothesis for one choice of θ₀ and θ₁. The dashed red segments are the residuals: the vertical distance between the prediction on the line and the real price. A better line is one that makes these segments, overall, shorter. Moving θ₀ raises or lowers the whole line; changing θ₁ rotates it around its point on the vertical axis.

The cost function: mean squared error

We must summarise the m errors in a single number. Adding them up does not work: positive and negative errors cancel. We square them and take the mean:

J(θ₀, θ₁) = (1 / 2m) · Σ_{i=1..m} ( h(x_i) − y_i )²

This is the mean squared error (MSE) with a factor ½ for convenience: when you differentiate the square a 2 appears that the ½ cancels. The factor does not change where the minimum is. But why the square and not the absolute value? Two reasons. The first is practical: the square is differentiable everywhere and produces a smooth function, convenient for optimisation. The second is deeper and you derive it below.

Why precisely the square: the Gaussian likelihood

Suppose the data arise like this: y_i = θ₀ + θ₁x_i + ε_i, where the noise ε_i is Gaussian with mean 0 and variance σ², independent from one example to the next. Then the density of y_i given x_i is

p(y_i | x_i; θ) = 1/√(2πσ²) · exp( − (y_i − h(x_i))² / (2σ²) )

The likelihood of the whole dataset is the product of the densities; taking the logarithm (which turns products into sums and does not move the maximum):

log L(θ) = Σ_i log p(y_i | x_i; θ) = −m · log√(2πσ²) − (1 / 2σ²) · Σ_i (y_i − h(x_i))²

The first term does not depend on θ. Maximising the likelihood is therefore equivalent to minimising Σ(y_i − h(x_i))²: least squares is the maximum-likelihood estimate under Gaussian noise. This also tells you when the square is a bad choice: if the noise has heavy tails (frequent outliers), the Gaussian assumption is false and the square overweights anomalous values.

The cost surface and the closed-form solution

J is a function of two variables (θ₀, θ₁): a bowl-shaped surface, a paraboloid. It has a single minimum, no deceptive local minima: J is convex. To find it, you set the partial derivatives to zero. In matrix notation: X is the m × 2 matrix with a column of ones (for the intercept) and one with the x values; θ is the vector (θ₀, θ₁).

J(θ) = (1/2m) · (Xθ − y)ᵀ(Xθ − y)

∇J(θ) = (1/m) · Xᵀ(Xθ − y) = 0       →       XᵀX θ = Xᵀy       →       θ = (XᵀX)⁻¹ Xᵀ y

This is the normal equation. For the single-feature case it simplifies into two formulas you can compute by hand:

θ₁ = Σ(x_i − x̄)(y_i − ȳ) / Σ(x_i − x̄)²  =  cov(x, y) / var(x)
θ₀ = ȳ − θ₁ · x̄

The second says something simple: the optimal line always passes through the mean point (x̄, ȳ).

A worked example by hand with three points

Before relying on the computer, do the calculation with pen and paper on three tiny listings: (x, y) = (1, 2), (2, 3), (3, 5). The means are x̄ = 2 and ȳ = 10/3. Let us compute the two ingredients of the formula:

Σ(x − x̄)(y − ȳ) = (−1)(2 − 10/3) + (0)(3 − 10/3) + (1)(5 − 10/3) = 4/3 + 0 + 5/3 = 3
Σ(x − x̄)²       = 1 + 0 + 1 = 2

θ₁ = 3 / 2 = 1.5          θ₀ = 10/3 − 1.5 · 2 = 1/3

The line is h(x) = 1/3 + 1.5x. The predictions are 1.833 / 3.333 / 4.833, so the residuals are −0.167 / +0.333 / −0.167 and the cost is J = (1/6)(0.0278 + 0.1111 + 0.0278) = 0.0278. Note one property: the residuals sum to zero. This is no accident. If the sum were not zero, you could raise or lower the whole line by a constant and reduce the cost, so you would not be at the minimum. Let us verify with the code:

import numpy as np

xs = np.array([1.0, 2.0, 3.0])
ys = np.array([2.0, 3.0, 5.0])

theta1 = ((xs - xs.mean()) * (ys - ys.mean())).sum() / ((xs - xs.mean()) ** 2).sum()
theta0 = ys.mean() - theta1 * xs.mean()
residui = (theta0 + theta1 * xs) - ys

print("theta0 = %.4f (1/3 = %.4f)   theta1 = %.4f" % (theta0, 1 / 3, theta1))
print("residuals:", np.round(residui, 4), "  sum of residuals: %.1e" % residui.sum())
print("cost J = %.4f" % ((residui ** 2).sum() / (2 * len(xs))))

theta0 = 0.3333 (1/3 = 0.3333)   theta1 = 1.5000
residuals: [-0.1667  0.3333 -0.1667]   sum of residuals: 8.9e-16
cost J = 0.0278

The code repeats exactly the steps done by hand: two means, two sums, a ratio. residui (Italian for "residuals", kept as the variable name) is computed as prediction minus true value, the same sign as the definition e_i = h(x_i) − y_i. The sum of the residuals is zero up to rounding errors (of the order of 10⁻¹⁶).

load  x (t)     2     4     6     8
fuel  y (l)    30    34    42    46

Checking on more data, in four ways

We generate data from the relation y = 3 + 2x with Gaussian noise, and compute the line in four independent ways:

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
m = 100
x = rng.uniform(0, 10, m)
y = 3 + 2 * x + rng.normal(0, 2, m)                 # truth: theta0 = 3, theta1 = 2, sigma = 2

# 1) formulas for one feature
t1 = ((x - x.mean()) * (y - y.mean())).sum() / ((x - x.mean()) ** 2).sum()
t0 = y.mean() - t1 * x.mean()

# 2) normal equation in matrix form
X = np.column_stack([np.ones(m), x])                # shape (m, 2): column of ones + x
theta = np.linalg.solve(X.T @ X, X.T @ y)           # better than inverting the matrix

# 3) numpy.polyfit   4) scikit-learn
pf = np.polyfit(x, y, 1)[::-1]
lr = LinearRegression().fit(x.reshape(-1, 1), y)

print("formulas:         theta0 = %.4f  theta1 = %.4f" % (t0, t1))
print("normal equation:  theta0 = %.4f  theta1 = %.4f" % tuple(theta))
print("polyfit:          theta0 = %.4f  theta1 = %.4f" % tuple(pf))
print("scikit-learn:     theta0 = %.4f  theta1 = %.4f" % (lr.intercept_, lr.coef_[0]))
print("shape of X:", X.shape, " shape of theta:", theta.shape)

formulas:         theta0 = 2.9851  theta1 = 1.9739
normal equation:  theta0 = 2.9851  theta1 = 1.9739
polyfit:          theta0 = 2.9851  theta1 = 1.9739
scikit-learn:     theta0 = 2.9851  theta1 = 1.9739
shape of X: (100, 2)  shape of theta: (2,)

What the code does, in order. rng.uniform(0, 10, m) generates 100 values of x; y is built by applying the true line (3 + 2x) and adding Gaussian noise with standard deviation 2: the data are synthetic, so we know the right answer and can see how close the methods get. The first method applies the two covariance formulas. The second builds the matrix X with a column of ones (so the intercept θ₀ becomes a coefficient like the others: x₀ = 1) and solves XᵀXθ = Xᵀy. The third and the fourth are library routines; polyfit returns the coefficients from the highest degree to the lowest, which is why the array is reversed with [::-1].

The four methods agree. The estimated values are close to the true ones (3 and 2) but not equal: with 100 points and noise, the estimate has an uncertainty. Note also np.linalg.solve instead of inv: solving the system is numerically more stable than inverting the matrix and multiplying.

The cost surface, numerically

To see that the solution really is the bottom of the bowl, we compute J on a grid of (θ₀, θ₁) values and search for the minimum by brute force:

def J(t0, t1):
    return ((t0 + t1 * x - y) ** 2).sum() / (2 * m)

g0 = np.linspace(0, 6, 241); g1 = np.linspace(1, 3, 241)
valori = np.array([[J(a, b) for b in g1] for a in g0])
i, j = np.unravel_index(valori.argmin(), valori.shape)

print("minimum on the grid: theta0 = %.3f  theta1 = %.3f   J = %.4f" % (g0[i], g1[j], valori[i, j]))
print("exact solution:      theta0 = %.3f  theta1 = %.3f   J = %.4f" % (theta[0], theta[1], J(*theta)))
print("J at nearby points:  J(t0+0.5, t1) = %.3f   J(t0, t1+0.2) = %.3f" % (J(theta[0] + 0.5, theta[1]), J(theta[0], theta[1] + 0.2)))

minimum on the grid: theta0 = 2.975  theta1 = 1.975   J = 1.8801
exact solution:      theta0 = 2.985  theta1 = 1.974   J = 1.8801
J at nearby points:  J(t0+0.5, t1) = 2.005   J(t0, t1+0.2) = 2.665

The code defines J as a function of (θ₀, θ₁), evaluates it on a grid of 241 × 241 pairs and takes the position of the smallest value with argmin/unravel_index. The grid finds 2.975 and 1.975, one step from the exact solution (2.985, 1.974), and J is 1.8801 in both. Moving θ₀ by +0.5 the cost rises to 2.005; moving θ₁ by +0.2 it rises to 2.665: the surface is steeper in the direction of θ₁, because an error on the slope is amplified as x grows.

Moving away from the minimum in any direction raises J: it is the behaviour of a bowl. This convexity is what will let you, in the next lesson, descend towards the minimum with gradient descent without risking stopping in the wrong place. For linear regression the closed-form solution exists, but with many features (or with nonlinear models) it does not exist or costs too much: this is where the need for iterative optimisation comes from.

Why it really is a minimum

Setting the gradient to zero finds a stationary point: it could be a maximum or a saddle. To make sure it is a minimum you look at the second derivative (the Hessian matrix). Differentiating ∇J = (1/m)Xᵀ(Xθ − y) again with respect to θ gives H = (1/m)XᵀX. For any vector v, vᵀHv = (1/m)·‖Xv‖² ≥ 0: a sum of squares is never negative. A Hessian with this property (positive semidefinite) describes a bowl without saddles or maxima: the function is convex and every stationary point is a global minimum. If the columns of X are linearly independent, the inequality is strict and the minimum is unique.

Measuring quality: MSE, RMSE, MAE, R²

MSE  = (1/m) Σ (y_i − ŷ_i)²                    same unit as y, squared
RMSE = √MSE                                    same unit as y
MAE  = (1/m) Σ |y_i − ŷ_i|                     less sensitive to outliers
R²   = 1 − Σ (y_i − ŷ_i)² / Σ (y_i − ȳ)²       share of variance explained

R² compares your model with the dumbest possible model, the one that always predicts the mean ȳ: R² = 1 is perfect, R² = 0 means "no better than the mean", a negative R² means worse than the mean (possible on the test set).

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

pred = X @ theta
mse = ((y - pred) ** 2).mean()
r2 = 1 - ((y - pred) ** 2).sum() / ((y - y.mean()) ** 2).sum()
print("by hand:   MSE %.3f  RMSE %.3f  MAE %.3f  R2 %.4f" % (mse, np.sqrt(mse), np.abs(y - pred).mean(), r2))
print("sklearn:   MSE %.3f  RMSE %.3f  MAE %.3f  R2 %.4f" % (
    mean_squared_error(y, pred), mean_squared_error(y, pred) ** 0.5, mean_absolute_error(y, pred), r2_score(y, pred)))
print("the true noise has sigma = 2, i.e. variance 4")

by hand:   MSE 3.760  RMSE 1.939  MAE 1.559  R2 0.9049
sklearn:   MSE 3.760  RMSE 1.939  MAE 1.559  R2 0.9049
the true noise has sigma = 2, i.e. variance 4

The two calculations coincide (MSE 3.760, RMSE 1.939, MAE 1.559, R² 0.9049). The MSE is close to the variance of the noise, which is 4 and which no model can eliminate. It is even a little below: the model has "absorbed" a small part of the noise of the 100 points on which it was trained. It is a foretaste of overfitting (the lesson on overfitting and regularisation will return to this): a low MSE on the training set does not prove that the model is better. It becomes so only if computed on data not used for training. The RMSE (1.94) has the same unit as the target and reads as the "typical error"; the MAE (1.56) is the mean of the absolute values and is less sensitive to large errors. R² = 0.905 means that the line explains about 90% of the variance of y.

The trap: same numbers, different data

In 1973 Francis Anscombe published four small datasets built on purpose to have the same summary statistics. They have the same mean, the same variance, the same correlation and therefore the same regression line. But they tell very different stories. Here are the data and the calculation:

x123 = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5], dtype=float)
x4   = np.array([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8], dtype=float)
ys = [np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68]),
      np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74]),
      np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73]),
      np.array([6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])]
for k, (xx, yy) in enumerate(zip([x123, x123, x123, x4], ys), 1):
    b1, b0 = np.polyfit(xx, yy, 1)
    r = np.corrcoef(xx, yy)[0, 1]
    print("dataset %d: theta0 = %.2f  theta1 = %.3f  correlation = %.3f  R2 = %.3f" % (k, b0, b1, r, r ** 2))

dataset 1: theta0 = 3.00  theta1 = 0.500  correlation = 0.816  R2 = 0.667
dataset 2: theta0 = 3.00  theta1 = 0.500  correlation = 0.816  R2 = 0.666
dataset 3: theta0 = 3.00  theta1 = 0.500  correlation = 0.816  R2 = 0.666
dataset 4: theta0 = 3.00  theta1 = 0.500  correlation = 0.817  R2 = 0.667

The code uses a trick: zip pairs each list of x with its y, keeping in mind that the first three datasets share the same x and the fourth does not. For each one np.polyfit(xx, yy, 1) returns slope and intercept (in order from the highest degree, which is why b1, b0) and np.corrcoef the correlation, whose square is R² for a single feature. The output is the surprising fact: the four rows are almost indistinguishable, and the differences in the last digit (0.666 against 0.667) are only rounding.

The four datasets give practically the same line and the same R² (about 0.67). And yet: the first is a reasonable cloud around the line; the second is a curve (a parabola) that the line approximates badly; the third is a perfect line with a single outlier that drags its slope; the fourth has all points with the same x except one, which alone determines the line. No summary statistic distinguishes these cases. Always plot the data and the residuals before trusting a number.

Outliers and the square

Because the errors are squared, a distant point weighs enormously: a residual of 30 contributes nine hundred times more to the cost than a residual of 1. To see it unambiguously we reuse the dataset from before (x, y, 100 points with true line 3 + 2x) and add just one, at x = 9.5 with y = −30, a value that the true line would put around 22. Then we compare the two fits with np.polyfit.

x_out = np.append(x, 9.5); y_out = np.append(y, -30.0)      # one point with y completely out of scale
b1_out, b0_out = np.polyfit(x_out, y_out, 1)
print("without outlier:  theta0 = %.3f  theta1 = %.3f" % (theta[0], theta[1]))
print("with 1 outlier:   theta0 = %.3f  theta1 = %.3f   (on %d points)" % (b0_out, b1_out, len(x_out)))

without outlier:  theta0 = 2.985  theta1 = 1.974
with 1 outlier:   theta0 = 3.689  theta1 = 1.754   (on 101 points)

A single point out of 101 moves the line noticeably: the slope goes from 1.974 to 1.754 and the intercept from 2.985 to 3.689, a change to which the other 100 points did not contribute in any way. Remedies: find and understand the outliers (are they measurement errors or real cases?), use a robust loss such as Huber, or quantile regression.

load  x (t)     2     4     6     8    10
fuel  y (l)    30    34    42    46    20

The bridge to neural networks

Look at the hypothesis with new eyes: h(x) = θ₁x + θ₀ is a neuron with one input, a weight (θ₁), a bias (θ₀) and no activation function. With more inputs it becomes w·x + b. The cost function you have just derived is the same one you will use to train networks that predict numbers. In the rest of the module only two things will change: the model will become more complex (more layers, nonlinear activations) and the closed-form solution will give way to gradient descent.

Trade-offs

Linear regression has three advantages: the closed-form solution is fast, every coefficient reads as an effect, and the low variance makes it an excellent baseline to beat before moving to a more complex model. The price is that it represents only relationships that are linear in the features, it is sensitive to outliers, and the inversion of XᵀX becomes unstable when the features are correlated with one another (multicollinearity) or too numerous. The scikit-learn documentation describes the problem like this: if some columns of X are nearly linearly dependent, the matrix becomes nearly singular and the least-squares estimate very sensitive to random errors in the target.

Real pitfalls

• High R² = good model: false (Anscombe). R² always grows when you add features, even useless ones, on the training set.
• Extrapolating: the line continues to infinity, reality does not. A price model for houses between 40 and 150 m² does not hold for 1000 m².
• Correlation and causation: the coefficient describes an association in the data, not a causal effect.
• Shape mismatch: scikit-learn wants X with shape (m, d) even for a single feature: x.reshape(-1, 1). Passing a vector (m,) gives an error, and a column vector (m, 1) in place of y (m,) can produce silent (m, m) broadcasts in hand calculations.
• Leakage: a feature built using the target (for instance "mean price of the neighbourhood" computed including the house itself) makes the training error collapse and lies about everything else.

The five pitfalls have a common root: a number computed on the training set (R², a coefficient, an error) is read as if it described the world. The sections that follow measure three of them: the uncertainty of the coefficients, the residuals that tell what R² hides, and extrapolation. The other two, correlation versus causation and leakage, are not solved with a calculation but with knowledge of how the data were collected.

How far to trust the coefficients: standard error and confidence interval

So far we have found the best line and nothing more. But the data are a sample: with another sample of the same phenomenon we would have obtained slightly different coefficients. How different? The answer is needed to say whether a slope is "really positive" or could be zero. For linear regression a closed formula exists. If the noise has variance σ², the covariance matrix of the estimated coefficients is σ² (XᵀX)⁻¹, and σ² is estimated from the residuals: σ̂² = Σ eᵢ² / (n − p), where n is the number of points and p the number of coefficients (the "− p" accounts for the degrees of freedom used to estimate the coefficients). The standard error of each coefficient is the square root of the corresponding diagonal element. A 95% confidence interval is the coefficient plus or minus about two standard errors (more precisely the quantile of Student's t distribution with n − p degrees of freedom, 1.99 here).

We verify the formula in two independent ways: with the bootstrap (resampling the points with replacement and redoing the regression 2,000 times), and by really simulating many experiments to see how often the interval contains the true value.

rng_ci = np.random.default_rng(0)
n_ci = 80
x_s = rng_ci.uniform(0, 10, n_ci)
y_s = 2.0 + 1.5 * x_s + rng_ci.normal(0, 2, n_ci)                         # truth: intercept 2, slope 1.5, noise 2

A = np.column_stack([np.ones(n_ci), x_s])
beta = np.linalg.solve(A.T @ A, A.T @ y_s)
residuals = y_s - A @ beta
sigma2 = residuals @ residuals / (n_ci - 2)                               # estimate of the noise variance
std_err = np.sqrt(np.diag(sigma2 * np.linalg.inv(A.T @ A)))
t_crit = 1.99                                                             # 97.5% quantile of t with 78 degrees of freedom
print("intercept %.3f (std error %.3f)   slope %.3f (std error %.3f)   estimated sigma %.3f (true 2)" % (beta[0], std_err[0], beta[1], std_err[1], np.sqrt(sigma2)))
print("95%% interval for the slope (formula):   [%.3f, %.3f]" % (beta[1] - t_crit * std_err[1], beta[1] + t_crit * std_err[1]))

boot_slopes = []
for _ in range(2000):
    idx = rng_ci.integers(0, n_ci, n_ci)
    boot_slopes.append(np.linalg.lstsq(A[idx], y_s[idx], rcond=None)[0][1])
print("95%% interval for the slope (bootstrap): [%.3f, %.3f]   (std of the slopes %.3f)" % (*np.percentile(boot_slopes, [2.5, 97.5]), np.std(boot_slopes)))

intercept 2.365 (std error 0.424)   slope 1.424 (std error 0.072)   estimated sigma 1.886 (true 2)
95% interval for the slope (formula):   [1.282, 1.567]
95% interval for the slope (bootstrap): [1.274, 1.578]   (std of the slopes 0.078)

Formula and bootstrap agree on this sample, but a single sample does not tell you whether "95%" is a promise kept. To find out you must repeat the whole experiment many times, with new data but the same truth, and count how often the interval contains the true slope 1.5.

hits = 0
for r in range(1000):                                                     # 1000 independent experiments with the same truth
    gen = np.random.default_rng(1000 + r)
    xx = gen.uniform(0, 10, n_ci); yy = 2 + 1.5 * xx + gen.normal(0, 2, n_ci)
    AA = np.column_stack([np.ones(n_ci), xx]); bb = np.linalg.solve(AA.T @ AA, AA.T @ yy)
    ee = yy - AA @ bb; se = np.sqrt(ee @ ee / (n_ci - 2) * np.linalg.inv(AA.T @ AA)[1, 1])
    hits += bb[1] - t_crit * se <= 1.5 <= bb[1] + t_crit * se
print("the 95%% interval contained the true slope in %.1f%% of cases" % (hits / 10))

the 95% interval contained the true slope in 94.7% of cases

What the code does. It generates 80 points from a line with intercept 2, slope 1.5 and noise of standard deviation 2, solves the normal equation by hand and computes the standard error with the formula σ̂² (AᵀA)⁻¹. The estimated slope is 1.424 with standard error 0.072; the noise deviation estimated from the residuals is 1.886 (the true one is 2: a sample of 80 points does not recover it exactly). The 95% interval with the formula is [1.282; 1.567], the bootstrap one [1.274; 1.578]: two independent methods that agree, and that both contain the true value 1.5.

The second block repeats the whole experiment 1,000 times with different samples and counts how often the interval contains the true slope: 94.7%, close to the promised 95%. This is the operational meaning of "95% confidence interval": not "the slope has a 95% probability of lying here", but "if you repeated the experiment many times, 95% of the intervals built this way would contain the true value". Remember, though, that the 94.7% coverage holds because the data were generated respecting the assumptions (Gaussian, homoscedastic noise, a truly linear relationship). With real data the assumptions are approximate, and the real coverage will be lower: the bootstrap is more robust to the assumptions about the noise distribution, but not to the one about the shape of the relationship.

Looking at the residuals: when the line is wrong in a systematic way

The standard-error formulas assume that the noise has the same variance everywhere (homoscedasticity). If it is not so (heteroscedasticity: the error grows with the value, as often happens with prices, incomes, counts) the coefficients remain reasonable but the intervals are wrong. The way to notice is to look at the residuals, never only at R². We generate data in which the noise deviation grows with x and measure the dispersion of the residuals in four bands.

x_h = rng_ci.uniform(0, 10, 400)
y_h = 2 + 1.5 * x_h + rng_ci.normal(0, 1, 400) * (0.3 + 0.5 * x_h)         # the noise deviation grows with x
A_h = np.column_stack([np.ones(400), x_h])
res_h = y_h - A_h @ np.linalg.solve(A_h.T @ A_h, A_h.T @ y_h)
print("x band          standard deviation of the residuals")
for lo, hi in ((0, 2.5), (2.5, 5), (5, 7.5), (7.5, 10)):
    mask = (x_h >= lo) & (x_h < hi)
    print("%4.1f - %4.1f     %.3f" % (lo, hi, res_h[mask].std()))
print("correlation between |residual| and x: %.2f" % np.corrcoef(np.abs(res_h), x_h)[0, 1])

x band          standard deviation of the residuals
 0.0 -  2.5     0.927
 2.5 -  5.0     2.084
 5.0 -  7.5     3.081
 7.5 - 10.0     4.819
correlation between |residual| and x: 0.58

What the code does. It generates 400 points on a line whose noise has deviation 0.3 + 0.5·x, therefore increasing, fits the line and splits the residuals into four bands of x. The standard deviation of the residuals is 0.927 in the first band, 2.084 in the second, 3.081 in the third and 4.819 in the last: it grows fivefold from left to right, and the correlation between the absolute value of the residual and x is 0.58. A single summary number like R² would not show it; the plot of the residuals against x (funnel-shaped) or this table would.

The consequences are practical: the standard error computed with a single σ̂² describes badly a model in which the true error is 0.9 in one zone and 4.8 in another; the prediction intervals are too wide on the left and too narrow on the right; and least squares gives the same weight to all points, including the noisiest ones, instead of weighting them less. The remedies are three: transform the response variable (the logarithm turns a proportional error into a constant one), use weighted least squares (weights inversely proportional to the variance), or compute standard errors robust to heteroscedasticity. None of the three changes the line dramatically: what changes is your knowledge of how far to trust it.

Outliers: least squares against robust regressions

In the section on outliers we saw that squaring the error gives a lot of weight to distant points. Estimators exist built on purpose to resist: Huber regression uses the square for small errors and the absolute value for large ones (so a distant point weighs linearly, not quadratically); Theil–Sen, in the classic one-dimensional version, takes the median of the slopes of all pairs of points (scikit-learn uses a generalisation of the median to more dimensions); RANSAC repeatedly fits a line to small random subsets and keeps the one with the most nearby points. We corrupt 6 points out of 80 by adding 40 to their y value, and compare with the line on the clean data (intercept 2.365, slope 1.424).

from sklearn.linear_model import HuberRegressor, RANSACRegressor, TheilSenRegressor

y_corrupt = y_s.copy()
y_corrupt[:6] += 40                                                       # 6 points out of 80 become outliers
print("estimator       intercept    slope       (truth: 2.0 and 1.5; on clean data 2.365 and 1.424)")
for name, model in (("OLS", LinearRegression()), ("Huber", HuberRegressor()), ("RANSAC", RANSACRegressor(random_state=0)), ("Theil-Sen", TheilSenRegressor(random_state=0))):
    model.fit(x_s[:, None], y_corrupt)
    fitted = model.estimator_ if name == "RANSAC" else model
    print("%-15s  %8.3f    %8.3f" % (name, np.ravel(fitted.intercept_)[0], np.ravel(fitted.coef_)[0]))

estimator       intercept    slope       (truth: 2.0 and 1.5; on clean data 2.365 and 1.424)
OLS                 6.519       1.199
Huber               2.879       1.392
RANSAC              3.460       1.225
Theil-Sen           2.730       1.400

What the code does. It takes the same 80 points, adds 40 to y for the first 6 (7.5% of the data) and fits four estimators. Ordinary least squares (OLS) is dragged: the slope falls from 1.424 to 1.199 and the intercept rises from 2.365 to 6.519, an error of more than four units on a parameter that is worth about 2. Huber (slope 1.392, intercept 2.879) and Theil–Sen (1.400 and 2.730) stay close to the clean values. RANSAC, with the default parameters, does better than OLS but not well (1.225 and 3.460): it is an honest result to report, because RANSAC is sensitive to the choice of the "closeness" threshold and here we did not tune it.

Do not conclude that "Huber is better than OLS": on clean data OLS is the most efficient estimator, and a robust one costs a little precision when the outliers are absent. The choice depends on what you know about the data: if the gross errors are measurement errors not to be chased, a robust estimator is more prudent; if they are real and important events (the huge transaction that really matters), chasing them with OLS is correct, and a robust estimator would hide them. The first step is never choosing the estimator: it is going to look at the six points.

Extrapolating: the model does not know what lies outside the data

The last pitfall in the list, extrapolating, is measured with an experiment in which we know the truth. We generate data from a square root (a curve that flattens), observe them only for x between 0 and 5, and fit a line. Then we ask the line what happens at x = 10 and x = 20, where the data do not reach, and compare with the true value.

x_e = rng_ci.uniform(0, 5, 100)
y_e = 3 + 4 * np.sqrt(x_e) + rng_ci.normal(0, 0.3, 100)                    # truth: a square root, not a line
line = LinearRegression().fit(x_e[:, None], y_e)
print("x     line prediction         true value (noise-free)     error")
for x_query in (2.5, 5, 10, 20):
    truth = 3 + 4 * np.sqrt(x_query)
    prediction = line.predict([[x_query]])[0]
    print("%-5g %-22.2f  %-26.2f  %+.2f" % (x_query, prediction, truth, prediction - truth))

x     line prediction         true value (noise-free)     error
2.5   8.98                    9.32                        -0.35
5     12.45                   11.94                       +0.51
10    19.40                   15.65                       +3.75
20    33.29                   20.89                       +12.40

What the code does. It generates data from a relationship that is not a line but a square root (3 + 4·√x), observed only for x between 0 and 5, and fits a line. Inside the observed range the line is fine: at x = 2.5 it is off by −0.35 and at x = 5 by +0.51, errors of the same order as the noise (0.3). Outside the range the model has no way of knowing that the curve bends: at x = 10 the error is +3.75, at x = 20 it is +12.40 (33.29 predicted against 20.89 true, 59% more). Inside the observed range the fit looks good: no statistic computed on the training data (residuals, R², standard errors of the coefficients) warns you about what happens outside. The standard error measures the uncertainty due to the sample within the assumed model, not the error due to the wrong model.

This is why prediction outside the range of the training data (extrapolation) is the classic source of disasters in applications: the price of a 1,000 m² house, the demand in a month never seen. Practical rule: check whether the new inputs fall inside the envelope of the training data and, if they do not, treat it as a question the model cannot answer, not as a prediction with a wider margin of error.

Three exercises

Exercise 1. Repeat the coverage simulation with the same truth but with n = 10 points instead of 80: still using 1.99 in place of the correct quantile of the t with 8 degrees of freedom (2.31), how much does the coverage fall? What does it teach about using 2 "by eye"? Exercise 2. In the outlier experiment, raise the corrupted points from 6 to 20 (25%): which estimators still hold? Exercise 3. For the heteroscedastic data, transform y with the logarithm (after making the values positive) and look at the residuals in the four bands again: does the dispersion become uniform? What is the price of this transformation for the interpretation of the coefficients?

To practise on a new domain, the fuel consumption of a river ferry, there is the file esercizi.py (Python 3.12+ and NumPy 2.x; its comments and messages are in Italian): the least-squares line written from scratch in closed form, the four metrics with a case of negative R², and the leave-one-out analysis that finds the wrong row among thirteen on its own. The file runs with python3 esercizi.py even with the exercises not completed and tells you which ones are missing: launched as it is, it prints "Esercizi completati: 0 su 3" ("Exercises completed: 0 of 3").

In the next lesson, gradient descent, you tackle the case in which the closed form is not enough: how to find the minimum of J by descending along its slope.

Previous lesson in the course: Dimensionality reduction: PCA from scratch and t-SNE. Other lessons in the series that link to this one: Statistics: Estimating, Comparing, Deciding, Reinforcement Learning: The Agent That Learns by Trial and Error (and the Exploration-Exploitation Dilemma).