03 - Regressione Logistica e classificazione binaria

Un ospedale vuole stimare, da trenta misure di un nodulo, se è benigno o maligno. La risposta che serve non è «numero 0,3» come uscirebbe da una retta, ma una probabilità compresa fra 0 e 1, a cui poi qualcuno applicherà una soglia decisa in base al costo degli errori. La regressione logistica fa esattamente questo: prende la combinazione lineare delle misure del capitolo 02 e la schiaccia in un intervallo [0, 1] con una funzione a S, la sigmoide. Nonostante il nome è un classificatore. Qui la ricaviamo, la addestriamo in due modi che raggiungono lo stesso minimo, la confrontiamo con scikit-learn, e vediamo come leggerne le probabilità quando le classi sono sbilanciate o più di due.

Gli esempi usano scikit-learn 1.9.0, numpy 2.4.6 e scipy 1.17.1, con i dataset Breast Cancer Wisconsin (569 noduli, 30 misure) e Digits (1797 immagini 8x8 di cifre) inclusi nella libreria, più dati sintetici con seme fisso. Le uscite sono identiche in due esecuzioni.

Dalla retta alla probabilità

Sia z = w·x + b il punteggio lineare di un esempio. La sigmoide σ(z) = 1/(1 + e^−z) trasforma qualsiasi numero reale in un valore fra 0 e 1: σ(0) = 0,5, valori molto positivi danno probabilità vicine a 1, molto negativi vicine a 0. Il modello dice P(y = 1 | x) = σ(z). Il significato dei pesi si legge meglio in termini di odds, il rapporto p/(1 − p) fra la probabilità di un evento e quella del suo contrario: un 75% di probabilità equivale a odds di 3 a 1. Invertendo la sigmoide si ottiene z = ln(p/(1 − p)), il logaritmo degli odds: il modello è lineare nei log-odds, e aumentare di una unità una misura moltiplica gli odds per e^w. Il manuale di James, Witten, Hastie e Tibshirani (An Introduction to Statistical Learning, paragrafo 4.3) sviluppa la stessa lettura.

La frontiera di decisione, cioè il luogo in cui il modello è indeciso (p = 0,5), è w·x + b = 0: un iperpiano. La regressione logistica separa le classi con superfici piatte; se la frontiera vera è curva servono misure trasformate o un altro modello, come quelli dei capitoli 06 e 07.

La perdita: l'entropia incrociata

Se ogni esempio è una moneta con probabilità p_i = σ(z_i), la verosimiglianza dell'insieme è il prodotto delle probabilità assegnate alle etichette osservate, e conviene minimizzare l'opposto del suo logaritmo: L = −(1/n) Σ [y_i ln p_i + (1 − y_i) ln(1 − p_i)], detta entropia incrociata o log loss. Il gradiente rispetto ai pesi ha una forma semplicissima, (1/n) Xᵀ(p − y), la stessa della regressione lineare con p al posto della previsione. L'Hessiana è (1/n) XᵀWX con W = diag(p_i(1 − p_i)), semidefinita positiva: la perdita è convessa, quindi ha un solo minimo e ogni metodo che scende lo trova.

Esercizio svolto: una previsione, una perdita

Un esempio ha punteggio z = 0,8. Gli odds sono e^0,8 = 2,2255 e la probabilità σ(0,8) = 0,69. Se l'etichetta vera è 1 la perdita è −ln 0,69 = 0,3711; se è 0 la perdita è −ln 0,31 = 1,1711. La stessa previsione costa più del triplo quando sbaglia. Il gradiente rispetto a z è p − y: −0,31 se y = 1 (spinge z in su) e +0,69 se y = 0 (lo spinge in giù, con più forza, perché l'errore è più grande).

import numpy as np
from scipy.special import expit

