Algorithms From Scratch

How To Implement Simple Linear Regression From Scratch With Python

Suresh Madhra·May 11, 2020·10 min read

Linear regression is one of the oldest and most widely-used prediction methods in statistics. It's also one of the best algorithms to build from scratch: the maths is tractable, the code is short, and the mental model transfers directly to more advanced techniques.

In this walkthrough we implement simple linear regression— a single input variable predicting a single output — in pure Python. Every snippet below runs live in your browser via Pyodide. Hit Run and edit anything you like.

The model

Simple linear regression fits a straight line through the data:

y = b0 + b1 * x

We need to estimate two coefficients from the training data — the intercept b0 and the slope b1. The closed-form solution uses the mean, variance and covariance of the inputs and outputs.

1. Mean and variance

Start with two helpers: the arithmetic mean and the sum of squared differences from the mean. We use a tiny five-row contrived dataset so you can verify the numbers by hand.

mean_variance.py

# Estimate mean and variance of x and y
def mean(values):
    return sum(values) / float(len(values))

def variance(values, m):
    return sum((v - m) ** 2 for v in values)

dataset = [[1, 1], [2, 3], [4, 3], [3, 2], [5, 5]]
x = [row[0] for row in dataset]
y = [row[1] for row in dataset]
mean_x, mean_y = mean(x), mean(y)
var_x, var_y = variance(x, mean_x), variance(y, mean_y)
print('x stats: mean=%.3f variance=%.3f' % (mean_x, var_x))
print('y stats: mean=%.3f variance=%.3f' % (mean_y, var_y))

2. Covariance

Covariance tells us how two variables change together. It generalises variance from one variable to two, and it's the numerator in our slope estimate.

covariance.py

# Covariance between x and y
def mean(values):
    return sum(values) / float(len(values))

def covariance(x, mean_x, y, mean_y):
    covar = 0.0
    for i in range(len(x)):
        covar += (x[i] - mean_x) * (y[i] - mean_y)
    return covar

dataset = [[1, 1], [2, 3], [4, 3], [3, 2], [5, 5]]
x = [row[0] for row in dataset]
y = [row[1] for row in dataset]
mean_x, mean_y = mean(x), mean(y)
covar = covariance(x, mean_x, y, mean_y)
print('Covariance: %.3f' % covar)

3. Estimating coefficients

With mean, variance and covariance in hand, the slope and intercept fall out directly:

b1 = covariance(x, y) / variance(x)
b0 = mean(y) - b1 * mean(x)

coefficients.py

# Estimate b0 and b1 from the training data
def mean(values):
    return sum(values) / float(len(values))

def variance(values, m):
    return sum((v - m) ** 2 for v in values)

def covariance(x, mx, y, my):
    return sum((x[i] - mx) * (y[i] - my) for i in range(len(x)))

def coefficients(dataset):
    x = [row[0] for row in dataset]
    y = [row[1] for row in dataset]
    mx, my = mean(x), mean(y)
    b1 = covariance(x, mx, y, my) / variance(x, mx)
    b0 = my - b1 * mx
    return b0, b1

dataset = [[1, 1], [2, 3], [4, 3], [3, 2], [5, 5]]
b0, b1 = coefficients(dataset)
print('Coefficients: B0=%.3f, B1=%.3f' % (b0, b1))

4. Fit, predict, evaluate

Wrap it all up: fit the model on the training rows, predict the test rows, then score with root mean squared error (RMSE).

toy_model.py

# Fit and predict on the contrived dataset
from math import sqrt

def mean(values): return sum(values) / float(len(values))
def variance(values, m): return sum((v - m) ** 2 for v in values)
def covariance(x, mx, y, my): return sum((x[i] - mx) * (y[i] - my) for i in range(len(x)))

def coefficients(dataset):
    x = [r[0] for r in dataset]; y = [r[1] for r in dataset]
    mx, my = mean(x), mean(y)
    b1 = covariance(x, mx, y, my) / variance(x, mx)
    return my - b1 * mx, b1

def simple_linear_regression(train, test):
    b0, b1 = coefficients(train)
    return [b0 + b1 * row[0] for row in test]

def rmse_metric(actual, predicted):
    err = sum((predicted[i] - actual[i]) ** 2 for i in range(len(actual)))
    return sqrt(err / float(len(actual)))

def evaluate(dataset, algorithm):
    test_set = [[row[0], None] for row in dataset]
    predicted = algorithm(dataset, test_set)
    actual = [row[-1] for row in dataset]
    return rmse_metric(actual, predicted), predicted

dataset = [[1, 1], [2, 3], [4, 3], [3, 2], [5, 5]]
rmse, preds = evaluate(dataset, simple_linear_regression)
print('Predictions:', ['%.2f' % p for p in preds])
print('RMSE: %.3f' % rmse)

# Fit and predict on the contrived dataset
from math import sqrt

def mean(values): return sum(values) / float(len(values))
def variance(values, m): return sum((v - m) ** 2 for v in values)
def covariance(x, mx, y, my): return sum((x[i] - mx) * (y[i] - my) for i in range(len(x)))

def coefficients(dataset):
    x = [r[0] for r in dataset]; y = [r[1] for r in dataset]
    mx, my = mean(x), mean(y)
    b1 = covariance(x, mx, y, my) / variance(x, mx)
    return my - b1 * mx, b1

def simple_linear_regression(train, test):
    b0, b1 = coefficients(train)
    return [b0 + b1 * row[0] for row in test]

