42 - Regressione multivariata, polinomiale e normalizzazione delle feature

Il prezzo di una casa non dipende solo dai metri quadri: contano anche il numero di stanze, l'età dell'edificio, il piano. Immagina il prezzo come l'uscita di una console di mixaggio con tanti cursori: ogni feature è un cursore, e ogni cursore ha una sua «sensibilità», cioè quanto sposta il prezzo quando lo alzi di una tacca. Il modello lineare multivariato è la somma di tutti i cursori, ciascuno moltiplicato per la propria sensibilità. In questa lezione estendi la regressione lineare a molte feature, scopri cosa significano i coefficienti (e quando non significano quello che pensi), vedi perché la normalizzazione delle feature non è un dettaglio, e poi usi un trucco elegante, la regressione polinomiale, per far imparare curve a un modello che resta lineare.

Dal numero al vettore: il modello multivariato

Con n feature x₁, ..., x_n l'ipotesi diventa:

h(x) = θ₀ + θ₁x₁ + θ₂x₂ + ... + θ_n x_n = θ · x      (con x₀ = 1)

Aggiungendo alla riga di ogni esempio un valore fittizio x₀ = 1, l'intercetta diventa un coefficiente come gli altri e l'intera previsione è un prodotto scalare (somma dei prodotti termine a termine) tra il vettore dei parametri θ e il vettore delle feature. Per tutti gli m esempi insieme si scrive h = Xθ, dove X ha forma (m, n+1) (m righe, una per esempio; n+1 colonne: l'uno e le n feature) e θ ha forma (n+1,). Il risultato h ha forma (m,): una previsione per esempio. Controllare le forme a ogni passaggio è l'abitudine che evita metà degli errori.

Il bello è che nulla cambia nella matematica vista nelle lezioni precedenti. La funzione di costo è la stessa, J = (1/2m)‖Xθ − y‖², il gradiente è lo stesso, ∇J = (1/m)Xᵀ(Xθ − y), e l'equazione normale è la stessa, θ = (XᵀX)⁻¹Xᵀy. Cambia solo la dimensione di θ.

Un caso concreto con dati sintetici

Costruiamo un mercato immobiliare artificiale dove conosciamo la verità: il prezzo vale 20 000 + 1 500 per metro quadro + 8 000 per stanza − 600 per anno di età, più un rumore di 15 000 euro di deviazione standard. Le stanze crescono con i metri quadri (case più grandi hanno più stanze): questo introduce una correlazione realistica tra feature. I numeri sono inventati per l'esercizio, non provengono da alcun mercato reale.

import numpy as np

rng = np.random.default_rng(0)
m = 300
area = rng.uniform(30, 200, m)
stanze = np.round(area / 35 + rng.normal(0, 0.6, m)).clip(1, 8)      # correlate con l'area
eta = rng.uniform(0, 60, m)
prezzo = 20000 + 1500 * area + 8000 * stanze - 600 * eta + rng.normal(0, 15000, m)

X = np.column_stack([np.ones(m), area, stanze, eta])                  # forma (300, 4)
theta = np.linalg.solve(X.T @ X, X.T @ prezzo)

print("forma di X:", X.shape, " forma di theta:", theta.shape)
print("coefficienti stimati [intercetta, area, stanze, eta]:", np.round(theta, 1))
print("coefficienti veri                                  : [20000, 1500, 8000, -600]")
print("correlazione tra area e stanze: %.3f" % np.corrcoef(area, stanze)[0, 1])

forma di X: (300, 4)  forma di theta: (4,)
coefficienti stimati [intercetta, area, stanze, eta]: [19831.5  1489.2  7907.1  -569.6]
coefficienti veri                                  : [20000, 1500, 8000, -600]
correlazione tra area e stanze: 0.913

Cosa fa il codice. Le prime righe generano le feature e il target applicando la formula vera più rumore gaussiano. np.column_stack affianca la colonna di uni e le tre feature ottenendo la matrice (300, 4). La soluzione è la stessa equazione normale della lezione precedente. I coefficienti stimati sono vicini ai veri ma non uguali, e lo scarto relativo più grande è sull’età (−570 contro −600, il 5%).

Come si interpreta un coefficiente: θ_area ≈ 1 489 significa «a parità di stanze ed età, un metro quadro in più è associato a circa 1 489 euro in più». La clausola «a parità di tutte le altre» è essenziale: ogni coefficiente è l'effetto di quella feature con le altre tenute ferme. Se le feature sono correlate (qui area e stanze lo sono al 91%), quella clausola diventa fragile, e lo vediamo tra poco.

La normalizzazione delle feature

Le feature hanno unità e scale diverse: l'area va da 30 a 200, le stanze da 1 a 8, l'età da 0 a 60. Come hai visto nella lezione sulla discesa del gradiente, scale diverse rendono la superficie del costo allungata e la discesa lenta o instabile. Le due trasformazioni più comuni:

Standardizzazione (z-score):   x' = (x − μ) / σ           media 0, deviazione standard 1
Min-max:                       x' = (x − min) / (max − min)   valori in [0, 1]

La standardizzazione è la scelta di default: non è distorta da pochi valori estremi quanto il min-max (un solo outlier cambia max e comprime tutti gli altri valori). Ricorda la regola di igiene: μ e σ (o min e max) si calcolano sul training set e si applicano identici al test set; calcolarli sull'intero dataset è data leakage.

mu, sd = X[:, 1:].mean(axis=0), X[:, 1:].std(axis=0)
Z = np.column_stack([np.ones(m), (X[:, 1:] - mu) / sd])              # feature standardizzate (la colonna di uni resta)
theta_z = np.linalg.solve(Z.T @ Z, Z.T @ prezzo)

print("numero di condizione di X^T X (grezze):          %.0f" % np.linalg.cond(X.T @ X))
print("numero di condizione di Z^T Z (standardizzate):  %.1f" % np.linalg.cond(Z.T @ Z))
print("coefficienti standardizzati [intercetta, area, stanze, eta]:", np.round(theta_z, 0))
print("in unità originali (theta_z[j]/sd[j]):", np.round(theta_z[1:] / sd, 1), " intercetta:", np.round(theta_z[0] - (theta_z[1:] * mu / sd).sum(), 1))

numero di condizione di X^T X (grezze):          179421
numero di condizione di Z^T Z (standardizzate):  22.0
coefficienti standardizzati [intercetta, area, stanze, eta]: [211544.  75330.  12258.  -9719.]
in unità originali (theta_z[j]/sd[j]): [1489.2 7907.1 -569.6]  intercetta: 19831.5

Cosa fa il codice. La standardizzazione è applicata colonna per colonna (axis=0) alle sole feature, mantenendo la colonna di uni. Il numero di condizione di una matrice misura quanto è «sbilanciata» (rapporto tra il suo autovalore più grande e il più piccolo): un valore enorme significa che alcune direzioni sono molto più ripide di altre, quindi il problema è numericamente delicato e la discesa del gradiente lenta. Con le feature grezze vale circa 179 000; dopo la standardizzazione scende a 22.

I coefficienti ottenuti sulle feature standardizzate hanno un vantaggio: sono confrontabili. Il coefficiente di «area» (circa 75 000) indica di quanto cambia il prezzo per una deviazione standard di area; quello dell'età (circa −9 700) per una deviazione standard di età. Senza standardizzare, 1 489 e −570 non si possono confrontare, perché i metri quadri e gli anni sono unità diverse. Riconvertendo (ultima riga) si ritrovano i coefficienti in unità originali: la trasformazione non cambia il modello, solo la sua rappresentazione.

La trappola della multicollinearità

Due feature sono collineari se una è (quasi) una combinazione lineare dell'altra. Il modello non riesce a distinguere i loro contributi: può attribuire l'effetto all'una o all'altra, in proporzioni diverse a seconda di piccoli cambiamenti nei dati. La previsione resta buona, ma i coefficienti diventano instabili e non interpretabili. Per misurarlo ricampioniamo i dati con il bootstrap (estrarre m righe con reinserimento) 300 volte, rifacciamo la regressione ogni volta e guardiamo la deviazione standard di ciascun coefficiente. Poi aggiungiamo una feature quasi duplicata, «area in un'altra unità» (l'area più un piccolo rumore), e ripetiamo:

def bootstrap_coefficienti(X, y, B=300, seed=1):
    r = np.random.default_rng(seed)
    out = []
    for _ in range(B):
        i = r.integers(0, len(y), len(y))
        out.append(np.linalg.lstsq(X[i], y[i], rcond=None)[0])
    return np.array(out)

b_normale = bootstrap_coefficienti(X, prezzo)
X_dup = np.column_stack([X, area + np.random.default_rng(2).normal(0, 0.5, m)])   # area quasi duplicata
b_dup = bootstrap_coefficienti(X_dup, prezzo)

print("modello normale   media  :", np.round(b_normale.mean(axis=0), 0))
print("modello normale   dev.std:", np.round(b_normale.std(axis=0), 0))
print("con area duplicata media :", np.round(b_dup.mean(axis=0), 0))
print("con area duplicata dev.std:", np.round(b_dup.std(axis=0), 0))
print("somma dei due coefficienti d'area: %.0f (nel modello normale: %.0f)" % (b_dup.mean(axis=0)[1] + b_dup.mean(axis=0)[4], b_normale.mean(axis=0)[1]))

modello normale   media  : [19654.  1491.  7855.  -567.]
modello normale   dev.std: [2709.   49. 1542.   60.]
con area duplicata media : [19846.   758.  7884.  -567.   730.]
con area duplicata dev.std: [2755. 1835. 1532.   60. 1837.]
somma dei due coefficienti d'area: 1489 (nel modello normale: 1491)