z = 0.8                                   # punteggio lineare w.x + b di un esempio
p = expit(z)
print("odds = e^z:", round(float(np.exp(z)), 4), " p = sigmoide(z):", round(float(p), 4))
print("perdita se l'etichetta vera e' 1:", round(float(-np.log(p)), 4))
print("perdita se l'etichetta vera e' 0:", round(float(-np.log(1 - p)), 4))
print("gradiente rispetto a z (p - y), y=1 e y=0:", round(float(p - 1), 4), round(float(p - 0), 4))

# Output:
# odds = e^z: 2.2255  p = sigmoide(z): 0.69
# perdita se l'etichetta vera e' 1: 0.3711
# perdita se l'etichetta vera e' 0: 1.1711
# gradiente rispetto a z (p - y), y=1 e y=0: -0.31 0.69

I valori stampati confermano i calcoli a mano. La perdita cresce senza limite quando il modello è sicuro e sbaglia (−ln p tende a infinito per p → 0): per questo l'entropia incrociata punisce la sicurezza sbagliata molto più dell'errore moderato, e per questo nei programmi si calcola con np.logaddexp(0, z) − y·z invece che con log(sigmoide), che per z molto negativi produrrebbe log(0).

Newton-Raphson e discesa del gradiente

Il listato implementa la perdita con penalità L2, il suo gradiente e la sua Hessiana, controlla il gradiente con le differenze finite, e minimizza in due modi. La discesa del gradiente usa un passo fisso di 0,5. Il metodo di Newton-Raphson (nella statistica detto IRLS, iteratively reweighted least squares) a ogni passo risolve H·δ = g e sposta i pesi di δ: usa la curvatura per scegliere la lunghezza del passo, e vicino al minimo converge in pochissime iterazioni. La penalità è scritta nella convenzione di scikit-learn: nella documentazione, la funzione minimizzata è la media della perdita più r(W)/(S·C), con r = ½||W||² e S uguale al numero di esempi quando non ci sono pesi. Con C = 1 ciò equivale alla mia lam = 1/(C·n), applicata ai pesi e non all'intercetta.

import numpy as np
from scipy.special import expit
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

X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=42, stratify=y)
sc = StandardScaler().fit(Xtr)
A = np.c_[np.ones(len(Xtr)), sc.transform(Xtr)]
n, d = A.shape
C = 1.0
lam = 1 / (C * n)                    # stessa penalita' di sklearn: media della perdita + (lam/2)||w||^2, intercetta esclusa

def perdita(w):
    z = A @ w
    ll = np.mean(np.logaddexp(0, z) - ytr * z)           # cross-entropy, stabile per z grandi
    return ll + 0.5 * lam * np.sum(w[1:] ** 2)

def grad(w):
    g = A.T @ (expit(A @ w) - ytr) / n
    g[1:] += lam * w[1:]
    return g

def hess(w):
    p = expit(A @ w)
    H = (A * (p * (1 - p))[:, None]).T @ A / n
    H[1:, 1:] += lam * np.eye(d - 1)
    return H

# 1) controllo del gradiente con le differenze finite
rng = np.random.default_rng(0)
w = 0.1 * rng.standard_normal(d)
num = np.array([(perdita(w + 1e-6 * e) - perdita(w - 1e-6 * e)) / 2e-6 for e in np.eye(d)])
print("gradiente analitico = differenze finite:", np.allclose(num, grad(w), atol=1e-8))

# 2) Newton-Raphson (IRLS)
w = np.zeros(d)
for it in range(1, 30):
    passo = np.linalg.solve(hess(w), grad(w))
    w = w - passo
    if np.linalg.norm(passo) < 1e-10:
        break
print("Newton: iterazioni =", it, " perdita =", round(perdita(w), 6))

# 3) discesa del gradiente con passo fisso
wg = np.zeros(d)
for k in range(1, 20001):
    wg = wg - 0.5 * grad(wg)
    if np.linalg.norm(grad(wg)) < 1e-6:
        break
print("discesa del gradiente: iterazioni =", k, " perdita =", round(perdita(wg), 6))

