ML · Week 1 · Note 2 of 15

Linear regression

Linear regression predicts a number as a weighted sum of features. It is the simplest model that learns, the one interviewers ask about most, and the pattern that every bigger model repeats inside itself. This note builds it from a line through three points to the matrix formula, solves it two ways, and shows how to read and distrust the result.

Why this note matters

“Explain linear regression” is the most common opening question in ML interviews, because a good answer touches everything: the model, the loss, how it is optimized, what the coefficients mean, what can go wrong. It is also the model you should fit first on any numeric problem, because it trains in milliseconds and gives you a baseline that a fancier model must beat. And every layer of a neural network ? (week 5) is literally this model with many outputs, so the shapes and the update rule you learn here are reused in week 5 without change.

Questions that are really this note:

• “Explain linear regression, how it is trained, and its assumptions.”
• “Why do we minimize squared error rather than absolute error?”
• “When would you use the closed-form solution and when gradient descent?”
• “What does R² mean? Can it be negative?”
• “The coefficient for a feature is huge. How do you interpret it, and when should you not?”
• “Is a polynomial fit still linear regression?”

In plain words

The model: a line, then a hyperplane

With one feature:

With $d$ features $\mathbf{x} = (x_1, \dots, x_d)$:

The parameters are the weight vector $\mathbf{w}$ and the intercept $b$, $d + 1$ numbers in all. The model is called linear because the prediction is a weighted sum of the inputs: double a feature’s value and its contribution exactly doubles, and the effect of one feature does not depend on the value of another.

A small trick makes the formulas cleaner and is used in every textbook: pretend there is an extra feature $x_0$ that is always equal to 1, and call the intercept $w_0$. Then $\hat{y} = w_0 \cdot 1 + w_1 x_1 + \dots = \boldsymbol{\theta}\cdot\mathbf{x}$ with no separate $b$. For a whole dataset, stack the examples as rows of a table (a matrix) $X$ of shape $n \times (d+1)$ whose first column is all ones, and every prediction at once is one matrix-vector product: every row of $X$ dotted with $\boldsymbol{\theta}$, the “many dot products at once” idea from the end of the vectors note. (The matrices note, next in the reading order, does the bookkeeping of shapes properly.)

The loss: why squared errors

The mean squared error over $n$ examples is

The second form uses the L2 norm: the MSE is the squared Euclidean distance between the vector of true values and the vector of predictions, divided by $n$. The quantities $r_i = y_i - \hat{y}_i$ are the residuals, and the loss is their mean square. Four reasons squared error is the default, each of which an interviewer may ask for:

• It is smooth with a single minimum. As a function of the parameters, MSE is a parabola in each one and a bowl overall (Fig. 2). There is exactly one bottom, no local traps, and the slope exists everywhere (the technical word is convex ?). Absolute error has a sharp corner at zero, which makes both the formula and gradient descent messier.
• It has a closed-form solution. “Find where the bowl is flat” can be written as equations you can solve directly. The next section does it.
• It is what “most likely” gives you under bell-curve noise. If the truth is a line plus random wobble shaped like a bell curve, then the line that makes the observed prices most probable is exactly the least-squares line. The maximum likelihood note in week 3 derives this; it is the deep reason behind the choice.
• It penalizes large errors heavily. An error of 10 costs 100; two errors of 5 cost 50. If big misses are much worse than small ones in your application, this is right. If your data has outliers you do not trust, it is wrong, and you switch to absolute error or a hybrid; the loss functions note next week compares them.

Solving it exactly: least squares

For one feature, write $\bar{x}$ and $\bar{y}$ for the means. The least-squares slope and intercept are:

The numerator of $w$ is how $x$ and $y$ vary together (the covariance, up to a constant) and the denominator is how $x$ varies on its own (its variance). The intercept formula says the line passes through the point of means $(\bar{x}, \bar{y})$. Work the houses:

So the best line is $\hat{y} = 180x + 23.3$, with predictions 203.3, 293.3, 383.3 and MSE $(3.33^2 + 6.67^2 + 3.33^2)/3 = 22.2$. No search was needed; the formula went straight to the bottom of the bowl. The reason the formula is true is that the derivative of the loss is zero at the minimum; the gradient note next week gives you the tools to derive it yourself in ten lines.

