43 - Logistic regression: from a line to a probability

A doctor does not say "the patient has a 100% or a 0% probability of being ill": they reason in degrees of confidence, and only at the end decide whether to intervene. Logistic regression does the same. It takes a score, built as a weighted sum of the features, and turns it into a probability between 0 and 1; then, if a decision is needed, the probability is compared with a threshold. Despite the name it is a classification algorithm, and it matters for two reasons: it is the most widely used linear classifier in practice, and it is the neuron with which you will build neural networks in the next lessons. Understanding it deeply, with its cost function and its gradient, means already holding the heart of how a network is trained.

Why linear regression is not enough

With two classes coded 0 and 1 you might want to use linear regression and say "positive if the prediction exceeds 0.5". Two problems, which we can measure. The first: the line is unbounded, so it produces "probabilities" below 0 and above 1. The second: a very distant but easily classified point shifts the line and with it the decision threshold. Let us check on two groups of points separated along one axis:

A note on the code: the Python blocks in this article are meant to be run in order, in a single script or notebook. Each one reuses the imports and variables defined in the previous blocks. All of them were run twice (Python 3.14.4, NumPy 2.3.5, scikit-learn 1.8.0) with identical output; with scikit-learn 1.8.0 the lines with C=np.inf print a UserWarning on stderr («Setting penalty=None will ignore the C and l1_ratio parameters»), harmless, which the output reported below does not include.

import numpy as np
from sklearn.linear_model import LinearRegression, LogisticRegression

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(2, 1, 50), rng.normal(5, 1, 50)])        # 50 points of class 0, 50 of class 1
y = np.concatenate([np.zeros(50), np.ones(50)])

def linear_threshold(x, y):
    lr = LinearRegression().fit(x[:, None], y)
    return lr, -(lr.intercept_ - 0.5) / lr.coef_[0]                      # x where the line equals 0.5

lr, s0 = linear_threshold(x, y)
print("line predictions at x=-2 and x=9: %s  (they are not probabilities)" % np.round(lr.predict([[-2], [9]]), 2))

x_extra, y_extra = np.append(x, 40.0), np.append(y, 1.0)                  # a single class-1 point, very far away
_, s1 = linear_threshold(x_extra, y_extra)
log0 = LogisticRegression(C=np.inf).fit(x[:, None], y)
log1 = LogisticRegression(C=np.inf).fit(x_extra[:, None], y_extra)
print("decision threshold, linear regression: %.3f -> %.3f after adding the distant point" % (s0, s1))
print("decision threshold, logistic regression: %.4f -> %.4f" % (-log0.intercept_[0] / log0.coef_[0, 0], -log1.intercept_[0] / log1.coef_[0, 0]))

line predictions at x=-2 and x=9: [-0.84  1.8 ]  (they are not probabilities)
decision threshold, linear regression: 3.581 -> 3.853 after adding the distant point
decision threshold, logistic regression: 3.5713 -> 3.5713

What the code does. It generates two clouds on a single dimension (class 0 centred at 2, class 1 at 5). For linear regression it computes the threshold as the value of x where the line equals 0.5. Then it adds a single class-1 point at x = 40, perfectly consistent (it is on the right side). For linear regression the threshold moves from 3.58 to 3.85 just because of a point that should have changed nothing, and the predictions at x = −2 and x = 9 are −0.84 and 1.80: they are not probabilities. With logistic regression the threshold stays at 3.571 in both cases. We need a function that "squashes" any score into the interval (0, 1).

The sigmoid function

The sigmoid (or logistic function) turns any real number z into a value between 0 and 1:

σ(z) = 1 / (1 + e^(−z))

What it does: for very negative z, e^(−z) is huge and σ(z) is close to 0; for very positive z, e^(−z) is close to 0 and σ(z) is close to 1; for z = 0 it is exactly 0.5. The shape is a smooth "S". Two properties will be needed shortly:

σ(−z) = 1 − σ(z)                (symmetry)
σ'(z) = σ(z) · (1 − σ(z))       (derivative)

