In the left panel the classes are two overlapping clouds, and the contours for 0.1, 0.5 and 0.9 are parallel lines: the probability depends only on the distance from the boundary along . In the middle, one class encloses the other. Every straight line does roughly as badly as guessing, and the model's test accuracy is 0.52. In the right panel the model is the same, but the features include , and . A circle is linear in those features, so the model can draw it.

Feature engineering is how logistic regression handles non-linearity: polynomial terms, interactions, logs, splines and bins. The cost is the same as in linear regression: more features mean more ways to overfit, and the need for regularisation grows with them. When you find yourself engineering many interactions by hand, that is the point to consider tree ensembles or neural networks, which learn them. A neural network classifier can be seen as logistic regression whose features are learned by the layers before it.

A real problem: breast tumour diagnosis

The deliveries were invented to be easy to compute by hand. Now for a real problem.

The Wisconsin Diagnostic Breast Cancer dataset (Wolberg, Mangasarian, Street and Street, UCI repository) was introduced by Street, Wolberg and Mangasarian (1993). Each of its 569 rows is a breast mass sampled by fine-needle aspiration. An image of the sample was digitised and the cell nuclei outlined, and ten properties of the nuclei were measured: radius, texture, perimeter, area, smoothness, compactness, concavity, concave points, symmetry and fractal dimension. Each property appears three times: its mean across the nuclei in the image, its standard error (named ... error), and the mean of the three largest values (named worst ...). That gives 30 features. The label is the diagnosis: 212 malignant, 357 benign.

One quirk matters. scikit-learn's load_breast_cancer codes malignant as 0 and benign as 1. Every metric below that depends on a "positive class" would then be about detecting benign tumours. The code flips the labels so that 1 means malignant.

The experiment:

• Split: 80% training (455 patients), 20% test (114), stratified so both sets have the same malignant fraction, random_state=42. The test set is not touched until the end.
• Model: standardise every feature using training-set statistics, then LogisticRegression(C=1.0), scikit-learn's defaults, inside a Pipeline so the scaler is fitted only on training data.
• Baseline: a DummyClassifier that predicts the training-set malignant rate, 37.4%, for every patient.
• Metrics: accuracy, log loss, Brier score, ROC-AUC and PR-AUC on the test set, explained in the evaluation section.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.dummy import DummyClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (accuracy_score, average_precision_score, brier_score_loss,
                             confusion_matrix, log_loss, roc_auc_score)
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

data = load_breast_cancer(as_frame=True)
X = data.data
y = 1 - data.target            # sklearn codes malignant as 0; make malignant the positive class
print(X.shape, "malignant:", int(y.sum()), "benign:", int((1 - y).sum()))
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)
model = make_pipeline(StandardScaler(), LogisticRegression(C=1.0))
model.fit(X_train, y_train)
baseline = DummyClassifier(strategy="prior").fit(X_train, y_train)
for name, est in [("baseline (prior)", baseline), ("logistic regression", model)]:
    proba = est.predict_proba(X_test)[:, 1]
    pred = est.predict(X_test)
    print(f"{name:>20}: accuracy={accuracy_score(y_test, pred):.3f}  "
          f"log loss={log_loss(y_test, proba):.3f}  Brier={brier_score_loss(y_test, proba):.3f}  "
          f"ROC-AUC={roc_auc_score(y_test, proba):.3f}  PR-AUC={average_precision_score(y_test, proba):.3f}")
print(confusion_matrix(y_test, model.predict(X_test)))

(569, 30) malignant: 212 benign: 357
    baseline (prior): accuracy=0.632  log loss=0.658  Brier=0.233  ROC-AUC=0.500  PR-AUC=0.368
 logistic regression: accuracy=0.965  log loss=0.077  Brier=0.021  ROC-AUC=0.996  PR-AUC=0.994
[[71  1]
 [ 3 39]]

Run with NumPy 2.5, scikit-learn 1.9 and Python 3.12.

The model gets 110 of 114 test patients right. In the confusion matrix (rows are the true class, benign then malignant; columns the prediction), the errors are one benign tumour flagged as malignant and three malignant tumours predicted benign. Those three are the errors that matter most, and the evaluation section is largely about them.

From scratch, and what the library adds