For many features, the same “find where the bowl is flat” step produces the normal equation. Two matrix moves appear in it that the matrices note (next in the reading order) makes precise, so here is the one-line version of each: $X^\top$, the transpose, is the table flipped so its rows become columns; and a superscript $-1$ on a square matrix means its inverse, the matrix that undoes it, the way multiplying by $1/5$ undoes multiplying by 5 (in code, np.linalg.inv). With $X$ the $n \times (d+1)$ data matrix (ones column included) and $\mathbf{y}$ the $n$ targets:

Read it with the shape rule. $X^\top X$ is $(d+1) \times (d+1)$, a small square matrix however many examples you have; $X^\top\mathbf{y}$ is a $(d+1)$-vector; the inverse times the vector is the $(d+1)$ parameters. For the houses, $d = 1$, so it is a 2 × 2 inverse, and the code shows it giving the same 23.33 and 180. Three practical notes:

• Cost. Forming $X^\top X$ is about $n d^2$ operations and inverting it about $d^3$. With a million rows and a hundred features that is instant. With a hundred thousand features (text with one feature per word) the inverse alone is $10^{15}$ operations, and you use gradient descent instead.
• When it breaks. $X^\top X$ has no inverse if some feature is an exact combination of others (a column of temperatures in Celsius next to one in Fahrenheit, or more features than examples). The model is then not unique. Adding a penalty on big weights (regularization ?, week 3) or using lstsq, which returns the shortest of the tied solutions, handles it.
• In code, do not invert. np.linalg.lstsq or solve is more stable and faster than inv. The explicit inverse is for the whiteboard.

Solving it by steps: gradient descent

“Which way is downhill” is the gradient ?, the subject of week 2. For MSE on a line it comes out as two simple expressions, which you can take on trust for one week (the symbol $\partial$ just means “the slope of the loss when only this one parameter moves”):

Read them: each is an average of the errors, weighted by the input for $w$ and unweighted for $b$. If predictions are too high on average, $\partial L/\partial b$ is positive and the update below pushes $b$ down. The update rule, with a learning rate $\eta$ (“eta”) that sets the step size:

One step by hand, starting from $w = 200, b = 0$ (loss 133.3) with $\eta = 0.1$:

• Predictions 200, 300, 400; errors $\hat{y} - y$ are 0, 0, +20.
• $\partial L/\partial w = \frac{2}{3}(0 \cdot 1 + 0 \cdot 1.5 + 20 \cdot 2) = 26.67$; $\partial L/\partial b = \frac{2}{3}(20) = 13.33$.
• New $w = 200 - 0.1 \times 26.67 = 197.33$; new $b = 0 - 0.1 \times 13.33 = -1.33$.
• New predictions 196, 294.7, 393.3; errors −4, −5.3, +13.3; new loss $(16 + 28.4 + 177.8)/3 = 74.1$. Down from 133.3 in one step.

Keep going and the pair drifts to $(180, 23.3)$, the closed-form answer, as the code shows. Two things to notice now and study properly next week: the learning rate matters (too small and it crawls, too large and it overshoots and diverges), and the two parameters move at very different speeds here because $x$ and the constant 1 have different scales. That second problem is why the feature scaling note exists.

Reading the result: coefficients and R²

Coefficients. $w_j$ is the change in the prediction for a one-unit increase in feature $j$, holding the other features fixed. For the houses, $w = 180$ means “180 thousand more per extra thousand square feet”. Three cautions:

• Units. Measure size in square feet instead of thousands and the coefficient becomes 0.18. Same model, same predictions, coefficient a thousand times smaller. To compare the importance of features, fit on standardized ? features (each scaled to mean 0, standard deviation 1); then each coefficient is the effect of one standard deviation, in the target’s units, and they are comparable.
• Correlated features. When two features move together (house size and number of rooms), the model can split the credit between them almost arbitrarily. Individual coefficients become unstable and can even flip sign, while predictions stay fine. This is multicollinearity; it is a problem for interpretation, not for prediction, and the fix is to drop or combine features, or to add a penalty on big weights (week 3).
• Not causation. The coefficient says what the data associates with what. Ice-cream sales have a positive coefficient for drowning deaths; the hidden feature is summer. Do not say “increasing X by one unit causes Y to rise by $w$” unless the data came from an experiment.

R². The coefficient of determination compares the model’s residuals with those of the dumbest baseline, always predicting the mean $\bar{y}$:

For the houses, $\text{SS}_\text{res} = 66.7$ and $\text{SS}_\text{tot} = 93.33^2 + 6.67^2 + 86.67^2 = 16{,}267$, so $R^2 = 1 - 66.7/16{,}267 = 0.996$: the line explains 99.6% of the variation in price. In words, $R^2$ is “the fraction of the target’s variance the model accounts for”. Four facts about it:

• On the training data, $R^2 = 1$ is a perfect fit and $R^2 = 0$ is no better than predicting the mean.
• On test data it can be negative: the model is doing worse than the test set’s mean. Negative test $R^2$ is a loud sign of overfitting or of a data shift.
• Adding any feature never lowers training $R^2$, even a random one, so it rewards bloat. Adjusted R² penalizes the feature count; better still, judge on held-out data.
• A low $R^2$ is not necessarily a bad model: if the target is genuinely noisy, 0.3 may be the best achievable and still valuable. Compare with a baseline and the domain, not with 1.

For reporting, also quote the RMSE, the square root of the MSE, because it is in the target’s own units: an RMSE of 4.7 here means “typically off by about 4.7 thousand”. Interviewers like a candidate who reports both a relative score and an absolute one.

The residual plot (residuals against the prediction or against each feature) is worth more than any single number. A pattern in it means the model is systematically wrong somewhere: a U shape means missing curvature, a fan shape means the noise grows with the prediction, a few far-off points mean outliers. The famous Anscombe’s quartet is four datasets with the same line, the same $R^2$ and completely different pictures; look at it once and you will never trust a summary number alone again.

What linear regression assumes

Two interview-grade remarks. First, the assumptions are needed for interpretation and inference far more than for prediction: a linear model with mildly non-normal, mildly heteroscedastic errors still predicts fine. Second, the assumption people forget is independence, and it is the one most often violated in practice: multiple rows from the same customer, repeated measurements, or time series. It does not bias the line much, but it makes you far more confident than you should be, and it is why train/test splits must respect groups.

Curves from a linear model

Common feature transformations, each still fitted by ordinary least squares:

• Polynomial features. Add $x^2, x^3, \dots$. Powerful and dangerous: high degrees overfit violently (the wiggly curve from the last note). Degree 2 or 3 is usually the limit.
• Log transforms. Use $\log x$ for features that span orders of magnitude (income, population), or predict $\log y$ when the target is positive and skewed (prices, counts). A coefficient on a log feature reads as “effect per percent change”, which is often the natural scale.
• Interactions. Add $x_1 x_2$ to let the effect of one feature depend on another: the value of an extra bedroom depends on the size of the house.
• One-hot encoding. Turn a category (city) into several 0/1 features so the model can give each city its own offset. The next note covers how.

This is also the bridge to neural networks ? (week 5). A network is “linear regression on features it invented itself”: its inner layers ? learn the transformations, and the last layer is exactly the weighted sum on this page. Classical feature engineering and deep learning are two answers to the same question of which features to hand the final line.

Code

Type these out. The first two reproduce the hand calculations exactly; the next two use real data; the last one draws the residual plot. Every name that is new in a block is explained in the yellow box right above it, and every library name inside the code is a link to its entry on the reference pages.

Closed form, three ways

import numpy as np

x = np.array([1.0, 1.5, 2.0])          # sizes (thousand sq ft)
y = np.array([200.0, 300.0, 380.0])    # prices (thousands)

# Simple linear regression: the two formulas from the closed-form section.
w = np.sum((x - x.mean()) * (y - y.mean())) / np.sum((x - x.mean()) ** 2)
b = y.mean() - w * x.mean()
print(w, b)                            # 180.0 23.333

# The same answer from the normal equation, which also works for many
# features.
ones = np.ones_like(x)                 # [1. 1. 1.]
# the data table with a ones column: shape (3, 2)
X = np.column_stack([ones, x])
theta = np.linalg.inv(X.T @ X) @ X.T @ y
# b first, because the ones column is first
print(theta)
# prints [ 23.333 180. ]

# In practice never invert by hand; lstsq is more stable and copes with
# redundant columns.
theta2, residuals, rank, singular_values = np.linalg.lstsq(X, y, rcond=None)
print(theta2)                          # [ 23.333 180. ]  the same answer

predictions = X @ theta                # every prediction at once
mse = np.mean((predictions - y) ** 2)
# the smallest MSE any line can reach on these points
print(mse)
# prints 22.222