The derivative is proved with the chain rule. Write σ(z) = (1 + e^(−z))⁻¹ and differentiate "from the outside in": σ'(z) = −(1 + e^(−z))⁻² · (−e^(−z)) = e^(−z) / (1 + e^(−z))². Now split the ratio into two factors: [1 / (1 + e^(−z))] · [e^(−z) / (1 + e^(−z))] = σ(z) · [1 − σ(z)], where the last equality follows from e^(−z)/(1 + e^(−z)) = 1 − 1/(1 + e^(−z)). The derivative is maximal at z = 0 (it equals 0.25) and tends to zero in the tails: a fact with important consequences in deep networks.

The model: score, probability, boundary

The model builds the score as in linear regression, z = θ·x (with x₀ = 1 for the intercept), and passes it through the sigmoid:

p = P(y = 1 | x) = σ(θ · x) = 1 / (1 + e^(−θ·x))

You decide "class 1" if p ≥ 0.5. Since σ(z) ≥ 0.5 if and only if z ≥ 0, the decision rule is simply θ·x ≥ 0: the decision boundary is the hyperplane θ·x = 0, linear. The model is therefore a linear classifier that in addition provides a graded probability: the farther a point is from the boundary, the closer p gets to 0 or 1.

On the left the sigmoid turns the score into a probability, with p = 0.5 at z = 0. On the right, in the feature plane, the boundary θ·x = 0 is a line; the weight vector (the orange arrow) is perpendicular to it and points towards class 1.

Interpretation: odds and log-odds

The odds of an event with probability p are p/(1 − p): if p = 0.75, the odds are 3 to 1. Inverting the sigmoid shows that θ·x is the logarithm of the odds (log-odds or logit):

log( p / (1 − p) ) = θ · x

So the model is linear in the logit. This gives the coefficients an interpretation: increasing the feature x_j by 1 (the others held fixed) increases the log-odds by θ_j, i.e. multiplies the odds by e^(θ_j). A coefficient of 0.7 doubles the odds (e^0.7 ≈ 2). Careful: it does not mean "increases the probability by 0.7": the relation between features and probability is non-linear, and the effect on the probability depends on where you are on the curve.

The cost function: cross-entropy, derived from the likelihood

As for linear regression, we do not choose the cost at random: we derive it from a probabilistic assumption. Each label is a coin toss with probability of heads p_i = σ(θ·x_i). The probability of observing the label y_i ∈ {0, 1} is written compactly as:

P(y_i | x_i; θ) = p_i^(y_i) · (1 − p_i)^(1 − y_i)

(for y_i = 1 only p_i remains, for y_i = 0 only 1 − p_i). The likelihood of the whole dataset is the product over all examples; the logarithm turns it into a sum. Changing the sign and taking the mean, so as to minimize, gives the log-loss or binary cross-entropy:

J(θ) = −(1/m) Σ_i [ y_i · log(p_i) + (1 − y_i) · log(1 − p_i) ]

It has an intuitive reading: it severely penalizes confident errors. If y = 1 and the model gives p = 0.01, the cost of that example is −log(0.01) ≈ 4.6; if p = 0.4 it is −log(0.4) ≈ 0.92; if p = 0.99 it is ≈ 0.01. A model that says "I am 99% sure" and is wrong pays a lot.

Why not squared error? It can be used, but combined with the sigmoid it produces a non-convex cost surface, with flat areas where the gradient is tiny: the descent struggles to get started. Cross-entropy, instead, is convex in the parameters for logistic regression, and is backed by a probabilistic derivation.

The gradient: the same as linear regression

Let us compute ∂J/∂θ with the chain rule. For a single example, call ℓ = −[y log p + (1 − y) log(1 − p)], with p = σ(z) and z = θ·x. Three derivatives, from the outside in:

∂ℓ/∂p = −y/p + (1 − y)/(1 − p)
∂p/∂z = p (1 − p)                      (derivative of the sigmoid)
∂z/∂θ_j = x_j

∂ℓ/∂z = [ −y/p + (1 − y)/(1 − p) ] · p(1 − p)
      = −y(1 − p) + (1 − y)p
      = p − y

The working: multiplying, the term −y/p by p(1 − p) gives −y(1 − p); the term (1 − y)/(1 − p) by p(1 − p) gives (1 − y)p; expanding −y + yp + p − yp = p − y. A wonderful result: the derivative of the cost with respect to the score is simply prediction minus label. Multiplying by ∂z/∂θ_j = x_j and averaging over the examples:

∂J/∂θ_j = (1/m) Σ_i (p_i − y_i) · x_ij        →        ∇J = (1/m) Xᵀ(p − y)

It is identical in form to the gradient of linear regression, with the prediction replaced by p = σ(Xθ). This "coincidence" is no accident: it is a property of generalized linear models with the right link function, and you will find it again in the last layer of a neural network.

Hand-written implementation

We apply gradient descent, with standardized features. The dataset is two-dimensional so we can reason about the boundary:

from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X2, y2 = make_classification(n_samples=600, n_features=2, n_redundant=0, n_informative=2,
                             n_clusters_per_class=1, class_sep=1.0, flip_y=0.03, random_state=3)
Xtr, Xte, ytr, yte = train_test_split(X2, y2, test_size=0.3, stratify=y2, random_state=0)

sc = StandardScaler().fit(Xtr)                                            # mean and std from the training set ONLY
A = np.column_stack([np.ones(len(Xtr)), sc.transform(Xtr)])               # (420, 3): ones + 2 features
B = np.column_stack([np.ones(len(Xte)), sc.transform(Xte)])               # (180, 3)

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-np.clip(z, -500, 500)))                  # clip: avoids exp overflow for very negative z

def logloss(p, y):
    p = np.clip(p, 1e-12, 1 - 1e-12)                                      # clip: never computes log(0)
    return -(y * np.log(p) + (1 - y) * np.log(1 - p)).mean()

def train_logistic(A, y, alpha=0.5, iterations=2000):
    theta = np.zeros(A.shape[1]); history = []
    for _ in range(iterations):
        p = sigmoid(A @ theta)
        theta = theta - alpha * A.T @ (p - y) / len(y)                    # gradient: (1/m) X^T (p - y)
        history.append(logloss(sigmoid(A @ theta), y))
    return theta, history

theta, history = train_logistic(A, ytr)
sk = LogisticRegression(C=np.inf).fit(A[:, 1:], ytr)

print("theta by hand:   ", np.round(theta, 4))
print("theta sklearn:   ", np.round(np.r_[sk.intercept_, sk.coef_[0]], 4))
print("log-loss with theta = 0: %.4f (ln 2 = %.4f);  after 100 iterations %.4f;  final %.4f" % (logloss(sigmoid(A @ np.zeros(3)), ytr), np.log(2), history[99], history[-1]))
print("test accuracy: by hand %.3f   sklearn %.3f" % (((sigmoid(B @ theta) >= 0.5) == yte).mean(), sk.score(B[:, 1:], yte)))

theta by hand:    [ 0.387  -3.7381  0.3253]
theta sklearn:    [ 0.3869 -3.7379  0.3253]
log-loss with theta = 0: 0.6931 (ln 2 = 0.6931);  after 100 iterations 0.2553;  final 0.2511
test accuracy: by hand 0.939   sklearn 0.939

What the code does. sigmoid uses np.clip to avoid overflow in np.exp with very negative scores (with SciPy there is scipy.special.expit, numerically stable, documentation). logloss bounds the probabilities to [10⁻¹², 1 − 10⁻¹²] so as never to compute log(0), which would be −∞. The training loop is the one of gradient descent for linear regression, with a single change: the "prediction" is sigmoid(A @ theta). At θ = 0 all probabilities equal 0.5 and the initial log-loss is ln 2 = 0.6931 (the model says "I don't know" about everything): a useful sanity check for every binary classifier. The coefficients match scikit-learn's (C=np.inf switches off regularization, which we discuss in the next lesson; the documentation of LogisticRegression says that C=np.inf is equivalent to the old penalty=None, deprecated since version 1.8), and so does the accuracy.

On the standardized features θ ≈ (0.387, −3.738, 0.325): the first feature has a strong negative weight, the second a weak one. In standardized units, increasing the first feature by one standard deviation multiplies the odds by e^(−3.74) ≈ 0.024, i.e. reduces them about 40-fold.

More than two classes: softmax

With K classes the model produces K scores z_k = θ_k·x (one per class) and turns them into probabilities with the softmax function:

p_k = e^(z_k) / Σ_j e^(z_j)