# 4) confronto con scikit-learn (stessa C, stessa penalita' L2)
sk = LogisticRegression(C=C, max_iter=5000, tol=1e-10).fit(sc.transform(Xtr), ytr)
w_sk = np.r_[sk.intercept_, sk.coef_.ravel()]
print("Newton = sklearn (tolleranza 1e-5):", np.allclose(w, w_sk, atol=1e-5))
print("gradient descent = Newton (tolleranza 1e-3):", np.allclose(wg, w, atol=1e-3))
At = np.c_[np.ones(len(Xte)), sc.transform(Xte)]
print("accuratezza test:", round(float(np.mean((expit(At @ w) > 0.5) == yte)), 4))

# 5) senza penalita' (C=inf) il dataset intero e' linearmente separabile: i pesi non hanno un massimo
Xs = StandardScaler().fit_transform(X)
for c in (1.0, np.inf):
    m = LogisticRegression(C=c, max_iter=5000, tol=1e-10).fit(Xs, y)
    print(f"C={c}: peso massimo in valore assoluto = {np.abs(m.coef_).max():.1f}, accuratezza train = {m.score(Xs, y):.3f}")

# Output:
# gradiente analitico = differenze finite: True
# Newton: iterazioni = 10  perdita = 0.068116
# discesa del gradiente: iterazioni = 5963  perdita = 0.068116
# Newton = sklearn (tolleranza 1e-5): True
# gradient descent = Newton (tolleranza 1e-3): True
# accuratezza test: 0.986
# C=1.0: peso massimo in valore assoluto = 1.3, accuratezza train = 0.988
# C=inf: peso massimo in valore assoluto = 64176.2, accuratezza train = 1.000

Il gradiente analitico coincide con quello numerico. Newton raggiunge il minimo (perdita 0,068116) in 10 iterazioni; la discesa del gradiente ne impiega 5963 per arrivare alla stessa perdita, perché non conosce la curvatura. Entrambi coincidono con LogisticRegression(C=1.0) entro le tolleranze indicate. Newton costa di più a ogni passo, perché risolve un sistema d × d, e non conviene quando le misure sono migliaia; scikit-learn offre per questo lbfgs (il default), che stima la curvatura dai gradienti passati.

L'accuratezza sul test è 0,986, cioè 141 casi giusti su 143. Togliere la penalità (C=np.inf; l'argomento penalty è deprecato nella versione installata e al suo posto si usano C e l1_ratio) cambia il comportamento: nelle ultime righe del listato, sull'intero dataset, i dati risultano separabili da un iperpiano, il peso massimo supera 60000 (valore che dipende dalla tolleranza del solutore: non ha un limite) e l'accuratezza sul training è 1,000, mentre con C=1.0 il peso massimo resta sotto 1,5. Quando le classi sono separabili la perdita non ha un minimo finito, perché basta far crescere i pesi per rendere le probabilità sempre più nette; la penalità serve anche a questo.

Dalla probabilità alla decisione

Il modello restituisce probabilità; la decisione richiede una soglia. I termini sono quelli della matrice di confusione: veri positivi (TP), falsi positivi (FP), falsi negativi (FN), veri negativi (TN). La precisione è TP/(TP+FP), la quota di allarmi giusti; la recall è TP/(TP+FN), la quota di positivi trovati. La ROC-AUC e la PR-AUC riassumono la qualità del ranking su tutte le soglie, e il log loss e il punteggio di Brier (errore quadratico medio delle probabilità) misurano le probabilità in sé. Sul dataset ora la classe positiva è il tumore maligno (53 casi su 143 nel test).

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

X, y = load_breast_cancer(return_X_y=True)
y = 1 - y                                   # nel dataset 0 = maligno: lo rendiamo la classe "positiva"
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=42, stratify=y)
modello = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)).fit(Xtr, ytr)
p = modello.predict_proba(Xte)[:, 1]
print("ROC-AUC:", round(roc_auc_score(yte, p), 4), " PR-AUC:", round(average_precision_score(yte, p), 4))
print("log loss:", round(log_loss(yte, p), 4), " Brier:", round(brier_score_loss(yte, p), 4))
for soglia in (0.5, 0.1):
    pred = (p >= soglia).astype(int)
    tn, fp, fn, tp = confusion_matrix(yte, pred).ravel()
    print(f"soglia {soglia}: TP={tp} FP={fp} FN={fn} TN={tn}  precisione={precision_score(yte, pred):.3f}  recall={recall_score(yte, pred):.3f}")