The from-scratch Newton solver above, with an L2 penalty added that leaves the intercept alone, fits the same model:

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))
def fit_logistic_newton(X, y, C=1.0, tol=1e-10, max_iter=50):
    """Minimise  sum_i log_loss_i + ||w||^2 / (2C)  by Newton's method.
X must not contain a column of ones; the intercept is handled separately
and is not penalised, matching scikit-learn.
"""
n, d = X.shape
A = np.column_stack([np.ones(n), X])
penalty = np.full(d + 1, 1.0 / C)
penalty[0] = 0.0                                  # no penalty on the intercept
theta = np.zeros(d + 1)
for it in range(1, max_iter + 1):
    p = sigmoid(A @ theta)
    grad = A.T @ (p - y) + penalty * theta
    hess = A.T @ (A * (p * (1 - p))[:, None]) + np.diag(penalty)
    step = np.linalg.solve(hess, grad)
    theta -= step
    if np.linalg.norm(step) < tol:
        break
return theta[0], theta[1:], it

data = load_breast_cancer(as_frame=True)
X, y = data.data, 1 - data.target
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)
scaler = StandardScaler().fit(X_train)
Z_train = scaler.transform(X_train)
b, w, iters = fit_logistic_newton(Z_train, y_train.to_numpy(), C=1.0)
sk = LogisticRegression(C=1.0, tol=1e-10, max_iter=1000).fit(Z_train, y_train)
print(f"Newton converged in {iters} iterations")
print("max |coef difference| :", np.abs(w - sk.coef_[0]).max())
print("intercept difference  :", abs(b - sk.intercept_[0]))

Newton converged in 10 iterations
max |coef difference| : 1.2118984182851733e-06
intercept difference  : 6.782782258074604e-08

Agreement to about , the limit of L-BFGS's stopping tolerance. Same objective, same optimum. What LogisticRegression adds:

• A choice of solvers, documented in the user guide. lbfgs (default) is a quasi-Newton method and a good general choice. newton-cholesky is the exact Newton method above, fast when there are many more examples than features. liblinear is a coordinate-descent library (Fan et al., 2008), good for small sparse problems, but it is binary-only and it penalises the intercept, which the others don't. saga is a stochastic method for very large data and the only solver that supports elastic net.
• Convergence diagnostics. A ConvergenceWarning means the solver hit max_iter before meeting its tolerance. On this dataset, fitting without the StandardScaler triggers it: the features range from under 0.001 (fractal dimension error) to over 4,000 (worst area), the problem is badly conditioned, and L-BFGS runs out of its default 100 iterations. The resulting test accuracy drops to 0.921. Scale first, as for any gradient-based method.
• Class weights, sample weights, multiclass support and the estimator interface, so the model fits into pipelines and cross-validation.

What it doesn't do is report standard errors or p-values. For inference on coefficients, statsmodels' Logit fits the unpenalised maximum-likelihood model and reports them; on this dataset it would run straight into the separation problem.

Choosing the regularisation strength

C=1.0 is a default, not a decision. Like ridge regression's , it is a hyperparameter: the training loss always prefers no penalty, so it must be chosen on data the model didn't fit, by cross-validation on the training set. The natural score is log loss, since that is what the model is supposed to get right.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss
from sklearn.model_selection import StratifiedKFold, cross_val_score, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

data = load_breast_cancer(as_frame=True)
X, y = data.data, 1 - data.target
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
for C in [0.01, 0.1, 1.0, 10.0, 100.0, np.inf]:
    pipe = make_pipeline(StandardScaler(), LogisticRegression(C=C, max_iter=10_000))
    scores = -cross_val_score(pipe, X_train, y_train, cv=cv, scoring="neg_log_loss")
    pipe.fit(X_train, y_train)
    train_acc = pipe.score(X_train, y_train)
    max_coef = np.abs(pipe[-1].coef_).max()
    print(f"C={C:>6}: CV log loss = {scores.mean():.3f} ± {scores.std():.3f}, "
          f"train accuracy = {train_acc:.3f}, max |w| = {max_coef:7.2f}")
Only now, for the two extremes, look at the test set.
for C in [1.0, np.inf]:
    pipe = make_pipeline(StandardScaler(), LogisticRegression(C=C, max_iter=10_000))
    proba = pipe.fit(X_train, y_train).predict_proba(X_test)[:, 1]
    print(f"test, C={C}: accuracy = {pipe.score(X_test, y_test):.3f}, "
          f"log loss = {log_loss(y_test, proba):.3f}")