The probabilities are positive and sum to 1. With K = 2 softmax coincides with the sigmoid (the difference of the two scores plays the role of z). The loss is the categorical cross-entropy, J = −(1/m) Σ_i log p_{i, y_i}, and its gradient again has the "prediction minus label" form: ∇_W J = (1/m) Xᵀ(P − Y), where Y is the matrix of labels in one-hot encoding (one row per example, with 1 in the column of the true class and 0 elsewhere). On Iris (3 species):

from sklearn.datasets import load_iris

Xi, yi = load_iris(return_X_y=True)
Xa, Xb, ya, yb = train_test_split(Xi, yi, test_size=0.3, stratify=yi, random_state=0)
s2 = StandardScaler().fit(Xa)
Aa = np.column_stack([np.ones(len(Xa)), s2.transform(Xa)])                # (105, 5)
Ab = np.column_stack([np.ones(len(Xb)), s2.transform(Xb)])                # (45, 5)

def softmax(Z):
    Z = Z - Z.max(axis=1, keepdims=True)                                  # subtracting the max does not change the result and avoids overflow
    E = np.exp(Z)
    return E / E.sum(axis=1, keepdims=True)

Y = np.eye(3)[ya]                                                         # one-hot: (105, 3)
W = np.zeros((5, 3))                                                      # one weight vector per class
for _ in range(3000):
    P = softmax(Aa @ W)                                                   # (105, 3)
    W = W - 0.3 * Aa.T @ (P - Y) / len(ya)                                # gradient: (1/m) X^T (P - Y)

lib = LogisticRegression(C=np.inf, max_iter=2000).fit(Aa[:, 1:], ya)
print("shape of W:", W.shape, " rows of P sum to:", np.round(softmax(Ab @ W).sum(axis=1)[:3], 3))
print("training cross-entropy: %.4f" % -np.log(softmax(Aa @ W)[np.arange(len(ya)), ya]).mean())
print("test accuracy: by hand %.3f   sklearn %.3f" % ((softmax(Ab @ W).argmax(axis=1) == yb).mean(), lib.score(Ab[:, 1:], yb)))

shape of W: (5, 3)  rows of P sum to: [1. 1. 1.]
training cross-entropy: 0.0646
test accuracy: by hand 1.000   sklearn 1.000

What the code does. np.eye(3)[ya] builds the one-hot matrix by choosing the row of the identity corresponding to the label. softmax subtracts the maximum of each row before exponentiating: mathematically it does not change the result (numerator and denominator are multiplied by the same factor) but it keeps e^(z) from overflowing. W has shape (5, 3): 5 coefficients (intercept + 4 features) for each of the 3 classes, so Aa @ W has shape (105, 3), one score per example and per class. The gradient is the same line as before, with P and Y matrices instead of p and y vectors. Shapes are the key: (5, 105) @ (105, 3) = (5, 3), the same shape as W.

The perfect separation trap

If the classes are perfectly separable by a linear boundary, the log-loss can always decrease by making the weights larger: the larger they are, the closer the probabilities get to 0 and 1, the more the cost falls towards zero. There is no finite minimum: training never converges. Let us check:

xs = np.concatenate([rng.normal(-3, 0.5, 50), rng.normal(3, 0.5, 50)])    # two well-separated groups
ys = np.concatenate([np.zeros(50), np.ones(50)])
As = np.column_stack([np.ones(100), xs])
for it in (10, 100, 1000, 10000):
    t, st = train_logistic(As, ys, alpha=0.5, iterations=it)
    print("after %5d iterations: theta = %s   log-loss = %.5f" % (it, np.round(t, 3), st[-1]))

after    10 iterations: theta = [0.006 1.336]   log-loss = 0.02277
after   100 iterations: theta = [0.025 2.132]   log-loss = 0.00288
after  1000 iterations: theta = [0.068 3.01 ]   log-loss = 0.00033
after 10000 iterations: theta = [0.145 3.941]   log-loss = 0.00004

The weight of the feature keeps growing, slowly but without limit (about 1.3 → 2.1 → 3.0 → 3.9 as the iterations are multiplied by ten), while the log-loss tends to zero. In practice, with nearly separable classes, weights become huge, and probabilities extreme and poorly calibrated. It is one of the main reasons for regularization, which penalizes large weights: the topic of the next lesson.

