How to Implement Linear Regression From Scratch in Python

Description

Multivariate Linear Regression

Linear regression is a technique for predicting a real value. Confusingly, problems where a real value is to be predicted are called regression problems. Linear regression uses a straight line — or, in more than two dimensions, a plane or hyperplane — to model the relationship between input and output values.

Each input attribute (x) is weighted using a coefficient (b), and the goal of the learning algorithm is to discover a set of coefficients that result in good predictions (y). Coefficients can be found using stochastic gradient descent.

Stochastic Gradient Descent

Gradient descent is the process of minimizing a function by following the gradients of the cost function. In machine learning, we can use a technique that evaluates and updates the coefficients every iteration — called stochastic gradient descent (SGD) — to minimize the error of a model on our training data.

Each training instance is shown to the model one at a time. The model makes a prediction, the error is calculated, and the model is updated to reduce the error for the next prediction. This process repeats for a fixed number of iterations.

where b is the coefficient being optimized, learning_rate is a step size you configure (e.g. 0.01), error is the prediction error attributed to the weight, and x is the input value.

Wine Quality Dataset

After we develop our linear regression algorithm, we will use it to model the wine quality dataset. This dataset comprises 4,898 white wines with measurements like acidity and pH. The goal is to use these objective measures to predict wine quality on a scale between 0 and 10.

Each attribute has different units and scales, so the dataset must be normalized to the range 0–1. By predicting the mean value (Zero-Rule Algorithm) on the normalized dataset, a baseline root mean squared error (RMSE) of 0.148 can be achieved.

Tutorial

This tutorial is broken down into three parts, giving you the foundation you need to implement and apply linear regression with SGD on your own predictive modeling problems:

• Making predictions.
• Estimating coefficients.
• Wine quality prediction.

1. Making Predictions

The first step is a function that can make predictions. We'll need this both for evaluating candidate coefficient values during SGD and after the model is finalized, when we want to make predictions on test data or new data.

Below is a predict() function that returns an output value for a row given a set of coefficients. The first coefficient is always the intercept (also called the bias or b0) — it is standalone and not tied to a specific input value.

Hit Run below — the first click loads the Python runtime in your browser (~5 MB), subsequent runs are instant. Edit the code and re-run to experiment.

def predict(row, coefficients):
    yhat = coefficients[0]
    for i in range(len(row) - 1):
        yhat += coefficients[i + 1] * row[i]
    return yhat

dataset = [[1, 1], [2, 3], [4, 3], [3, 2], [5, 5]]
coef = [0.4, 0.8]
for row in dataset:
    yhat = predict(row, coef)
    print("Expected=%.3f, Predicted=%.3f" % (row[-1], yhat))

2. Estimating Coefficients

Now we estimate the coefficient values for our training data using SGD. SGD requires two parameters:

• Learning rate: how much each coefficient is corrected each time it is updated.
• Epochs: how many times to run through the training data while updating the coefficients.

Three loops are required — over each epoch, over each row in the training data, and over each coefficient and update it for a row in an epoch. Coefficients are updated based on the error the model made:

We track the sum of squared error each epoch so we can watch it fall. Click Run — error should drop from ~3.0 to ~2.5 over 50 epochs.

def predict(row, coefficients):
    yhat = coefficients[0]
    for i in range(len(row) - 1):
        yhat += coefficients[i + 1] * row[i]
    return yhat

def coefficients_sgd(train, l_rate, n_epoch):
    coef = [0.0 for i in range(len(train[0]))]
    for epoch in range(n_epoch):
        sum_error = 0
        for row in train:
            yhat = predict(row, coef)
            error = yhat - row[-1]
            sum_error += error ** 2
            coef[0] = coef[0] - l_rate * error
            for i in range(len(row) - 1):
                coef[i + 1] = coef[i + 1] - l_rate * error * row[i]
        print('>epoch=%d, lrate=%.3f, error=%.3f' % (epoch, l_rate, sum_error))
    return coef

dataset = [[1, 1], [2, 3], [4, 3], [3, 2], [5, 5]]
l_rate = 0.001
n_epoch = 50
coef = coefficients_sgd(dataset, l_rate, n_epoch)
print('Coefficients:', coef)