C=  0.01: CV log loss = 0.193 ± 0.017, train accuracy = 0.954, max |w| =    0.24
C=   0.1: CV log loss = 0.103 ± 0.016, train accuracy = 0.980, max |w| =    0.61
C=   1.0: CV log loss = 0.083 ± 0.031, train accuracy = 0.987, max |w| =    1.43
C=  10.0: CV log loss = 0.126 ± 0.081, train accuracy = 0.989, max |w| =    3.39
C= 100.0: CV log loss = 0.227 ± 0.144, train accuracy = 0.993, max |w| =    8.20
C=   inf: CV log loss = 1.021 ± 0.321, train accuracy = 1.000, max |w| =  576.54
test, C=1.0: accuracy = 0.965, log loss = 0.077
test, C=inf: accuracy = 0.939, log loss = 1.969

The whole bias-variance trade-off is in this table.

• Training accuracy rises monotonically as the penalty weakens, reaching 1.000 without a penalty: the training set is separable, and the largest weight reaches 577 before L-BFGS gives up.
• Cross-validated log loss is U-shaped. A strong penalty () shrinks the weights so much that the probabilities are timid: underfitting. A weak penalty lets the weights grow to chase the training labels: overfitting. The minimum is near here, so the default happens to be a good choice on this dataset; that is luck, not a rule.
• The unpenalised model is a disaster on log loss and only mildly worse on accuracy. Its test accuracy, 0.939, looks respectable. Its test log loss, 1.969, is 25 times worse than the regularised model's, and worse than the baseline's 0.658. It is confidently wrong: a few test patients it has never seen get probabilities like for the diagnosis they actually have. Accuracy hides this; log loss doesn't.

LogisticRegressionCV automates this search efficiently, reusing each fit as the starting point for the next value of C. As always, the scaler lives inside the pipeline so it is refitted in every fold; scaling the full dataset before cross-validating leaks the held-out folds' statistics into training. The linear regression article's section on preprocessing without leaking covers the rules.

Evaluating a classifier

A regression model is wrong by some amount. A classifier can be wrong in two different directions, and its probabilities can be good or bad separately from its decisions. So classification needs more than one metric, and it is worth being clear about what each one measures.

Is it better than guessing?

Accuracy, the fraction of correct predictions, is the obvious first number, and on its own it is close to useless. The baseline that never predicts malignant gets 63.2% here, because 63.2% of patients are benign. On a problem with 1% positives, a model that always says "negative" scores 99%. Always report accuracy next to the majority-class rate, and prefer the metrics below.

Are the probabilities good?

Two metrics score the probabilities directly, before any threshold.

Log loss, the training loss, evaluated on test data: 0.077 for the model against 0.658 for the baseline. It is the strictest metric here, and the one that exposed the unpenalised model.

The Brier score, introduced by Glenn Brier in 1950 for weather forecasts, is the mean squared error between probability and label, : 0.021 against 0.233. It is the squared error we rejected as a training loss; as an evaluation metric it is useful because it is bounded and less dominated by a single confident mistake than log loss.

Both are proper scoring rules: they are minimised in expectation by reporting the true probability, so they reward honesty. Accuracy is not; it cannot tell a 51% prediction from a 99% one.

Are the probabilities calibrated?

A model is calibrated if, among the patients it gives a 20% probability, about 20% turn out to be malignant. Calibration is what lets you treat the output as a probability in a decision rule, and it is checked by binning predictions and comparing each bin's mean prediction with the observed rate. With only 114 test patients the bins would be tiny, so this check uses out-of-fold predictions on the training set: each patient is scored by a cross-validation model that did not see it.

import numpy as np
from sklearn.calibration import calibration_curve
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_predict, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

data = load_breast_cancer(as_frame=True)
X, y = data.data, 1 - data.target
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)
Out-of-fold probabilities: each training patient is scored by a model
that did not see it. 455 predictions instead of the test set's 114.
pipe = make_pipeline(StandardScaler(), LogisticRegression(C=1.0))
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
oof = cross_val_predict(pipe, X_train, y_train, cv=cv, method="predict_proba")[:, 1]
observed, predicted = calibration_curve(y_train, oof, n_bins=5, strategy="quantile")
for p, o in zip(predicted, observed):
    print(f"mean predicted {p:.3f}  ->  observed malignant rate {o:.3f}")
