Logistic regression

The linear classifiers note ended with a table of what each surrogate loss estimates. The hinge recovers only the side of the boundary; the logistic loss recovers the log-odds $\log\frac{\eta}{1-\eta}$, from which the probability $\eta$ itself follows. Many decisions need that probability, not just the side: an advertiser compares the chance of a click with its price, a clinic weighs a patient’s risk against the cost of a test, and the next note chooses thresholds from costs. Logistic regression is the linear model built for this. This note derives it from a probability model of the label, fits it by maximum likelihood, works out its gradient and its failure on separable data, and asks when its probabilities can be taken at face value.

Notation. In this note the labels are $y_i \in \{0, 1\}$, which suits the Bernoulli distribution; the margin form of the previous note uses $\pm 1$ labels, $2y_i - 1$. The features $x_i \in \mathbb R^d$ include a constant 1, so the bias is one of the weights $w \in \mathbb R^d$. The score (or logit) is $z_i = \langle w, x_i\rangle$, the sigmoid is $\sigma(z) = 1/(1 + e^{-z})$, and the model’s probability is $p_i = \sigma(z_i)$. $X$ is the $n \times d$ matrix with rows $x_i$, $y$ and $p$ are the vectors of labels and probabilities, and $\eta(x) = P(y = 1 \mid x)$ is the true probability.

Why the log-odds?#

A linear function $\langle w, x\rangle$ takes every real value; a probability lives in $[0, 1]$. Modelling $p$ itself as linear produces “probabilities” of $1.3$ and $-0.2$. The odds $p/(1 - p)$ live in $(0, \infty)$ and the log-odds in $(-\infty, \infty)$, so the log-odds is a quantity that a linear model can describe without leaving its range:

$$\log\frac{P(y = 1 \mid x)}{1 - P(y = 1 \mid x)} = \langle w, x\rangle \quad\Longleftrightarrow\quad P(y = 1 \mid x) = \frac{1}{1 + e^{-\langle w, x\rangle}} = \sigma(\langle w, x\rangle) .$$

The sigmoid is just the inverse of the log-odds. Two of its properties are used below: $\sigma(-z) = 1 - \sigma(z)$, and $\sigma'(z) = \sigma(z)\,(1 - \sigma(z))$.

Why linear log-odds are natural. By Bayes’ formula in odds form, the posterior log-odds are the prior log-odds plus the log likelihood ratio $\log\frac{p(x \mid y = 1)}{p(x \mid y = 0)}$. When the two classes are Gaussian with a common covariance $\Sigma$ and means $\mu_1, \mu_0$, the quadratic terms of the two log-densities cancel and

$$\log\frac{p(x \mid y = 1)}{p(x \mid y = 0)} = (\mu_1 - \mu_0)^\top \Sigma^{-1} x + \text{const},$$

which is linear in $x$ (Bishop, 2006, Section 4.2.1; the Gaussian itself is in the multivariate Gaussian note). The same happens for many other class-conditional families. Logistic regression keeps the conclusion, linear log-odds, and drops the assumption: it models $P(y \mid x)$ directly and never describes how $x$ is distributed.

Worked example. Classes $N(1, 1)$ and $N(-1, 1)$ on a line, equally likely. The log likelihood ratio is $\big[(x + 1)^2 - (x - 1)^2\big]/2 = 2x$, so the log-odds are exactly $2x$ and at $x = 0.5$ the probability of class 1 is $\sigma(1) = 0.731$. If class 1 is rarer, with prior $0.2$, the prior log-odds $\log(0.2/0.8) = -1.386$ are added: $\sigma(1 - 1.386) = 0.405$. A change of base rate shifts every log-odds by the same constant, the rule the probability course derived for deployment under a different prevalence.

Reading the weights. Increasing feature $j$ by one unit adds $w_j$ to the log-odds, which multiplies the odds by $e^{w_j}$, whatever the other features are: a weight of $0.7$ doubles the odds ($e^{0.7} = 2.01$). This is the odds ratio that medical and social-science papers report. It describes the model, and it is causal only if the data allow a causal reading.

Maximum likelihood gives the log-loss#

