40 - Regressione lineare: ipotesi, errore quadratico medio e funzione di costo

Vuoi stimare il prezzo di una casa dai suoi metri quadri. Hai venti annunci, li metti su un foglio con i metri quadri in orizzontale e il prezzo in verticale, e ti viene naturale prendere un righello e tracciare la retta che «passa meglio» in mezzo ai punti. Hai appena fatto, a occhio, una regressione lineare. Quello che manca per farla fare a un computer è una definizione precisa di «passa meglio»: una funzione di costo, un numero che misura quanto una retta è sbagliata. Tutto il resto del modulo, dalla regressione logistica alle reti neurali, è costruito sulla stessa idea: scegli un modello con dei parametri, definisci un costo, cerca i parametri che lo minimizzano.

Questa lezione è il primo mattone. Derivi la funzione di costo da un'ipotesi sul rumore, trovi la retta ottima con la matematica, la verifichi con il codice, e scopri perché i numeri riassuntivi (come R²) possono ingannare.

Le parole che useremo (e che non daremo più per scontate)

Un esempio (o campione, o riga) è una singola osservazione, per esempio un annuncio di casa: ne abbiamo m. La feature è un valore in ingresso che descrive l'esempio (i metri quadri) e si indica con x; il target (o etichetta) è il valore che vogliamo prevedere (il prezzo), si indica con y, ed essendo un numero ci porta nel territorio della regressione. I parametri sono i numeri che il modello impara (θ₀ e θ₁): non li scegli tu, li cerca l'algoritmo. Il modello (o ipotesi) è la regola che trasforma le feature in una previsione, resa specifica da quei parametri.

L'ultimo termine regge tutto il resto: la funzione di costo è un numero, calcolato su tutti i dati, che dice quanto il modello sbaglia; più è basso, meglio è. Addestrare significa trovare i parametri che lo rendono il più basso possibile. Tieni a mente la distinzione fra i tre livelli, perché tornerà in ogni lezione: i dati sono dati, i parametri sono ciò che cambia, il costo è il metro con cui si decide che cosa cambiare.

L'ipotesi: una retta

Con una sola feature x (metri quadri) e un target y (prezzo), l'ipotesi è una funzione lineare con due parametri:

h(x) = θ₀ + θ₁ · x

θ₁ è la pendenza (quanto cresce il prezzo per ogni metro quadro in più), θ₀ l'intercetta (il valore a x = 0). Per ogni coppia (θ₀, θ₁) hai una retta diversa. Imparare significa scegliere la coppia giusta. Per ogni esempio i, l'errore (o residuo) è la differenza tra previsione e valore vero:

e_i = h(x_i) − y_i

Come si legge il disegno: ogni punto blu è un annuncio (metri quadri in orizzontale, prezzo in verticale). La retta nera è l'ipotesi per una scelta di θ₀ e θ₁. I segmenti rossi tratteggiati sono i residui: la distanza verticale tra la previsione sulla retta e il prezzo reale. Una retta migliore è una retta che rende questi segmenti, nel complesso, più corti. Spostare θ₀ fa salire o scendere tutta la retta; cambiare θ₁ la fa ruotare attorno al suo punto sull'asse verticale.

La funzione di costo: errore quadratico medio

Dobbiamo riassumere gli m errori in un solo numero. Sommarli non funziona: errori positivi e negativi si annullano. Si eleva al quadrato e si fa la media:

J(θ₀, θ₁) = (1 / 2m) · Σ_{i=1..m} ( h(x_i) − y_i )²

Questo è l'errore quadratico medio (MSE) con un fattore ½ per comodità: quando derivi il quadrato compare un 2 che il ½ cancella. Il fattore non cambia dove sta il minimo. Ma perché il quadrato e non il valore assoluto? Due ragioni. La prima è pratica: il quadrato è derivabile ovunque e produce una funzione liscia, comoda per l'ottimizzazione. La seconda è più profonda e la derivi qui sotto.

Perché proprio il quadrato: la verosimiglianza gaussiana

Supponi che i dati nascano così: y_i = θ₀ + θ₁x_i + ε_i, dove il rumore ε_i è gaussiano con media 0 e varianza σ², indipendente da un esempio all'altro. Allora la densità di y_i dato x_i è

p(y_i | x_i; θ) = 1/√(2πσ²) · exp( − (y_i − h(x_i))² / (2σ²) )

La verosimiglianza di tutto il dataset è il prodotto delle densità; prendendo il logaritmo (che trasforma i prodotti in somme e non sposta il massimo):

log L(θ) = Σ_i log p(y_i | x_i; θ) = −m · log√(2πσ²) − (1 / 2σ²) · Σ_i (y_i − h(x_i))²

Il primo termine non dipende da θ. Massimizzare la verosimiglianza equivale quindi a minimizzare Σ(y_i − h(x_i))²: i minimi quadrati sono la stima di massima verosimiglianza sotto rumore gaussiano. Questo ti dice anche quando il quadrato è una cattiva scelta: se il rumore ha code pesanti (outlier frequenti), l'ipotesi gaussiana è falsa e il quadrato sovrappesa i valori anomali.

