Esercizio - Regressione lineare a tratti con k-means

Questa pagina non ha ancora la versione ripasso: qui sotto c'è il testo completo.

Testo (1° appello 2025, parte 2). In un reattore chimico, per ogni lotto si conoscono temperatura (°C) e tempo (minuti) di reazione; l'obiettivo è prevedere la concentrazione finale del prodotto (%). Sono forniti i dati storici (900 lotti). Si sa che il reattore ha operato in 3 punti di lavoro diversi, e che in ciascuno la relazione è lineare ma con coefficienti diversi; il punto di lavoro non ha un'etichetta.

1. Costruire un modello lineare globale $f(x)=\beta_0+\beta_1x_1+\beta_2x_2$ con OLS e valutarne correttamente l'errore di generalizzazione.
2. Costruire il modello a tratti $g(x)=f_{k(x)}(x)$: raggruppare i dati con il k-means ($K=3$) per trovare i punti di lavoro, adattare un modello lineare per ogni cluster, combinarli.
3. Confrontare i due modelli con la tecnica di valutazione più adatta.

Teoria usata: Clustering e k-meansIl clustering raggruppa osservazioni simili senza etichette, come preprocessing (un modello per ogni cluster) o come obiettivo (segmentazione clienti, organizzazione di documenti). K-means: si sceglie $K$, si inizializzano $K$ centroidi, si alterna assegnazione di ogni punto al centroide più vicino e aggiornamento di ogni centroide alla media dei suoi punti, fino a convergenza; minimizza $\mathrm{MSE}{\text{within}}=\frac1N\sum_k\sum{x_i\in C_k}|x_i-\mu_k|^2$ ma solo fino a un minimo locale, quindi dipende dall'inizializzazione. Il numero di cluster si sceglie col metodo del gomito (la dispersione cala sempre, si cerca dove rallenta) o con la gap statistic $\mathrm{Gap}(K)=E[\log W_K^{ref}]-\log W_K$ (si prende il più piccolo $K$ con $\mathrm{Gap}(K)\ge\mathrm{Gap}(K+1)-s_{K+1}$). Il clustering gerarchico agglomerativo parte da un cluster per punto e fonde i due più vicini (linkage single, complete, average, Ward) costruendo un dendrogramma; quello divisivo parte da un solo cluster. Programma di Telecomunicazioni: clustering.Clustering e k-means → (k-means, modelli locali per cluster); Regressione lineareNell'apprendimento supervisionato si impara una funzione $F(x)$ dagli esempi $(x,y)$: regressione se $y$ è continua, classificazione se è categorica. Il modello lineare è $F_\beta(x)=\beta_0+\beta_1x_1+\dots+\beta_px_p=X\beta$ (con una colonna di uni per $\beta_0$) e i parametri si scelgono minimizzando l'errore quadratico medio $\mathrm{MSE}=\frac1n\sum_i(y_i-F_\beta(x_i))^2$, funzione convessa dei parametri. Annullando il gradiente di $J(\beta)=|y-X\beta|^2$ si ottengono le equazioni normali $X^TX\beta=X^Ty$ e la soluzione dei minimi quadrati ordinari $\beta=(X^TX)^{-1}X^Ty$. Il coefficiente di determinazione $R^2=1-SS_{res}/SS_{tot}$ misura la qualità del fit (0 = come la media, negativo = peggio della media). Un modello va valutato su un test set mai usato per addestrare: l'errore sul training è ottimistico e un polinomio di grado alto lo azzera senza generalizzare.Regressione lineare → (OLS, MSE, $R^2$); la valutazione è quella di Overfitting, ridge regression e cross-validationUna buona prestazione sul training non basta: serve stimare quella su dati nuovi. La cross-validation (K-fold: $k$ parti, ciascuna a turno come test, errore medio; Monte Carlo: $k$ divisioni casuali con quota di test $q$; leave-one-out se $k=n$) evita di dipendere da una sola divisione casuale. L'errore atteso si scompone in $\text{bias}^2+\text{varianza}+\sigma^2$: i modelli semplici fanno underfitting (bias alto), quelli complessi overfitting (varianza alta). La regolarizzazione aggiunge alla perdita una penalità: la ridge regression minimizza $|y-X\beta|^2+\lambda\sum_{j\ge1}\beta_j^2$ e ha soluzione $\beta=(X^TX+\lambda\tilde I)^{-1}X^Ty$ (l'intercetta non si penalizza, le feature si standardizzano): riduce i coefficienti, rende l'inversa stabile con feature collineari, e $\lambda$ è un iperparametro scelto con la validazione (cross-validation annidata per non contaminare il test).Overfitting, ridge regression e cross-validation →.