Each label is a Bernoulli variable with success probability $p_i = \sigma(\langle w, x_i\rangle)$, and the labels are independent given the features. The likelihood of the training labels is $\prod_i p_i^{y_i} (1 - p_i)^{1 - y_i}$, and its negative logarithm is

$$\mathcal L(w) = -\sum_{i=1}^n \Big[ y_i \log p_i + (1 - y_i) \log(1 - p_i) \Big] ,$$

the log-loss, or binary cross-entropy. That cross-entropy losses are negative log-likelihoods of a categorical model is derived in the maximum likelihood note; the information-theoretic reading, the surprise of the observed labels under the model’s beliefs, is in the entropy note of the generative models course.

It is the logistic loss. For $y_i = 1$ the term is $-\log\sigma(z_i) = \log(1 + e^{-z_i})$. For $y_i = 0$ it is $-\log(1 - \sigma(z_i)) = -\log\sigma(-z_i) = \log(1 + e^{z_i})$. With the margin $M_i = (2y_i - 1)\, z_i$ both cases read

$$\mathcal L(w) = \sum_{i=1}^n \log\big(1 + e^{-M_i}\big) ,$$

the logistic surrogate of the previous note, with the natural logarithm. So the surrogate chosen there for its shape is exactly maximum likelihood under the model of linear log-odds.

Worked example. One feature, $w = 1$, no bias, three points: $x = (2, -1, 0.5)$ with labels $y = (1, 0, 0)$. The probabilities are $p = \sigma(x) = (0.881, 0.269, 0.622)$. The losses are $-\log 0.881 = 0.127$ for the confident correct point, $-\log(1 - 0.269) = 0.313$ for the moderately correct one, and $-\log(1 - 0.622) = 0.974$ for the third point, which the model places on the wrong side; the total is $1.414$. A model that said $p = 0.99$ for the third point would pay $-\log 0.01 = 4.6$ for it: the log-loss punishes confident mistakes without limit.

The gradient: residuals times inputs#

Differentiate one term with respect to its score. For $y_i = 1$, $\frac{d}{dz}\big[-\log\sigma(z)\big] = -(1 - \sigma(z))$; for $y_i = 0$, $\frac{d}{dz}\big[-\log(1 - \sigma(z))\big] = \sigma(z)$; both cases equal $\sigma(z_i) - y_i$. By the chain rule, with $\partial z_i / \partial w = x_i$,

$$\nabla \mathcal L(w) = \sum_{i=1}^n \big(\sigma(z_i) - y_i\big)\, x_i = X^\top (p - y) .$$

The form is the one of least squares, $X^\top(Xw - y)$: prediction minus target, times the input. A point pulls the weights in proportion to how wrong its probability is, and a correctly and confidently classified point, with $p_i \approx y_i$, hardly pulls at all. In the worked example $p - y = (-0.119, 0.269, 0.622)$, the gradient with respect to $w$ is $-0.119 \cdot 2 + 0.269 \cdot (-1) + 0.622 \cdot 0.5 = -0.196$, and with respect to a bias it would be $\sum_i (p_i - y_i) = 0.772$: the predictions are too high on average (their mean is $0.591$, the share of positives $1/3$), and a gradient step lowers the bias.

Convexity. Differentiating once more, with $\sigma' = \sigma(1 - \sigma)$, gives the Hessian

$$\nabla^2 \mathcal L(w) = \sum_{i=1}^n p_i (1 - p_i)\, x_i x_i^\top = X^\top S X, \qquad S = \operatorname{diag}\big(p_i (1 - p_i)\big) .$$

For any vector $v$, $v^\top X^\top S X v = \sum_i p_i(1 - p_i) \langle x_i, v\rangle^2 \ge 0$, so the log-loss is convex, and strictly convex when $X$ has full column rank (every $p_i$ is strictly between 0 and 1 at finite $w$). Every local minimum is global.

No closed form, but Newton’s method. Setting the gradient to zero gives $X^\top(\sigma(Xw) - y) = 0$, equations that are nonlinear in $w$, unlike the normal equations of least squares. They are solved iteratively. Gradient descent works; with a few hundred features, Newton’s method converges in a handful of steps:

$$w \leftarrow w - (X^\top S X)^{-1} X^\top (p - y) .$$