La superficie del costo e la soluzione in forma chiusa

J è una funzione di due variabili (θ₀, θ₁): una superficie a forma di conca, un paraboloide. Ha un solo minimo, nessun minimo locale ingannevole: J è convessa. Per trovarlo si impongono a zero le derivate parziali. Con la notazione matriciale: X è la matrice m × 2 con una colonna di uni (per l'intercetta) e una con gli x; θ il vettore (θ₀, θ₁).

J(θ) = (1/2m) · (Xθ − y)ᵀ(Xθ − y)

∇J(θ) = (1/m) · Xᵀ(Xθ − y) = 0       →       XᵀX θ = Xᵀy       →       θ = (XᵀX)⁻¹ Xᵀ y

Questa è l'equazione normale. Per il caso di una sola feature si semplifica in due formule che puoi calcolare a mano:

θ₁ = Σ(x_i − x̄)(y_i − ȳ) / Σ(x_i − x̄)²  =  cov(x, y) / var(x)
θ₀ = ȳ − θ₁ · x̄

La seconda dice una cosa semplice: la retta ottima passa sempre per il punto medio (x̄, ȳ).

Un esempio a mano con tre punti

Prima di affidarti al computer, fai il calcolo con carta e penna su tre annunci minuscoli: (x, y) = (1, 2), (2, 3), (3, 5). Le medie sono x̄ = 2 e ȳ = 10/3. Calcoliamo i due ingredienti della formula:

Σ(x − x̄)(y − ȳ) = (−1)(2 − 10/3) + (0)(3 − 10/3) + (1)(5 − 10/3) = 4/3 + 0 + 5/3 = 3
Σ(x − x̄)²       = 1 + 0 + 1 = 2

θ₁ = 3 / 2 = 1,5          θ₀ = 10/3 − 1,5 · 2 = 1/3

La retta è h(x) = 1/3 + 1,5x. Le previsioni sono 1,833 / 3,333 / 4,833, quindi i residui sono −0,167 / +0,333 / −0,167 e il costo è J = (1/6)(0,0278 + 0,1111 + 0,0278) = 0,0278. Nota una proprietà: i residui sommano a zero. Non è un caso. Se la somma non fosse zero, potresti alzare o abbassare tutta la retta di una costante e ridurre il costo, quindi non saresti al minimo. Verifichiamo con il codice:

import numpy as np

xs = np.array([1.0, 2.0, 3.0])
ys = np.array([2.0, 3.0, 5.0])

theta1 = ((xs - xs.mean()) * (ys - ys.mean())).sum() / ((xs - xs.mean()) ** 2).sum()
theta0 = ys.mean() - theta1 * xs.mean()
residui = (theta0 + theta1 * xs) - ys

print("theta0 = %.4f (1/3 = %.4f)   theta1 = %.4f" % (theta0, 1 / 3, theta1))
print("residui:", np.round(residui, 4), "  somma dei residui: %.1e" % residui.sum())
print("costo J = %.4f" % ((residui ** 2).sum() / (2 * len(xs))))

theta0 = 0.3333 (1/3 = 0.3333)   theta1 = 1.5000
residui: [-0.1667  0.3333 -0.1667]   somma dei residui: 8.9e-16
costo J = 0.0278

Il codice ripete esattamente i passaggi fatti a mano: due medie, due somme, un rapporto. residui è calcolato come previsione meno valore vero, lo stesso segno della definizione e_i = h(x_i) − y_i. La somma dei residui vale zero a meno di errori di arrotondamento (dell'ordine di 10⁻¹⁶).

carico  x (t)    2     4     6     8
gasolio y (l)   30    34    42    46

Verifica su più dati, in quattro modi

Generiamo dati dalla relazione y = 3 + 2x con rumore gaussiano, e calcoliamo la retta in quattro modi indipendenti:

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
m = 100
x = rng.uniform(0, 10, m)
y = 3 + 2 * x + rng.normal(0, 2, m)                 # verità: theta0 = 3, theta1 = 2, sigma = 2

# 1) formule per una feature
t1 = ((x - x.mean()) * (y - y.mean())).sum() / ((x - x.mean()) ** 2).sum()
t0 = y.mean() - t1 * x.mean()

# 2) equazione normale in forma matriciale
X = np.column_stack([np.ones(m), x])                # shape (m, 2): colonna di uni + x
theta = np.linalg.solve(X.T @ X, X.T @ y)           # meglio di invertire la matrice

# 3) numpy.polyfit   4) scikit-learn
pf = np.polyfit(x, y, 1)[::-1]
lr = LinearRegression().fit(x.reshape(-1, 1), y)

print("formule:          theta0 = %.4f  theta1 = %.4f" % (t0, t1))
print("equazione normale: theta0 = %.4f  theta1 = %.4f" % tuple(theta))
print("polyfit:          theta0 = %.4f  theta1 = %.4f" % tuple(pf))
print("scikit-learn:     theta0 = %.4f  theta1 = %.4f" % (lr.intercept_, lr.coef_[0]))
print("forma di X:", X.shape, " forma di theta:", theta.shape)

formule:          theta0 = 2.9851  theta1 = 1.9739
equazione normale: theta0 = 2.9851  theta1 = 1.9739
polyfit:          theta0 = 2.9851  theta1 = 1.9739
scikit-learn:     theta0 = 2.9851  theta1 = 1.9739
forma di X: (100, 2)  forma di theta: (2,)

Cosa fa il codice, in ordine. rng.uniform(0, 10, m) genera 100 valori di x; y è costruita applicando la retta vera (3 + 2x) e aggiungendo rumore gaussiano con deviazione standard 2: i dati sono sintetici, quindi conosciamo la risposta giusta e possiamo vedere quanto ci arrivano vicini i metodi. Il primo metodo applica le due formule della covarianza. Il secondo costruisce la matrice X con una colonna di uni (così l'intercetta θ₀ diventa un coefficiente come gli altri: x₀ = 1) e risolve XᵀXθ = Xᵀy. Il terzo e il quarto sono le routine di libreria; polyfit restituisce i coefficienti dal grado più alto al più basso, per questo si rovescia l'array con [::-1].

I quattro metodi concordano. I valori stimati sono vicini a quelli veri (3 e 2) ma non uguali: con 100 punti e rumore, la stima ha un'incertezza. Nota anche np.linalg.solve invece di inv: risolvere il sistema è più stabile numericamente che invertire la matrice e moltiplicare.

La superficie del costo, numericamente

Per vedere che la soluzione è davvero il fondo della conca, calcoliamo J su una griglia di valori (θ₀, θ₁) e cerchiamo il minimo con la forza bruta:

def J(t0, t1):
    return ((t0 + t1 * x - y) ** 2).sum() / (2 * m)

g0 = np.linspace(0, 6, 241); g1 = np.linspace(1, 3, 241)
valori = np.array([[J(a, b) for b in g1] for a in g0])
i, j = np.unravel_index(valori.argmin(), valori.shape)

print("minimo sulla griglia: theta0 = %.3f  theta1 = %.3f   J = %.4f" % (g0[i], g1[j], valori[i, j]))
print("soluzione esatta:     theta0 = %.3f  theta1 = %.3f   J = %.4f" % (theta[0], theta[1], J(*theta)))
print("J nei punti vicini:  J(t0+0.5, t1) = %.3f   J(t0, t1+0.2) = %.3f" % (J(theta[0] + 0.5, theta[1]), J(theta[0], theta[1] + 0.2)))

minimo sulla griglia: theta0 = 2.975  theta1 = 1.975   J = 1.8801
soluzione esatta:     theta0 = 2.985  theta1 = 1.974   J = 1.8801
J nei punti vicini:  J(t0+0.5, t1) = 2.005   J(t0, t1+0.2) = 2.665

Il codice definisce J come funzione di (θ₀, θ₁), la valuta su una griglia di 241 × 241 coppie e prende la posizione del valore più piccolo con argmin/unravel_index. La griglia trova 2.975 e 1.975, a un passo dalla soluzione esatta (2.985, 1.974), e J vale 1.8801 in entrambi. Spostando θ₀ di +0,5 il costo sale a 2.005; spostando θ₁ di +0,2 sale a 2.665: la superficie è più ripida nella direzione di θ₁, perché un errore sulla pendenza si amplifica al crescere di x.

Allontanarsi dal minimo in qualunque direzione fa salire J: è il comportamento di una conca. Questa convessità è ciò che permetterà, nella prossima lezione, di scendere verso il minimo con la discesa del gradiente senza rischiare di fermarsi nel posto sbagliato. Per la regressione lineare esiste la soluzione in forma chiusa, ma con molte feature (o con modelli non lineari) non esiste o costa troppo: da qui nasce la necessità dell'ottimizzazione iterativa.

Perché è davvero un minimo

Annullare il gradiente trova un punto stazionario: potrebbe essere un massimo o una sella. Per assicurarsi che sia un minimo si guarda la derivata seconda (la matrice Hessiana). Derivando ancora ∇J = (1/m)Xᵀ(Xθ − y) rispetto a θ si ottiene H = (1/m)XᵀX. Per qualunque vettore v, vᵀHv = (1/m)·‖Xv‖² ≥ 0: una somma di quadrati non è mai negativa. Una Hessiana con questa proprietà (semidefinita positiva) descrive una conca senza selle né massimi: la funzione è convessa e ogni punto stazionario è un minimo globale. Se le colonne di X sono linearmente indipendenti, la disuguaglianza è stretta e il minimo è unico.

Misurare la qualità: MSE, RMSE, MAE, R²

MSE  = (1/m) Σ (y_i − ŷ_i)²                    stessa unità di y, al quadrato
RMSE = √MSE                                    stessa unità di y
MAE  = (1/m) Σ |y_i − ŷ_i|                     meno sensibile agli outlier
R²   = 1 − Σ (y_i − ŷ_i)² / Σ (y_i − ȳ)²       quota di varianza spiegata

R² confronta il tuo modello con il modello più stupido possibile, quello che predice sempre la media ȳ: R² = 1 è perfetto, R² = 0 vuol dire «non meglio della media», R² negativo vuol dire peggio della media (possibile sul test set).

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

pred = X @ theta
mse = ((y - pred) ** 2).mean()
r2 = 1 - ((y - pred) ** 2).sum() / ((y - y.mean()) ** 2).sum()
print("a mano:    MSE %.3f  RMSE %.3f  MAE %.3f  R2 %.4f" % (mse, np.sqrt(mse), np.abs(y - pred).mean(), r2))
print("sklearn:   MSE %.3f  RMSE %.3f  MAE %.3f  R2 %.4f" % (
    mean_squared_error(y, pred), mean_squared_error(y, pred) ** 0.5, mean_absolute_error(y, pred), r2_score(y, pred)))
print("il rumore vero ha sigma = 2, cioè varianza 4")

a mano:    MSE 3.760  RMSE 1.939  MAE 1.559  R2 0.9049
sklearn:   MSE 3.760  RMSE 1.939  MAE 1.559  R2 0.9049
il rumore vero ha sigma = 2, cioè varianza 4

I due calcoli coincidono (MSE 3.760, RMSE 1.939, MAE 1.559, R² 0.9049). L'MSE è vicino alla varianza del rumore, che vale 4 e che nessun modello può eliminare. È addirittura un po' sotto: il modello ha «assorbito» una piccola parte del rumore dei 100 punti su cui è stato addestrato. È un assaggio dell'overfitting (la lezione su overfitting e regolarizzazione tornerà su questo): un MSE basso sul training non prova che il modello sia migliore. Lo diventa solo se calcolato su dati non usati per l'addestramento. L'RMSE (1.94) ha la stessa unità del target e si legge come «errore tipico»; il MAE (1.56) è la media dei valori assoluti ed è meno sensibile agli errori grandi. R² = 0.905 significa che la retta spiega circa il 90% della varianza di y.

La trappola: numeri uguali, dati diversi

Nel 1973 Francis Anscombe pubblicò quattro piccoli dataset costruiti apposta per avere le stesse statistiche riassuntive. Hanno la stessa media, la stessa varianza, la stessa correlazione e quindi la stessa retta di regressione. Ma raccontano storie molto diverse. Ecco i dati e il calcolo:

x123 = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5], dtype=float)
x4   = np.array([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8], dtype=float)
ys = [np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68]),
      np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74]),
      np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73]),
      np.array([6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])]
