43 - Regressione logistica: dalla retta alla probabilità

Un medico non dice «il paziente ha il 100% o lo 0% di probabilità di essere malato»: ragiona per gradi di fiducia, e solo alla fine decide se intervenire. La regressione logistica fa lo stesso. Prende un punteggio, costruito come somma pesata delle feature, e lo trasforma in una probabilità tra 0 e 1; poi, se serve una decisione, si confronta la probabilità con una soglia. Nonostante il nome è un algoritmo di classificazione, ed è importante per due motivi: è il classificatore lineare più usato in pratica, e è il neurone con cui costruirai le reti neurali nelle prossime lezioni. Capirla a fondo, con la sua funzione di costo e il suo gradiente, significa avere già in mano il cuore dell'addestramento di una rete.

Perché non basta la regressione lineare

Con due classi codificate 0 e 1 ti verrebbe da usare la regressione lineare e dire «positivo se la previsione supera 0,5». Due problemi, che possiamo misurare. Il primo: la retta non è limitata, quindi produce «probabilità» sotto 0 e sopra 1. Il secondo: un punto molto lontano, ma facilmente classificabile, sposta la retta e con essa la soglia di decisione. Verifichiamolo su due gruppi di punti separati lungo un asse:

Nota sul codice: i blocchi Python di questo articolo vanno eseguiti in ordine, in un unico script o notebook. Ciascuno riusa gli import e le variabili definiti nei blocchi precedenti. Tutti sono stati eseguiti due volte (Python 3.14.4, NumPy 2.3.5, scikit-learn 1.8.0) con uscita identica; con scikit-learn 1.8.0 le righe con C=np.inf stampano su stderr un UserWarning («Setting penalty=None will ignore the C and l1_ratio parameters»), innocuo, che l'uscita riportata qui sotto non 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 punti di classe 0, 50 di classe 1
y = np.concatenate([np.zeros(50), np.ones(50)])

def soglia_lineare(x, y):
    lr = LinearRegression().fit(x[:, None], y)
    return lr, -(lr.intercept_ - 0.5) / lr.coef_[0]                      # x dove la retta vale 0.5

lr, s0 = soglia_lineare(x, y)
print("previsioni della retta in x=-2 e x=9: %s  (non sono probabilità)" % np.round(lr.predict([[-2], [9]]), 2))