balanced = make_pipeline(StandardScaler(), LogisticRegression(C=1.0, class_weight="balanced"))
balanced.fit(X_train, y_train)
plain = pipe.fit(X_train, y_train)
print(f"intercept: plain = {plain[-1].intercept_[0]:.3f}, "
      f"balanced = {balanced[-1].intercept_[0]:.3f}, "
      f"log(n_benign / n_malignant) = {np.log((1 - y_train).sum() / y_train.sum()):.3f}")

mean predicted 0.000  ->  observed malignant rate 0.000
mean predicted 0.002  ->  observed malignant rate 0.011
mean predicted 0.047  ->  observed malignant rate 0.033
mean predicted 0.811  ->  observed malignant rate 0.824
mean predicted 1.000  ->  observed malignant rate 1.000
intercept: plain = -0.243, balanced = 0.167, log(n_benign / n_malignant) = 0.517

(The last line is used in the class imbalance section.) Each bin holds 91 patients. Predicted and observed rates agree closely: when the model says about 81%, about 82% are malignant. This is typical of logistic regression. Because it is trained with a proper scoring rule, it tends to produce well-calibrated probabilities when the model is roughly right and regularisation is not extreme. Niculescu-Mizil and Caruana (2005) compared ten learning methods and found logistic regression among the best calibrated out of the box, while boosted trees, SVMs and naive Bayes needed recalibration. It is not guaranteed: a misspecified model, heavy regularisation or resampled training data all distort calibration. scikit-learn's calibration guide covers checking and fixing it.

The five bins also show something about this dataset: the predictions are bunched near 0 and 1. Most tumours are easy. The interesting patients are the handful in between.

Choosing a threshold

To act, we need a threshold $t$: predict malignant when . Each choice of $t$ gives a confusion matrix, and the matrix gives the error-specific metrics.

• Recall (sensitivity, true positive rate) : of the malignant tumours, what fraction did we catch?
• Precision : of the tumours we flagged, what fraction were malignant?
• Specificity : of the benign tumours, what fraction did we clear? One minus specificity is the false positive rate.
• F1 is the harmonic mean of precision and recall: a single number that is high only if both are.

Lowering the threshold flags more patients: recall goes up, precision usually goes down. Move it yourself. Every dot below is a real test patient, at the probability the model above gave it.

At a few thresholds:

At the default 0.5, three cancers are missed. Their probabilities were 0.49, 0.31 and 0.06. Lowering the threshold to 0.3 catches two of them at no extra cost on this test set. Catching the third requires going down to 0.05, at the price of eleven false alarms out of 72 benign patients.

Which is right? That is not a statistical question. If the probabilities are calibrated and the costs of the two errors are (a missed cancer) and (an unnecessary biopsy), predicting positive is worth it when the expected cost of saying "benign", , exceeds the expected cost of saying "malignant", . Solving for $p$:

Equal costs give , which is where the default comes from. If a missed cancer is judged 20 times worse than a false alarm, . The model supplies probabilities; the threshold encodes values. Keep the two separate, and choose the threshold on validation data, not on the test set: picking it by looking at the table above, as we just did for illustration, would make the test results optimistic.

[!WARNING] Don't tune the threshold to maximise accuracy or F1 on the test set and then report those numbers. The threshold is a hyperparameter like any other. Choose it with cross-validation or a separate validation split, and touch the test set once. scikit-learn's TunedThresholdClassifierCV does this inside cross-validation.

Ranking quality: ROC and PR curves

The ROC curve plots the true positive rate against the false positive rate as the threshold sweeps from 1 to 0; it is the curve in the explorer. A model that guesses lies on the diagonal. The area under it, ROC-AUC, has a direct reading: the probability that a randomly chosen malignant case gets a higher score than a randomly chosen benign one (Hanley and McNeil, 1982). Here it is 0.996: the model's ranking is nearly perfect, and the misses at 0.5 are a threshold problem as much as a model problem. ROC-AUC ignores calibration entirely; any increasing transformation of the scores gives the same curve.

When positives are rare, ROC-AUC can look excellent while precision is poor, because the false positive rate divides by the large number of negatives. The precision-recall curve and its summary, average precision (PR-AUC, average_precision_score), are more informative then; Saito and Rehmsmeier (2015) make the case with examples. A guessing model's PR-AUC equals the positive rate, which is why the baseline scored 0.368.

Class imbalance

Imbalance changes logistic regression less than people expect. The model's intercept absorbs the base rate: with fewer positives, $b$ is lower and every probability is lower, which is correct, because positives are less likely.

