Linear regression

Lesson 46 of 62, about 6 min

Linear regression fits a straight line (or, with several features, a flat plane) through data so that predictions are as close as possible to the real values. It is the first supervised algorithm to learn because everything else builds on its three ideas: a model with parameters, a cost that measures how wrong it is, and a method for making the cost smaller.

The model

With one feature, the prediction is a line:

\hat{y} = w x + b

$x$ is the input (flat size)
$\hat{y}$ (y hat) is the prediction (rent)
$w$ is the weight or slope: how much $\hat{y}$ changes when $x$ goes up by 1
$b$ is the bias or intercept: the prediction when $x = 0$

With several features the pattern is the same, one weight per feature:

\hat{y} = w_1 x_1 + w_2 x_2 + \dots + w_n x_n + b

Training means finding the values of $w$ and $b$ that fit the data best.

The cost function

"Best" needs a number. For each training example we take the error $\hat{y}_i - y_i$ , square it (so positive and negative errors do not cancel, and large errors are punished more), and average over all $m$ examples. This is the mean squared error:

J(w, b) = \frac{1}{m} \sum_{i=1}^{m} \left(\hat{y}_i - y_i\right)^2

Training minimises $J$ . Because $J$ is a smooth bowl shaped function of $w$ and $b$ , it has a single lowest point.

Computing the cost by hand

python

import numpy as np

x = np.array([450, 600, 750, 900, 1200], dtype=float)
y = np.array([9000, 12000, 15500, 18000, 24500], dtype=float)

def cost(w, b):
    predictions = w * x + b
    return ((predictions - y) ** 2).mean()

for w in (10, 15, 20, 25):
    print(f"w={w}: cost={cost(w, 0):,.0f}")

A weight near 20 gives the smallest cost. The algorithm's job is to find that automatically.

Two ways to minimise the cost

The normal equation solves for the optimum in one shot using linear algebra. It is exact and fast for small feature counts. scikit-learn's LinearRegression uses a variant of this.
Gradient descent starts from a guess and walks downhill in small steps. It scales to huge datasets and to models where no closed form exists, which is why it has its own lesson next.

The closed form, for the curious, with $X$ the feature matrix including a column of ones for the bias:

\theta = (X^T X)^{-1} X^T y

python

import numpy as np

x = np.array([450, 600, 750, 900, 1200], dtype=float)
y = np.array([9000, 12000, 15500, 18000, 24500], dtype=float)

X = np.column_stack([x, np.ones_like(x)])
w, b = np.linalg.inv(X.T @ X) @ X.T @ y
print(f"w={w:.3f}, b={b:.1f}")

With scikit-learn

python

import numpy as np
from sklearn.linear_model import LinearRegression

X = np.array([[450], [600], [750], [900], [1200]])
y = np.array([9000, 12000, 15500, 18000, 24500])

model = LinearRegression().fit(X, y)
print("weight:", model.coef_.round(3))
print("bias:  ", round(model.intercept_, 1))
print("predict 1000 sq ft:", model.predict([[1000]]).round(0))
print("R squared:", round(model.score(X, y), 4))

coef_ holds the weights, intercept_ the bias. score returns $R^2$ , the fraction of the variance in $y$ that the model explains (1.0 is perfect).

Several features

python

import numpy as np
from sklearn.linear_model import LinearRegression

# size in sq ft, bedrooms, distance to metro in km
X = np.array([[450, 1, 2.0], [600, 1, 1.0], [750, 2, 3.5], [900, 2, 0.5], [1200, 3, 4.0], [1000, 2, 1.5]])
y = np.array([9000, 13000, 14500, 21000, 23000, 20500])

model = LinearRegression().fit(X, y)
for name, w in zip(["size", "bedrooms", "distance"], model.coef_):
    print(f"{name:10s} {w:9.2f}")
print("predict:", model.predict([[800, 2, 1.0]]).round(0))

Each weight reads as "holding the other features fixed, one more unit of this feature changes the prediction by this much". The negative weight on distance says flats further from the metro rent for less.

Assumptions and limits

The relationship is roughly linear. Curves need polynomial features or a different model.
Errors have similar spread across the range of $x$ .
Features are not near duplicates of each other (multicollinearity makes weights unstable).
Outliers pull the line hard, because errors are squared.

Polynomial features

A line can still fit a curve if you add powers of the features. The model is linear in the parameters, so everything above still applies.

python

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

x = np.linspace(0, 5, 20).reshape(-1, 1)
y = 2 * x.ravel() ** 2 - 3 * x.ravel() + 1

model = make_pipeline(PolynomialFeatures(degree=2), LinearRegression()).fit(x, y)
print(model.predict([[6]]).round(2), "expected", 2 * 36 - 18 + 1)

Practice

Fit a line to hours studied [1, 2, 3, 4, 5, 6] and marks [35, 48, 55, 66, 72, 85]. Print the weight and bias and interpret them in a sentence.
Compute the mean squared error of your model on its training data by hand and compare with sklearn.metrics.mean_squared_error.
Add an outlier (hours 7, marks 20) and see what happens to the slope.