Logistic regression

What this chapter covers

Back in the threshold chapter we drew one vertical line through one feature and called it a classifier. It worked, it was honest, and it left seven of the eight columns on the floor. Logistic regression is the natural next move: keep the straight boundary, but let every feature vote, and instead of a hard yes/no at the line, hand back a probability. It's the model I reach for first on any linear classification problem, and most of the time it's the one I ship.

We build it from scratch in NumPy — the sigmoid, the linear score it wraps, the cross-entropy loss, the gradient, and the gradient descent loop that ties them together. Then we let scikit-learn fit the same model and check that the numbers line up. Along the way you get to watch the training happen: a decision boundary sliding and rotating across a 2-D dataset while the loss curve drops underneath it, every frame one real gradient step. No closed form this time, no argmax over a grid — an actual optimizer walking downhill.

The results run on the same Pima diabetes set from week 1, so we get a clean callback. One feature and a threshold scored 73% back then. Eight features and a sigmoid should do better. Let's see by how much.

A bit of history

The curve at the center of all this is older than statistics as a field. Pierre François Verhulst wrote down the logistic function in 1838 to model population growth — an S-shaped curve that starts slow, accelerates, then saturates as it runs out of room. He needed something that lived between two bounds and moved smoothly from one to the other, which is exactly the shape you want for a probability.

The classification story starts a century later. Joseph Berkson, a physicist and statistician at the Mayo Clinic, coined the word "logit" in 1944 — log-odds, the thing the model is actually linear in — and spent years arguing it was a better tool than the then-dominant probit for bioassay work. He was mostly right, and the name stuck. Then in 1958 David Cox laid out logistic regression as a general method for binary outcomes, and Nelder and Wedderburn folded it into the generalized linear model framework in 1972. That lineage is why the model feels so settled: it isn't a machine learning trick that happened to work, it's a statistical model with sixty years of theory under it. When people say logistic regression is the default linear classifier, that's the weight behind the word default.

The intuition

Here's a small 2-D dataset, two classes, some overlap. Two features on the axes, color for the true label. This is the toy set we'll train on so you can actually see the boundary move.

A straight line can mostly separate these two blobs, the same way one cutoff mostly separated glucose. The difference is the line now lives in two dimensions, so it can tilt — it doesn't have to sit perpendicular to an axis. In eight dimensions it's a hyperplane you can't picture, but the idea is identical: a flat boundary, oriented by a weight per feature.

The other new idea is what happens near the line. The threshold classifier was brutal — one side sick, the other side healthy, no in-between. Logistic regression softens that. Far from the boundary on the positive side it says "almost certainly class 1"; far on the other side, "almost certainly class 0"; right at the line, "I genuinely don't know, call it a coin flip." That softness is the whole point. It gives you a probability you can threshold wherever you want, and it gives training something smooth to optimize.

The math

Start with the linear score. For a feature vector $x_i$ with $D$ entries, a weight vector $w$, and a bias $b$:

That $z_i$ is the log-odds — the logit. It runs from $-\infty$ to $+\infty$, which is no good as a probability, so we squash it with the sigmoid:

The sigmoid is the S-curve Verhulst wrote down. It maps every real number into $(0, 1)$, crosses $0.5$ at $z = 0$, and flattens toward the ends. The model's prediction is that squash applied to the score:

To classify, threshold the probability at $0.5$ — which, because the sigmoid crosses $0.5$ exactly at $z = 0$, is the same as asking whether $z_i$ is positive. The boundary is the flat surface $w^{\top} x + b = 0$.

Now training. We need a loss that's small when $p_i$ is confident and correct, large when it's confident and wrong. That's cross-entropy, also called log loss:

When the truth is $y_i = 1$ only the first term survives, and it's $-\log p_i$ — zero if you said $1.0$, blowing up as you approach $0$. When $y_i = 0$ the second term does the same in reverse. A confident correct call is nearly free; a confident wrong call is very expensive.

There's no closed form for the $w$ and $b$ that minimize this, so we descend the gradient. And here's the payoff for pairing the sigmoid with cross-entropy — the gradient collapses to something clean:

The error is just $p_i - y_i$, the gap between the predicted probability and the truth, and the gradient is that error weighted by the features. The sigmoid's derivative and the log's derivative cancel each other perfectly; nothing ugly survives. Gradient descent then repeats one line — $w \leftarrow w - \eta\, \partial\mathcal{L}/\partial w$ — until the loss stops dropping.