x_extra, y_extra = np.append(x, 40.0), np.append(y, 1.0)                  # un solo punto di classe 1, molto lontano
_, s1 = soglia_lineare(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("soglia di decisione, regressione lineare: %.3f -> %.3f dopo aver aggiunto il punto lontano" % (s0, s1))
print("soglia di decisione, regressione logistica: %.4f -> %.4f" % (-log0.intercept_[0] / log0.coef_[0, 0], -log1.intercept_[0] / log1.coef_[0, 0]))

previsioni della retta in x=-2 e x=9: [-0.84  1.8 ]  (non sono probabilità)
soglia di decisione, regressione lineare: 3.581 -> 3.853 dopo aver aggiunto il punto lontano
soglia di decisione, regressione logistica: 3.5713 -> 3.5713

Cosa fa il codice. Genera due nuvole su una sola dimensione (classe 0 centrata in 2, classe 1 in 5). Per la regressione lineare calcola la soglia come il valore di x in cui la retta vale 0,5. Poi aggiunge un solo punto di classe 1 a x = 40, perfettamente coerente (è dal lato giusto). Per la regressione lineare la soglia si sposta da 3,58 a 3,85 solo per un punto che non avrebbe dovuto cambiare nulla, e le previsioni in x = −2 e x = 9 sono −0,84 e 1,80: non sono probabilità. Con la regressione logistica la soglia resta a 3,571 in entrambi i casi. Serve una funzione che «schiacci» qualsiasi punteggio nell'intervallo (0, 1).

La funzione sigmoide

La sigmoide (o funzione logistica) trasforma un numero reale z qualunque in un valore tra 0 e 1:

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

Cosa fa: per z molto negativo e^(−z) è enorme e σ(z) è vicino a 0; per z molto positivo e^(−z) è vicino a 0 e σ(z) è vicino a 1; per z = 0 vale esattamente 0,5. La forma è una «S» liscia. Due proprietà ci serviranno tra poco:

σ(−z) = 1 − σ(z)                (simmetria)
σ'(z) = σ(z) · (1 − σ(z))       (derivata)

La derivata si dimostra con la regola della catena. Scrivi σ(z) = (1 + e^(−z))⁻¹ e deriva «da fuori verso dentro»: σ'(z) = −(1 + e^(−z))⁻² · (−e^(−z)) = e^(−z) / (1 + e^(−z))². Ora spezza il rapporto in due fattori: [1 / (1 + e^(−z))] · [e^(−z) / (1 + e^(−z))] = σ(z) · [1 − σ(z)], dove l'ultima uguaglianza segue da e^(−z)/(1 + e^(−z)) = 1 − 1/(1 + e^(−z)). La derivata è massima in z = 0 (vale 0,25) e tende a zero nelle code: un fatto che avrà conseguenze importanti nelle reti profonde.

Il modello: punteggio, probabilità, frontiera

Il modello costruisce il punteggio come nella regressione lineare, z = θ·x (con x₀ = 1 per l'intercetta), e lo passa alla sigmoide:

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

Si decide «classe 1» se p ≥ 0,5. Poiché σ(z) ≥ 0,5 se e solo se z ≥ 0, la regola di decisione è semplicemente θ·x ≥ 0: la frontiera di decisione è l'iperpiano θ·x = 0, lineare. Il modello è quindi un classificatore lineare che in più fornisce una probabilità graduata: più un punto è lontano dalla frontiera, più p si avvicina a 0 o a 1.

A sinistra la sigmoide trasforma il punteggio in probabilità, con p = 0,5 in z = 0. A destra, nel piano delle feature, la frontiera θ·x = 0 è una retta; il vettore dei pesi (la freccia arancione) è perpendicolare ad essa e punta verso la classe 1.

Interpretazione: odds e log-odds

Le odds (quote) di un evento di probabilità p sono p/(1 − p): se p = 0,75, le odds sono 3 a 1. Invertendo la sigmoide si ottiene che θ·x è il logaritmo delle odds (log-odds o logit):

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

Quindi il modello è lineare nel logit. Questo dà un'interpretazione ai coefficienti: aumentare di 1 la feature x_j (a parità delle altre) aumenta il logaritmo delle odds di θ_j, cioè moltiplica le odds per e^(θ_j). Un coefficiente di 0,7 raddoppia le odds (e^0,7 ≈ 2). Attenzione: non significa «aumenta la probabilità di 0,7»: la relazione tra feature e probabilità è non lineare, e l'effetto sulla probabilità dipende da dove ti trovi sulla curva.

La funzione di costo: la cross-entropia, derivata dalla verosimiglianza

Come per la regressione lineare, non scegliamo il costo a caso: lo ricaviamo da un'ipotesi probabilistica. Ogni etichetta è un lancio di moneta con probabilità di testa p_i = σ(θ·x_i). La probabilità di osservare l'etichetta y_i ∈ {0, 1} si scrive in forma compatta così:

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

(per y_i = 1 resta p_i, per y_i = 0 resta 1 − p_i). La verosimiglianza di tutto il dataset è il prodotto su tutti gli esempi; il logaritmo la trasforma in una somma. Cambiando segno e facendo la media, per poter minimizzare, si ottiene la log-loss o cross-entropia binaria:

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

Ha un'interpretazione intuitiva: penalizza in modo severo gli errori sicuri. Se y = 1 e il modello dà p = 0,01, il costo di quell'esempio è −log(0,01) ≈ 4,6; se p = 0,4 è −log(0,4) ≈ 0,92; se p = 0,99 è ≈ 0,01. Un modello che dice «sono sicuro al 99%» e sbaglia paga moltissimo.

Perché non l'errore quadratico? Si può usare, ma combinato con la sigmoide produce una superficie di costo non convessa, con zone piatte dove il gradiente è minuscolo: la discesa fatica a partire. La cross-entropia, invece, è convessa nei parametri per la regressione logistica, ed è sostenuta da una derivazione probabilistica.

Il gradiente: lo stesso della regressione lineare

Calcoliamo ∂J/∂θ con la regola della catena. Per un singolo esempio, chiamiamo ℓ = −[y log p + (1 − y) log(1 − p)], con p = σ(z) e z = θ·x. Tre derivate, da fuori verso dentro:

∂ℓ/∂p = −y/p + (1 − y)/(1 − p)
∂p/∂z = p (1 − p)                      (derivata della sigmoide)
∂z/∂θ_j = x_j

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

Lo svolgimento: moltiplicando, il termine −y/p per p(1 − p) dà −y(1 − p); il termine (1 − y)/(1 − p) per p(1 − p) dà (1 − y)p; sviluppando −y + yp + p − yp = p − y. Un risultato meraviglioso: la derivata del costo rispetto al punteggio è semplicemente previsione meno etichetta. Moltiplicando per ∂z/∂θ_j = x_j e mediando sugli esempi:

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

È identica nella forma al gradiente della regressione lineare, con la previsione sostituita da p = σ(Xθ). Questa «coincidenza» non è un caso: è una proprietà dei modelli lineari generalizzati con la funzione di collegamento giusta, e la ritroverai nell'ultimo strato di una rete neurale.

Implementazione a mano

Applichiamo la discesa del gradiente, con feature standardizzate. Il dataset è bidimensionale così da poter ragionare sulla frontiera:

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)                                            # media e std SOLO dal training
A = np.column_stack([np.ones(len(Xtr)), sc.transform(Xtr)])               # (420, 3): uni + 2 feature
B = np.column_stack([np.ones(len(Xte)), sc.transform(Xte)])               # (180, 3)

def sigmoide(z):
    return 1.0 / (1.0 + np.exp(-np.clip(z, -500, 500)))                  # clip: evita overflow di exp per z molto negativo

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

def addestra_logistica(A, y, alpha=0.5, iterazioni=2000):
    theta = np.zeros(A.shape[1]); storia = []
    for _ in range(iterazioni):
        p = sigmoide(A @ theta)
        theta = theta - alpha * A.T @ (p - y) / len(y)                    # gradiente: (1/m) X^T (p - y)
        storia.append(logloss(sigmoide(A @ theta), y))
    return theta, storia

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

print("theta a mano:    ", np.round(theta, 4))
print("theta sklearn:   ", np.round(np.r_[sk.intercept_, sk.coef_[0]], 4))
print("log-loss con theta = 0: %.4f (ln 2 = %.4f);  dopo 100 iterazioni %.4f;  finale %.4f" % (logloss(sigmoide(A @ np.zeros(3)), ytr), np.log(2), storia[99], storia[-1]))
print("accuratezza sul test: a mano %.3f   sklearn %.3f" % (((sigmoide(B @ theta) >= 0.5) == yte).mean(), sk.score(B[:, 1:], yte)))

theta a mano:     [ 0.387  -3.7381  0.3253]
theta sklearn:    [ 0.3869 -3.7379  0.3253]
log-loss con theta = 0: 0.6931 (ln 2 = 0.6931);  dopo 100 iterazioni 0.2553;  finale 0.2511
accuratezza sul test: a mano 0.939   sklearn 0.939

Cosa fa il codice. sigmoide usa np.clip per evitare overflow in np.exp con punteggi molto negativi (con SciPy esiste scipy.special.expit, numericamente stabile, documentazione). logloss limita le probabilità a [10⁻¹², 1 − 10⁻¹²] per non calcolare mai log(0), che sarebbe −∞. Il ciclo di addestramento è quello della discesa del gradiente per la regressione lineare, con una sola modifica: la «previsione» è sigmoide(A @ theta). A θ = 0 tutte le probabilità valgono 0,5 e la log-loss iniziale è ln 2 = 0,6931 (il modello dice «non so» su tutto): è un controllo di sanità utile per ogni classificatore binario. I coefficienti coincidono con quelli di scikit-learn (C=np.inf disattiva la regolarizzazione, di cui parliamo nella prossima lezione; la documentazione di LogisticRegression dice che C=np.inf equivale al vecchio penalty=None, deprecato dalla versione 1.8), e anche l'accuratezza.

Sulle feature standardizzate θ ≈ (0,387, −3,738, 0,325): la prima feature ha un peso forte e negativo, la seconda debole. In unità standardizzate, aumentare di una deviazione standard la prima feature moltiplica le odds di e^(−3,74) ≈ 0,024, cioè le riduce di circa 40 volte.

Più di due classi: softmax

Con K classi il modello produce K punteggi z_k = θ_k·x (uno per classe) e li trasforma in probabilità con la funzione softmax:

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

Le probabilità sono positive e sommano a 1. Con K = 2 la softmax coincide con la sigmoide (la differenza dei due punteggi gioca il ruolo di z). La perdita è la cross-entropia categorica, J = −(1/m) Σ_i log p_{i, y_i}, e il suo gradiente ha di nuovo la forma «previsione meno etichetta»: ∇_W J = (1/m) Xᵀ(P − Y), dove Y è la matrice delle etichette in codifica one-hot (una riga per esempio, con 1 nella colonna della classe vera e 0 altrove). Su Iris (3 specie):

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)                                  # sottrarre il massimo non cambia il risultato ed evita 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))                                                      # un vettore di pesi per classe
for _ in range(3000):
    P = softmax(Aa @ W)                                                   # (105, 3)
    W = W - 0.3 * Aa.T @ (P - Y) / len(ya)                                # gradiente: (1/m) X^T (P - Y)