Gradient descent reaching the same answer

import numpy as np

x = np.array([1.0, 1.5, 2.0])
y = np.array([200.0, 300.0, 380.0])
w = 200.0      # starting slope
b = 0.0        # starting intercept
lr = 0.1       # learning rate: how big each step is

for step in range(2000):
    predictions = w * x + b
    errors = predictions - y
    # slope of the loss in the w direction
    grad_w = (2 / len(x)) * np.sum(errors * x)
    # slope of the loss in the b direction
    grad_b = (2 / len(x)) * np.sum(errors)
    w = w - lr * grad_w                            # step downhill
    b = b - lr * grad_b
    if step in (0, 1, 10, 100, 1999):
        loss = np.mean((w * x + b - y) ** 2)
        print(f"step {step:>4}: w = {w:8.3f}   b = {b:8.3f}   loss = {loss:8.3f}")

# step    0: w =  197.333   b =   -1.333   loss =   74.074
# step    1: w =  196.356   b =   -1.600   loss =   66.966
# step   10: w =  194.592   b =    0.306   loss =   59.009
# step  100: w =  185.960   b =   13.928   loss =   28.359
# step 1999: w =  180.000   b =   23.333   loss =   22.222   <- the same answer as the closed form

Many features: R² and RMSE on real data

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

data = fetch_california_housing()
X = data.data            # 20,640 districts, 8 features each
y = data.target          # median house value, in units of $100,000
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=0)

model = LinearRegression()
model.fit(X_train, y_train)
predictions = model.predict(X_test)

# typical miss, in $100k
rmse = np.sqrt(mean_squared_error(y_test, predictions))
r2 = r2_score(y_test, predictions)
# test RMSE: 0.727   test R²: 0.594
print(f"test RMSE: {rmse:.3f}   test R²: {r2:.3f}")
# A typical miss of about 3,000, and the line removes about 59% of the
# "guess the average" error.

Comparing coefficients fairly

The raw coefficients are in mixed units (rooms, years, degrees of latitude), so they cannot be compared. The fix from the coefficients section is to standardize ? the features first, so that each coefficient means “the effect of one standard deviation”.

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression

data = fetch_california_housing()
X_train, X_test, y_train, y_test = train_test_split(
    data.data, data.target, test_size=0.2, random_state=0)

# Rescale every feature to mean 0 and spread 1, using the training data
# only.
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
model = LinearRegression().fit(X_train_scaled, y_train)

# Pair each feature name with its coefficient, then sort by size, ignoring
# the sign.
pairs = list(zip(data.feature_names, model.coef_))
pairs = sorted(pairs, key=lambda pair: abs(pair[1]), reverse=True)
for name, coef in pairs:
    print(f"{name:>12}: {coef:+.3f}")

#     Latitude: -0.900
#    Longitude: -0.871
#       MedInc: +0.826     median income: the biggest effect that is not geography
#    AveBedrms: +0.290
#     AveRooms: -0.249     opposite sign to AveBedrms: the two move together (multicollinearity)
#     HouseAge: +0.117
#     AveOccup: -0.031
#   Population: -0.009
# Latitude and Longitude are large and opposite together: they describe
# location (the coast) between them.

The residual plot catching a curved relationship

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
x = np.linspace(0, 10, 60)
# the true relationship is curved
y = 0.3 * x ** 2 + rng.normal(0, 1, 60)

# scikit-learn wants a table: shape (60, 1)
X = x[:, None]
model = LinearRegression().fit(X, y)         # but we fit a straight line
line = model.predict(X)
# what the line missed, point by point
residuals = y - line

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(x, y, s=12)
axes[0].plot(x, line, color="red")
axes[0].set_title("the straight-line fit")
axes[1].scatter(x, residuals, s=12)
axes[1].axhline(0, color="gray")
axes[1].set_title("residuals: a U shape means missing curvature")
plt.show()
# left end: the line is too low
print(residuals[:3].round(2))
# prints [5.24 4.47 4.75]
# middle: the line is too high
print(residuals[28:31].round(2))
# prints [-3.09 -2.44 -3.67]
# right end: too low again
print(residuals[-3:].round(2))
# prints [4.3  2.89 3.9 ]
# Positive, negative, positive: the U shape. Fix: add x**2 as a feature (the
# curves section).