The common reflex is to "fix" imbalance by weighting the minority class up (class_weight="balanced") or resampling. For logistic regression, the main effect of reweighting is to shift the intercept. Weighting positives by $a$ and negatives by $c$ multiplies the odds the model reports by roughly . With balanced weights that factor is , a log-odds shift of on this training set. The calibration script above shows the intercept moving from to $0.167$, a shift of $0.41$: close to the theory, which holds exactly only for a correctly specified, unpenalised model. The ranking of patients barely changes. The probabilities, though, are no longer calibrated: they describe a world where the classes are balanced, which is not the one the model will be used in.

So for logistic regression, the usual advice is:

• Keep the probabilities honest and move the threshold. Train without reweighting, check calibration, and pick $t$ from the costs as above. This gives the same decisions reweighting would, without corrupting the probabilities.
• Use metrics that respect imbalance: PR-AUC, recall at a fixed precision, or expected cost. Not accuracy.
• Reweight when the optimiser needs it, for example with extreme imbalance and few positives, where the regularised fit may barely move off the base rate. If you do, recalibrate afterwards, or shift the intercept back by .

Reading the coefficients

A large part of logistic regression's appeal is that its weights can be read. With standardised features, is the factor by which the odds of malignancy are multiplied for a one-standard-deviation increase in feature $j$, holding the other features fixed.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

data = load_breast_cancer(as_frame=True)
X, y = data.data, 1 - data.target
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42
)
model = make_pipeline(StandardScaler(), LogisticRegression(C=1.0)).fit(X_train, y_train)
coef = model[-1].coef_[0]
for j in np.argsort(-np.abs(coef))[:8]:
    print(f"{X.columns[j]:>20}: w = {coef[j]:6.3f}, odds ratio per SD = {np.exp(coef[j]):.2f}")
print(X_train[["mean radius", "mean perimeter", "mean area"]].corr().round(3))
Refit on 200 bootstrap resamples of the training set: how stable are the weights?
rng = np.random.default_rng(0)
boot = []
for _ in range(200):
    idx = rng.integers(0, len(X_train), len(X_train))
    refit = make_pipeline(StandardScaler(), LogisticRegression(C=1.0))
    boot.append(refit.fit(X_train.iloc[idx], y_train.iloc[idx])[-1].coef_[0])
lo, hi = np.percentile(np.array(boot), [2.5, 97.5], axis=0)
print(f"weights whose 95% bootstrap interval contains 0: {((lo < 0) & (hi > 0)).sum()} of 30")

       worst texture: w =  1.434, odds ratio per SD = 4.20
        radius error: w =  1.233, odds ratio per SD = 3.43
      worst symmetry: w =  1.061, odds ratio per SD = 2.89
 mean concave points: w =  0.953, odds ratio per SD = 2.59
     worst concavity: w =  0.911, odds ratio per SD = 2.49
          area error: w =  0.909, odds ratio per SD = 2.48
   compactness error: w = -0.907, odds ratio per SD = 0.40
          worst area: w =  0.900, odds ratio per SD = 2.46
                mean radius  mean perimeter  mean area
mean radius           1.000           0.998      0.987
mean perimeter        0.998           1.000      0.986
mean area             0.987           0.986      1.000
weights whose 95% bootstrap interval contains 0: 12 of 30

Read naively: each standard deviation of worst texture multiplies the odds of malignancy by 4.2, and compactness error is protective. Be careful with both readings.

• The features are heavily correlated. Radius, perimeter and area are three measurements of the same nucleus size, correlated at 0.99 or more, and the same holds for their error and worst versions. As in the linear regression article's multicollinearity section, the data determine the combined effect of correlated features well and how to split it between them poorly. Across 200 bootstrap refits, 12 of the 30 weights have intervals that include zero: the data don't even determine their sign. compactness error is not one of them (it is negative in every refit), but a stable negative weight on one of a cluster of correlated shape features is more plausibly a correction to its partners' positive weights than a protective effect. The model only promises that the combination predicts well.
• The weights here are shrunk. With C=1.0 they are pulled toward zero, and correlated features share the shrinkage. Odds ratios from a penalised model are biased toward 1.
• Odds ratios are not probability changes. A factor of 4.2 on the odds moves a 1% probability to 4%, and a 50% probability to 81%.
• Association, not causation. These are cell measurements made after the fact. Nothing about the model says changing a nucleus's texture changes anything.