lib = LogisticRegression(C=np.inf, max_iter=2000).fit(Aa[:, 1:], ya)
print("forma di W:", W.shape, " righe di P sommano a:", np.round(softmax(Ab @ W).sum(axis=1)[:3], 3))
print("cross-entropia sul training: %.4f" % -np.log(softmax(Aa @ W)[np.arange(len(ya)), ya]).mean())
print("accuratezza sul test: a mano %.3f   sklearn %.3f" % ((softmax(Ab @ W).argmax(axis=1) == yb).mean(), lib.score(Ab[:, 1:], yb)))

forma di W: (5, 3)  righe di P sommano a: [1. 1. 1.]
cross-entropia sul training: 0.0646
accuratezza sul test: a mano 1.000   sklearn 1.000

Cosa fa il codice. np.eye(3)[ya] costruisce la matrice one-hot scegliendo la riga dell'identità corrispondente all'etichetta. softmax sottrae il massimo di ogni riga prima di esponenziare: matematicamente non cambia il risultato (numeratore e denominatore si moltiplicano per lo stesso fattore) ma evita che e^(z) vada in overflow. W ha forma (5, 3): 5 coefficienti (intercetta + 4 feature) per ciascuna delle 3 classi, quindi Aa @ W ha forma (105, 3), un punteggio per esempio e per classe. Il gradiente è la stessa riga di prima, con P e Y matrici al posto di p e y vettori. Le forme sono la chiave: (5, 105) @ (105, 3) = (5, 3), la stessa forma di W.