for k, (xx, yy) in enumerate(zip([x123, x123, x123, x4], ys), 1):
    b1, b0 = np.polyfit(xx, yy, 1)
    r = np.corrcoef(xx, yy)[0, 1]
    print("dataset %d: theta0 = %.2f  theta1 = %.3f  correlazione = %.3f  R2 = %.3f" % (k, b0, b1, r, r ** 2))

dataset 1: theta0 = 3.00  theta1 = 0.500  correlazione = 0.816  R2 = 0.667
dataset 2: theta0 = 3.00  theta1 = 0.500  correlazione = 0.816  R2 = 0.666
dataset 3: theta0 = 3.00  theta1 = 0.500  correlazione = 0.816  R2 = 0.666
dataset 4: theta0 = 3.00  theta1 = 0.500  correlazione = 0.817  R2 = 0.667

Il codice usa un trucco: zip accoppia ogni lista di x con il suo y, tenendo presente che i primi tre dataset condividono gli stessi x e il quarto no. Per ciascuno np.polyfit(xx, yy, 1) restituisce pendenza e intercetta (nell'ordine dal grado più alto, per questo b1, b0) e np.corrcoef la correlazione, il cui quadrato è R² per una sola feature. L'uscita è il fatto sorprendente: le quattro righe sono quasi indistinguibili, e le differenze sull'ultima cifra (0.666 contro 0.667) sono solo arrotondamenti.

I quattro dataset danno praticamente la stessa retta e lo stesso R² (circa 0.67). Eppure: il primo è una nuvola ragionevole attorno alla retta; il secondo è una curva (una parabola) che la retta approssima male; il terzo è una retta perfetta con un solo outlier che ne spinge la pendenza; il quarto ha tutti i punti con lo stesso x tranne uno, che da solo determina la retta. Nessuna statistica riassuntiva distingue questi casi. Disegna sempre i dati e i residui prima di fidarti di un numero.

Gli outlier e il quadrato

Poiché gli errori sono elevati al quadrato, un punto lontano pesa moltissimo: un residuo di 30 contribuisce al costo novecento volte più di un residuo di 1. Per vederlo senza ambiguità riusiamo il dataset di prima (x, y, 100 punti con retta vera 3 + 2x) e ne aggiungiamo uno solo, a x = 9,5 con y = −30, un valore che la retta vera darebbe intorno a 22. Poi confrontiamo i due adattamenti con np.polyfit.

x_out = np.append(x, 9.5); y_out = np.append(y, -30.0)      # un punto con y completamente fuori scala
b1_out, b0_out = np.polyfit(x_out, y_out, 1)
print("senza outlier:  theta0 = %.3f  theta1 = %.3f" % (theta[0], theta[1]))
print("con 1 outlier:  theta0 = %.3f  theta1 = %.3f   (su %d punti)" % (b0_out, b1_out, len(x_out)))

senza outlier:  theta0 = 2.985  theta1 = 1.974
con 1 outlier:  theta0 = 3.689  theta1 = 1.754   (su 101 punti)

Un solo punto su 101 sposta sensibilmente la retta: la pendenza passa da 1.974 a 1.754 e l'intercetta da 2.985 a 3.689, un cambiamento a cui il resto dei 100 punti non ha contribuito in alcun modo. Rimedi: individuare e capire gli outlier (sono errori di misura o casi reali?), usare una perdita robusta come la Huber, o la regressione quantile.

carico  x (t)    2     4     6     8    10
gasolio y (l)   30    34    42    46    20

Il ponte verso le reti neurali

Guarda l'ipotesi con occhi nuovi: h(x) = θ₁x + θ₀ è un neurone con un ingresso, un peso (θ₁), un bias (θ₀) e nessuna funzione di attivazione. Con più ingressi diventa w·x + b. La funzione di costo che hai appena derivato è la stessa che userai per addestrare reti che predicono numeri. Nel resto del modulo cambieranno solo due cose: il modello diventerà più complesso (più strati, attivazioni non lineari) e la soluzione in forma chiusa lascerà il posto alla discesa del gradiente.

Trade-off

La regressione lineare ha tre vantaggi: la soluzione in forma chiusa è veloce, ogni coefficiente si legge come un effetto, e la bassa varianza ne fa una baseline eccellente da battere prima di passare a un modello più complesso. Il prezzo è che rappresenta solo relazioni lineari nelle feature, è sensibile agli outlier, e l'inversione di XᵀX diventa instabile quando le feature sono correlate fra loro (multicollinearità) o troppo numerose. La documentazione di scikit-learn descrive il problema così: se alcune colonne di X sono quasi linearmente dipendenti, la matrice diventa quasi singolare e la stima dei minimi quadrati molto sensibile agli errori casuali del target.

Trappole reali

• R² alto = modello buono: falso (Anscombe). R² cresce sempre aggiungendo feature, anche inutili, sul training.
• Estrapolare: la retta continua all'infinito, la realtà no. Un modello di prezzo per case tra 40 e 150 m² non vale per 1000 m².
• Correlazione e causalità: il coefficiente descrive un'associazione nei dati, non un effetto causale.
• Shape mismatch: scikit-learn vuole X con forma (m, d) anche per una sola feature: x.reshape(-1, 1). Passare un vettore (m,) dà errore, e un vettore colonna (m, 1) al posto di y (m,) può produrre broadcast silenziosi (m, m) nei calcoli a mano.
• Leakage: una feature costruita usando il target (per esempio «prezzo medio del quartiere» calcolato includendo la casa stessa) fa crollare l'errore sul training e mente sul resto.

Le cinque trappole hanno una radice comune: un numero calcolato sul training (R², un coefficiente, un errore) si legge come se descrivesse il mondo. Le sezioni che seguono ne misurano tre: l'incertezza dei coefficienti, i residui che raccontano ciò che R² nasconde, e l'estrapolazione. Le altre due, correlazione contro causalità e leakage, non si risolvono con un calcolo ma con la conoscenza di come i dati sono stati raccolti.

Quanto fidarsi dei coefficienti: errore standard e intervallo di confidenza

Finora abbiamo trovato la retta migliore e basta. Ma i dati sono un campione: con un altro campione dello stesso fenomeno avremmo ottenuto coefficienti leggermente diversi. Quanto diversi? La risposta serve per dire se una pendenza è «davvero positiva» o se potrebbe essere zero. Per la regressione lineare esiste una formula chiusa. Se il rumore ha varianza σ², la matrice di covarianza dei coefficienti stimati è σ² (XᵀX)⁻¹, e σ² si stima dai residui: σ̂² = Σ eᵢ² / (n − p), dove n è il numero di punti e p il numero di coefficienti (il «− p» tiene conto dei gradi di libertà usati per stimare i coefficienti). L'errore standard di ciascun coefficiente è la radice del rispettivo elemento diagonale. Un intervallo di confidenza al 95% è il coefficiente più o meno circa due errori standard (più precisamente il quantile della distribuzione t di Student con n − p gradi di libertà, 1,99 qui).

Verifichiamo la formula in due modi indipendenti: con il bootstrap (ricampionare i punti con reinserimento e rifare la regressione 2 000 volte), e simulando davvero molti esperimenti per vedere con che frequenza l'intervallo contiene il valore vero.

rng_ci = np.random.default_rng(0)
n_ci = 80
x_s = rng_ci.uniform(0, 10, n_ci)
y_s = 2.0 + 1.5 * x_s + rng_ci.normal(0, 2, n_ci)                         # verità: intercetta 2, pendenza 1,5, rumore 2

A = np.column_stack([np.ones(n_ci), x_s])
beta = np.linalg.solve(A.T @ A, A.T @ y_s)
residuals = y_s - A @ beta
sigma2 = residuals @ residuals / (n_ci - 2)                               # stima della varianza del rumore
std_err = np.sqrt(np.diag(sigma2 * np.linalg.inv(A.T @ A)))
t_crit = 1.99                                                             # quantile 97,5% di t con 78 gradi di libertà
print("intercetta %.3f (errore std %.3f)   pendenza %.3f (errore std %.3f)   sigma stimato %.3f (vero 2)" % (beta[0], std_err[0], beta[1], std_err[1], np.sqrt(sigma2)))
print("intervallo al 95%% per la pendenza con la formula: [%.3f, %.3f]" % (beta[1] - t_crit * std_err[1], beta[1] + t_crit * std_err[1]))

boot_slopes = []
for _ in range(2000):
    idx = rng_ci.integers(0, n_ci, n_ci)
    boot_slopes.append(np.linalg.lstsq(A[idx], y_s[idx], rcond=None)[0][1])
print("intervallo al 95%% con il bootstrap:       [%.3f, %.3f]   (deviazione std delle pendenze %.3f)" % (*np.percentile(boot_slopes, [2.5, 97.5]), np.std(boot_slopes)))

intercetta 2.365 (errore std 0.424)   pendenza 1.424 (errore std 0.072)   sigma stimato 1.886 (vero 2)
intervallo al 95% per la pendenza con la formula: [1.282, 1.567]
intervallo al 95% con il bootstrap:       [1.274, 1.578]   (deviazione std delle pendenze 0.078)

Formula e bootstrap concordano su questo campione, ma un solo campione non dice se «95%» sia una promessa mantenuta. Per saperlo bisogna ripetere l'intero esperimento molte volte, con dati nuovi ma stessa verità, e contare quante volte l'intervallo contiene la pendenza vera 1,5.

hits = 0
for r in range(1000):                                                     # 1000 esperimenti indipendenti con la stessa verità
    gen = np.random.default_rng(1000 + r)
    xx = gen.uniform(0, 10, n_ci); yy = 2 + 1.5 * xx + gen.normal(0, 2, n_ci)
    AA = np.column_stack([np.ones(n_ci), xx]); bb = np.linalg.solve(AA.T @ AA, AA.T @ yy)
    ee = yy - AA @ bb; se = np.sqrt(ee @ ee / (n_ci - 2) * np.linalg.inv(AA.T @ AA)[1, 1])
    hits += bb[1] - t_crit * se <= 1.5 <= bb[1] + t_crit * se
print("l'intervallo al 95%% ha contenuto la pendenza vera in %.1f%% dei casi" % (hits / 10))

l'intervallo al 95% ha contenuto la pendenza vera in 94.7% dei casi

Cosa fa il codice. Genera 80 punti da una retta con intercetta 2, pendenza 1,5 e rumore di deviazione 2, risolve l'equazione normale a mano e calcola l'errore standard con la formula σ̂² (AᵀA)⁻¹. La pendenza stimata è 1,424 con errore standard 0,072; la deviazione del rumore stimata dai residui è 1,886 (la vera è 2: un campione di 80 punti non la ritrova esattamente). L'intervallo al 95% con la formula è [1,282; 1,567], quello del bootstrap [1,274; 1,578]: due metodi indipendenti che concordano, e che contengono entrambi il valore vero 1,5.

Il secondo blocco ripete l'intero esperimento 1 000 volte con campioni diversi e conta quante volte l'intervallo contiene la pendenza vera: 94,7%, vicino al 95% promesso. Questo è il significato operativo di «intervallo di confidenza al 95%»: non «la pendenza ha il 95% di probabilità di stare qui», ma «se ripetessi l'esperimento molte volte, il 95% degli intervalli costruiti così conterrebbe il valore vero». Ricorda però che la copertura del 94,7% vale perché i dati erano stati generati rispettando le ipotesi (rumore gaussiano, omoschedastico, relazione davvero lineare). Con dati reali le ipotesi sono approssimate, e la copertura reale sarà inferiore: il bootstrap è più robusto alle ipotesi sulla distribuzione del rumore, ma non a quella sulla forma della relazione.

Guardare i residui: quando la retta sbaglia in modo sistematico

Le formule dell'errore standard assumono che il rumore abbia la stessa varianza ovunque (omoschedasticità). Se non è vero (eteroschedasticità: l'errore cresce con il valore, come accade spesso con prezzi, redditi, conteggi) i coefficienti restano ragionevoli ma gli intervalli sono sbagliati. Il modo per accorgersene è guardare i residui, mai solo R². Generiamo dati in cui la deviazione del rumore cresce con x e misuriamo la dispersione dei residui in quattro fasce.