def rmse_metric(actual, predicted):
    err = sum((predicted[i] - actual[i]) ** 2 for i in range(len(actual)))
    return sqrt(err / float(len(actual)))

def evaluate(dataset, algorithm):
    test_set = [[row[0], None] for row in dataset]
    predicted = algorithm(dataset, test_set)
    actual = [row[-1] for row in dataset]
    return rmse_metric(actual, predicted), predicted

dataset = [[1, 1], [2, 3], [4, 3], [3, 2], [5, 5]]
rmse, preds = evaluate(dataset, simple_linear_regression)
print('Predictions:', ['%.2f' % p for p in preds])
print('RMSE: %.3f' % rmse)

5. Real data: Swedish Auto Insurance

The Swedish Auto Insurance dataset relates the number of claims to the total payment for those claims (in thousands of Swedish Kronor). 63 rows, one input, one output — a perfect fit for simple linear regression. The runner below fetches the CSV over the network, splits 60/40 into train and test, and reports the fitted model and the test RMSE.

insurance.py

# Simple linear regression on the Swedish Auto Insurance dataset.
# Fetched over the network inside your browser via Pyodide.
from random import seed, randrange
from csv import reader
from math import sqrt
from io import StringIO
from pyodide.http import pyfetch

DATA_URL = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/insurance.csv"

async def load_csv(url):
    resp = await pyfetch(url)
    text = await resp.string()
    rows = [r for r in reader(StringIO(text)) if r]
    # Skip header if present
    try:
        float(rows[0][0])
    except ValueError:
        rows = rows[1:]
    return [[float(c) for c in r] for r in rows]

def train_test_split(dataset, split=0.6):
    train, copy = [], list(dataset)
    n_train = int(len(dataset) * split)
    while len(train) < n_train:
        train.append(copy.pop(randrange(len(copy))))
    return train, copy

def mean(values): return sum(values) / float(len(values))
def variance(values, m): return sum((v - m) ** 2 for v in values)
def covariance(x, mx, y, my): return sum((x[i] - mx) * (y[i] - my) for i in range(len(x)))

def coefficients(dataset):
    x = [r[0] for r in dataset]; y = [r[1] for r in dataset]
    mx, my = mean(x), mean(y)
    b1 = covariance(x, mx, y, my) / variance(x, mx)
    return my - b1 * mx, b1

def simple_linear_regression(train, test):
    b0, b1 = coefficients(train)
    return [b0 + b1 * row[0] for row in test]

def rmse_metric(actual, predicted):
    err = sum((predicted[i] - actual[i]) ** 2 for i in range(len(actual)))
    return sqrt(err / float(len(actual)))

seed(1)
print("Fetching Swedish Auto Insurance dataset…")
dataset = await load_csv(DATA_URL)
print("Loaded", len(dataset), "rows.")
train, test = train_test_split(dataset, 0.6)
predicted = simple_linear_regression(train, [[row[0], None] for row in test])
actual = [row[-1] for row in test]
b0, b1 = coefficients(train)
print('Model: y = %.3f + %.3f * x' % (b0, b1))
print('RMSE: %.3f' % rmse_metric(actual, predicted))

# Simple linear regression on the Swedish Auto Insurance dataset.
# Fetched over the network inside your browser via Pyodide.
from random import seed, randrange
from csv import reader
from math import sqrt
from io import StringIO
from pyodide.http import pyfetch

DATA_URL = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/insurance.csv"

async def load_csv(url):
    resp = await pyfetch(url)
    text = await resp.string()
    rows = [r for r in reader(StringIO(text)) if r]
    # Skip header if present
    try:
        float(rows[0][0])
    except ValueError:
        rows = rows[1:]
    return [[float(c) for c in r] for r in rows]

def train_test_split(dataset, split=0.6):
    train, copy = [], list(dataset)
    n_train = int(len(dataset) * split)
    while len(train) < n_train:
        train.append(copy.pop(randrange(len(copy))))
    return train, copy

def coefficients(dataset):
    x = [r[0] for r in dataset]; y = [r[1] for r in dataset]
    mx, my = mean(x), mean(y)
    b1 = covariance(x, mx, y, my) / variance(x, mx)
    return my - b1 * mx, b1

def simple_linear_regression(train, test):
    b0, b1 = coefficients(train)
    return [b0 + b1 * row[0] for row in test]

def rmse_metric(actual, predicted):
    err = sum((predicted[i] - actual[i]) ** 2 for i in range(len(actual)))
    return sqrt(err / float(len(actual)))

seed(1)
print("Fetching Swedish Auto Insurance dataset…")
dataset = await load_csv(DATA_URL)
print("Loaded", len(dataset), "rows.")
train, test = train_test_split(dataset, 0.6)
predicted = simple_linear_regression(train, [[row[0], None] for row in test])
actual = [row[-1] for row in test]
b0, b1 = coefficients(train)
print('Model: y = %.3f + %.3f * x' % (b0, b1))
print('RMSE: %.3f' % rmse_metric(actual, predicted))

Expect an RMSE around 35 thousand Kronor — a solid baseline that any more sophisticated model has to beat.

Extensions

Try multiple random train/test splits and average the RMSE.
Extend the algorithm to multivariate linear regression using the normal equation.
Add a baseline (predict the mean of y) and confirm your linear model beats it.

Adapted from Jason Brownlee's original tutorial. Python execution is powered by Pyodide, running entirely in your browser — no server required.

← Back to blog