That cancellation is also the answer to a question people ask early: why not just use squared error, like linear regression? Two reasons. Squared error on top of a sigmoid is non-convex, so descent can stall in local dips. And it punishes confident wrong answers only mildly, because the sigmoid flattens the gradient right where you most need a strong correction. Cross-entropy is convex here and keeps pushing hard on confident mistakes. Use the loss the model was built for.

What it's good at, what it isn't

The upside is a package that's hard to beat for the price. You get calibrated probabilities, not just labels — when the model says 0.7 it means it, which matters enormously the moment a downstream decision depends on how likely, not just whether. You get weights you can read: each $w_j$ is the change in log-odds per unit of feature $j$, so the model tells you which measurements move the prediction and in which direction. It trains fast, it barely overfits with a little regularization, and it's convex, so you get the same answer every time. As a strong baseline on tabular data it's the number your gradient-boosted whatever has to actually beat, not just tie.

The limit is the same one the threshold had, lifted by one dimension: the boundary is still flat. Logistic regression draws a hyperplane, and if the true boundary curves or the classes wrap around each other, a plane can't follow it. You can buy some of that back by hand — add polynomial or interaction features and the plane bends in the original space — but at that point you're doing the feature engineering the tree-based models do for free. It also takes the linear score at face value, so wildly scaled features and strong outliers can drag the weights around. Reach past it when the signal is genuinely nonlinear and you'd rather the model find the interactions than you.

The data

Same patients as week 1: the Pima Indians Diabetes set, 768 women, eight measurements each — glucose, BMI, age, blood pressure, and so on — and a label for whether they were diagnosed with diabetes. Last time we used one column. This time all eight go in.

One preprocessing step matters enough to call out. These features live on completely different scales: glucose runs past 120, the diabetes pedigree score hovers around 0.5, pregnancy counts are single digits. Gradient descent with a single learning rate can't serve ranges that different — the step that's right for glucose is a thousand times too big for the pedigree score. So we standardize first: subtract each column's mean, divide by its standard deviation, using statistics from the training set only and applying the same shift to the test set. After that every feature is mean 0, unit variance, and one learning rate works for all of them. It also makes the weights comparable, which we'll cash in at the end.

Build it, one function at a time

Seven small pieces, bottom up, in the order you'd write them at a terminal. Start with the squash, because everything rides on it.

def sigmoid(z):
    """Squash any real number into (0, 1).

    Written in the numerically stable branch form so a large negative z
    never overflows exp(). For z >= 0 use 1 / (1 + e^-z); for z < 0 use
    e^z / (1 + e^z). Both are the same function, algebraically.
    """
    z = np.asarray(z, dtype=float)
    out = np.empty_like(z)
    pos = z >= 0
    out[pos] = 1.0 / (1.0 + np.exp(-z[pos]))
    ez = np.exp(z[~pos])
    out[~pos] = ez / (1.0 + ez)
    return out

That branchy form is deliberate. The textbook $1/(1+e^{-z})$ overflows for large negative $z$ — exp of a big positive number is inf — so we use the algebraically identical $e^{z}/(1+e^{z})$ whenever $z$ is negative. Same curve, no NaNs. It's the kind of thing that never bites you on toy data and always bites you on real data.

Next the linear score, the log-odds before squashing:

def logit(X, w, b):
    """The linear part: z = Xw + b.

    X is (N, D), w is (D,), b is a scalar bias. Returns one score z per
    row — the log-odds of class 1 before the sigmoid squashes it.
    """
    return X @ w + b

Wrap the two together and you have the model's probability output:

def predict_proba(X, w, b):
    """P(y = 1 | x) for every row: sigmoid of the linear score."""
    return sigmoid(logit(X, w, b))

To turn probabilities into labels, threshold at 0.5:

def predict_proba(X, w, b):
    """P(y = 1 | x) for every row: sigmoid of the linear score."""
    return sigmoid(logit(X, w, b))

Note the threshold is a parameter, not a magic constant. 0.5 is the natural choice, but it's the same dial as the cutoff from week 1 — slide it to trade missed positives against false alarms when one mistake costs more than the other.

Now the loss. This is what training minimizes:

def cross_entropy(X, y, w, b, eps=1e-12):
    """Mean binary cross-entropy (log loss) over the dataset.

    For each row, -log(p) when the truth is 1 and -log(1 - p) when it's 0.
    A confident, correct probability costs almost nothing; a confident,
    wrong one costs a lot. Clipping keeps log() away from 0.
    """
    p = predict_proba(X, w, b)
    p = np.clip(p, eps, 1.0 - eps)
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))

The clip keeps log away from zero when the model gets cocky and predicts exactly 0 or 1. Now the gradient — the direction that makes the loss worse, which we'll step against:

def gradient(X, y, w, b):
    """Gradient of the mean cross-entropy w.r.t. w and b.

    The algebra collapses to something clean: the error is just
    (p - y), and the gradient is that error dotted with the features.
    dw = X^T (p - y) / N, db = mean(p - y). No sigmoid derivative left
    over — cross-entropy and the sigmoid were built for each other.
    """
    p = predict_proba(X, w, b)
    err = p - y                      # shape (N,)
    dw = X.T @ err / len(y)          # shape (D,)
    db = float(np.mean(err))         # scalar
    return dw, db

Look how little is there. The whole gradient is (p - y) dotted with the features. All the calculus from the math section boils down to one subtraction and one matrix multiply, and that's not a coincidence — it's the sigmoid and cross-entropy canceling, the same cancellation that made squared error the wrong choice.

Last piece, the training loop. Start the weights at zero and step downhill:

def fit(X, y, lr=0.1, n_iter=300, record=False):
    """Train by batch gradient descent: step downhill on the log loss.

    Start the weights at zero, and on every iteration compute the
    gradient over the whole training set and take one step against it.
    With `record=True` we also return the per-iteration (w, b, loss)
    history so the chapter can animate the descent.
    """
    N, D = X.shape
    w = np.zeros(D)
    b = 0.0
    history = []
    for t in range(n_iter):
        if record:
            history.append((w.copy(), b, cross_entropy(X, y, w, b)))
        dw, db = gradient(X, y, w, b)
        w -= lr * dw
        b -= lr * db
    if record:
        history.append((w.copy(), b, cross_entropy(X, y, w, b)))
        return w, b, history
    return w, b

Batch gradient descent, the plain version: every iteration uses the whole training set to compute one gradient and take one step. The record flag stashes the weights and loss at each iteration so we can replay the descent in a minute. That replay is the whole reason the loop looks the way it does.

Watch it work

This is the training itself, running on the 2-D toy set so every point is visible. The top panel is the data; the orange line is the current decision boundary, the place where the model's probability is exactly 0.5. A point gets a red ring the moment the boundary puts it on the wrong side. The bottom panel traces the cross-entropy loss, one dot per iteration.

Press play. Every frame is one real gradient step out of fit — the same weights the code produced, no smoothing. Reset and run it again as many times as you like.

Watch the first move. The weights start at zero, so there's no line at all and the model calls everything one class — half the points wear a red ring. Then the gradient kicks in and a boundary appears, crude at first, cutting the plane roughly in half. From there it rotates and slides, and the red rings peel off in batches as it finds the orientation that separates the blobs. The loss curve underneath tells the same story from the other side: a steep drop early while the big mistakes get fixed, then a long flat glide as it fine-tunes. By the end the boundary has settled into the gap between the classes and the loss has nearly bottomed out, landing at 96% accuracy on this toy set. That flattening is what convergence looks like — the gradient shrinking toward zero because there's little left to fix.

The one thing thresholding couldn't do is right there in the tilt of the line. It isn't vertical, it isn't horizontal, it's angled — because both features carry signal and the weights found the mix. That angle is the extra dimension of freedom logistic regression buys you over a single cutoff.

If you want to see the squashing function on its own, here it is: the sigmoid across the range of scores the model produces. Flat near the ends, steep through the middle, crossing 0.5 at zero.

The steep middle band is the zone of uncertainty near the boundary; the flat tails are where the model is confident. A point far from the line has a large $|z|$ and sits out on a tail, which is why the model can be so sure about it.

The full implementation

The whole file, top to bottom — the model, the loss, the gradient, the loop, plus the standardize and accuracy helpers the traces use. This is the code the animation above actually ran:

"""Logistic regression, built from scratch.

A linear model that outputs a probability. Push the features through a
weighted sum, squash that sum to (0, 1) with the sigmoid, and call it
P(class = 1). Training is gradient descent on the cross-entropy loss — no
closed form, no library, just the gradient and a step size.

Pure NumPy (+ pandas for loading). Every function below appears in the
chapter one step at a time; the `# region:` markers are what the book's
include directives pull in.
"""

import numpy as np
import pandas as pd


# region: sigmoid
def sigmoid(z):
    """Squash any real number into (0, 1).

    Written in the numerically stable branch form so a large negative z
    never overflows exp(). For z >= 0 use 1 / (1 + e^-z); for z < 0 use
    e^z / (1 + e^z). Both are the same function, algebraically.
    """
    z = np.asarray(z, dtype=float)
    out = np.empty_like(z)
    pos = z >= 0
    out[pos] = 1.0 / (1.0 + np.exp(-z[pos]))
    ez = np.exp(z[~pos])
    out[~pos] = ez / (1.0 + ez)
    return out
# endregion


# region: logit
def logit(X, w, b):
    """The linear part: z = Xw + b.

    X is (N, D), w is (D,), b is a scalar bias. Returns one score z per
    row — the log-odds of class 1 before the sigmoid squashes it.
    """
    return X @ w + b
# endregion


# region: predict_proba
def predict_proba(X, w, b):
    """P(y = 1 | x) for every row: sigmoid of the linear score."""
    return sigmoid(logit(X, w, b))
# endregion


# region: predict
def predict(X, w, b, threshold=0.5):
    """Turn probabilities into 0/1 labels at a decision threshold.

    0.5 is the natural cutoff — predict class 1 when it's more likely than
    not — but the threshold is a knob you can move to trade the two kinds
    of mistake, exactly like the cutoff in week 1.
    """
    return (predict_proba(X, w, b) >= threshold).astype(int)
# endregion


# region: cross_entropy
def cross_entropy(X, y, w, b, eps=1e-12):
    """Mean binary cross-entropy (log loss) over the dataset.

    For each row, -log(p) when the truth is 1 and -log(1 - p) when it's 0.
    A confident, correct probability costs almost nothing; a confident,
    wrong one costs a lot. Clipping keeps log() away from 0.
    """
    p = predict_proba(X, w, b)
    p = np.clip(p, eps, 1.0 - eps)
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))
# endregion


# region: gradient
def gradient(X, y, w, b):
    """Gradient of the mean cross-entropy w.r.t. w and b.

    The algebra collapses to something clean: the error is just
    (p - y), and the gradient is that error dotted with the features.
    dw = X^T (p - y) / N, db = mean(p - y). No sigmoid derivative left
    over — cross-entropy and the sigmoid were built for each other.
    """
    p = predict_proba(X, w, b)
    err = p - y                      # shape (N,)
    dw = X.T @ err / len(y)          # shape (D,)
    db = float(np.mean(err))         # scalar
    return dw, db
# endregion


# region: fit
def fit(X, y, lr=0.1, n_iter=300, record=False):
    """Train by batch gradient descent: step downhill on the log loss.

    Start the weights at zero, and on every iteration compute the
    gradient over the whole training set and take one step against it.
    With `record=True` we also return the per-iteration (w, b, loss)
    history so the chapter can animate the descent.
    """
    N, D = X.shape
    w = np.zeros(D)
    b = 0.0
    history = []
    for t in range(n_iter):
        if record:
            history.append((w.copy(), b, cross_entropy(X, y, w, b)))
        dw, db = gradient(X, y, w, b)
        w -= lr * dw
        b -= lr * db
    if record:
        history.append((w.copy(), b, cross_entropy(X, y, w, b)))
        return w, b, history
    return w, b
# endregion


def accuracy(X, y, w, b, threshold=0.5):
    """Fraction of rows the model labels correctly at a threshold."""
    return float((predict(X, w, b, threshold) == y).mean())


def standardize(X, mean=None, std=None):
    """Center and scale each column to mean 0, unit variance.

    Gradient descent on raw features with wildly different ranges crawls;
    standardizing puts every feature on the same footing so one learning
    rate works for all of them. Returns the transformed X plus the mean
    and std, so the exact same shift can be applied to test data.
    """
    if mean is None:
        mean = X.mean(axis=0)
    if std is None:
        std = X.std(axis=0)
        std = np.where(std == 0, 1.0, std)
    return (X - mean) / std, mean, std