La trappola della separazione perfetta

Se le classi sono perfettamente separabili da una frontiera lineare, la log-loss può sempre diminuire rendendo i pesi più grandi: più sono grandi, più le probabilità si avvicinano a 0 e 1, più il costo scende verso zero. Non esiste un minimo finito: l'addestramento non converge mai. Verifichiamolo:

xs = np.concatenate([rng.normal(-3, 0.5, 50), rng.normal(3, 0.5, 50)])    # due gruppi ben separati
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 = addestra_logistica(As, ys, alpha=0.5, iterazioni=it)
    print("dopo %5d iterazioni: theta = %s   log-loss = %.5f" % (it, np.round(t, 3), st[-1]))

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

Il peso della feature continua a crescere, in modo lento ma senza limite (circa 1,3 → 2,1 → 3,0 → 3,9 al moltiplicarsi per dieci delle iterazioni), mentre la log-loss tende a zero. Nella pratica, con classi quasi separabili, i pesi diventano enormi, le probabilità estreme e mal calibrate. È una delle ragioni principali della regolarizzazione, che penalizza i pesi grandi: argomento della prossima lezione.

Trade-off e trappole

• Pro: semplice, veloce, probabilità graduate, coefficienti interpretabili (log-odds), buona baseline, convesso (un solo minimo).
• Contro: frontiera lineare (non cattura interazioni a meno di costruire feature); sensibile a feature correlate; richiede feature scalate per la discesa del gradiente; le probabilità sono ragionevoli solo se il modello è ben specificato.
• Soglia 0,5 non sacra: con classi sbilanciate o costi asimmetrici scegli la soglia dalla curva precision-recall (vedi il processo di un progetto di ML: split, leakage e metriche). Usare class_weight="balanced" è un'altra leva.
• Accuratezza su classi sbilanciate: mente, come spiega la lezione su split, leakage e metriche. Guarda precision, recall e log-loss.
• MSE come costo per la classificazione: funziona male con la sigmoide (gradiente quasi nullo quando il modello è sicuro e sbaglia).
• Interpretare un coefficiente come variazione di probabilità: è una variazione del log-odds.
• Overflow e log(0): usa funzioni stabili (expit, clip, softmax con sottrazione del massimo).
• Shape mismatch nella softmax: Z.max(axis=1) senza keepdims=True produce un vettore (m,) che non si sottrae a una matrice (m, K) nel modo atteso (broadcast sbagliato o errore).