Trade-offs and traps

• Pros: simple, fast, graded probabilities, interpretable coefficients (log-odds), good baseline, convex (a single minimum).
• Cons: linear boundary (it does not capture interactions unless you build features); sensitive to correlated features; requires scaled features for gradient descent; probabilities are reasonable only if the model is well specified.
• The 0.5 threshold is not sacred: with imbalanced classes or asymmetric costs choose the threshold from the precision-recall curve (see the process of an ML project: splits, leakage and metrics). Using class_weight="balanced" is another lever.
• Accuracy on imbalanced classes: it lies, as the lesson on splits, leakage and metrics explains. Look at precision, recall and log-loss.
• MSE as a classification cost: it works badly with the sigmoid (almost null gradient when the model is confident and wrong).
• Interpreting a coefficient as a change in probability: it is a change in log-odds.
• Overflow and log(0): use stable functions (expit, clip, softmax with max subtraction).
• Shape mismatch in softmax: Z.max(axis=1) without keepdims=True produces an (m,) vector that does not subtract from an (m, K) matrix as expected (wrong broadcast or error).

The eight items fall into three families. The first two (pros and cons) tell you when to choose the model: as a linear baseline with interpretable probabilities. The threshold, the accuracy on imbalanced classes and the wrong cost concern evaluation: the model can be right and the way of judging it wrong. The last three (log-odds read as probabilities, overflow, the shape of softmax) are implementation errors that raise no exception but produce wrong numbers, and this is why the experiments above always compare the hand-written result with scikit-learn's.

Reading the coefficients: odds ratio with a confidence interval

In logistic regression the coefficient w_j is the change in the log-odds for a unit increase of feature j. Exponentiating makes it more readable: exp(w_j) is the odds ratio, the factor by which the odds (p / (1 − p)) are multiplied when the feature increases by one unit, all else equal. An odds ratio of 2 means "doubled odds", one of 0.5 "halved odds", one of 1 "no effect". To know how far to trust an odds ratio we need an interval; here we obtain it with the bootstrap, resampling the points and refitting the model. We simulate 2,000 patients with two standardized variables (age and dose), whose truth is log-odds = −1.0 + 0.8·age − 1.2·dose.

import warnings; warnings.filterwarnings("ignore")
from sklearn.datasets import make_circles
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score

rng_l = np.random.default_rng(0)
n_l = 2000
age, dose = rng_l.normal(0, 1, n_l), rng_l.normal(0, 1, n_l)
prob = 1 / (1 + np.exp(-(-1.0 + 0.8 * age - 1.2 * dose)))
outcome = (rng_l.uniform(size=n_l) < prob).astype(int)
X_l = np.column_stack([age, dose])

fit = LogisticRegression(C=np.inf).fit(X_l, outcome)                     # no regularization
boot = np.array([LogisticRegression(C=np.inf).fit(X_l[i], outcome[i]).coef_[0] for i in (np.random.default_rng(s).integers(0, n_l, n_l) for s in range(300))])
low, high = np.exp(np.percentile(boot, [2.5, 97.5], axis=0))
print("outcome frequency: %.3f" % outcome.mean())
for name, j, true_value in (("age", 0, 0.8), ("dose", 1, -1.2)):
    print("%-5s coefficient %+.3f (%s %+.1f)   odds ratio %.3f   95%% interval: [%.3f, %.3f]   true odds ratio %.3f" % (name, fit.coef_[0, j], "true", true_value, np.exp(fit.coef_[0, j]), low[j], high[j], np.exp(true_value)))

outcome frequency: 0.325
age   coefficient +0.792 (true +0.8)   odds ratio 2.207   95% interval: [1.977, 2.505]   true odds ratio 2.226
dose  coefficient -1.215 (true -1.2)   odds ratio 0.297   95% interval: [0.255, 0.337]   true odds ratio 0.301

What the code does. It generates 2,000 patients with a binary outcome (frequency 0.325), fits a logistic regression without regularization and then repeats it on 300 bootstrap resamples. Age has coefficient +0.792 (true +0.8), i.e. an odds ratio of exp(0.792) = 2.207: for each extra standard deviation of age the odds of the outcome more than double. Dose has coefficient −1.215 (true −1.2), odds ratio 0.297: each extra standard deviation divides the odds by about 3.4. The 95% bootstrap intervals are [1.977, 2.505] for age and [0.255, 0.337] for dose, and both contain the true values (2.226 and 0.301).