Each step is a weighted least-squares problem with weights $S$ and a working response $Xw + S^{-1}(y - p)$, which is why it is called iteratively reweighted least squares (Bishop, 2006, Section 4.3.3). On the breast cancer data with two standardized features, 25 Newton steps reach a gradient of $10^{-15}$ and the same weights as scikit-learn (Code section).

What the optimum guarantees. The bias is the weight of the constant feature, so the gradient equation for it reads $\sum_i (p_i - y_i) = 0$: at the maximum likelihood fit, the average predicted probability equals the share of positives in the training set. The same holds within any group marked by a binary feature of the model. On the breast cancer fit, the predicted probabilities of malignancy add up to $212.0$, the number of malignant tumours.

Figure

Logistic regression: the boundary and the contours of probability

The figure shows the geometry. The probability depends on $x$ only through the score $z = \langle w, x\rangle$, so its contours are the hyperplanes $z = \operatorname{logit}(p)$, all parallel to the boundary. Moving along $w$ by a distance $1/\lVert w\rVert$ changes the log-odds by one; the larger $\lVert w\rVert$, the steeper the transition from “surely 0” to “surely 1”. The width of the band between $p = 0.1$ and $p = 0.9$, $2\log 9/\lVert w\rVert$, is a direct picture of how confident the model is.

Separable data: the weights run away#

If some $w_0$ classifies the training set perfectly, every margin $M_i = (2y_i - 1)\langle w_0, x_i\rangle$ is positive, and scaling $w_0$ by $c > 1$ multiplies them all by $c$. Each term $\log(1 + e^{-cM_i})$ decreases towards zero as $c$ grows, while $\mathcal L(w) > 0$ for every finite $w$. The infimum, zero, is not attained: the maximum likelihood estimate does not exist. It exists exactly when the classes overlap, neither completely nor quasi-completely separated (Albert and Anderson, 1984).

Gradient descent on separable data keeps going. Its direction settles down, towards the maximum-margin separator of the support vector machines note (Soudry et al., 2018), while its norm grows without bound, only logarithmically in the number of steps, and all predicted probabilities drift to 0 and 1. On four points $x = -2, -1, 1, 2$ with the two right ones positive, gradient descent with step $0.5$ from zero has slope $5.55$ after $10^3$ steps, $7.83$ after $10^4$ and $10.13$ after $10^5$: about $\log(t/4)$, so every tenfold increase in steps adds $\log 10 = 2.30$ to the weight and divides the loss, about $2/t$, by ten.

Separation is not exotic. As the previous note showed, random labels on $n$ points in $d$ dimensions are separable with probability $1/2$ when $n = 2(d + 1)$, and real data are far from random. The breast cancer data with all 30 standardized features are linearly separable: as the penalty is relaxed, the fitted weight norm goes from $3.8$ to $1{,}891$ and the training log-loss from $0.053$ to $0.004$ (Code section).

The penalty fixes it. Minimizing

$$\mathcal L(w) + \frac{\lambda}{2}\lVert w\rVert^2$$

adds $\lambda I$ to the Hessian, so the objective is strictly convex and grows without bound in every direction: a unique, finite minimizer exists for any data. The penalty is a Gaussian prior $w \sim N(0, \lambda^{-1} I)$, and the penalized fit is its MAP estimate; the bias is usually left unpenalized. Larger $\lambda$ keeps $\lVert w\rVert$ small and the probability contours wide apart. scikit-learn’s LogisticRegression applies this penalty by default, with C = 1 the inverse of its strength; C = np.inf gives plain maximum likelihood.

Calibration: do the probabilities mean what they say?#

A model is calibrated if, among the examples to which it assigns probability $p$, a share $p$ is positive:

$$P\big(y = 1 \mid \hat p(x) = p\big) = p \qquad \text{for every } p .$$

A weather forecaster who says “70% chance of rain” is calibrated if it rains on 70% of such days. Calibration is checked with a reliability diagram: split the test predictions into bins (ten equal-width bins are common), and plot the mean prediction in each bin against the share of positives in it. A calibrated model’s points lie on the diagonal. The expected calibration error summarizes the distance,

$$\mathrm{ECE} = \sum_{b} \frac{n_b}{n} \,\big\lvert \bar y_b - \bar p_b \big\rvert ,$$

