How To Implement Logistic Regression

Logistic regression is the go-to linear classification algorithm for two-class problems. It is easy to implement, easy to understand and gets great results on a wide variety of problems — even when the assumptions of the method are violated.

In this tutorial you will discover how to implement logistic regression with stochastic gradient descent from scratch with Python — and, right here in your browser, run every snippet interactively.

Description

Logistic regression

Logistic regression is named for the logistic function at the core of the method. It uses a linear combination of the inputs, squashed by the sigmoid, to model the probability of the positive class:

Stochastic gradient descent

Gradient descent is an iterative optimiser: at each step, evaluate the slope of the loss with respect to every parameter and nudge each parameter a small amount in the opposite direction of that slope. Repeat until the loss stops falling. Batch gradient descent averages the gradient over the entire training set before each update; stochastic gradient descent uses the gradient from a single example, so it takes many more, noisier steps per epoch — but each step is cheap and, in practice, the noise helps the optimiser escape flat regions and shallow local minima.

For a single training example the log-loss gradient of logistic regression works out to a neat closed form. Combined with the learning rate η it gives the per-step update rule:

The picture below shows the same idea on a one-dimensional convex loss. Starting from a poor guess, each iteration takes a step proportional to the local slope — long strides where the surface is steep, tiny ones as the parameter approaches the minimum where the gradient vanishes.

Two knobs govern the behaviour. The learning rate η scales every step: too small and training crawls; too large and the updates overshoot the minimum and can diverge. The number of epochs sets how many full passes we make through the training data — each pass gives every example one chance to nudge the coefficients.

1. Making predictions

The first step is a predict() function that takes a row and a set of coefficients and returns a probability. The first coefficient is the intercept (b0).

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.

# Make a prediction
from math import exp

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

dataset = [
    [2.7810836, 2.550537003, 0],
    [1.465489372, 2.362125076, 0],
    [3.396561688, 4.400293529, 0],
    [1.38807019, 1.850220317, 0],
    [3.06407232, 3.005305973, 0],
    [7.627531214, 2.759262235, 1],
    [5.332441248, 2.088626775, 1],
    [6.922596716, 1.77106367, 1],
    [8.675418651, -0.242068655, 1],
    [7.673756466, 3.508563011, 1],
]
coef = [-0.406605464, 0.852573316, -1.104746259]
for row in dataset:
    yhat = predict(row, coef)
    print("Expected=%.3f, Predicted=%.3f [%d]" % (row[-1], yhat, round(yhat)))

2. Estimating coefficients

Now we estimate the coefficients with SGD. Three nested loops: over epochs, over rows in the training set, and over coefficients. We track the sum of squared error each epoch so we can watch it fall.

from math import exp

