Python linear regression from scratch

Linear regression is the smallest useful model: one input, one output, and a straight line. Building it by hand shows what fitting means, how a prediction is made, and why a low error on the training points is not the same as a model that will generalize. You will also compare the line with the mean using R² and inspect an outlier's leverage.

A line is a prediction rule

The slope says how much the prediction changes when the input increases by one. The intercept is the prediction when the input is zero. Together they turn an input x into a predicted output.

def predict(x, slope, intercept):
    return slope * x + intercept

print(predict(3, 2.0, 1.0))

Two numbers define the line

Fit with centered sums

For one input, the least-squares slope is the covariance of x and y divided by the variance of x. Centering both lists around their means keeps the arithmetic small. The intercept then passes the line through the mean point.

xs = [1, 2, 3]
ys = [2, 4, 6]
mean_x = sum(xs) / len(xs)
mean_y = sum(ys) / len(ys)
numerator = sum((x - mean_x) * (y - mean_y) for x, y in zip(xs, ys))
denominator = sum((x - mean_x) ** 2 for x in xs)
slope = numerator / denominator
intercept = mean_y - slope * mean_x
print(round(slope, 2), round(intercept, 2))

The line passes through the means

Predict, then measure error

Mean squared error averages the squared distance between each real value and its prediction. Squaring punishes large misses more than small ones, and the result is in squared output units. Report the metric on data the fit did not use whenever possible.

xs = [1, 2, 3]
ys = [2, 4, 5]
mean_x = sum(xs) / len(xs)
mean_y = sum(ys) / len(ys)
slope = sum((x - mean_x) * (y - mean_y) for x, y in zip(xs, ys)) / sum((x - mean_x) ** 2 for x in xs)
intercept = mean_y - slope * mean_x
mse = sum((y - (slope * x + intercept)) ** 2 for x, y in zip(xs, ys)) / len(xs)
print(round(slope * 4 + intercept, 2))
print(round(mse, 2))

Prediction and error are different questions

Guard the input range

A line can produce a number for any x, but that does not make the answer trustworthy outside the observed range. Record the training range and return an explicit warning or fallback instead of silently extrapolating.

def within_range(x, xs):
    return min(xs) <= x <= max(xs)

print(within_range(4, [1, 2, 3]))
print(within_range(2, [1, 2, 3]))

A number is not a guarantee

R² compares the line with the mean

R² is one minus the residual sum of squares divided by the total sum of squares around the mean. A value near one means the line explains much of the variation; zero means it is no better than always predicting the mean. It can be negative on new data when the line is worse than that baseline.

xs = [1, 2, 3, 4]
ys = [2, 4, 5, 8]
mean_y = sum(ys) / len(ys)
mean_x = sum(xs) / len(xs)
slope = sum((x - mean_x) * (y - mean_y) for x, y in zip(xs, ys)) / sum((x - mean_x) ** 2 for x in xs)
intercept = mean_y - slope * mean_x
ss_res = sum((y - (slope * x + intercept)) ** 2 for x, y in zip(xs, ys))
ss_tot = sum((y - mean_y) ** 2 for y in ys)
print(round(1 - ss_res / ss_tot, 2))

Improvement over a constant prediction

One outlier can move the line

Least squares minimizes squared residuals, so a point far from the pattern can pull the line toward itself. Compare the fit with and without the point, inspect residuals, and report the decision. Do not silently delete inconvenient observations.

def fit(xs, ys):
    mx = sum(xs) / len(xs)
    my = sum(ys) / len(ys)
    slope = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
    return slope, my - slope * mx

for ys in ([2, 4, 6, 8], [2, 4, 6, 20]):
    slope, intercept = fit([1, 2, 3, 4], ys)
    print(round(slope, 2), round(intercept, 2))

Sensitivity is part of the result

Try it yourself

Change the values, predict the output, then run the code. The Python runtime is downloaded the first time you press Run, so the first run can take a few seconds.