x_h = rng_ci.uniform(0, 10, 400)
y_h = 2 + 1.5 * x_h + rng_ci.normal(0, 1, 400) * (0.3 + 0.5 * x_h)         # la deviazione del rumore cresce con x
A_h = np.column_stack([np.ones(400), x_h])
res_h = y_h - A_h @ np.linalg.solve(A_h.T @ A_h, A_h.T @ y_h)
print("fascia di x     deviazione standard dei residui")
for lo, hi in ((0, 2.5), (2.5, 5), (5, 7.5), (7.5, 10)):
    mask = (x_h >= lo) & (x_h < hi)
    print("%4.1f - %4.1f     %.3f" % (lo, hi, res_h[mask].std()))
print("correlazione tra |residuo| e x: %.2f" % np.corrcoef(np.abs(res_h), x_h)[0, 1])

fascia di x     deviazione standard dei residui
 0.0 -  2.5     0.927
 2.5 -  5.0     2.084
 5.0 -  7.5     3.081
 7.5 - 10.0     4.819
correlazione tra |residuo| e x: 0.58

Cosa fa il codice. Genera 400 punti su una retta il cui rumore ha deviazione 0,3 + 0,5·x, quindi crescente, adatta la retta e divide i residui in quattro fasce di x. La deviazione standard dei residui è 0,927 nella prima fascia, 2,084 nella seconda, 3,081 nella terza e 4,819 nell'ultima: aumenta di cinque volte da sinistra a destra, e la correlazione tra il valore assoluto del residuo e x è 0,58. Un unico numero di sintesi come R² non lo mostrerebbe; il grafico dei residui contro x (a forma di imbuto) o questa tabella sì.