Two interpretation pitfalls. First: the odds ratio is not a relative risk. When the outcome is rare (small probabilities) odds and probability almost coincide and the two readings are close; here the outcome concerns a third of the cases and they are not, so "doubled odds" does not mean "doubled probability". Second: the coefficients comparable with each other are those of features on the same scale (we standardized them); if age were in years and dose in milligrams, the two numbers could not be compared.

Finally, an odds ratio describes an association conditional on the other variables in the model, not a causal effect: saying so requires a study design, not a coefficient.

How the library really trains: Newton's method (IRLS)

In the lesson on gradient descent we saw that Newton's method corrects the gradient with the inverse of the Hessian. For logistic regression the cost function is convex but not quadratic, so one step is not enough; a few steps are, though. The gradient of the log-likelihood is Xᵀ(y − p) and the Hessian (sign changed) is Xᵀ W X, with W the diagonal matrix with elements p_i (1 − p_i). The update w ← w + (XᵀWX)⁻¹ Xᵀ(y − p) is called IRLS (iteratively reweighted least squares), because each step is a weighted least-squares regression where the weights depend on the current probabilities. Let us compare it with plain gradient descent.

X_b = np.column_stack([np.ones(n_l), X_l])
w_newton = np.zeros(3)
print("iteration   log-likelihood   largest step")
for iteration in range(7):
    p = 1 / (1 + np.exp(-X_b @ w_newton))
    log_lik = np.sum(outcome * np.log(p) + (1 - outcome) * np.log(1 - p))
    weights = p * (1 - p)
    step = np.linalg.solve(X_b.T @ (X_b * weights[:, None]), X_b.T @ (outcome - p))
    w_newton = w_newton + step
    print("%-11d   %-19.4f   %.2e" % (iteration, log_lik, np.abs(step).max()))
print("Newton weights: %s     scikit-learn weights: %s" % (np.round(w_newton, 4), np.round(np.r_[fit.intercept_, fit.coef_[0]], 4)))

w_gd = np.zeros(3)
for iteration in range(1, 100001):
    w_gd = w_gd + 0.5 * X_b.T @ (outcome - 1 / (1 + np.exp(-X_b @ w_gd))) / n_l
    if np.abs(w_gd - w_newton).max() < 1e-6:
        break
print("gradient descent (learning rate 0.5): %d iterations to get within 1e-6" % iteration)

iteration   log-likelihood   largest step
0             -1386.2944            7.69e-01
1             -1009.3351            3.43e-01
2             -976.6892             9.69e-02
3             -975.1331             6.14e-03
4             -975.1277             2.26e-05
5             -975.1277             3.06e-10
6             -975.1277             1.40e-16
Newton weights: [-1.0035  0.7916 -1.2151]     scikit-learn weights: [-1.0038  0.7916 -1.2149]
gradient descent (learning rate 0.5): 352 iterations to get within 1e-6

What the code does. It implements Newton's method for logistic regression by hand. At iteration 0, with all weights at zero, every probability is 0.5 and the log-likelihood is −1,386.29 (which is exactly 2,000 · ln 0.5, the value of a model that guesses at random). Then it rises to −1,009.34, −976.69, −975.13 and stops at −975.1277 at iteration 4: from there the largest step falls to 10⁻¹⁰ and 10⁻¹⁶, i.e. rounding level.

This is Newton's quadratic convergence: the largest step goes from 0.77 to 0.34, 0.097, 0.0061, 2.3·10⁻⁵, 3·10⁻¹⁰, and at each iteration the number of correct digits roughly doubles. The weights obtained (−1.0035, 0.7916, −1.2151) coincide with scikit-learn's (−1.0038, 0.7916, −1.2149) up to differences of about 3·10⁻⁴ on the intercept and 2·10⁻⁴ on the dose weight, due to scikit-learn's stopping tolerance.

Gradient descent with learning rate 0.5, which uses only the gradient, needs 352 iterations to get within 10⁻⁶: almost sixty times more than Newton (6 iterations to fall below 10⁻⁶).