with $n_b$ examples, mean prediction $\bar p_b$ and share of positives $\bar y_b$ in bin $b$. The number depends on the bins and on the sample size. Each bin’s share of positives is noisy, so even perfectly calibrated predictions, with labels drawn as $y \sim \mathrm{Bernoulli}(p)$, have a positive ECE: on 10,000 examples it averages $0.010$ with ten bins and $0.022$ with fifty. A reliability diagram should therefore show the count of each bin, and the ECE should be read next to the log-loss and the Brier score below, which need no bins. Calibration alone is not quality: predicting the base rate for everyone is perfectly calibrated and useless. A good probabilistic classifier is calibrated and separates the classes, which the next note measures.

Why the log-loss aims at calibrated probabilities. At a single point with true probability $\eta$, a prediction $p$ has expected log-loss $-\eta\log p - (1 - \eta)\log(1 - p)$. Its derivative $-\eta/p + (1 - \eta)/(1 - p)$ vanishes at $p = \eta$: the expected loss is minimized by telling the truth. Losses with this property are called proper scoring rules (Gneiting and Raftery, 2007). The squared error $(y - p)^2$, the Brier score (Brier, 1950), is proper too; the absolute error $\lvert y - p\rvert$ is not, since its expected value $\eta(1 - p) + (1 - \eta)p$ is linear in $p$ and minimized at $0$ or $1$. Propriety is a statement about the best possible predictor. A fitted logistic regression is calibrated in the large on its training set, by the bias equation above, and calibrated everywhere only to the extent that the true log-odds are linear and the sample is large. Overfitting makes it overconfident; a strong penalty makes it underconfident.

Other models are not even aiming. A hinge-trained score is not a probability at all. Naive Bayes multiplies the evidence of correlated features as if it were independent, counts the same signal several times, and tends to push its probabilities towards 0 and 1. Max-margin methods and boosted trees tend to push them away from 0 and 1, giving S-shaped reliability curves (Niculescu-Mizil and Caruana, 2005). Modern deep networks trained to zero training error tend to be overconfident (Guo et al., 2017). These are observed tendencies of model families, not guarantees. In the Code section, naive Bayes on Gaussian classes with correlated features has ECE $0.236$ and log-loss $1.133$, while logistic regression, the correct model there, has $0.012$ and $0.633$. On these 10,000 test examples, ECE values near $0.01$ are at the level of perfectly calibrated predictions, so they cannot rank well-calibrated models against each other.

Recalibration. A score that ranks well but is miscalibrated can be repaired by a monotone map fitted on data not used for training (a held-out set or cross-validation, as in the validation note):

• Platt scaling (Platt, 1999) fits a one-feature logistic regression on the score: $\hat p = \sigma(a s + b)$, where $s$ is the model’s score or the logit of its probability, by maximum likelihood. To limit overfitting on small sets, Platt replaced the targets $1$ and $0$ by $(N_+ + 1)/(N_+ + 2)$ and $1/(N_- + 2)$, with $N_\pm$ the class counts. Temperature scaling of neural networks is the special case $b = 0$, a single factor on the logits (Guo et al., 2017).
• Isotonic regression (Zadrozny and Elkan, 2002) fits the best non-decreasing step function of the score in squared error, by the pool-adjacent-violators algorithm (Ayer et al., 1955): walk through the examples sorted by score and, whenever a block’s mean exceeds the next block’s, merge them. For six labels in score order, $0, 1, 0, 0, 1, 1$, the second label exceeds the third, the merged block $(1, 0)$ with mean $0.5$ exceeds the fourth label, and the result is $0, \tfrac13, \tfrac13, \tfrac13, 1, 1$. It can fix any monotone distortion, but it needs more data: with fewer than a few hundred to a thousand calibration examples Platt scaling usually does better (Niculescu-Mizil and Caruana, 2005).

Platt’s map is strictly increasing (for $a > 0$) and leaves the ranking, and therefore the ROC curve of the next note, unchanged. The isotonic map is only non-decreasing: it can merge neighbouring scores into ties, which changes the AUC little. Both statements are about one fixed score. CalibratedClassifierCV with cv=5, as in the Code section, fits the model and its map on each of five folds and averages the five calibrated outputs; an average of different models can rank the examples differently from any one of them. In the Code section the AUC of naive Bayes is $0.695$ before and after either recalibration.