For explaining individual predictions, the score decomposes exactly: , so each feature's contribution to a patient's log-odds is . That additivity in log-odds is what makes logistic regression easy to audit, and it is why it remains common in regulated settings such as credit scoring.

More than two classes: softmax

With classes, logistic regression generalises to multinomial logistic regression, also called softmax regression. Each class $k$ gets its own weight vector and bias, and so its own score . The softmax function turns the $K$ scores into a probability distribution:

Every output is positive and they sum to 1. Adding the same constant to every score leaves the probabilities unchanged, so only differences between scores matter. With , the probability of class 1 is : binary logistic regression, with the score being the difference of the two class scores.

The loss is again the negative log-likelihood, , the categorical cross-entropy. Writing the label as a one-hot vector and the predicted distribution as , the same cancellation as before gives

Prediction minus label again, now a vector. This is precisely the gradient that starts backpropagation in every neural-network classifier: the last layer of such a network is multinomial logistic regression, on features learned by the layers below.

scikit-learn fits the multinomial model by default for three or more classes (all solvers except liblinear). The fitted model is literally the formula above:

import numpy as np
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression

X, y = load_iris(return_X_y=True)          # 150 flowers, 4 features, 3 species
model = LogisticRegression(max_iter=1000).fit(X, y)
scores = X @ model.coef_.T + model.intercept_          # (150, 3) logits
scores -= scores.max(axis=1, keepdims=True)            # stability; softmax is shift-invariant
softmax = np.exp(scores) / np.exp(scores).sum(axis=1, keepdims=True)
print("coef_ shape:", model.coef_.shape)
print("max |softmax − predict_proba|:", np.abs(softmax - model.predict_proba(X)).max())
print("rows sum to 1:", np.allclose(softmax.sum(axis=1), 1.0))

coef_ shape: (3, 4)
max |softmax − predict_proba|: 0.0
rows sum to 1: True

The weight matrix is : one row of four weights per species. Subtracting each row's maximum score before exponentiating is the standard trick to avoid overflow; the shift invariance means it changes nothing else.

The alternative is one-vs-rest: train $K$ independent binary classifiers, "class $k$ or not", and pick the most confident. It is simple and parallel, but the $K$ probabilities don't sum to 1 and aren't jointly calibrated. scikit-learn provides it as OneVsRestClassifier; prefer the multinomial model unless you have a reason not to.

Common misconceptions

"Logistic regression is a regression method." The name is historical. It regresses the log-odds on the features, and it is used as a classifier. It is a generalised linear model: a linear predictor, a Bernoulli distribution for the response, and the logit as the link between them (Nelder and Wedderburn, 1972).

"Its output is a probability, so it is calibrated." Its output is in $[0, 1]$. Whether it means what it says depends on the model being roughly right, the regularisation, and whether the training data were reweighted or resampled. Check with a calibration curve.

"Classify at 0.5." 0.5 is optimal only when both errors cost the same and the probabilities are calibrated. The breast cancer example shows thresholds from 0.05 to 0.3 being defensible.

"Logistic regression can only draw straight boundaries." Its boundary is linear in the features you give it. With polynomial, interaction or spline features, it draws curves and circles.

"Separable data is the easy case." For unpenalised maximum likelihood it is the broken case: the weights go to infinity and the probabilities become meaningless. Regularise.

"Accuracy tells you how good the model is." Accuracy ignores the probabilities, hides the direction of errors, and is inflated by imbalance. The unpenalised model had 94% test accuracy and a test log loss worse than predicting the base rate.

"Imbalanced data needs resampling or class weights." For logistic regression these mostly shift the intercept and break calibration. Threshold choice usually does the same job better.

"You must scale features." Scaling doesn't change the unpenalised maximum-likelihood fit at all, only its conditioning. It does matter for the optimiser's convergence and, because the penalty depends on weight size, it changes the regularised solution. With scikit-learn's default penalty, you are in the second case.

When to use it

Logistic regression is a strong choice when:

• You need probabilities, not just labels, especially well-calibrated ones that feed a cost-based decision or another system.
• You need to explain or audit the model. Each prediction is a sum of per-feature contributions to the log-odds, and each weight has a stated meaning.
• Data are limited relative to the number of features. A regularised linear model is hard to beat with a few hundred examples.
• Training and serving must be cheap. A prediction is one dot product and one exponential; it scales to billions of sparse features with first-order solvers, which is why it was long the workhorse of click-through-rate prediction.
• You need a baseline. If a gradient-boosted ensemble can't clearly beat regularised logistic regression on log loss, the extra complexity may not be worth it.