Le conseguenze sono pratiche: l'errore standard calcolato con una σ̂² unica descrive male un modello in cui l'errore vero è 0,9 in una zona e 4,8 in un'altra; gli intervalli di previsione sono troppo larghi a sinistra e troppo stretti a destra; e i minimi quadrati danno lo stesso peso a tutti i punti, anche a quelli più rumorosi, anziché pesarli meno. I rimedi sono tre: trasformare la variabile risposta (il logaritmo trasforma un errore proporzionale in un errore costante), usare i minimi quadrati pesati (pesi inversamente proporzionali alla varianza), oppure calcolare errori standard robusti all'eteroschedasticità. Nessuno dei tre cambia la retta in modo drammatico: cambia la tua conoscenza di quanto fidarti.

Outlier: minimi quadrati contro regressioni robuste

Nella sezione sugli outlier abbiamo visto che il quadrato dell'errore dà molto peso ai punti lontani. Esistono stimatori costruiti apposta per resistere: la regressione di Huber usa il quadrato per errori piccoli e il valore assoluto per errori grandi (quindi un punto lontano pesa in modo lineare, non quadratico); Theil–Sen, nella versione classica a una dimensione, prende la mediana delle pendenze di tutte le coppie di punti (scikit-learn ne usa una generalizzazione della mediana a più dimensioni); RANSAC adatta ripetutamente una retta a piccoli sottoinsiemi casuali e tiene quella che ha più punti vicini. Corrompiamo 6 punti su 80 aggiungendo 40 al loro valore di y, e confrontiamo con la retta sui dati puliti (intercetta 2,365, pendenza 1,424).