1. Modello globale

Si usano le equazioni normali $\hat\beta=(X^TX)^{-1}X^Ty$ con la colonna di uni. Per valutare in modo corretto si divide in training e test (per esempio $80\%/20\%$, ripetuto 10 volte: Monte Carlo cross-validation) e si calcolano MSE e $R^2$ sul test, non sul training. Con le dimensioni dei dati (temperatura media $50$, tempo $31$, concentrazione media $58{,}3$ con deviazione standard $12{,}6$) i risultati sono:

$$\hat y=1{,}56+0{,}357\cdot\text{temperatura}+1{,}256\cdot\text{tempo},\qquad\mathrm{MSE}_{\text{test}}\approx78{,}8,\quad R^2\approx0{,}50.$$

Il modello spiega solo la metà della varianza: una sola retta non può rappresentare tre relazioni diverse.

2. Modello a tratti

K-means. Si raggruppano i punti $(\text{temperatura},\text{tempo})$ (le feature, non il target: i punti di lavoro dipendono dalle condizioni operative) in $K=3$ cluster. Il k-means è sensibile all'inizializzazione, quindi si fissa un seme o si ripete e si tiene il $W$ minimo: i tre centroidi trovati sono circa $(30;\,29)$, $(50;\,35)$ e $(70;\,29)$ (temperatura, tempo), con $305$, $290$ e $305$ lotti: i tre punti di lavoro alle temperature $30$, $50$ e $70$ °C.

Modelli locali. Per ogni cluster $k$ si adatta un OLS sui soli lotti del cluster: coefficienti $(\beta_0,\beta_{\text{temp}},\beta_{\text{tempo}})$:

• cluster a circa $30$ °C: $(19{,}21;\ 0{,}515;\ 0{,}310)$;
• cluster a circa $50$ °C: $(29{,}98;\ 0{,}610;\ 0{,}385)$;
• cluster a circa $70$ °C: $(107{,}5;\ -0{,}663;\ -0{,}107)$.

I coefficienti sono molto diversi (nel cluster a $70$ °C la concentrazione scende con la temperatura): ecco perché il modello globale fallisce.

Combinazione. $g(x)=f_{k(x)}(x)$ con $k(x)$ il centroide più vicino a $x$ (si assegna ogni nuovo punto al cluster più vicino con la distanza euclidea; è il passo di assegnazione del k-means).

3. Confronto con la validazione ripetuta

Per ogni ripetizione (10 volte, training $720$, test $180$): si assegna ogni punto al centroide più vicino, si adatta un OLS per ogni cluster sul training e si prevede sul test con il modello del cluster di appartenenza. Importante: nel confronto deve valere lo stesso test set per i due modelli, e l'errore del modello a tratti va calcolato sull'insieme completo del test (non si media gli errori per cluster come fa la soluzione del notebook, che dà lo stesso peso a cluster di dimensione diversa).

L'MSE si riduce di circa settanta volte e l'$R^2$ sale da $0{,}50$ a $0{,}99$: usare tre modelli lineari locali, uno per cluster, è molto meglio di un modello lineare globale. Va ricordato che (1) i centroidi sono stati calcolati sull'intero dataset (compreso il test): il k-means non usa il target, quindi non c'è fuga diretta di informazione, ma per un'applicazione vera conviene calcolarli sul solo training; (2) i modelli per cluster hanno meno dati ciascuno (circa $240$ lotti): se i dati fossero pochi il guadagno potrebbe svanire.

Codice

labels, centroids = kmeans(X, k=3)                           # sulle sole feature (o sul solo training)
for tr, te in repeated_splits(900, n=10, p_train=0.8):
    b = ols(X[tr], y[tr]); mse_global = mean((y[te] - predict(X[te], b)) ** 2)
    c_tr = nearest(X[tr], centroids); c_te = nearest(X[te], centroids)
    pred = np.zeros(len(te))
    for k in range(3):
        bk = ols(X[tr][c_tr == k], y[tr][c_tr == k])         # un modello lineare per cluster
        pred[c_te == k] = predict(X[te][c_te == k], bk)
    mse_piecewise = mean((y[te] - pred) ** 2)