Le otto voci si raggruppano in tre famiglie. Le prime due (pro e contro) dicono quando scegliere il modello: come baseline lineare con probabilità interpretabili. La soglia, l'accuratezza su classi sbilanciate e il costo sbagliato riguardano la valutazione: il modello può essere giusto e il modo di giudicarlo no. Le ultime tre (log-odds letti come probabilità, overflow, forma della softmax) sono errori di implementazione che non producono eccezioni ma numeri sbagliati, e per questo gli esperimenti sopra confrontano sempre il risultato a mano con quello di scikit-learn.

Leggere i coefficienti: odds ratio con un intervallo di confidenza

Nella regressione logistica il coefficiente w_j è la variazione del log-odds per un aumento unitario della feature j. Esponenziare lo rende più leggibile: exp(w_j) è l'odds ratio, cioè il fattore per cui si moltiplicano gli odds (p / (1 − p)) quando la feature aumenta di un'unità, a parità delle altre. Un odds ratio di 2 significa «odds doppi», uno di 0,5 «odds dimezzati», uno di 1 «nessun effetto». Per sapere quanto fidarsi di un odds ratio serve un intervallo; qui lo otteniamo con il bootstrap, ricampionando i punti e riadattando il modello. Simuliamo 2 000 pazienti con due variabili standardizzate (età e dose), la cui verità è log-odds = −1,0 + 0,8·età − 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)                     # nessuna regolarizzazione
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("frequenza dell'esito: %.3f" % outcome.mean())
for name, j, true_value in (("età", 0, 0.8), ("dose", 1, -1.2)):
    print("%-5s coefficiente %+.3f (%s %+.1f)   odds ratio %.3f   intervallo 95%%: [%.3f, %.3f]   odds ratio vero %.3f" % (name, fit.coef_[0, j], "vero", true_value, np.exp(fit.coef_[0, j]), low[j], high[j], np.exp(true_value)))

frequenza dell'esito: 0.325
età   coefficiente +0.792 (vero +0.8)   odds ratio 2.207   intervallo 95%: [1.977, 2.505]   odds ratio vero 2.226
dose  coefficiente -1.215 (vero -1.2)   odds ratio 0.297   intervallo 95%: [0.255, 0.337]   odds ratio vero 0.301