Each Newton step, however, costs more than a gradient descent step because it forms and solves a p × p system. With three weights it is irrelevant; with thousands of features scikit-learn moves to intermediate methods (such as L-BFGS, which approximates the Hessian from past gradients) and with a lot of data to stochastic methods.

Are the probabilities calibrated? And regularization makes them honest

A classifier with high accuracy can give unreliable probabilities. Calibration (definition in the scikit-learn guide) asks: among the cases to which the model assigns probability 0.8, are 80% really positive? It is checked by grouping the predictions into bands and comparing the mean predicted probability with the observed frequency. A correctly specified model on abundant data is typically well calibrated; one trained on little data with many features tends instead to be overconfident. We check both cases.

X_a, X_b_test, y_a, y_b = train_test_split(X_l, outcome, test_size=0.5, random_state=0)
model = LogisticRegression().fit(X_a, y_a)
predicted = model.predict_proba(X_b_test)[:, 1]
print("probability band   cases   predicted (mean)   observed (frequency)")
edges = np.linspace(0, 1, 6)
for lo, hi in zip(edges[:-1], edges[1:]):
    mask = (predicted >= lo) & (predicted < hi + (hi == 1))
    print("%.1f - %.1f              %-5d  %-17.3f  %.3f" % (lo, hi, mask.sum(), predicted[mask].mean(), y_b[mask].mean()))

X_wide = rng_l.normal(size=(300, 100))                                      # 300 examples, 100 features, only one informative and lots of noise
y_wide = (X_wide[:, 0] + rng_l.normal(0, 2, 300) > 0).astype(int)
Xw_a, Xw_b, yw_a, yw_b = train_test_split(X_wide, y_wide, test_size=0.5, random_state=0)
print("\nC (inverse regularization)   train accuracy   test accuracy   mean confidence   test log-loss")
for C in (np.inf, 1, 0.01):
    wide = LogisticRegression(C=C, max_iter=5000).fit(Xw_a, yw_a)
    p_test = wide.predict_proba(Xw_b)[:, 1]
    print("%-35s %-17.3f   %-15.3f   %-15.3f    %.3f" % (C, wide.score(Xw_a, yw_a), wide.score(Xw_b, yw_b), (np.abs(p_test - 0.5) * 2).mean(), -np.mean(yw_b * np.log(p_test + 1e-12) + (1 - yw_b) * np.log(1 - p_test + 1e-12))))

probability band   cases   predicted (mean)   observed (frequency)
0.0 - 0.2              373    0.103              0.094
0.2 - 0.4              268    0.297              0.250
0.4 - 0.6              185    0.501              0.546
0.6 - 0.8              122    0.700              0.648
0.8 - 1.0              52     0.864              0.942

C (inverse regularization)   train accuracy   test accuracy   mean confidence   test log-loss
inf                                 1.000               0.520             0.966              9.775
1                                   0.993               0.520             0.809              2.076
0.01                                0.833               0.560             0.166              0.699

What the code does. The first table divides the predicted probabilities on the test (1,000 cases) into five bands and compares the predicted mean with the observed frequency. The model is well calibrated overall: in the 0.0–0.2 band it predicts 0.103 and observes 0.094; in 0.4–0.6 it predicts 0.501 and observes 0.546; in 0.8–1.0 it predicts 0.864 and observes 0.942. The gaps reach 7.8 points in the last band, but there the cases are only 52, and the standard error of a frequency observed on 52 cases is about √(0.94·0.06/52) ≈ 3.2 points: the gap is about 2.4 standard errors, a borderline case. The lesson is that a calibration table with a few thousand cases and five bands can only be read roughly.

The second table is the serious warning: 300 examples, 100 features and only one really informative (with a lot of noise), 150 for training and 150 for testing. Without regularization (infinite C) the model has accuracy 1.000 on training and 0.520 on test (almost a coin flip), but a mean confidence of 0.966: it is very sure of answers that are almost random, and the test log-loss is 9.775 (a constant model that always says 0.5 would have 0.693).

With C = 1 things improve but remain ugly (log-loss 2.076, confidence 0.809). With C = 0.01 the confidence falls to 0.166 and the log-loss to 0.699: the probabilities are honest, but we are practically at the level of the constant model, and the accuracy of 0.560 on 150 cases is not distinguishable from chance.