def load_data(path="../data/diabetes.csv"):
    """Pima Indians Diabetes dataset: 768 patients, 8 features, Outcome."""
    return pd.read_csv(path)

The library version

You wouldn't hand-roll this in production, and once you understand it you don't need to. Scikit-learn's LogisticRegression is the same model — linear score, sigmoid, cross-entropy — and it hands you the fitted weights and bias the same way:

def sklearn_logreg(Xtr, ytr, Xte, yte):
    """Fit sklearn's LogisticRegression on standardized features.

    Returns the learned weight vector, bias, and test accuracy. We hand it
    the already-standardized arrays so it's judged on the same footing as
    our from-scratch model.
    """
    clf = LogisticRegression(max_iter=1000)
    clf.fit(Xtr, ytr)
    w = clf.coef_[0]
    b = float(clf.intercept_[0])
    acc = float(clf.score(Xte, yte))
    return w, b, acc

Two differences worth knowing. First, sklearn regularizes by default: it adds an L2 penalty on the weights, which nudges them toward zero and guards against overfitting. Our from-scratch version has no penalty at all, which is fine on this data but is the first thing I'd add for anything higher-dimensional. Second, it doesn't use plain gradient descent — the default solver is lbfgs, a quasi-Newton method that uses curvature to take much smarter steps, so it converges in a few dozen iterations where our fixed-step descent takes a couple thousand. Same destination, faster road. We feed it the already-standardized arrays so the comparison is honest — same features, same scaling, same split.

Scratch versus library

Split the Pima data 70/30 — 537 patients to train on, 231 held out — standardize on the training statistics, fit both models, and score on the held-out set. Here's the face-off, with the week-1 threshold classifier dropped in as the reference bar:

Our from-scratch logistic regression scores 0.7662 on the held-out set. Scikit-learn scores 0.7662 — the same number to four decimals, with a bias of -0.888 against sklearn's -0.881, which is as close as two different optimizers on the same convex loss ever need to get. The convexity is why: there's one minimum, and both roads reach it. That agreement is the point of the exercise. When your hand-written gradient lands on the library's answer, you've understood the model, not just called it.

Now the callback. The week-1 stump — one feature, one threshold on glucose — scores 0.7316 on the same split. Logistic regression, with all eight features and the same held-out patients, scores 0.7662. Three and a half points isn't a landslide, and that's honest too: glucose really is the dominant signal in this dataset, so a single-feature model was never going to be embarrassed. But the extra features do pull real weight. Standardized, the three heaviest are glucose at +1.07, BMI at +0.85, and pregnancy count at +0.39 — all pushing toward a diabetes prediction, in an order a clinician would nod at. That's the interpretability dividend: the model didn't just classify better, it told you why, in units you can rank.

The same two caveats from week 1 still apply. The classes are lopsided — about a third of these patients are positive — so a model that guessed "no diabetes" for everyone would already clear 65%, and accuracy alone flatters every model on the chart. And the held-out split is what keeps these numbers honest; training accuracy would read higher and mean less. Logistic regression doesn't escape those traps. It just gives you a better boundary and a probability to reason about once you're inside them.

Takeaways

Logistic regression is the linear classifier I default to, and I think you should too. It's the honest next rung above the threshold: keep the straight boundary, let every feature vote, get a calibrated probability instead of a bare label. On tabular data it's the baseline a heavier model has to genuinely beat before it's earned its place — and often it doesn't, which is exactly the thing you want to find out in two seconds of fit, not two weeks into a project.

The three things to carry forward. Probabilities, not just labels: the 0.7 means 0.7, which is what lets you set a threshold to fit the cost of each mistake instead of accepting 0.5 by default. Readable weights: standardize the features and the coefficients rank your signal for free, direction and all. And the log-loss lesson, which outlives this chapter — the sigmoid and cross-entropy were built as a pair, their gradients cancel to a clean (p - y), and that's why you don't reach for squared error on a classification problem. That same cross-entropy loss and that same clean gradient are the engine underneath the neural networks later in the course; a network is mostly this model stacked and bent until the boundary can finally curve. Which is the one thing logistic regression can't do, and the reason the course keeps going. But it's where linear classification grows up, and most days it's all you need.