Figure

Reliability diagram: an overconfident model before and after calibration

Several classes: the softmax#

With $K$ classes, each class gets its own weight vector $w_k$ and score $z_k = \langle w_k, x\rangle$, and the scores are turned into probabilities by the softmax:

$$P(y = k \mid x) = \frac{e^{z_k}}{\sum_{j=1}^K e^{z_j}} .$$

Maximum likelihood gives the cross-entropy loss $-\sum_i \log P(y_i \mid x_i)$, and its gradient has the same residual form, $X^\top(P - Y)$, with $P$ the $n \times K$ matrix of predicted probabilities and $Y$ the one-hot labels. Only differences of scores matter: adding the same vector to every $w_k$ changes nothing, so one class can be fixed at $w_K = 0$. For $K = 2$ the softmax probability of class 1 is $e^{z_1}/(e^{z_1} + e^{z_2}) = \sigma(z_1 - z_2)$, which is binary logistic regression with $w = w_1 - w_2$. The multiclass note develops this model and its alternatives.

Where logistic regression reappears#

• The first text classifier of the NLP course is multinomial logistic regression on word counts, the bag-of-words model; it reaches 0.918 accuracy on AG News and remains the baseline in the text classification note.
• Negative sampling in word2vec trains a logistic regression that tells observed (word, context) pairs from random ones, and a language model’s output layer is a softmax over the vocabulary trained with cross-entropy, as in the n-gram note.
• The discriminator of a GAN is a classifier trained with exactly this log-loss (GAN note), and noise-contrastive estimation learns an unnormalized density by logistic regression of data against noise (energy-based models note).

What changed#

	Linear classifier with a surrogate loss	Logistic regression
Output	a side of the hyperplane	a probability $\sigma(\langle w, x\rangle)$
Loss	chosen for its shape	derived: the negative Bernoulli log-likelihood
What the score means	depends on the loss	the log-odds; $e^{w_j}$ is an odds ratio
Fitting	gradient methods	the same, plus Newton’s method (IRLS) with $\nabla = X^\top(p - y)$
Separable data	many solutions	no maximum likelihood estimate; the penalty (a Gaussian prior) restores one
Probabilities	not available	calibrated if the model is right; checked by a reliability diagram, repaired by Platt or isotonic maps

Limitations#

• Linear log-odds. Interactions and nonlinear effects must be supplied as features; otherwise the model is miscalibrated where the true log-odds bend.
• Separation and near-separation. Unpenalized fits diverge or give huge, unstable coefficients; the penalty strength then matters and must be chosen by validation, with features on comparable scales.
• Coefficients are fragile to interpret. Correlated features share and trade weight, and an odds ratio is conditional on all the other features in the model.
• Calibration is not discrimination. A calibrated model can still separate the classes poorly, and calibration measured on one population breaks when the base rate or the feature distribution shifts.
• Rare classes. The probabilities can be right while the default threshold of 0.5 almost never fires; choosing the threshold is a separate decision.

What you should remember#

• Logistic regression models the log-odds as linear: $P(y = 1 \mid x) = \sigma(\langle w, x\rangle)$. Gaussian classes with a common covariance give exactly this form.
• Maximum likelihood under the Bernoulli model gives the log-loss, which is the logistic surrogate $\sum_i \log(1 + e^{-M_i})$.
• The gradient is $X^\top(p - y)$, the Hessian $X^\top S X$ with $S = \operatorname{diag}(p_i(1 - p_i))$; the loss is convex, has no closed-form minimizer, and Newton’s method is IRLS.
• At the optimum with a bias, the mean predicted probability equals the share of positives.
• On separable data the weights grow without bound and the probabilities go to 0 and 1; an L2 penalty, a Gaussian prior, gives a unique finite fit.
• The log-loss is a proper scoring rule. Check calibration with a reliability diagram that shows its bin counts, and with ECE next to the log-loss or Brier score; repair it with Platt scaling or isotonic regression fitted on held-out data. Platt scaling keeps the ranking; isotonic regression can only add ties.
• The softmax with cross-entropy is the $K$-class version; for $K = 2$ it reduces to the sigmoid.