3. Wine Quality Prediction

Finally we apply the algorithm to the wine quality dataset. We use k-fold cross-validation to estimate the performance of the learned model on unseen data — constructing and evaluating k models and reporting the mean error. Root mean squared error (RMSE) is used to evaluate each model.

The runner below fetches the CSV over the network from GitHub, normalises the columns to 0–1, runs 5-fold cross-validation, and reports the mean RMSE. Training takes a couple of minutes in the browser — feel free to reduce n_epoch in the source to speed it up.

# Linear Regression With Stochastic Gradient Descent for Wine Quality
# Fetches the CSV directly from GitHub inside the 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/winequality-white.csv"

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

def str_column_to_float(dataset, column):
    for row in dataset:
        row[column] = float(row[column].strip())

def dataset_minmax(dataset):
    return [[min(r[i] for r in dataset), max(r[i] for r in dataset)] for i in range(len(dataset[0]))]

def normalize_dataset(dataset, minmax):
    for row in dataset:
        for i in range(len(row)):
            row[i] = (row[i] - minmax[i][0]) / (minmax[i][1] - minmax[i][0])

def cross_validation_split(dataset, n_folds):
    split, copy = [], list(dataset)
    fold_size = int(len(dataset) / n_folds)
    for _ in range(n_folds):
        fold = []
        while len(fold) < fold_size:
            fold.append(copy.pop(randrange(len(copy))))
        split.append(fold)
    return split

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

def evaluate_algorithm(dataset, algorithm, n_folds, *args):
    folds = cross_validation_split(dataset, n_folds)
    scores = []
    for fold in folds:
        train_set = [r for f in folds if f is not fold for r in f]
        test_set = []
        for row in fold:
            rc = list(row); rc[-1] = None; test_set.append(rc)
        predicted = algorithm(train_set, test_set, *args)
        actual = [row[-1] for row in fold]
        scores.append(rmse_metric(actual, predicted))
    return scores

def predict(row, coefficients):
    yhat = coefficients[0]
    for i in range(len(row) - 1):
        yhat += coefficients[i + 1] * row[i]
    return yhat

def coefficients_sgd(train, l_rate, n_epoch):
    coef = [0.0 for i in range(len(train[0]))]
    for epoch in range(n_epoch):
        for row in train:
            yhat = predict(row, coef)
            error = yhat - row[-1]
            coef[0] = coef[0] - l_rate * error
            for i in range(len(row) - 1):
                coef[i + 1] = coef[i + 1] - l_rate * error * row[i]
    return coef

def linear_regression_sgd(train, test, l_rate, n_epoch):
    predictions = []
    coef = coefficients_sgd(train, l_rate, n_epoch)
    for row in test:
        predictions.append(predict(row, coef))
    return predictions

seed(1)
print("Fetching wine quality dataset…")
dataset = await load_csv_from_url(DATA_URL)
print("Loaded", len(dataset), "rows.")
for i in range(len(dataset[0])):
    str_column_to_float(dataset, i)
minmax = dataset_minmax(dataset)
normalize_dataset(dataset, minmax)
print("Training 5 folds x 50 epochs…")
scores = evaluate_algorithm(dataset, linear_regression_sgd, 5, 0.01, 50)
print('Scores: %s' % scores)
print('Mean RMSE: %.3f' % (sum(scores) / float(len(scores))))

Expect a mean RMSE around 0.126 — well below the 0.148 Zero-Rule baseline.

Extensions

• Tune the example. Adjust the learning rate, the number of epochs, or the data preparation to improve the score on the wine quality dataset.
• Batch stochastic gradient descent. Accumulate updates across each epoch and only update the coefficients in a batch at the end of the epoch.
• Additional regression problems. Apply the technique to other regression problems on the UCI machine learning repository.

Review

In this tutorial you discovered how to implement linear regression using stochastic gradient descent from scratch with Python. You learned how to make predictions for a multivariate linear regression problem, how to optimize a set of coefficients using SGD, and how to apply the technique to a real regression predictive modeling problem.