Where you will meet this again

• Week 2. The loss functions note compares squared error with absolute error and Huber; the gradient descent note derives the two slope formulas above and explains the learning rate properly.
• Week 3. Logistic regression is this model with the output squashed into a probability. Ridge and Lasso are this model with a penalty on big weights (regularization ?), and maximum likelihood explains where squared error comes from.
• Week 5, neural networks. Every layer ? is multi-output linear regression: a 768-to-3072 layer is 3,072 linear regressions sharing the same 768 inputs, computed as one matrix multiply. The “error times input” slope formula above is exactly what backpropagation ? computes for such a layer.
• Weeks 7 to 9, language models. The final step of a language model is a linear regression from its internal state onto 50,257 scores (logits ?), one per possible next token, which the scores-to-percentages recipe (softmax ?) turns into probabilities. Fine-tuning ? for classification replaces that last linear step with a new, small one. And “linear probes”, the standard tool for asking what a model knows internally, are linear regressions fitted from its internal numbers to a label.

The full interview answer

The question is simply “explain linear regression”. Here is a complete answer in about ninety seconds that leaves the interviewer nothing obvious to add.

If they push further, be ready for:

• “Why not absolute error?” Not differentiable at zero and no closed form; but it is robust to outliers, and Huber loss gets both.
• “What if $X^\top X$ is not invertible?” Collinear features or $d > n$. Use the pseudo-inverse (lstsq) or ridge regression, which adds $\lambda I$ and always inverts.
• “Why can test R² be negative?” The model’s residuals are larger than the test set’s own variance around its mean: worse than guessing the mean.
• “Two features are highly correlated. What happens?” Predictions are fine, individual coefficients become unstable and uninterpretable. Drop one, combine them, or regularize.
• “Derive the gradient.” $\partial L/\partial \mathbf{w} = \frac{2}{n}X^\top(X\theta - y)$; week 2 shows the steps.

Practice questions

Answer each one out loud before opening it.

Practice problems

Paper first, then check. Solutions are hidden under each problem.

By hand

Code

import numpy as np

def simple_fit(x, y):
    w = np.sum((x - x.mean()) * (y - y.mean())) / np.sum((x - x.mean()) ** 2)
    b = y.mean() - w * x.mean()
    return w, b

def normal_eq_fit(x, y):
    X = np.column_stack([np.ones_like(x), x])
    # solve the 2 x 2 system; never inv
    b, w = np.linalg.solve(X.T @ X, X.T @ y)
    return w, b

def gd_fit(x, y, lr=0.05, steps=20000):
    w = 0.0
    b = 0.0
    for step in range(steps):
        errors = w * x + b - y
        w = w - lr * 2 * np.mean(errors * x)
        b = b - lr * 2 * np.mean(errors)
    return w, b

rng = np.random.default_rng(0)
# 100 random x values between 0 and 5
x = rng.uniform(0, 5, 100)
# true line w = 3, b = -2, plus noise
y = 3 * x - 2 + rng.normal(0, 0.5, 100)

for name, fit in [("simple", simple_fit), ("normal equation", normal_eq_fit), ("gradient descent", gd_fit)]:
    w, b = fit(x, y)
    print(f"{name:>16}: w = {w:.4f}   b = {b:.4f}")

#           simple: w = 2.9869   b = -2.0037
#  normal equation: w = 2.9869   b = -2.0037
# gradient descent: w = 2.9869   b = -2.0037    all three agree, close to the true 3 and -2

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score

rng = np.random.default_rng(1)
x = rng.uniform(-2, 2, 40)
# the truth is a cubic, plus noise
y = x ** 3 - x + rng.normal(0, 0.5, 40)
# first 20 points for training
x_train = x[:20]
y_train = y[:20]
x_test = x[20:]                                  # last 20 for testing
y_test = y[20:]

degrees = list(range(1, 11))
train_r2 = []
test_r2 = []
for degree in degrees:
    # least squares on x, x², ..., x^degree
    coeffs = np.polyfit(x_train, y_train, degree)
    train_r2.append(r2_score(y_train, np.polyval(coeffs, x_train)))
    test_r2.append(r2_score(y_test, np.polyval(coeffs, x_test)))

best = degrees[np.argmax(test_r2)]
print("best degree on test:", best)
for degree, tr, te in zip(degrees, train_r2, test_r2):
    print(f"degree {degree:>2}: train R² {tr:.3f}   test R² {te:.3f}")