The signal, in this experiment, is too weak for the available features. It is a negative result to report as such: regularization does not create information, it only avoids pretending to have it.

When the linear boundary is not enough: concentric circles

X_c, y_c = make_circles(600, noise=0.1, factor=0.4, random_state=0)         # one circle inside another
linear = LogisticRegression()
quadratic = make_pipeline(PolynomialFeatures(2), LogisticRegression(C=10, max_iter=5000))
print("CV accuracy (5 folds)   linear boundary: %.3f   with second-degree features: %.3f" % (cross_val_score(linear, X_c, y_c, cv=5).mean(), cross_val_score(quadratic, X_c, y_c, cv=5).mean()))
fit_q = quadratic.fit(X_c, y_c)
print("features:      %s" % str(list(fit_q[0].get_feature_names_out())))
print("coefficients: %s" % str(np.round(fit_q[-1].coef_[0], 2)))

CV accuracy (5 folds)   linear boundary: 0.483   with second-degree features: 0.993
features:      ['1', 'x0', 'x1', 'x0^2', 'x0 x1', 'x1^2']
coefficients: [-6.000e-02  1.000e-02  1.400e-01 -1.250e+01 -6.900e-01 -1.249e+01]

What the code does. It generates 600 points on two concentric circles (one of small radius inside a larger one) and compares two models with cross-validation. Logistic regression on the two coordinates alone reaches an accuracy of 0.483: slightly worse than chance, because the boundary is a line and there is no line separating a circle from the one containing it (it is the same limit as XOR for the perceptron). With second-degree features (1, x₀, x₁, x₀², x₀x₁, x₁²) the accuracy rises to 0.993, with the same training algorithm: only the space in which the model is linear changes.

The coefficients tell why: the terms x₀² and x₁² are −12.50 and −12.49, practically equal, while the linear terms and the mixed term x₀x₁ are small (0.01, 0.14, −0.69). The model has learned a function of x₀² + x₁², the square of the distance from the origin, which is exactly the quantity that separates the two circles.

It is the bridge to neural networks: instead of choosing by hand the features to add, a network learns them from the data with the hidden layers. At the same time it shows the limit of the manual method: it works because we guessed that degree 2 was enough, which in high-dimensional problems cannot be done by eye.

Three exercises

Exercise 1. In the calibration bands experiment, reduce the sample to 200 patients: do the bands stay close to the diagonal? How does the variability between bands change and why does a band with 10 cases say nothing? Exercise 2. In the 100-feature experiment, use cross-validation to choose C on a grid from 10⁻³ to 10²: is the chosen C the one with the lowest log-loss or the one with the highest accuracy, and do the two choices coincide? Exercise 3. On the concentric circles try degree 3 and degree 6 with the same C: does the accuracy improve? Look at the coefficients to understand which terms the regression really finds useful.

Previous lesson in the course: Multivariate and polynomial regression, and feature scaling. Other lessons in the series that link to this one: Entropy, Cross-Entropy and KL Divergence: Why Loss Functions Are All the Same Thing, Naive Bayes: Bayes' theorem with a naive assumption, Boosting: AdaBoost, Gradient Boosting and a note on XGBoost, Dimensionality reduction: PCA from scratch and t-SNE.

Sources

• scikit-learn, LogisticRegression - the C and penalty parameters (C=np.inf as unpenalized regression, penalty deprecated since 1.8), consulted on 2026-11-12
• scikit-learn, Linear Models (Logistic regression section) - logistic regression as a linear model for classification, binary, one-vs-rest and multinomial, consulted on 2026-11-12
• scikit-learn, Probability calibration - the definition of a well-calibrated classifier used in the calibration section, consulted on 2026-11-12
• SciPy, scipy.special.expit - the numerically stable logistic sigmoid, consulted on 2026-11-12

The figures in the lesson are the author's measurements with a fixed seed. Historical reference not reopened during verification: D. R. Cox, The regression analysis of binary sequences, Journal of the Royal Statistical Society B, 1958, for the logistic model.

In the next lesson you see why a model with too much freedom memorizes noise, and how regularization reins it in, in overfitting and regularization: L2, L1 and the accuracy trap.