# Output:
# ROC-AUC: 0.9962  PR-AUC: 0.9943
# log loss: 0.0776  Brier: 0.0224
# soglia 0.5: TP=49 FP=1 FN=4 TN=89  precisione=0.980  recall=0.925
# soglia 0.1: TP=52 FP=6 FN=1 TN=84  precisione=0.897  recall=0.981

Con soglia 0,5 il modello trova 49 maligni su 53 (recall 0,925) con un solo falso allarme; abbassando la soglia a 0,1 i maligni mancati scendono a 1 (recall 0,981) e i falsi allarmi salgono a 6, la precisione scende a 0,897. Nessuna delle due soglie è «giusta»: dipende dal costo relativo di mancare un tumore e di mandare una persona a un esame in più. ROC-AUC 0,9962 e PR-AUC 0,9943 dicono che il ranking è eccellente, non quale soglia usare. Con 53 positivi nel test, un singolo caso sposta la recall di quasi due punti: queste cifre vanno lette come ordini di grandezza e non come classifica fra modelli. La famiglia completa delle metriche è nel capitolo 11.

Classi sbilanciate

Quando il 97% degli esempi è negativo, un modello che dice sempre «negativo» ha il 97% di accuratezza. Costruiamo un caso sintetico con il 3% di positivi e vediamo tre scelte: nessun intervento, abbassare la soglia, e il parametro class_weight="balanced", che moltiplica la perdita di ogni esempio per un peso inversamente proporzionale alla frequenza della sua classe.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import average_precision_score, precision_score, recall_score
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=5000, n_features=10, n_informative=4, weights=[0.97, 0.03],
                           class_sep=1.5, flip_y=0.0, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
print("quota positivi nel test:", round(yte.mean(), 3), " accuratezza del 'sempre 0':", round(1 - yte.mean(), 3))

def riga(nome, pred, p):
    print(f"{nome:34s} precisione={precision_score(yte, pred, zero_division=0):.3f} "
          f"recall={recall_score(yte, pred):.3f} PR-AUC={average_precision_score(yte, p):.3f}")

senza = LogisticRegression(max_iter=1000).fit(Xtr, ytr)
p0 = senza.predict_proba(Xte)[:, 1]
riga("senza pesi, soglia 0.5", (p0 >= 0.5).astype(int), p0)
riga("senza pesi, soglia 0.15", (p0 >= 0.15).astype(int), p0)
bil = LogisticRegression(class_weight="balanced", max_iter=1000).fit(Xtr, ytr)
p1 = bil.predict_proba(Xte)[:, 1]
riga("class_weight='balanced', soglia 0.5", (p1 >= 0.5).astype(int), p1)
print("probabilita' media prevista: senza pesi", round(p0.mean(), 3), " balanced", round(p1.mean(), 3))

# Output:
# quota positivi nel test: 0.03  accuratezza del 'sempre 0': 0.97
# senza pesi, soglia 0.5             precisione=0.583 recall=0.311 PR-AUC=0.354
# senza pesi, soglia 0.15            precisione=0.383 recall=0.511 PR-AUC=0.354
# class_weight='balanced', soglia 0.5 precisione=0.110 recall=0.733 PR-AUC=0.170
# probabilita' media prevista: senza pesi 0.035  balanced 0.295

Senza pesi e con soglia 0,5 il modello è prudente: precisione 0,583 e recall 0,311. Abbassare la soglia a 0,15 sposta l'equilibrio (recall 0,511, precisione 0,383) senza cambiare il modello, e lascia la PR-AUC a 0,354. Con class_weight="balanced" la recall sale a 0,733 ma la precisione crolla a 0,110, e la PR-AUC è 0,170, cioè il ranking è peggiorato; inoltre la probabilità media prevista passa da 0,035 (vicina alla frequenza vera, 0,03) a 0,295, quindi le probabilità non sono più probabilità.