Code#

The Gaussian log-odds, Newton’s method, separation, calibration of naive Bayes, pool adjacent violators and the two-class softmax:

import numpy as np
from scipy.special import expit as sigmoid, softmax
from scipy.stats import norm
from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from sklearn.calibration import CalibratedClassifierCV
from sklearn.isotonic import IsotonicRegression
from sklearn.metrics import log_loss, roc_auc_score

rng = np.random.default_rng(0)

# 1. Gaussian classes N(+1, 1) and N(-1, 1): the log-odds are exactly 2x; a prior of 0.2 adds log(0.25)
x0 = 0.5
log_odds = norm.logpdf(x0, 1) - norm.logpdf(x0, -1)
print(log_odds, round(sigmoid(log_odds), 4), round(sigmoid(log_odds + np.log(0.2 / 0.8)), 4))   # 1.0 0.7311 0.4046

# 2. maximum likelihood by Newton's method (IRLS) on two standardized breast cancer features
data = load_breast_cancer()
X = StandardScaler().fit_transform(data.data[:, :2])        # mean radius, mean texture
y = 1.0 - data.target                                       # 1 = malignant
Xb = np.c_[np.ones(len(X)), X]
w = np.zeros(3)
for _ in range(25):
    p = sigmoid(Xb @ w)
    grad = Xb.T @ (p - y)                                   # the gradient X^T (sigma - y)
    H = Xb.T @ (Xb * (p * (1 - p))[:, None])                # the Hessian X^T S X
    w -= np.linalg.solve(H, grad)
sk = LogisticRegression(C=np.inf, solver='newton-cholesky').fit(X, y)
print(w.round(4), np.r_[sk.intercept_, sk.coef_[0]].round(4))   # [-0.7076  3.722   0.9374] twice
print(f'{np.abs(grad).max():.1e}', round(p.sum(), 4), y.sum())  # 1.2e-15 212.0 212.0: mean prediction = share of positives

# 3. all 30 features: the training set is separable, and the unpenalized loss runs to zero
X30 = StandardScaler().fit_transform(data.data)
for C in (1.0, 1e2, 1e4, 1e6):
    m = LogisticRegression(C=C, max_iter=100_000).fit(X30, y)
    p = m.predict_proba(X30)[:, 1]
    print(C, round(np.linalg.norm(m.coef_), 1), f'{log_loss(y, p):.2g}', int((m.predict(X30) != y).sum()))
# C = 1: |w| 3.8, loss 0.053, 7 errors | 1e2: 21.4, 0.03, 4 | 1e4: 113.8, 0.022, 4 | 1e6: 1891.1, 0.0042, 0

# 4. gradient descent on four separable points: the weight grows like log(t/4), the loss falls like 2/t
Xs, ys = np.array([[1., -2], [1, -1], [1, 1], [1, 2]]), np.array([0, 0, 1, 1.])
w = np.zeros(2)
for t in range(1, 100_001):
    w -= 0.5 * Xs.T @ (sigmoid(Xs @ w) - ys) / 4
    if t in (1_000, 10_000, 100_000):
        print(t, round(w[1], 3), f'{np.mean(np.logaddexp(0, -(2 * ys - 1) * (Xs @ w))):.2e}')
# 1000 5.548 1.95e-03 | 10000 7.828 1.99e-04 | 100000 10.127 2.00e-05

# 5. calibration: Gaussian classes with correlated features; naive Bayes double-counts the correlated evidence
d, rho = 8, 0.7
L = np.linalg.cholesky(rho * np.ones((d, d)) + (1 - rho) * np.eye(d))
yc = rng.integers(0, 2, 20_000)
Xc = rng.normal(size=(20_000, d)) @ L.T + np.where(yc[:, None] == 1, 0.3, -0.3)
tr, te = slice(0, 10_000), slice(10_000, None)
def ece(p, y, bins=10):
    k = np.minimum((p * bins).astype(int), bins - 1)
    return sum((k == b).mean() * abs(p[k == b].mean() - y[k == b].mean()) for b in range(bins) if (k == b).any())