Cosa fa il codice. Genera 2 000 pazienti con esito binario (frequenza 0,325), adatta una regressione logistica senza regolarizzazione e poi la ripete su 300 ricampionamenti bootstrap. L'età ha coefficiente +0,792 (vero +0,8), cioè un odds ratio di exp(0,792) = 2,207: per ogni deviazione standard in più di età gli odds dell'esito più che raddoppiano. La dose ha coefficiente −1,215 (vero −1,2), odds ratio 0,297: ogni deviazione standard in più divide gli odds per circa 3,4. Gli intervalli bootstrap al 95% sono [1,977; 2,505] per l'età e [0,255; 0,337] per la dose, e contengono entrambi i valori veri (2,226 e 0,301).

Due insidie di interpretazione. Primo: l'odds ratio non è un rischio relativo. Quando l'esito è raro (probabilità piccole) odds e probabilità quasi coincidono e le due letture sono vicine; qui l'esito riguarda un terzo dei casi e non lo sono, quindi «odds doppi» non significa «probabilità doppia». Secondo: i coefficienti confrontabili tra loro sono quelli di feature sulla stessa scala (le abbiamo standardizzate); se età fosse in anni e dose in milligrammi, i due numeri non si potrebbero mettere a confronto.

Infine, un odds ratio descrive un'associazione condizionata alle altre variabili del modello, non un effetto causale: dirlo richiede un disegno dello studio, non un coefficiente.

Come la libreria addestra davvero: il metodo di Newton (IRLS)

Nella lezione sulla discesa del gradiente abbiamo visto che il metodo di Newton corregge il gradiente con l'inversa dell'Hessiana. Per la regressione logistica la funzione di costo è convessa ma non quadratica, quindi un passo non basta; tuttavia pochi passi sì. Il gradiente della log-verosimiglianza è Xᵀ(y − p) e l'Hessiana (cambiata di segno) è Xᵀ W X, con W la matrice diagonale con elementi p_i (1 − p_i). L'aggiornamento w ← w + (XᵀWX)⁻¹ Xᵀ(y − p) è detto IRLS (minimi quadrati pesati iterati), perché ogni passo è una regressione ai minimi quadrati pesata dove i pesi dipendono dalle probabilità correnti. Confrontiamolo con la discesa del gradiente semplice.

X_b = np.column_stack([np.ones(n_l), X_l])
w_newton = np.zeros(3)
print("iterazione   log-verosimiglianza   passo massimo")
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("pesi di Newton: %s     pesi di scikit-learn: %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("discesa del gradiente (learning rate 0,5): %d iterazioni per arrivare entro 1e-6" % iteration)

iterazione   log-verosimiglianza   passo massimo
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
pesi di Newton: [-1.0035  0.7916 -1.2151]     pesi di scikit-learn: [-1.0038  0.7916 -1.2149]
discesa del gradiente (learning rate 0,5): 352 iterazioni per arrivare entro 1e-6

Cosa fa il codice. Implementa a mano il metodo di Newton per la regressione logistica. All'iterazione 0, con tutti i pesi a zero, ogni probabilità vale 0,5 e la log-verosimiglianza è −1 386,29 (che è esattamente 2 000 · ln 0,5, il valore di un modello che tira a caso). Poi sale a −1 009,34, −976,69, −975,13 e si ferma a −975,1277 all'iterazione 4: da lì il passo massimo scende a 10⁻¹⁰ e 10⁻¹⁶, cioè all'arrotondamento.

Questa è la convergenza quadratica di Newton: il passo massimo passa da 0,77 a 0,34, 0,097, 0,0061, 2,3·10⁻⁵, 3·10⁻¹⁰, e a ogni iterazione il numero di cifre corrette circa raddoppia. I pesi ottenuti (−1,0035; 0,7916; −1,2151) coincidono con quelli di scikit-learn (−1,0038; 0,7916; −1,2149) fino a differenze di circa 3·10⁻⁴ sul termine noto e 2·10⁻⁴ sul peso della dose, dovute alla tolleranza di arresto di scikit-learn.

La discesa del gradiente con learning rate 0,5, che usa solo il gradiente, richiede 352 iterazioni per arrivare entro 10⁻⁶: quasi sessanta volte più di Newton (6 iterazioni per scendere sotto 10⁻⁶).

Ogni passo di Newton, tuttavia, costa più di uno di discesa del gradiente perché forma e risolve un sistema p × p. Con tre pesi è irrilevante; con migliaia di feature scikit-learn passa a metodi intermedi (come L-BFGS, che approssima l'Hessiana dai gradienti passati) e con molti dati a metodi stocastici.

Le probabilità sono calibrate? E la regolarizzazione le rende oneste

Un classificatore con alta accuratezza può dare probabilità inaffidabili. La calibrazione (definizione nella guida scikit-learn) chiede: tra i casi a cui il modello assegna probabilità 0,8, il 80% è davvero positivo? Si controlla raggruppando le previsioni in fasce e confrontando la probabilità media prevista con la frequenza osservata. Un modello correttamente specificato su dati abbondanti è tipicamente ben calibrato; uno addestrato su pochi dati con molte feature tende invece a essere troppo sicuro di sé. Verifichiamo entrambi i casi.

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("fascia di probabilità   casi   prevista (media)   osservata (frequenza)")
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 esempi, 100 feature, una sola informativa e molto rumore
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 (inverso della regolarizzazione)   accuratezza train   accuratezza test   confidenza media   log-loss test")
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))))