# best degree on test: 6
# degree  1: train R² 0.721   test R² 0.693
# degree  3: train R² 0.978   test R² 0.960    the truth is a cubic: the big jump happens here
# degree  6: train R² 0.980   test R² 0.960    no real gain over degree 3
# degree  8: train R² 0.988   test R² 0.286
# degree 10: train R² 0.992   test R² -1.836   worse than guessing the average
# Train R² only goes up. Test R² jumps at degree 3, stays flat, then
# collapses: overfitting.

plt.plot(degrees, train_r2, label="train R²")
plt.plot(degrees, test_r2, label="test R²")
plt.axvline(best, color="gray", linestyle="--")
plt.xlabel("degree")
plt.legend()
plt.show()

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression

data = fetch_california_housing()
X = data.data
y = data.target

# features in their own units
raw = LinearRegression().fit(X, y)

# learns each feature's mean and spread
scaler = StandardScaler().fit(X)
X_scaled = scaler.transform(X)
standardized = LinearRegression().fit(X_scaled, y)

same = np.allclose(raw.predict(X), standardized.predict(X_scaled))
print("same predictions:", same)
# prints True: same model, only the units differ

for name, c_raw, c_std in zip(data.feature_names, raw.coef_, standardized.coef_):
    print(f"{name:>12}  raw {c_raw:+9.4f}   standardized {c_std:+.3f}")
#       MedInc  raw   +0.4367   standardized +0.830
#     HouseAge  raw   +0.0094   standardized +0.119
#    AveBedrms  raw   +0.6451   standardized +0.306    looks biggest raw, only because of units
#   Population  raw   -0.0000   standardized -0.005
#     Latitude  raw   -0.4213   standardized -0.900    biggest once units are removed
#   (AveRooms, AveOccup and Longitude omitted here)

# scaler.scale_ holds each feature's spread (standard deviation).
print(np.allclose(standardized.coef_, raw.coef_ * scaler.scale_))
# prints True: standardized = raw x spread

import numpy as np

rng = np.random.default_rng(0)
c = rng.uniform(-10, 35, 200)                    # temperature in Celsius
# the same temperature in Fahrenheit
f = 1.8 * c + 32
y = 0.5 * c + 3 + rng.normal(0, 0.3, 200)
X = np.column_stack([np.ones(200), c, f])        # ones, Celsius, Fahrenheit

# independent columns: 2, not 3
print(np.linalg.matrix_rank(X))

try:
    theta = np.linalg.inv(X.T @ X) @ X.T @ y
    # runs, but the numbers are garbage
    print("inv gave:", theta.round(1))
    # prints [-2.5 -0.2  0.2]
except np.linalg.LinAlgError as err:
    print("inv failed:", err)

theta, residuals, rank, singular_values = np.linalg.lstsq(X, y, rcond=None)
# rank 2
print("lstsq:", theta.round(3), "rank", rank)
# prints [-0.016  0.329  0.094]

rmse = np.sqrt(np.mean((X @ theta - y) ** 2))
print("prediction RMSE:", round(rmse, 3))
# prints 0.306: predictions are fine (the noise was 0.3)
# lstsq splits the effect between c and f: 0.329 + 1.8 * 0.094 = 0.498, the
# true 0.5 per degree C.
# Lesson: collinearity ruins interpretation, not prediction. Drop f, or add
# a penalty (week 3).

Common mistakes

All sources

Everything linked above, in the order to use it, plus one reference.

1. StatQuest, “The Main Ideas of Fitting a Line to Data” (video, 9 min), then “Linear Regression, Clearly Explained” (27 min) and “R-squared, Clearly Explained” (11 min).
2. Google ML Crash Course, Linear regression (reading and interactive exercises).
3. MLU-Explain, Linear Regression (interactive).
4. Khan Academy, Describing relationships in quantitative data (exercises on least squares and residuals).
5. Anscombe’s quartet (Wikipedia), for why you plot residuals.
6. NumPy, linalg.lstsq and scikit-learn, Linear models (documentation).
7. James, Witten, Hastie and Tibshirani, An Introduction to Statistical Learning, chapter 3 (free PDF). The reference for assumptions and interpretation; and Andrew Ng’s Supervised Machine Learning, week 1, for the gradient-descent view.

One-minute revision