from sklearn.linear_model import HuberRegressor, RANSACRegressor, TheilSenRegressor

y_corrupt = y_s.copy()
y_corrupt[:6] += 40                                                       # 6 punti su 80 diventano outlier
print("stimatore       intercetta   pendenza    (verità: 2,0 e 1,5; sui dati puliti 2,365 e 1,424)")
for name, model in (("OLS", LinearRegression()), ("Huber", HuberRegressor()), ("RANSAC", RANSACRegressor(random_state=0)), ("Theil-Sen", TheilSenRegressor(random_state=0))):
    model.fit(x_s[:, None], y_corrupt)
    fitted = model.estimator_ if name == "RANSAC" else model
    print("%-15s  %8.3f    %8.3f" % (name, np.ravel(fitted.intercept_)[0], np.ravel(fitted.coef_)[0]))

stimatore       intercetta   pendenza    (verità: 2,0 e 1,5; sui dati puliti 2,365 e 1,424)
OLS                 6.519       1.199
Huber               2.879       1.392
RANSAC              3.460       1.225
Theil-Sen           2.730       1.400

Cosa fa il codice. Prende gli stessi 80 punti, aggiunge 40 a y per i primi 6 (il 7,5% dei dati) e adatta quattro stimatori. I minimi quadrati ordinari (OLS) vengono trascinati: la pendenza scende da 1,424 a 1,199 e l'intercetta sale da 2,365 a 6,519, un errore di oltre quattro unità su un parametro che vale circa 2. Huber (pendenza 1,392, intercetta 2,879) e Theil–Sen (1,400 e 2,730) restano vicini ai valori puliti. RANSAC, con i parametri predefiniti, fa meglio dell'OLS ma non bene (1,225 e 3,460): è un risultato onesto da riportare, perché RANSAC è sensibile alla scelta della soglia di «vicinanza» e qui non l'abbiamo regolata.