fascia di probabilità   casi   prevista (media)   osservata (frequenza)
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 (inverso della regolarizzazione)   accuratezza train   accuratezza test   confidenza media   log-loss test
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

Cosa fa il codice. La prima tabella divide le probabilità previste sul test (1 000 casi) in cinque fasce e confronta la media prevista con la frequenza osservata. Il modello è ben calibrato nel complesso: nella fascia 0,0–0,2 prevede 0,103 e osserva 0,094; nella 0,4–0,6 prevede 0,501 e osserva 0,546; nella 0,8–1,0 prevede 0,864 e osserva 0,942. Gli scarti arrivano a 7,8 punti nell'ultima fascia, ma lì i casi sono solo 52, e l'errore standard di una frequenza osservata su 52 casi è circa √(0,94·0,06/52) ≈ 3,2 punti: lo scarto è circa 2,4 errori standard, un caso limite. La lezione è che una tabella di calibrazione con poche migliaia di casi e cinque fasce può essere letta solo grossolanamente.

La seconda tabella è l'avvertimento serio: 300 esempi, 100 feature e una sola davvero informativa (con molto rumore), 150 per l'addestramento e 150 per il test. Senza regolarizzazione (C infinito) il modello ha accuratezza 1,000 sul training e 0,520 sul test (quasi un lancio di moneta), ma una confidenza media di 0,966: è sicurissimo di risposte che sono quasi casuali, e la log-loss sul test è 9,775 (un modello costante che dice sempre 0,5 avrebbe 0,693).

Con C = 1 le cose migliorano ma restano brutte (log-loss 2,076, confidenza 0,809). Con C = 0,01 la confidenza scende a 0,166 e la log-loss a 0,699: le probabilità sono oneste, ma siamo praticamente al livello del modello costante, e l'accuratezza di 0,560 su 150 casi non è distinguibile dal caso.

Il segnale, in questo esperimento, è troppo debole per le feature disponibili. È un risultato negativo da riportare così: la regolarizzazione non crea informazione, evita soltanto di fingere di averla.

Quando la frontiera lineare non basta: cerchi concentrici