def predict(row, coefficients):
    yhat = coefficients[0]
    for i in range(len(row) - 1):
        yhat += coefficients[i + 1] * row[i]
    return 1.0 / (1.0 + exp(-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 = row[-1] - yhat
            sum_error += error ** 2
            coef[0] = coef[0] + l_rate * error * yhat * (1.0 - yhat)
            for i in range(len(row) - 1):
                coef[i + 1] = coef[i + 1] + l_rate * error * yhat * (1.0 - yhat) * row[i]
        if epoch % 10 == 0 or epoch == n_epoch - 1:
            print('>epoch=%d, lrate=%.3f, error=%.3f' % (epoch, l_rate, sum_error))
    return coef

dataset = [
    [2.7810836, 2.550537003, 0],
    [1.465489372, 2.362125076, 0],
    [3.396561688, 4.400293529, 0],
    [1.38807019, 1.850220317, 0],
    [3.06407232, 3.005305973, 0],
    [7.627531214, 2.759262235, 1],
    [5.332441248, 2.088626775, 1],
    [6.922596716, 1.77106367, 1],
    [8.675418651, -0.242068655, 1],
    [7.673756466, 3.508563011, 1],
]
l_rate = 0.3
n_epoch = 100
coef = coefficients_sgd(dataset, l_rate, n_epoch)
print(coef)

You should see the error drop from ~3.0 to under 0.03 by the last epoch, with the final coefficients printed at the bottom.

How the code works

The whole training loop is just a direct translation of the update rule above. Reading coefficients_sgd from the top:

• coef = [0.0 for i in range(len(train[0]))] — one coefficient per input column plus one for the bias, all seeded at zero. With the sigmoid, a zero vector predicts 0.5 for every row, a neutral starting point.
• for epoch in range(n_epoch) — the outer loop is one full pass over the training data. sum_error is reset each epoch so we can watch the loss trend downward.
• yhat = predict(row, coef) then error = row[-1] - yhat — compute the current prediction and its signed residual against the true label row[-1].
• coef[0] = coef[0] + l_rate * error * yhat * (1.0 - yhat) — the bias update. The yhat * (1 - yhat) factor is the derivative of the sigmoid; it damps the update near confident predictions (where the sigmoid saturates) and amplifies it near 0.5 where the model is most uncertain.
• coef[i + 1] = coef[i + 1] + l_rate * error * yhat * (1.0 - yhat) * row[i] — the weight update for each input, identical to the bias update but scaled by the input value x_j. Large inputs produce larger corrections, which is exactly why the wine-quality and diabetes examples normalise every column to 0–1 first.

Notice there is no matrix algebra, no autograd, and no library call — the entire optimiser is a dozen lines of arithmetic that faithfully implement the equation shown above.

3. Diabetes prediction

Finally we apply the algorithm to the Pima Indians diabetes dataset — 768 rows, 8 features, one binary label. The runner below fetches the CSV over the network from GitHub, normalises the features, and evaluates the model with 5-fold cross-validation. Training takes ~10–20 seconds in the browser.

# Logistic Regression on the Pima Indians Diabetes Dataset
# Fetches the CSV directly from GitHub inside the browser via Pyodide.
from random import seed, randrange
from csv import reader
from math import exp
from io import StringIO
from pyodide.http import pyfetch

DATA_URL = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/pima-indians-diabetes.csv"

async def load_csv_from_url(url):
    resp = await pyfetch(url)
    text = await resp.string()
    return [row for row in reader(StringIO(text)) if row]

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 accuracy_metric(actual, predicted):
    correct = sum(1 for a, p in zip(actual, predicted) if a == p)
    return correct / float(len(actual)) * 100.0

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

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

def logistic_regression(train, test, l_rate, n_epoch):
    coef = coefficients_sgd(train, l_rate, n_epoch)
    return [round(predict(row, coef)) for row in test]

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(accuracy_metric(actual, predicted))
    return scores

seed(1)
print("Fetching dataset…")
dataset = await load_csv_from_url(DATA_URL)
for i in range(len(dataset[0])):
    str_column_to_float(dataset, i)
minmax = dataset_minmax(dataset)
normalize_dataset(dataset, minmax)
print("Loaded", len(dataset), "rows. Training…")

scores = evaluate_algorithm(dataset, logistic_regression, 5, 0.1, 100)
print('Scores: %s' % scores)
print('Mean Accuracy: %.3f%%' % (sum(scores) / float(len(scores))))

Expect a mean accuracy around 77% — well above the 65% Zero-Rule baseline.

Watching the coefficients learn

Rather than only printing a final accuracy, we can log sum_error and the full coef vector at the end of every epoch and plot them. The instrumented loop is a two-line change to coefficients_sgd:

history = []
for epoch in range(n_epoch):
    sum_error = 0.0
    for row in train:
        yhat = predict(row, coef)
        error = row[-1] - yhat
        sum_error += error ** 2
        coef[0] += l_rate * error * yhat * (1.0 - yhat)
        for i in range(len(row) - 1):
            coef[i + 1] += l_rate * error * yhat * (1.0 - yhat) * row[i]
    history.append((epoch, sum_error, list(coef)))   # <- record trajectory

Running that on 600 training rows of the diabetes dataset with η = 0.1 for 100 epochs produces the two curves below. The top chart is the sum-of-squared-error per epoch — a smooth, monotonic descent that flattens as the model approaches its best fit. The bottom chart tracks all nine coefficients as they migrate away from zero.

A few things pop out. Glucose and BMI climb fastest and finish with the largest positive weights — high values push the sigmoid toward a "diabetic" prediction, which matches clinical intuition. The bias drifts strongly negative to offset that, keeping the baseline probability low for an average patient. Blood pressure and skin thickness barely move: after normalisation their signal is weak, so the gradient rarely pushes them. This is exactly what "the model is learning" looks like at the level of individual parameters.

Extensions

• Tune the learning rate, epochs, or data preparation.
• Switch stochastic updates for batch updates accumulated across an epoch.
• Apply the technique to other binary classification problems from the UCI ML repository.