Look further when the relationship involves strong non-linearities or interactions you can't engineer by hand (tree ensembles, neural networks), or when the data are images, audio or text, where the useful features have to be learned. Even then, logistic regression usually survives as the final layer.

Summary

• Logistic regression models the log-odds of the positive class as a linear function of the features. Inverting the logit gives .
• Each weight adds to the log-odds, so is an odds ratio. The decision boundary is a hyperplane; the threshold that turns probabilities into actions is a separate choice.
• Maximum likelihood under a Bernoulli model gives the log loss. It punishes confident mistakes without limit, which keeps the probabilities honest.
• The gradient of log loss with respect to the score is , so : the linear regression gradient with a sigmoid inside. The Hessian is positive semi-definite, so the loss is convex.
• There is no closed form. Newton's method (IRLS) converges in a handful of iterations; gradient descent and L-BFGS scale further.
• On separable data the maximum-likelihood estimate doesn't exist and the weights diverge. An L2 penalty, which is a Gaussian prior, restores a finite, unique solution. scikit-learn applies one by default.
• Evaluate the probabilities (log loss, Brier, calibration), the ranking (ROC-AUC, PR-AUC) and the decisions (confusion matrix, precision, recall) separately. Choose the threshold from the costs of errors, on validation data.
• Softmax regression extends all of this to $K$ classes with the same gradient, and is the output layer of neural-network classifiers.

Further reading

Textbooks

• Gareth James, Daniela Witten, Trevor Hastie, Robert Tibshirani and Jonathan Taylor, An Introduction to Statistical Learning, chapter 4: logistic regression alongside discriminant analysis and naive Bayes, with Python labs. Free online.
• Trevor Hastie, Robert Tibshirani and Jerome Friedman, The Elements of Statistical Learning, §4.4: the full treatment, including IRLS, L1-regularised logistic regression, and the comparison with linear discriminant analysis.
• Andrew Ng and Tengyu Ma, Stanford CS229 lecture notes, chapters 2 and 3: logistic regression, Newton's method, softmax, and generalised linear models derived from the exponential family.
• Ian Goodfellow, Yoshua Bengio and Aaron Courville, Deep Learning, chapters 5 and 6: maximum likelihood, cross-entropy, and why sigmoid and softmax outputs pair with log-likelihood losses.

Papers

• Joseph Berkson (1944), Application of the Logistic Function to Bio-Assay, Journal of the American Statistical Association. Coined "logit".
• D. R. Cox (1958), The Regression Analysis of Binary Sequences, Journal of the Royal Statistical Society, Series B. Logistic regression in its modern form.
• J. A. Nelder and R. W. M. Wedderburn (1972), Generalized Linear Models, Journal of the Royal Statistical Society, Series A. The framework in which logistic regression is one member of a family.
• A. Albert and J. A. Anderson (1984), On the Existence of Maximum Likelihood Estimates in Logistic Regression Models, Biometrika. Separation, quasi-separation and overlap.
• David Firth (1993), Bias Reduction of Maximum Likelihood Estimates, Biometrika. A penalised likelihood that stays finite under separation.
• Dong C. Liu and Jorge Nocedal (1989), On the Limited Memory BFGS Method for Large Scale Optimization, Mathematical Programming. The algorithm behind scikit-learn's default solver.
• Alexandru Niculescu-Mizil and Rich Caruana (2005), Predicting Good Probabilities with Supervised Learning, ICML. Which learning methods produce calibrated probabilities.
• James A. Hanley and Barbara J. McNeil (1982), The Meaning and Use of the Area under a Receiver Operating Characteristic (ROC) Curve, Radiology.
• Takaya Saito and Marc Rehmsmeier (2015), The Precision-Recall Plot Is More Informative than the ROC Plot When Evaluating Binary Classifiers on Imbalanced Datasets, PLOS ONE.
• W. Nick Street, William H. Wolberg and O. L. Mangasarian (1993), Nuclear Feature Extraction for Breast Tumor Diagnosis, Proc. SPIE 1905. The source of the breast cancer dataset.

Documentation

• scikit-learn: logistic regression user guide, LogisticRegression, probability calibration, tuning the decision threshold and model evaluation metrics.
• Google's Machine Learning Crash Course: Logistic regression and Classification, a gentler interactive introduction to the same ideas.