X_c, y_c = make_circles(600, noise=0.1, factor=0.4, random_state=0)         # un cerchio dentro un altro
linear = LogisticRegression()
quadratic = make_pipeline(PolynomialFeatures(2), LogisticRegression(C=10, max_iter=5000))
print("accuratezza CV (5 fold)   frontiera lineare: %.3f   con feature di secondo grado: %.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("feature:      %s" % str(list(fit_q[0].get_feature_names_out())))
print("coefficienti: %s" % str(np.round(fit_q[-1].coef_[0], 2)))

accuratezza CV (5 fold)   frontiera lineare: 0.483   con feature di secondo grado: 0.993
feature:      ['1', 'x0', 'x1', 'x0^2', 'x0 x1', 'x1^2']
coefficienti: [-6.000e-02  1.000e-02  1.400e-01 -1.250e+01 -6.900e-01 -1.249e+01]

Cosa fa il codice. Genera 600 punti su due cerchi concentrici (uno di raggio piccolo dentro uno più grande) e confronta due modelli con la validazione incrociata. La regressione logistica sulle sole due coordinate raggiunge un'accuratezza di 0,483: leggermente peggiore del caso, perché la frontiera è una retta e non esiste nessuna retta che separi un cerchio da quello che lo contiene (è lo stesso limite del XOR per il perceptrone). Con le feature di secondo grado (1, x₀, x₁, x₀², x₀x₁, x₁²) l'accuratezza sale a 0,993, con lo stesso algoritmo di addestramento: cambia solo lo spazio in cui il modello è lineare.

I coefficienti raccontano il perché: i termini x₀² e x₁² valgono −12,50 e −12,49, praticamente uguali, mentre i termini lineari e il misto x₀x₁ sono piccoli (0,01; 0,14; −0,69). Il modello ha imparato una funzione di x₀² + x₁², cioè del quadrato della distanza dall'origine, che è esattamente la quantità che separa i due cerchi.

È il ponte verso le reti neurali: invece di sceglierti a mano le feature da aggiungere, una rete le impara dai dati con gli strati nascosti. Allo stesso tempo mostra il limite del metodo manuale: funziona perché abbiamo indovinato che il grado 2 bastava, cosa che in problemi ad alta dimensione non si può fare a occhio.

Tre esercizi

Esercizio 1. Nell'esperimento delle fasce di calibrazione, riduci il campione a 200 pazienti: le fasce restano vicine alla diagonale? Come cambia la variabilità tra fasce e perché una fascia con 10 casi non dice nulla? Esercizio 2. Nell'esperimento con 100 feature, usa la validazione incrociata per scegliere C su una griglia da 10⁻³ a 10²: il C scelto è quello con la log-loss più bassa o con l'accuratezza più alta, e le due scelte coincidono? Esercizio 3. Sui cerchi concentrici prova il grado 3 e il grado 6 con la stessa C: l'accuratezza migliora? Guarda i coefficienti per capire quali termini la regressione ritiene davvero utili.

Lezione precedente del percorso: Regressione multivariata, polinomiale e normalizzazione delle feature. Nella stessa serie rimandano a questa lezione anche: Entropia, entropia incrociata e divergenza KL: perché le funzioni di costo sono tutte la stessa cosa, Naive Bayes: il teorema di Bayes con un'ipotesi ingenua, Boosting: AdaBoost, Gradient Boosting e cenni su XGBoost, Riduzione di dimensionalità: PCA da zero e t-SNE.

Fonti

• scikit-learn, LogisticRegression - parametri C e penalty (C=np.inf come regressione non penalizzata, penalty deprecato dalla 1.8), consultata il 2026-11-12
• scikit-learn, Linear Models (sezione Logistic regression) - la regressione logistica come modello lineare per la classificazione, binaria, one-vs-rest e multinomiale, consultata il 2026-11-12
• scikit-learn, Probability calibration - la definizione di classificatore ben calibrato usata nella sezione sulla calibrazione, consultata il 2026-11-12
• SciPy, scipy.special.expit - la sigmoide logistica numericamente stabile, consultata il 2026-11-12

Le cifre della lezione sono misure dell'autore con seme fisso. Riferimento storico non riaperto durante la verifica: D. R. Cox, The regression analysis of binary sequences, Journal of the Royal Statistical Society B, 1958, per il modello logistico.

Nella prossima lezione vedi perché un modello con troppa libertà memorizza il rumore, e come la regolarizzazione lo mette a freno, in overfitting e regolarizzazione: L2, L1 e la trappola dell'accuratezza.