Il risultato di questo esperimento, che vale per questi dati sintetici e non per ogni problema, è che pesare le classi cambia il modello e non solo la soglia; se ti serve soltanto più recall, spostare la soglia è più economico e lascia le probabilità interpretabili. La documentazione di scikit-learn sulla calibrazione osserva che la regressione logistica restituisce di norma probabilità ben calibrate perché ottimizza la perdita logaritmica, e il peso delle classi rompe proprio quella proprietà.

Più di due classi: softmax

Con K classi si hanno K vettori di pesi e la probabilità della classe k è la softmax dei punteggi: exp(z_k) / Σ_l exp(z_l), positiva e con somma 1 per costruzione. L'alternativa è l'uno-contro-tutti (one-vs-rest): un classificatore binario per classe, che non garantisce probabilità con somma 1. Nella versione 1.9.0 l'argomento multi_class di LogisticRegression non esiste più (passandolo si ottiene un TypeError): il caso multiclasse è multinomiale per default, e l'uno-contro-tutti si ottiene avvolgendo il modello in OneVsRestClassifier.

import numpy as np
from scipy.special import softmax
from sklearn.datasets import load_digits
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.multiclass import OneVsRestClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_digits(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=0, stratify=y)
soft = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(Xtr, ytr)
ovr = make_pipeline(StandardScaler(), OneVsRestClassifier(LogisticRegression(max_iter=2000))).fit(Xtr, ytr)
print("cifre 8x8: classi, feature:", len(set(y)), X.shape[1])
print("softmax (multinomiale):", round(soft.score(Xte, yte), 4), " uno-contro-tutti:", round(ovr.score(Xte, yte), 4))

lr = soft[-1]
Z = lr.decision_function(soft[0].transform(Xte[:3]))
print("softmax a mano = predict_proba:", np.allclose(softmax(Z, axis=1), soft.predict_proba(Xte[:3])))
print("somma per riga:", softmax(Z, axis=1).sum(axis=1).round(6).tolist())

# Output:
# cifre 8x8: classi, feature: 10 64
# softmax (multinomiale): 0.9689  uno-contro-tutti: 0.9689
# softmax a mano = predict_proba: True
# somma per riga: [1.0, 1.0, 1.0]

Sulle cifre 8x8 con 10 classi i due schemi hanno la stessa accuratezza (0,9689): su questo problema e con questo seme non c'è differenza, e non la si può attribuire a nessuna regola generale. La softmax calcolata a mano da decision_function coincide con predict_proba, e le righe sommano a 1.

Quando non usare la regressione logistica

• La frontiera tra le classi non è piatta e non puoi costruire misure che la rendano tale.
• Hai molte interazioni fra misure che non conosci in anticipo.
• Il target è un numero e non una categoria: serve la regressione.

Nei primi due casi un albero o un insieme di alberi cattura la struttura senza che tu la specifichi, a costo di spiegabilità; nel terzo il problema è mal formulato. Dove invece conta spiegare ogni previsione, la regressione logistica ha il vantaggio che ogni peso è un moltiplicatore di odds che si può leggere e discutere.

Che cosa sai fare adesso

Sai ricavare la perdita e il gradiente della regressione logistica e verificarli, addestrarla con Newton e con la discesa del gradiente e riconoscere la convenzione di scikit-learn, scegliere una soglia dopo aver letto la matrice di confusione, diffidare dei pesi di classe quando ti servono probabilità, e usare la softmax per più classi. Il capitolo successivo spiega perché un modello può essere troppo semplice o troppo flessibile: 04 - Bias-Variance Tradeoff.

Nella stessa serie trovi anche: 10 - Reti Neurali Feedforward: da perceptron a deep learning e 12 - Deployment ML in produzione: serving e monitoring.