Non concluderne che «Huber è meglio di OLS»: sui dati puliti OLS è lo stimatore più efficiente, e uno robusto costa un po' di precisione quando gli outlier non ci sono. La scelta dipende da cosa sai dei dati: se gli errori grossi sono errori di misura da non inseguire, un robusto è più prudente; se sono eventi reali e importanti (la transazione enorme che conta davvero), inseguirli con OLS è corretto, e un robusto li nasconderebbe. Il primo passo non è mai scegliere lo stimatore: è andare a vedere i sei punti.

Estrapolare: il modello non sa cosa c'è fuori dai dati

L'ultima trappola della lista, estrapolare, si misura con un esperimento in cui conosciamo la verità. Generiamo dati da una radice quadrata (una curva che si appiattisce), li osserviamo solo per x fra 0 e 5, e adattiamo una retta. Poi chiediamo alla retta che cosa succede a x = 10 e a x = 20, dove i dati non arrivano, e confrontiamo con il valore vero.

x_e = rng_ci.uniform(0, 5, 100)
y_e = 3 + 4 * np.sqrt(x_e) + rng_ci.normal(0, 0.3, 100)                    # verità: una radice quadrata, non una retta
line = LinearRegression().fit(x_e[:, None], y_e)
print("x     previsione della retta   valore vero (senza rumore)   errore")
for x_query in (2.5, 5, 10, 20):
    truth = 3 + 4 * np.sqrt(x_query)
    prediction = line.predict([[x_query]])[0]
    print("%-5g %-22.2f  %-26.2f  %+.2f" % (x_query, prediction, truth, prediction - truth))