models = {'naive Bayes': GaussianNB(), '+ Platt': CalibratedClassifierCV(GaussianNB(), method='sigmoid', cv=5),
          '+ isotonic': CalibratedClassifierCV(GaussianNB(), method='isotonic', cv=5), 'logistic regression': LogisticRegression()}
for name, m in models.items():
    p = m.fit(Xc[tr], yc[tr]).predict_proba(Xc[te])[:, 1]
    print(name, round(ece(p, yc[te]), 3), round(log_loss(yc[te], p), 3), round(roc_auc_score(yc[te], p), 3))
# ECE, log-loss, AUC: naive Bayes 0.236 1.133 0.695 | + Platt 0.008 0.64 0.695 | + isotonic 0.01 0.634 0.695
#                     logistic regression 0.012 0.633 0.695

# 6. isotonic regression (pool adjacent violators) on six labels sorted by score
print(IsotonicRegression().fit_transform(np.arange(6), [0, 1, 0, 0, 1, 1]).round(3))   # [0. 0.333 0.333 0.333 1. 1.]

# 7. two-class softmax = sigmoid of the score difference
z = np.array([1.3, -0.4])
print(softmax(z)[0].round(4), sigmoid(z[0] - z[1]).round(4))    # 0.8455 0.8455

# 8. ECE of perfectly calibrated predictions (y ~ Bernoulli(p)) on 10,000 examples, 200 repetitions, 10 and 50 bins
p0 = rng.random((200, 10_000)); y0 = (rng.random((200, 10_000)) < p0).astype(float)
print([round(float(np.mean([ece(p, y, b) for p, y in zip(p0, y0)])), 3) for b in (10, 50)])   # [0.01, 0.022]

Next#

Logistic regression returns a probability, and a probability is not yet a decision. Somewhere a threshold turns it into “flag” or “pass”, and the right threshold depends on how rare the positives are and on what each kind of error costs. Judging the model needs a number too, and accuracy, the obvious one, says almost nothing when positives are rare: a model that never flags anything scores 99% when 1% are positive. The next note builds the vocabulary for these questions: the confusion matrix, precision and recall, the ROC and precision–recall curves, and thresholds chosen from costs.

References#

• D. R. Cox. The regression analysis of binary sequences. Journal of the Royal Statistical Society, Series B, 20(2), 215–232, 1958.
• J. Berkson. Application of the logistic function to bio-assay. Journal of the American Statistical Association, 39(227), 357–365, 1944.
• A. Albert, J. A. Anderson. On the existence of maximum likelihood estimates in logistic regression models. Biometrika, 71(1), 1–10, 1984.
• D. Soudry, E. Hoffer, M. S. Nacson, S. Gunasekar, N. Srebro. The implicit bias of gradient descent on separable data. Journal of Machine Learning Research, 19(70), 1–57, 2018.
• J. C. Platt. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. In A. Smola, P. Bartlett, B. Schölkopf, D. Schuurmans (eds.), Advances in Large Margin Classifiers, 61–74. MIT Press, 1999.
• B. Zadrozny, C. Elkan. Transforming classifier scores into accurate multiclass probability estimates. KDD, 694–699, 2002.
• M. Ayer, H. D. Brunk, G. M. Ewing, W. T. Reid, E. Silverman. An empirical distribution function for sampling with incomplete information. The Annals of Mathematical Statistics, 26(4), 641–647, 1955.
• A. Niculescu-Mizil, R. Caruana. Predicting good probabilities with supervised learning. ICML, 625–632, 2005.
• C. Guo, G. Pleiss, Y. Sun, K. Q. Weinberger. On calibration of modern neural networks. ICML, 2017.
• T. Gneiting, A. E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477), 359–378, 2007.
• G. W. Brier. Verification of forecasts expressed in terms of probability. Monthly Weather Review, 78(1), 1–3, 1950.
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Section 4.2 (generative models with linear log-odds), Sections 4.3.2–4.3.4 (logistic regression, IRLS, the multiclass case).
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Section 4.4 (logistic regression).
• E. Sokolov et al. Машинное обучение (HSE): lecture 5 (logistic regression and the log-loss), lecture 6 (multinomial logistic regression), seminar 6: calibration.