Cosa fa il codice. bootstrap_coefficienti estrae indici con reinserimento, rifà la regressione con np.linalg.lstsq (il risolutore ai minimi quadrati di NumPy, documentazione, che accetta anche matrici non di rango pieno, dove l'equazione normale si rompe) e raccoglie i coefficienti. Il confronto è eloquente: con la feature duplicata, la deviazione standard del coefficiente dell'area passa da 49 a 1 835, circa trentasette volte più incerta, e lo stesso vale per la copia (1 837).

Le medie si spartiscono il lavoro in modo arbitrario (758 e 730), ma la loro somma (1 489) resta vicina al valore del modello normale (1 491): il modello sa quanto pesa «l'area», non a chi dei due attribuirlo. Anche le stanze, che sono correlate all'area al 91%, hanno un coefficiente incerto (deviazione standard 1 542 su un valore di circa 7 900).

Rimedi: eliminare o combinare le feature ridondanti, usare la regolarizzazione (tema della lezione su overfitting e regolarizzazione), o ridurre la dimensionalità con la PCA. E regola di interpretazione: non leggere i coefficienti di feature fortemente correlate uno per uno.

Regressione polinomiale: curve con un modello lineare

Una relazione curva, come il prezzo che cresce più che proporzionalmente con l'area, non si adatta bene a una retta. Il trucco è sorprendente nella sua semplicità: se manca la curvatura, creala come nuova feature. Da una sola feature x costruisci x, x², x³, ..., x^d e le dai alla stessa regressione lineare:

h(x) = θ₀ + θ₁x + θ₂x² + θ₃x³ + ... + θ_d x^d

Il modello resta lineare nei parametri θ: la funzione è una curva in x, ma per l'algoritmo di apprendimento i «dati» sono semplicemente le colonne [x, x², x³, ...]. Tutto ciò che hai visto (equazione normale, gradiente, discesa) funziona senza modifiche. Il grado d è l'iperparametro che controlla la flessibilità: d = 1 è una retta, d alto una curva sinuosa.

Il disegno mostra la stessa nuvola di punti (una sinusoide con rumore) adattata con tre gradi. Il grado 1 è sottoadattato (underfitting): troppo rigido, non può rappresentare l'andamento. Il grado 4 segue la curva. Il grado 15 è sovraadattato (overfitting): passa vicino a ogni punto di training inseguendo il rumore, e fa oscillazioni selvagge dove i dati sono radi. Vediamolo con i numeri, usando una divisione training/test e implementando la matrice delle potenze a mano:

from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression

r = np.random.default_rng(1)
x = r.uniform(0, 6, 60)
y = np.sin(x) + r.normal(0, 0.25, 60)                                  # relazione curva + rumore (sigma 0.25, varianza 0.0625)
x_tr, x_te, y_tr, y_te = train_test_split(x, y, test_size=0.5, random_state=1)

def matrice_potenze(x, grado):
    return np.column_stack([x ** k for k in range(1, grado + 1)])       # colonne x, x^2, ..., x^grado

def adatta_a_mano(x, y, grado):
    P = matrice_potenze(x, grado)
    mu, sd = P.mean(axis=0), P.std(axis=0)
    A = np.column_stack([np.ones(len(x)), (P - mu) / sd])               # standardizzo le potenze
    theta = np.linalg.lstsq(A, y, rcond=None)[0]
    previsione = lambda xn: np.column_stack([np.ones(len(xn)), (matrice_potenze(xn, grado) - mu) / sd]) @ theta
    return theta, previsione

# controllo: a mano e con la pipeline di scikit-learn danno le stesse previsioni (grado 3)
_, prev_mano = adatta_a_mano(x_tr, y_tr, 3)
pipe3 = make_pipeline(PolynomialFeatures(3, include_bias=False), StandardScaler(), LinearRegression()).fit(x_tr[:, None], y_tr)
print("differenza massima tra a mano e sklearn (grado 3): %.2e" % np.abs(prev_mano(x_te) - pipe3.predict(x_te[:, None])).max())

print("grado   MSE train   MSE test")
for grado in (1, 2, 3, 4, 5, 7, 10, 15):
    _, prev = adatta_a_mano(x_tr, y_tr, grado)
    print("%5d   %9.4f   %8.4f" % (grado, ((prev(x_tr) - y_tr) ** 2).mean(), ((prev(x_te) - y_te) ** 2).mean()))

differenza massima tra a mano e sklearn (grado 3): 3.18e-14
grado   MSE train   MSE test
    1      0.1757     0.2357
    2      0.1670     0.2067
    3      0.0318     0.0358
    4      0.0318     0.0359
    5      0.0318     0.0358
    7      0.0271     0.0340
   10      0.0232     0.1026
   15      0.0196     3.4167

Cosa fa il codice. matrice_potenze costruisce le colonne x¹, ..., x^d. adatta_a_mano le standardizza (le potenze alte hanno scale enormi: 6¹⁵ è dell'ordine di 5·10¹¹, e senza standardizzazione la matrice è praticamente singolare), risolve ai minimi quadrati e restituisce una funzione di previsione che riapplica la stessa standardizzazione ai nuovi punti. Il confronto con la pipeline di scikit-learn (PolynomialFeatures fa le potenze, StandardScaler le scala) conferma le stesse previsioni fino all'arrotondamento.

Leggi la tabella con attenzione. L'errore sul training scende in modo monotono al crescere del grado (da 0,176 a 0,020): più flessibilità, più il modello si adatta ai dati che ha visto. L'errore sul test invece scende fino al grado 3 (da 0,236 a 0,036), resta piatto fino al grado 7 (0,034), e poi esplode: 0,103 al grado 10 e 3,42 al grado 15, oltre cento volte l’errore di training. Un confronto basato solo sull'errore di training sceglierebbe sempre il grado più alto; il test dice l'opposto, e la forma a «U» dell'errore di test è la firma dell'overfitting. Una cautela: il test set ha 30 punti, quindi differenze piccole (0,034 contro 0,036) non sono significative. Anche il fatto che l'errore di test sia sotto la varianza del rumore (0,0625) è un effetto del campione piccolo, non una violazione del limite teorico.

Estrapolare con i polinomi è pericoloso

modello9 = make_pipeline(PolynomialFeatures(9, include_bias=False), StandardScaler(), LinearRegression()).fit(x_tr[:, None], y_tr)
for xn in (6.5, 7.0, 8.0):
    print("x = %.1f   previsione grado 9: %10.2f   valore vero sin(x) = %6.3f" % (xn, modello9.predict([[xn]])[0], np.sin(xn)))
print("coefficiente massimo (in valore assoluto) del grado 9:", np.round(np.abs(modello9[-1].coef_).max(), 0))

x = 6.5   previsione grado 9:     -13.61   valore vero sin(x) =  0.215
x = 7.0   previsione grado 9:     -67.92   valore vero sin(x) =  0.657
x = 8.0   previsione grado 9:    -663.09   valore vero sin(x) =  0.989
coefficiente massimo (in valore assoluto) del grado 9: 10946.0

Cosa fa il codice: addestra un polinomio di grado 9 e lo interroga appena fuori dall'intervallo dei dati (0-6). Un polinomio di grado alto ha un termine dominante x⁹ che esplode fuori dall'intervallo di training: già a x = 6,5 la previsione è −13,6 contro un valore vero di 0,215; a x = 7 è −67,9 e a x = 8 è −663, mentre il seno non esce mai da [−1, 1]. Anche i coefficienti sono enormi (il massimo supera 10 000), un altro sintomo di instabilità. La retta del grado 1 estrapola in modo noioso ma limitato; il polinomio di grado alto no.

Trade-off

• Più feature: più informazione, ma più rischio di collinearità, di overfitting e di costo computazionale (la matrice XᵀX è (n+1)×(n+1)).
• Grado polinomiale: più grado = più flessibilità = più varianza. Si sceglie con la validazione incrociata, non con l'errore di training.
• Equazione normale o discesa del gradiente: la prima è esatta ma ha costo circa cubico nel numero di feature (inversione di una matrice n×n) e non funziona per modelli non lineari; la seconda scala a molte feature e a dataset enormi, ma richiede feature standardizzate e un learning rate scelto con cura.

Le tre voci sono lo stesso compromesso tra capacità del modello e affidabilità della stima, visto da tre lati. Più feature e più grado aumentano la flessibilità e, con essa, la varianza; l'equazione normale o la discesa non cambiano ciò che il modello può imparare, ma quanto costa arrivarci. La regola pratica che ne esce è procedere per aggiunte misurate: una feature o un grado in più alla volta, ciascuno giudicato dalla validazione e non dall'errore sul training.

Trappole reali

• Interpretare un coefficiente senza la clausola «a parità di altre feature», specie con feature correlate.
• Scaler addestrato sull'intero dataset (data leakage). Va dentro una Pipeline o calcolato solo sul training.
• Dimenticare di applicare lo stesso scaler al test e in produzione: un modello addestrato su dati standardizzati riceve dati grezzi e produce numeri assurdi senza alcun errore.
• Scegliere il grado guardando il training (o il test, usato più volte).
• Estrapolare con polinomi di grado alto.
• Variabili categoriche trattate come numeri (per esempio il quartiere codificato 1, 2, 3): inventano un ordine che non esiste; si usa la codifica one-hot (una colonna 0/1 per categoria; in scikit-learn la fa OneHotEncoder, che con il parametro drop toglie la categoria di riferimento).
• Shape mismatch: PolynomialFeatures e scikit-learn vogliono X di forma (m, n); un vettore (m,) va trasformato con x[:, None] o reshape(-1, 1).

Queste sette trappole hanno una radice comune: il modello non protesta mai. Un coefficiente mal interpretato, uno scaler dimenticato, un grado scelto sul training producono numeri plausibili e nessun errore. Per questo gli esperimenti di questa lezione guardano i coefficienti e la validazione e non soltanto l'R², che nell'esempio delle zone differiva di meno di mezzo punto tra la codifica giusta e quella sbagliata.

Variabili categoriche: perché il codice numerico inganna e come funziona il one-hot

Un modello lineare moltiplica ogni feature per un coefficiente: capisce solo quantità. Una variabile come la «zona» di una casa (centro, semicentro, periferia) è una categoria, non una quantità. La tentazione è di codificarla 0, 1, 2. Questo impone però una struttura che i dati non hanno: dice che il passaggio da zona 0 a zona 1 vale esattamente quanto quello da zona 1 a zona 2, e che la zona 2 è «il doppio» della zona 1. La soluzione è il one-hot encoding: una colonna 0/1 per ogni categoria («la casa è in questa zona?»). Il modello impara allora un contributo separato per ciascuna. C'è una sottigliezza: se tieni tutte le colonne e anche l'intercetta, le colonne sommate valgono sempre 1, cioè esattamente la colonna dell'intercetta, e questa dipendenza esatta è multicollinearità perfetta (la trappola vista prima): per questo se ne toglie una (la «categoria di riferimento»). Facciamo l'esperimento con un prezzo generato come 50 + 2·metri quadri + effetto zona (0, 30 e 80) + rumore di deviazione 15.

from sklearn.model_selection import KFold, cross_val_score
from sklearn.linear_model import Ridge

rng_p = np.random.default_rng(0)
n_p = 300
zone = rng_p.integers(0, 3, n_p)
size = rng_p.uniform(40, 140, n_p)
zone_effect = np.array([0, 30, 80])[zone]
price = 50 + 2.0 * size + zone_effect + rng_p.normal(0, 15, n_p)

coded = LinearRegression().fit(np.column_stack([size, zone]), price)
one_hot = np.eye(3)[zone]
dummies = LinearRegression().fit(np.column_stack([size, one_hot[:, 1:]]), price)

print("codifica numerica 0/1/2:  coefficiente metri quadri %.3f   effetto per ogni «gradino» di zona %.2f   R² %.4f" % (coded.coef_[0], coded.coef_[1], coded.score(np.column_stack([size, zone]), price)))
print("one-hot (zona 0 di riferimento): metri quadri %.3f   zona 1 %.2f   zona 2 %.2f   intercetta %.2f   R² %.4f" % (dummies.coef_[0], dummies.coef_[1], dummies.coef_[2], dummies.intercept_, dummies.score(np.column_stack([size, one_hot[:, 1:]]), price)))
print("valori veri: metri quadri 2,000   zona 1 30   zona 2 80   intercetta 50")

full = LinearRegression().fit(np.column_stack([size, one_hot]), price)
rank = np.linalg.matrix_rank(np.column_stack([np.ones(n_p), size, one_hot]))
print("\ntutte e tre le colonne + intercetta: rango %d su 5 colonne; coefficienti delle tre zone %s (nessuno ha un significato da solo)" % (rank, np.round(full.coef_[1:], 2)))

codifica numerica 0/1/2:  coefficiente metri quadri 1.987   effetto per ogni «gradino» di zona 40.06   R² 0.9483
one-hot (zona 0 di riferimento): metri quadri 1.989   zona 1 30.98   zona 2 79.44   intercetta 50.58   R² 0.9520
valori veri: metri quadri 2,000   zona 1 30   zona 2 80   intercetta 50

tutte e tre le colonne + intercetta: rango 4 su 5 colonne; coefficienti delle tre zone [-36.8   -5.83  42.63] (nessuno ha un significato da solo)

Cosa fa il codice. Genera 300 case con tre zone e adatta due modelli. Con la codifica numerica 0/1/2 il coefficiente dei metri quadri è 1,987 (vero: 2) ma l'effetto della zona viene stimato come un unico valore, 40,06 per ogni «gradino»: significa che la zona 1 vale 40,06 più della zona 0 (vero: 30) e che la zona 2 ne vale 80,12 (vero: 80). Il modello deve far stare due effetti diversi su una retta, e sbaglia la zona intermedia di circa 10 unità. Il one-hot, invece, recupera 30,98 e 79,44 (veri: 30 e 80) con un'intercetta di 50,58 (vera: 50).

Guarda l'R²: 0,9483 contro 0,9520, una differenza di soli 0,4 punti. L'R² non ti avrebbe mai segnalato l'errore: la codifica numerica «funziona» in senso statistico, ma attribuisce alla zona intermedia un prezzo sbagliato di quasi un terzo. È un esempio di errore che si vede nei coefficienti e non nelle metriche aggregate.

L'ultima riga mostra la trappola delle variabili fittizie: con tutte e tre le colonne più l'intercetta il rango è 4 su 5 colonne, il sistema è singolare, e i coefficienti restituiti (−36,8; −5,83; 42,63) sono una delle infinite soluzioni equivalenti: nessuno ha un significato da solo, ma le loro differenze sì (−5,83 − (−36,8) = 30,97, la stessa stima della zona 1 rispetto alla zona 0).

Interazioni: quando l'effetto di una feature dipende da un'altra

Un modello lineare additivo assume che l'effetto di ogni feature sia lo stesso qualunque sia il valore delle altre. A volte non è così: l'effetto dei metri quadri sul prezzo può essere diverso in centro e in periferia. Si cattura aggiungendo come feature il prodotto delle due (termine di interazione): il modello resta lineare nei coefficienti, ma la pendenza rispetto a x ora dipende dal gruppo. Generiamo dati con due gruppi che hanno una pendenza diversa (+0,5 per il gruppo 0, e 0,5 − 1,5 = −1 per il gruppo 1) e confrontiamo i due modelli con la validazione incrociata.

x_i = rng_p.uniform(0, 10, n_p)
group = rng_p.integers(0, 2, n_p)
y_i = 1 + 0.5 * x_i + group * (2 - 1.5 * x_i) + rng_p.normal(0, 1, n_p)    # il gruppo 1 ha intercetta +2 e pendenza −1,5 in più

folds = KFold(5, shuffle=True, random_state=0)
additive = np.column_stack([x_i, group])
with_interaction = np.column_stack([x_i, group, x_i * group])
print("R² in validazione incrociata, modello additivo:           %.3f" % cross_val_score(LinearRegression(), additive, y_i, cv=folds, scoring="r2").mean())
print("R² in validazione incrociata, con il termine x·gruppo:    %.3f" % cross_val_score(LinearRegression(), with_interaction, y_i, cv=folds, scoring="r2").mean())
fit_i = LinearRegression().fit(with_interaction, y_i)
print("coefficienti con interazione: intercetta %.3f   x %.3f   gruppo %.3f   x·gruppo %.3f   (veri: 1, 0,5, 2, −1,5)" % (fit_i.intercept_, *fit_i.coef_))

R² in validazione incrociata, modello additivo:           0.585
R² in validazione incrociata, con il termine x·gruppo:    0.913
coefficienti con interazione: intercetta 1.146   x 0.486   gruppo 1.693   x·gruppo -1.450   (veri: 1, 0,5, 2, −1,5)

Cosa fa il codice. Confronta, con la validazione incrociata a 5 fold, un modello additivo (x e gruppo) con uno che aggiunge il prodotto x·gruppo. L'R² in validazione passa da 0,585 a 0,913: il modello additivo non può rappresentare due rette con pendenze diverse e ne adatta una media, mentre l'interazione lo permette. I coefficienti dell'interazione (intercetta 1,146; x 0,486; gruppo 1,693; x·gruppo −1,450) sono vicini a quelli veri (1; 0,5; 2; −1,5).

Per leggerli: la pendenza di x nel gruppo 0 è 0,486; nel gruppo 1 è 0,486 − 1,450 = −0,964; ed ecco il senso del termine di interazione, è la differenza tra le pendenze dei due gruppi.

Attenzione a due cose. Primo, con un'interazione il coefficiente di «x» da solo non è più «l'effetto di x», ma l'effetto di x quando il gruppo è 0: se non lo tieni presente leggerai male il modello. Secondo, le interazioni sono molte (tra p feature ce ne sono p(p − 1)/2 a coppie), e ognuna aggiunta senza una ragione è un grado di libertà in più per l'overfitting. Aggiungile guidato dalla conoscenza del dominio o da un'esplorazione dei residui, e verifica sempre il guadagno in validazione, come qui.

Scegliere il grado del polinomio: addestramento, validazione e verità

x_p = rng_p.uniform(-3, 3, 60)
y_p = np.sin(x_p) + 0.3 * rng_p.normal(size=60)
x_grid = np.linspace(-3, 3, 400)                                            # griglia fitta per misurare l'errore rispetto alla curva vera sin(x)

print("grado   MSE su training   MSE in validazione incrociata   MSE rispetto alla curva vera")
for degree in (1, 2, 3, 5, 7, 9, 12, 15):
    pipe = make_pipeline(PolynomialFeatures(degree), StandardScaler(), LinearRegression())
    cv_mse = -cross_val_score(pipe, x_p[:, None], y_p, cv=folds, scoring="neg_mean_squared_error").mean()
    pipe.fit(x_p[:, None], y_p)
    print("%-6d   %-14.3f   %-28.3f   %.3f" % (degree, np.mean((pipe.predict(x_p[:, None]) - y_p) ** 2), cv_mse, np.mean((pipe.predict(x_grid[:, None]) - np.sin(x_grid)) ** 2)))

grado   MSE su training   MSE in validazione incrociata   MSE rispetto alla curva vera
1        0.225            0.242                          0.176
2        0.215            0.261                          0.182
3        0.082            0.093                          0.014
5        0.080            0.104                          0.013
7        0.076            0.103                          0.018
9        0.071            0.135                          0.028
12       0.064            0.499                          0.045
15       0.058            2.559                          0.069

Cosa fa il codice. Genera 60 punti dalla curva sin(x) con rumore 0,3 (varianza 0,09), adatta polinomi di grado crescente e riporta tre errori: sul training, in validazione incrociata (5 fold) e rispetto alla curva vera su una griglia fitta, quantità che nella pratica non hai ma che qui possiamo calcolare perché conosciamo la verità. L'errore sul training scende sempre, da 0,225 con il grado 1 a 0,058 con il grado 15: più flessibilità, più adattamento ai punti, per costruzione. L'errore in validazione ha un minimo a grado 3 (0,093), sale lentamente a 0,104 con grado 5 e 0,103 con grado 7, poi esplode: 0,135 con grado 9, 0,499 con grado 12, 2,559 con grado 15. Lo scarto tra training e validazione è l'overfitting che si apre.

L'errore rispetto alla curva vera dà la stessa indicazione (minimo 0,013–0,014 per i gradi 3 e 5), ma con un'importante differenza: a grado 15 è solo 0,069, non 2,559. La validazione incrociata esagera il guasto, e ci sono due spiegazioni plausibili, che qui non ho isolato con un esperimento: ogni fold di validazione contiene punti ai bordi del range che il polinomio, addestrato sugli altri fold, deve estrapolare (dove i polinomi di grado alto sono pessimi, come visto), e la stima di MSE è dominata da pochi punti con errori enormi. Non sarebbe un errore del metodo, ma un'indicazione che con pochi dati e modelli molto flessibili la stima di CV è ad alta varianza.

La conclusione operativa rimane valida: il grado giusto è 3, scelto senza guardare la curva vera, e i gradi tra 3 e 7 sono quasi equivalenti; scegli il più semplice dentro la zona piatta.

Misurare la multicollinearità: il VIF e il suo costo sui coefficienti

Il fattore di inflazione della varianza (VIF) di una feature è 1 / (1 − R²_j), dove R²_j è l'R² della regressione di quella feature su tutte le altre. Se la feature è prevedibile dalle altre, R²_j è vicino a 1 e il VIF esplode: il coefficiente ha una varianza pari a VIF volte quella che avrebbe se la feature fosse indipendente. Creiamo tre feature: x₁, x₂ ≈ x₁ con un piccolo rumore (deviazione 0,1), e x₃ indipendente, con bersaglio y = x₁ + x₃ + rumore.

x1 = rng_p.normal(size=500)
x2 = x1 + rng_p.normal(0, 0.1, 500)                                         # quasi una copia di x1
x3 = rng_p.normal(size=500)
X_v = np.column_stack([x1, x2, x3])
y_v = x1 + x3 + rng_p.normal(0, 1, 500)

for j in range(3):
    others = np.delete(X_v, j, axis=1)
    r2_j = LinearRegression().fit(others, X_v[:, j]).score(others, X_v[:, j])
    print("VIF della feature x%d: %.1f" % (j + 1, 1 / (1 - r2_j)))

boot = np.array([LinearRegression().fit(X_v[i], y_v[i]).coef_ for i in (np.random.default_rng(s).integers(0, 500, 500) for s in range(200))])
print("coefficienti (media su 200 ricampionamenti):   %s   (verità: 1, 0, 1)" % np.round(boot.mean(axis=0), 3))
print("deviazione standard dei coefficienti:           %s" % np.round(boot.std(axis=0), 3))
print("deviazione standard della somma dei coefficienti di x1 e x2: %.3f" % (boot[:, 0] + boot[:, 1]).std())

VIF della feature x1: 112.8
VIF della feature x2: 112.8
VIF della feature x3: 1.0
coefficienti (media su 200 ricampionamenti):   [ 1.048 -0.072  1.032]   (verità: 1, 0, 1)
deviazione standard dei coefficienti:           [0.481 0.477 0.043]
deviazione standard della somma dei coefficienti di x1 e x2: 0.046

Cosa fa il codice. Per ciascuna feature calcola l'R² della regressione sulle altre due e ne ricava il VIF. x₁ e x₂ hanno un VIF di 112,8 ciascuna (la loro correlazione è praticamente totale), x₃ ha 1,0, il minimo possibile: è indipendente dalle altre. Come regola empirica (convenzione diffusa, non un teorema) si considera preoccupante un VIF sopra 5–10; qui siamo oltre venti volte la soglia più severa di 5.

Il costo si legge nell'ultima parte, che ricampiona i dati 200 volte e riadatta il modello. I coefficienti medi sono 1,048, −0,072 e 1,032 (verità: 1, 0, 1): in media il modello non sbaglia, ed è per questo che la multicollinearità è insidiosa, perché le previsioni sono buone.

Ma la deviazione standard dei coefficienti di x₁ e x₂ è 0,481 e 0,477, mentre per x₃ è 0,043: undici volte più grande, e √112,8 ≈ 10,6, esattamente come prevede la teoria (la deviazione standard cresce con la radice del VIF). La riga aggiunta alla fine lo conferma: la deviazione standard della somma dei due coefficienti è 0,046, dieci volte più piccola di quella di ciascuno (circa 0,48). Il modello sa di quanto dipende y dalla coppia, ma non può distribuirlo tra le due, e in un singolo campione i due coefficienti si compensano a vicenda. Se ti interessa solo prevedere, puoi ignorare il problema; se interpreti i coefficienti («quanto pesa x₁?»), le conclusioni sono inaffidabili, e conviene togliere una delle due feature, combinarle, o usare la regolarizzazione, come in overfitting e regolarizzazione.

Tre esercizi

Esercizio 1. Nell'esperimento delle zone, cambia la categoria di riferimento (usa la zona 2 invece della zona 0): come cambiano intercetta e coefficienti, e che cosa resta invariato? Perché una scelta arbitraria di riferimento non cambia le previsioni? Esercizio 2. Sostituisci la regressione lineare con Ridge (alpha = 1) nell'esperimento VIF: la deviazione standard dei coefficienti di x₁ e x₂ cambia? Quale compromesso paghi? Esercizio 3. Ripeti la selezione del grado con 600 punti invece di 60: il grado migliore in validazione incrociata resta 3? Che cosa succede alla distanza tra errore su training e in validazione con gradi alti?

Nella stessa serie rimandano a questa lezione anche: Algebra lineare 1: vettori, matrici, sistemi.

Fonti

• scikit-learn, Linear Models - la descrizione della multicollinearità (stima dei minimi quadrati sensibile all'errore sul target) e la sezione sulla regressione polinomiale come modello lineare su funzioni non lineari dei dati, consultata il 2026-11-12
• scikit-learn, Preprocessing data - la standardizzazione (togliere la media e dividere per la deviazione standard), consultata il 2026-11-12
• scikit-learn, OneHotEncoder - la codifica one-hot e il parametro drop, consultata il 2026-11-12
• NumPy, numpy.linalg.lstsq - il risolutore ai minimi quadrati usato nel bootstrap, consultata il 2026-11-12

Le cifre di tutta la lezione sono misure dell'autore: Python 3.14.4, NumPy 2.3.5, scikit-learn 1.8.0, ogni listato eseguito due volte con seme fisso e uscita identica. Le soglie del VIF (5-10) sono una convenzione, non una norma.

Nella prossima lezione la stessa macchina serve a prevedere una categoria invece di un numero: regressione logistica, dalla retta alla probabilità.