x     previsione della retta   valore vero (senza rumore)   errore
2.5   8.98                    9.32                        -0.35
5     12.45                   11.94                       +0.51
10    19.40                   15.65                       +3.75
20    33.29                   20.89                       +12.40

Cosa fa il codice. Genera dati da una relazione che non è una retta ma una radice quadrata (3 + 4·√x), osservata solo per x tra 0 e 5, e adatta una retta. Dentro l'intervallo osservato la retta va bene: a x = 2,5 sbaglia di −0,35 e a x = 5 di +0,51, errori dello stesso ordine del rumore (0,3). Fuori dall'intervallo il modello non ha modo di sapere che la curva si piega: a x = 10 l'errore è +3,75, a x = 20 è +12,40 (33,29 predetto contro 20,89 vero, il 59% in più). Dentro l'intervallo osservato l'adattamento sembra buono: nessuna statistica calcolata sui dati di training (residui, R², errori standard dei coefficienti) ti avverte di ciò che succede fuori. L'errore standard misura l'incertezza dovuta al campione dentro il modello assunto, non l'errore dovuto al modello sbagliato.

Per questo la previsione fuori dal range dei dati di training (estrapolazione) è la fonte classica dei disastri nelle applicazioni: il prezzo di una casa di 1 000 m², la domanda in un mese mai visto. Regola pratica: controlla se i nuovi input cadono dentro l'involucro dei dati di training e, se non cadono, trattala come una domanda a cui il modello non può rispondere, non come una previsione con una barra d'errore più larga.

Tre esercizi

Esercizio 1. Ripeti la simulazione di copertura con la stessa verità ma con n = 10 punti invece di 80: usando ancora 1,99 al posto del quantile corretto della t con 8 gradi di libertà (2,31), quanto scende la copertura? Che cosa insegna sull'uso del 2 «a occhio»? Esercizio 2. Nell'esperimento degli outlier, aumenta da 6 a 20 i punti corrotti (25%): quali stimatori reggono ancora? Esercizio 3. Per i dati eteroschedastici, trasforma y con il logaritmo (dopo aver reso positivi i valori) e riguarda i residui nelle quattro fasce: la dispersione diventa uniforme? Qual è il prezzo di questa trasformazione sull'interpretazione dei coefficienti?

Per esercitarti su un dominio nuovo, il consumo di gasolio di un traghetto fluviale, c'è il file esercizi.py (Python 3.12+ e NumPy 2.x): la retta dei minimi quadrati scritta da zero in forma chiusa, le quattro metriche con un caso di R² negativo, e l'analisi lascia-fuori-uno che trova da sola la riga sbagliata fra tredici. Il file si esegue con python3 esercizi.py anche con gli esercizi non completati e ti dice quali mancano: lanciato così com'è stampa «Esercizi completati: 0 su 3».

Nella prossima lezione, la discesa del gradiente, affronti il caso in cui la forma chiusa non basta: come trovare il minimo di J scendendo lungo la sua pendenza.

Lezione precedente del percorso: Riduzione di dimensionalità: PCA da zero e t-SNE. Nella stessa serie rimandano a questa lezione anche: Statistica: stimare, confrontare, decidere, Reinforcement learning: l'agente che impara per tentativi (e il dilemma esplorazione-sfruttamento).