Lab: linear regression in NumPy

Time to put the whole section together. In this lab you’ll write, from scratch, the same kind of model that real estate sites use as a baseline: multiple linear regression trained with gradient descent, all in vectorised NumPy. No ML library, no shortcuts, every line yours.

The cells on this page share one Python session, like a notebook, so run them in order. If you reload the page, use Run all above to catch up.

1. Meet the data

A hidden setup cell has loaded 120 house sales from a made-up town. Each house has four features, and the target is the sale price in thousands of dollars.

print("X shape:", X.shape, "  y shape:", y.shape)
print("first three houses:\n", X[:3], "\nprices:", y[:3])
print()
for j, name in enumerate(feature_names):
    col = X[:, j]
    print(f"{name:<25} min {col.min():8.1f}   max {col.max():8.1f}   mean {col.mean():8.1f}")
print(f"{'price ($1000s)':<25} min {y.min():8.1f}   max {y.max():8.1f}   mean {y.mean():8.1f}")

So $m = 120$ and $n = 4$. Notice the ranges: size is in the thousands, bedrooms in single digits. You know what that means for gradient descent.

2. Scale the features

Z-score every column, and keep mu and sigma: you’ll need them to price new houses later.

mu = X.mean(axis=0)
sigma = X.std(axis=0)
X_norm = (X - mu) / sigma
print("scaled ranges:")
for j, name in enumerate(feature_names):
    print(f"  {name:<25} {X_norm[:, j].min():6.2f} to {X_norm[:, j].max():5.2f}")

3. Cost and gradient

These are the vectorised versions from earlier lessons. Read them once more and match each line to the maths: $J = \frac{1}{2m}\sum_i (\hat{y}^{(i)} - y^{(i)})^2$ and $\nabla_{\vec{w}} J = \frac{1}{m} X^T(\hat{\vec{y}} - \vec{y})$.

def compute_cost(X, y, w, b):
    err = X @ w + b - y
    return (err ** 2).mean() / 2

def compute_gradient(X, y, w, b):
    err = X @ w + b - y
    return X.T @ err / len(y), err.mean()

print("cost at w = 0, b = 0:", round(compute_cost(X_norm, y, np.zeros(4), 0.0), 1))

4. Write gradient descent

This is the heart of the lab. Write the training loop yourself, then check it.

import numpy as np

rng = np.random.default_rng(7)
m = 50
X_norm = rng.normal(size=(m, 3))
y = X_norm @ np.array([30.0, -12.0, 5.0]) + 200 + rng.normal(0, 3, m)

def compute_cost(X, y, w, b):
    err = X @ w + b - y
    return (err ** 2).mean() / 2

def compute_gradient(X, y, w, b):
    err = X @ w + b - y
    return X.T @ err / len(y), err.mean()

def gradient_descent(X, y, w, b, alpha, iters):
    history = []
    for _ in range(iters):
        # your code here: get the gradient, update w and b, record the cost
        pass
    return w, b, history

w, b, history = gradient_descent(X_norm, y, np.zeros(3), 0.0, alpha=0.1, iters=300)
print(w, b, history[-1] if history else None)

When your exercise passes, the cell below uses the same algorithm on the house data. (It defines its own copy, so it works even if you skipped ahead.)

def gradient_descent(X, y, w, b, alpha, iters):
    history = []
    for _ in range(iters):
        dj_dw, dj_db = compute_gradient(X, y, w, b)
        w = w - alpha * dj_dw
        b = b - alpha * dj_db
        history.append(compute_cost(X, y, w, b))
    return w, b, history

w, b, history = gradient_descent(X_norm, y, np.zeros(4), 0.0, alpha=0.1, iters=500)
print("final cost:", round(history[-1], 2))
print("weights (scaled units):", np.round(w, 2), "  b:", round(b, 2))
bpl.plot(range(1, len(history) + 1), history, kind="line", title="Learning curve, alpha = 0.1")

Check the learning curve: does it fall on every step and level off? Then try the α ladder from the last lesson: change alpha to 0.03, 0.3, 1.0 and 2.0 and rerun. Which values converge within 500 steps, and which blow up?

5. Check against a direct solver

For linear regression there’s also a formula that jumps straight to the best weights, solving a small system of equations (called the normal equation). NumPy’s np.linalg.lstsq does this. It’s a great way to check a gradient descent implementation, but it only works for linear regression and gets slow when $n$ is very large, which is why gradient descent is the tool that scales to every model in the rest of this course.

A = np.column_stack([X_norm, np.ones(len(y))])     # a column of ones plays the role of b
exact = np.linalg.lstsq(A, y, rcond=None)[0]
print("gradient descent:", np.round(np.append(w, b), 3))
print("direct solver:   ", np.round(exact, 3))
print("match?", np.allclose(np.append(w, b), exact, atol=1e-2))

6. What did the model learn?

The weights above are in scaled units (“dollars per standard deviation of size”). To read them in everyday units, undo the scaling. Since $x_{\text{scaled}} = (x - \mu)/\sigma$, a weight $w_j$ on the scaled feature equals $w_j / \sigma_j$ on the original one, and the intercept absorbs the means:

w_orig = w / sigma
b_orig = b - np.sum(w * mu / sigma)
for name, wj in zip(feature_names, w_orig):
    print(f"{name:<25} {wj:+8.3f}  thousand dollars per unit")
print(f"{'intercept b':<25} {b_orig:+8.1f}")

Read them as sentences: each extra square foot adds about $160, each year of age takes off about $1,300, each kilometre from the centre about $3,900. The town was simulated with effects of $160, $1,200 and $3,500, so from 120 noisy examples the model has landed close to the truth, though not exactly on it: the noise in the prices nudges every estimate a little. The bedrooms weight (about $7,600 against a true $9,000) is the shakiest, because bedrooms and size move together and the data can’t fully separate their effects.

7. Price a new house

import numpy as np

mu = np.array([2050.0, 3.1, 30.0, 15.5])        # training means
sigma = np.array([830.0, 1.3, 17.0, 8.4])       # training standard deviations
w = np.array([133.0, 8.0, -20.4, -29.4])        # learned weights (scaled units)
b = 322.0

def predict_price(x_raw):
    # your code here: scale with the TRAINING mu and sigma, then apply the model
    return None

my_price = None
print(my_price)

Challenges

1. No scaling. Train on raw X instead of X_norm. Find the largest α that doesn’t diverge. How far is the cost from the scaled result after 500 steps?
2. Engineer a feature. Add a column X[:, 0] / X[:, 1] (square feet per bedroom) before scaling. Does the cost go down? Is the new weight large or close to zero?
3. Fewer examples. Train on only the first 20 houses. How much do the weights in original units change? Which one moves most?

Section quiz