38 - Clustering: K-Means e DBSCAN senza etichette

Hai una scatola di mille fotografie senza didascalie e devi ordinarle in pile. Nessuno ti ha detto quali pile esistano: «mare», «montagna», «feste» emergono guardando quali foto si somigliano. È il clustering: raggruppare esempi simili senza che nessuno abbia fornito le etichette. È apprendimento non supervisionato, e porta con sé una difficoltà che nella classificazione non c'era: non esiste una risposta corretta contro cui misurarsi. Il risultato dipende dall'idea di «somiglianza» che scegli, e un algoritmo ti darà sempre dei gruppi, anche se nei dati non ce ne sono.

Vediamo due algoritmi con filosofie opposte. K-Means cerca K gruppi compatti attorno a dei centri. DBSCAN cerca regioni dense separate da regioni vuote, e sa riconoscere il rumore.

K-Means: l'obiettivo

Si fissa il numero di gruppi K. Ogni gruppo C_k ha un centro μ_k (centroide). L'algoritmo cerca assegnazioni e centri che minimizzino l'inerzia, la somma dei quadrati delle distanze di ogni punto dal centro del proprio gruppo:

J = Σ_{k=1..K} Σ_{x ∈ C_k} ‖x − μ_k‖²

Il ciclo è lo schema da tenere a mente: assegna, aggiorna, ripeti. Un centroide è il «punto medio» di un gruppo, la media aritmetica delle coordinate dei suoi punti; non è necessariamente uno dei punti dei dati.

L'algoritmo di Lloyd e perché converge

1. Scegli K centri iniziali.
2. Assegnazione: ogni punto va al centro più vicino.
3. Aggiornamento: ogni centro diventa la media dei punti assegnati.
4. Ripeti 2-3 finché le assegnazioni non cambiano più.

I passi 2 e 3 sono un'alternanza, non due fasi indipendenti: il primo fissa i centri e sceglie le assegnazioni, il secondo fissa le assegnazioni e sceglie i centri. Ognuno minimizza J rispetto a una sola delle due incognite, ed è questo che permette la dimostrazione di convergenza che segue. Il passo 1 è l'unico non deterministico, e come vedrai conta più di quanto sembri.

Perché la media? Fissate le assegnazioni, vogliamo il μ che minimizza Σ_x ‖x − μ‖² in un gruppo. Derivando rispetto a μ e ponendo a zero:

∂/∂μ Σ_x ‖x − μ‖² = −2 Σ_x (x − μ) = 0      →      μ = (1/|C|) Σ_x x

Ogni passo non può aumentare J: l'assegnazione sceglie, per ogni punto, il centro che minimizza la sua distanza (dunque J non cresce a centri fissi), e l'aggiornamento sceglie il centro che minimizza J a assegnazioni fisse. J non è mai in aumento e può assumere un numero finito di valori, quindi l'algoritmo termina. Ma termina in un minimo locale, non per forza globale: dipende dai centri iniziali.

Inizializzazione: k-means++

Partire con centri scelti a caso può dare soluzioni pessime (due centri nello stesso gruppo vero, un gruppo vero senza centro). Il k-means++ sceglie il primo centro a caso, poi ogni successivo tra i punti con probabilità proporzionale a D(x)², dove D(x) è la distanza dal centro più vicino già scelto: i centri si allontanano tra loro.

Implementazione a mano

import numpy as np
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score, silhouette_score

X, y_vero = make_blobs(n_samples=600, centers=4, cluster_std=[1.0, 1.5, 0.5, 1.0], random_state=7)

def dist2(X, C):
    return ((X[:, None, :] - C[None, :, :]) ** 2).sum(axis=2)      # (n, K)

def init_pp(X, k, rng):
    centri = [X[rng.integers(len(X))]]
    for _ in range(k - 1):
        d = dist2(X, np.array(centri)).min(axis=1)                   # distanza al centro più vicino
        centri.append(X[rng.choice(len(X), p=d / d.sum())])          # prob. proporzionale a D(x)^2
    return np.array(centri)

def kmeans(X, k, init="pp", seed=0, max_iter=100):
    rng = np.random.default_rng(seed)
    C = init_pp(X, k, rng) if init == "pp" else X[rng.choice(len(X), k, replace=False)]
    for _ in range(max_iter):
        etichette = dist2(X, C).argmin(axis=1)                        # passo di assegnazione
        nuovi = np.array([X[etichette == j].mean(axis=0) if (etichette == j).any() else C[j]
                          for j in range(k)])                         # passo di aggiornamento
        if np.allclose(nuovi, C):
            break
        C = nuovi
    inerzia = dist2(X, C).min(axis=1).sum()
    return C, etichette, inerzia

C, lab, J = kmeans(X, 4)
sk = KMeans(n_clusters=4, n_init=10, random_state=0).fit(X)
print("inerzia a mano: %.1f    sklearn: %.1f" % (J, sk.inertia_))
print("accordo con i gruppi veri (ARI) a mano: %.3f   sklearn: %.3f" % (
    adjusted_rand_score(y_vero, lab), adjusted_rand_score(y_vero, sk.labels_)))

inerzia a mano: 1295.9    sklearn: 1295.9
accordo con i gruppi veri (ARI) a mano: 0.996   sklearn: 0.996

Nel listato dist2 calcola le distanze al quadrato tra tutti i punti e tutti i centri con il broadcasting, restituendo una matrice (n, K). init_pp implementa k-means++: sceglie il primo centro a caso, poi per ogni centro successivo calcola per ogni punto la distanza al centro più vicino già scelto e sceglie un punto con probabilità proporzionale a quel valore (rng.choice(..., p=d/d.sum())). kmeans alterna i due passi: argmin per l'assegnazione, la media del gruppo per l'aggiornamento (se un gruppo resta vuoto il centro non si sposta); si ferma quando i centri non cambiano (np.allclose). Ritorna anche l'inerzia finale. I dati sono 4 gruppi gaussiani con deviazioni diverse (0,5-1,5), e y_vero contiene i gruppi che abbiamo generato.

L'ARI (Adjusted Rand Index) misura l'accordo tra due partizioni: 1 se coincidono, circa 0 se sono indipendenti. Lo possiamo usare qui solo perché abbiamo generato noi i dati e conosciamo i gruppi veri: nella realtà non li conosci, ed è esattamente il problema del clustering.

Un solo avvio non basta: minimi locali

for init in ("casuale", "pp"):
    inerzie = [kmeans(X, 4, init="pp" if init == "pp" else "random", seed=s)[2] for s in range(30)]
    print("init %-8s  30 avvii: inerzia migliore %.1f   peggiore %.1f   avvii nel minimo migliore: %d/30" % (
        init, min(inerzie), max(inerzie), sum(abs(v - min(inerzie)) < 1.0 for v in inerzie)))

init casuale   30 avvii: inerzia migliore 1295.9   peggiore 14265.5   avvii nel minimo migliore: 19/30
init pp        30 avvii: inerzia migliore 1295.9   peggiore 1295.9   avvii nel minimo migliore: 30/30

Il listato ripete l'algoritmo con 30 semi diversi, con inizializzazione casuale (X[rng.choice(len(X), k, replace=False)]: k punti presi a caso come centri) oppure k-means++, e conta quante volte si raggiunge l'inerzia migliore (con tolleranza di 1.0).

Con centri iniziali casuali, su 30 avvii solo 19 raggiungono l'inerzia migliore (1295.9) e il peggiore finisce a 14265.5, undici volte peggio. Con k-means++ tutti i 30 avvii arrivano al minimo migliore su questo dataset. Avvii diversi finiscono quindi in soluzioni diverse. La regola pratica è quella che applica scikit-learn: lancia l'algoritmo più volte (n_init) con k-means++ e tieni la soluzione con inerzia minore.

centro      : (0,0)   (1,0)   (0,1)
università : (10,0)  (11,0)  (10,1)
stazione FS : (20,0)  (21,0)  (20,1)

Come si sceglie K?

L'inerzia scende sempre aumentando K (con K = n vale zero), quindi non si può minimizzare: si cerca il «gomito», il punto oltre il quale i guadagni diventano marginali. Il gomito spesso è ambiguo. Un criterio più informativo è il coefficiente di silhouette. Per un punto i, sia a(i) la distanza media dagli altri punti del proprio gruppo e b(i) la distanza media dal gruppo più vicino tra gli altri:

s(i) = ( b(i) − a(i) ) / max( a(i), b(i) )         valori in [−1, +1]

Vicino a +1: il punto è ben dentro il suo gruppo e lontano dagli altri. Vicino a 0: è al confine. Negativo: probabilmente è nel gruppo sbagliato. La silhouette del clustering è la media su tutti i punti. Qui sotto la calcoliamo a mano e la usiamo, insieme all'inerzia, per scegliere K.

def silhouette_a_mano(X, lab):
    D = np.sqrt(dist2(X, X))
    s = np.zeros(len(X))
    for i in range(len(X)):
        stesso = (lab == lab[i]); stesso[i] = False
        a = D[i, stesso].mean()
        b = min(D[i, lab == j].mean() for j in np.unique(lab) if j != lab[i])
        s[i] = (b - a) / max(a, b)
    return s.mean()

print("silhouette a mano (K=4): %.4f   sklearn: %.4f" % (silhouette_a_mano(X, lab), silhouette_score(X, lab)))
for k in range(2, 9):
    m = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    print("K=%d  inerzia %8.1f   silhouette %.3f" % (k, m.inertia_, silhouette_score(X, m.labels_)))

silhouette a mano (K=4): 0.7982   sklearn: 0.7982
K=2  inerzia  18575.6   silhouette 0.580
K=3  inerzia   5430.6   silhouette 0.761
K=4  inerzia   1295.9   silhouette 0.798
K=5  inerzia   1078.1   silhouette 0.705
K=6  inerzia    921.1   silhouette 0.712
K=7  inerzia    811.3   silhouette 0.619
K=8  inerzia    710.7   silhouette 0.488

La funzione silhouette_a_mano calcola la matrice completa delle distanze (radice di dist2), poi per ogni punto i calcola a come media delle distanze dai punti dello stesso gruppo (escluso sé stesso con stesso[i] = False) e b come il minimo, tra gli altri gruppi, della distanza media da quel gruppo. Il ciclo finale riaddestra K-Means per K da 2 a 8 con 10 avvii e stampa inerzia e silhouette.

La silhouette ha il massimo a K = 4 (0.798), il valore con cui abbiamo generato i dati. L'inerzia scende sempre (18575, 5430, 1296, 1078, ...); qui il gomito a K = 4 è visibile perché i gruppi sono ben separati, ma su dati reali raramente è così netto. La silhouette a mano coincide con quella di scikit-learn (0.7982). Nessuna delle due metriche è un oracolo: sono strumenti per formulare un'ipotesi, che poi controlli con la conoscenza del dominio (i gruppi hanno senso per chi conosce i dati?).

Dove K-Means sbaglia: forme non sferiche e scale diverse

K-Means assume implicitamente gruppi sferici, di dimensione simile e separati dai centri: assegna al centro più vicino, quindi i confini tra gruppi sono rette. Due lune intrecciate non hanno questa forma. Allo stesso modo, se una feature ha scala molto più grande delle altre, la distanza euclidea è dominata da quella (come in KNN).

from sklearn.datasets import make_moons

Xm, ym = make_moons(n_samples=500, noise=0.07, random_state=0)
km = KMeans(n_clusters=2, n_init=10, random_state=0).fit(Xm)
print("K-Means sulle due lune: ARI = %.3f   (1.0 = gruppi veri recuperati)" % adjusted_rand_score(ym, km.labels_))

K-Means sulle due lune: ARI = 0.245   (1.0 = gruppi veri recuperati)

Il listato genera 500 punti su due semicerchi con poco rumore (0,07) e lancia K-Means con K = 2; confronta il risultato con i gruppi veri tramite ARI.

L'ARI è 0.245, lontanissimo da 1: il K-Means taglia le lune a metà con una retta. Serve un algoritmo che ragioni sulla densità, non sulla distanza da un centro.

DBSCAN: clustering per densità

DBSCAN (Density-Based Spatial Clustering of Applications with Noise) ha due parametri: un raggio ε (eps) e un numero minimo di punti min_samples. Classifica ogni punto in uno di tre modi:

• Punto core: nel raggio ε attorno a sé (compreso sé stesso) ci sono almeno min_samples punti. È in una regione densa.
• Punto di bordo: non è core, ma sta nel raggio ε di un punto core.
• Rumore: né core né di bordo. Riceve l'etichetta −1.

Le tre categorie hanno un ordine logico: prima si decide chi è core guardando i soli conteggi, poi si assegnano i punti di bordo guardando i core, e quello che resta è rumore. Secondo la guida di scikit-learn un cluster è un insieme di core vicini fra loro più i non-core vicini a un core. Nel nostro codice min_samples conta anche il punto stesso, come fa il parametro di DBSCAN.

Un cluster è l'insieme dei punti core connessi per densità (raggiungibili l'uno dall'altro passando per punti core vicini a catena), più i punti di bordo agganciati. L'algoritmo parte da un punto core non visitato, esplora a catena i vicini dei core (una visita in ampiezza), e quando la catena si esaurisce chiude il cluster:

from sklearn.cluster import DBSCAN

def dbscan(X, eps, min_samples):
    n = len(X)
    D = np.sqrt(dist2(X, X))
    vicini = [np.where(D[i] <= eps)[0] for i in range(n)]
    core = np.array([len(v) >= min_samples for v in vicini])
    lab = np.full(n, -1)                       # -1 = rumore, finché non viene agganciato
    c = 0
    for i in range(n):
        if lab[i] != -1 or not core[i]:
            continue
        lab[i] = c
        coda = [i]
        while coda:
            p = coda.pop()
            if not core[p]:
                continue                       # i punti di bordo si agganciano ma non espandono
            for q in vicini[p]:
                if lab[q] == -1:
                    lab[q] = c
                    coda.append(q)
        c += 1
    return lab

lab_db = dbscan(Xm, eps=0.2, min_samples=5)
lib_db = DBSCAN(eps=0.2, min_samples=5).fit(Xm).labels_
print("cluster trovati (a mano): %d, punti di rumore: %d" % (len(set(lab_db)) - (1 if -1 in lab_db else 0), (lab_db == -1).sum()))
print("ARI DBSCAN a mano sulle lune:  %.3f" % adjusted_rand_score(ym, lab_db))
print("ARI DBSCAN sklearn:            %.3f" % adjusted_rand_score(ym, lib_db))
print("accordo tra le due implementazioni (ARI): %.3f" % adjusted_rand_score(lab_db, lib_db))

cluster trovati (a mano): 2, punti di rumore: 0
ARI DBSCAN a mano sulle lune:  1.000
ARI DBSCAN sklearn:            1.000
accordo tra le due implementazioni (ARI): 1.000

Nel listato vicini[i] contiene gli indici dei punti entro il raggio ε da i (sé stesso incluso); core è un vettore booleano. Il ciclo esterno scorre i punti e, per ogni punto core non ancora assegnato, apre un nuovo cluster c e usa una pila (coda) per esplorare: estrae un punto, e se è core aggiunge al cluster tutti i suoi vicini non ancora assegnati. I punti di bordo vengono etichettati ma, non essendo core, non estendono la ricerca. Quelli mai raggiunti restano a −1 (rumore). Il confronto finale con DBSCAN di scikit-learn dà ARI 1.000: stesse partizioni.

DBSCAN segue la forma delle lune perché ragiona per connessione, non per vicinanza a un centro. Non serve fissare il numero di gruppi: emerge dai dati.

La trappola di DBSCAN: scegliere ε

I risultati dipendono moltissimo da ε. Nell'esperimento qui sotto, sulle lune, il 99° percentile della k-distanza è 0.170: un raggio di 0.2 copre quasi tutti i punti. Con ε = 0.05 troviamo 30 frammenti e 253 punti di rumore, con 0.10 tre cluster, con 0.20 i due veri (ARI 1.000), con 0.40 tutto si fonde in un cluster solo (ARI 0).

In generale, ε troppo piccolo frammenta tutto o lo declassa a rumore, ε troppo grande fonde tutto in un cluster solo. Una guida pratica è il grafico delle k-distanze: per ogni punto calcola la distanza dal suo min_samples-esimo vicino, ordina i valori in modo crescente e cerca il «gomito», dove i punti passano da zone dense a isolate.

kd = np.sort(np.sort(np.sqrt(dist2(Xm, Xm)), axis=1)[:, 5])          # distanza dal 5° vicino (indice 0 = se stesso)
print("percentili della k-distanza (k=5):", {p: round(float(np.percentile(kd, p)), 3) for p in (50, 90, 95, 99)})
for eps in (0.05, 0.1, 0.2, 0.4, 1.0):
    l = dbscan(Xm, eps=eps, min_samples=5)
    print("eps=%.2f  cluster %d  rumore %3d   ARI %.3f" % (eps, len(set(l)) - (1 if -1 in l else 0), (l == -1).sum(), adjusted_rand_score(ym, l)))

percentili della k-distanza (k=5): {50: 0.071, 90: 0.11, 95: 0.128, 99: 0.17}
eps=0.05  cluster 30  rumore 253   ARI 0.020
eps=0.10  cluster 3  rumore  18   ARI 0.706
eps=0.20  cluster 2  rumore   0   ARI 1.000
eps=0.40  cluster 1  rumore   0   ARI 0.000
eps=1.00  cluster 1  rumore   0   ARI 0.000

La riga della k-distanza ordina per riga la matrice delle distanze e prende la colonna 5 (la distanza dal quinto vicino, essendo la colonna 0 il punto stesso), poi riordina i valori; si stampano alcuni percentili. Nel ciclo si prova DBSCAN con cinque valori di ε confrontando numero di cluster, punti di rumore e ARI.

Altre insidie: ε è un unico raggio per tutto il dataset, quindi DBSCAN fatica quando i gruppi hanno densità diverse; è basato su distanze, quindi richiede feature scalate; e in alta dimensione le distanze si concentrano (la maledizione della dimensionalità descritta in KNN) e il concetto di densità si indebolisce.

Trade-off tra i due

• K-Means: veloce, scala a milioni di punti, centri interpretabili (profili medi), ma richiede K, assume gruppi sferici, non gestisce il rumore (ogni punto va a un gruppo), dipende dall'inizializzazione.
• DBSCAN: forme arbitrarie, riconosce il rumore, non richiede K, ma è sensibile a ε e min_samples, non gestisce densità diverse, costa di più su grandi dataset (la ricerca dei vicini; con indici spaziali è accettabile in bassa dimensione).

I due elenchi sono speculari: K-Means chiede il numero di gruppi e ne impone la forma, DBSCAN chiede la densità e lascia libera la forma. Nessuno dei due vince in generale. Scegli K-Means se i gruppi sono compatti e ti serve un profilo medio per gruppo, DBSCAN se la forma è irregolare o ti interessa sapere quali punti non appartengono a nulla.

Trappole reali

• Prendere i cluster per verità: l'algoritmo trova gruppi anche nel rumore uniforme. Chiediti sempre se i gruppi hanno senso e sono stabili (rilancia con altri semi e su sottocampioni).
• Dimenticare lo scaling: età in anni e reddito in euro nello stesso K-Means = il reddito decide tutto.
• Feature categoriche codificate come numeri: K-Means sulle medie non ha senso su categorie. Servono varianti (K-Modes) o altre rappresentazioni.
• Usare ARI o accuratezza in produzione: servono etichette vere, che per definizione non hai. Usa silhouette, stabilità e validazione sul dominio.
• Leakage in fase di preparazione: se clusterizzi per creare una feature da dare a un modello supervisionato, il clustering va addestrato solo sul training (dentro la pipeline).

Le prime tre trappole sono errori di preparazione, le ultime due errori di valutazione. In tutte il clustering gira e restituisce gruppi, perché un algoritmo non supervisionato non ha modo di dirti che la risposta è insensata. La sezione che segue misura due di questi errori con numeri.

La scala delle feature e le feature inutili: due modi di rovinare un clustering

Nel clustering, dove non esistono etichette da controllare, gli errori di preparazione dei dati sono ancora più insidiosi che nel caso supervisionato: nessuna metrica di test ti avverte che qualcosa non va. Mostriamo due esperimenti. Nel primo, ai tre gruppi ben separati di prima si aggiunge una colonna di rumore puro con deviazione standard 50 (cioè una grandezza «grande» e senza informazione, come un identificativo numerico o un importo non normalizzato). Nel secondo, si cambia soltanto l'unità di misura di una delle due feature vere (da piedi a centimetri, un fattore 30,48), senza aggiungere nulla.

from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler

rng2 = np.random.default_rng(0)
X_true, y_true = make_blobs(n_samples=600, centers=3, cluster_std=1.0, random_state=4)
X_noisy = np.hstack([X_true, rng2.normal(0, 50, (600, 1))])                      # una colonna di rumore su scala molto grande

def evaluate(name, A):
    km = KMeans(3, n_init=10, random_state=0).fit(A)
    print("%-38s  ARI %.3f    silhouette %.3f" % (name, adjusted_rand_score(y_true, km.labels_), silhouette_score(A, km.labels_)))

evaluate("solo le due feature vere", X_true)
evaluate("più una colonna di rumore (dev. 50)", X_noisy)
evaluate("... dopo la standardizzazione", StandardScaler().fit_transform(X_noisy))

print()
evaluate("feature vere, unità originali", X_true)
evaluate("prima feature moltiplicata per 30,48", np.hstack([X_true[:, :1] * 30.48, X_true[:, 1:]]))

solo le due feature vere                ARI 0.932    silhouette 0.613
più una colonna di rumore (dev. 50)     ARI 0.001    silhouette 0.508
... dopo la standardizzazione           ARI 0.498    silhouette 0.439

feature vere, unità originali           ARI 0.932    silhouette 0.613
prima feature moltiplicata per 30,48    ARI 0.492    silhouette 0.626

La funzione evaluate esegue K-Means con K = 3 e riporta due numeri: l'ARI con i gruppi veri (che qui conosciamo perché i dati sono sintetici) e la silhouette (che nella pratica è l'unica cosa che avresti). Con le sole feature vere l'ARI è 0,932 e la silhouette 0,613. Aggiungendo la colonna di rumore l'ARI crolla a 0,001: il clustering è casuale, perché la distanza euclidea è dominata dalla colonna con deviazione 50 e K-Means separa i punti in base al rumore. Standardizzare aiuta, ma non basta: l'ARI risale solo a 0,498, perché la colonna di rumore ha ancora lo stesso peso delle due colonne informative (un terzo ciascuna) e continua a disturbare.

Nota la silhouette: 0,508 per il clustering rovinato dal rumore, un valore che sembrerebbe «buono», perché i gruppi sono ben separati lungo la dimensione di rumore. La silhouette misura la compattezza e la separazione geometrica dei gruppi trovati, non il loro significato: non ti dice se i gruppi corrispondono a qualcosa di reale.

Il secondo esperimento mostra che anche solo cambiare l'unità di misura di una feature, senza aggiungere nulla, sposta l'ARI da 0,932 a 0,492 (la silhouette, intanto, sale a 0,626): una scala diversa pesa di più la feature nella distanza e cambia i gruppi. Il clustering non è invariante rispetto alle unità, a differenza di molti modelli supervisionati. Le conseguenze pratiche sono tre: standardizza sempre (o scegli le scale con consapevolezza); seleziona le feature (togli quelle senza informazione e gli identificativi); e non fidarti di una metrica interna come la silhouette per dire che un clustering «è giusto».

Gaussian Mixture e l'algoritmo EM: K-Means con incertezza e forme

K-Means assegna ogni punto a un solo gruppo («assegnazione dura») e assume gruppi sferici di dimensione simile. Un modello a miscela di gaussiane (GMM) generalizza: i dati sono generati da K gaussiane, ciascuna con una propria media μ_k, varianza (o matrice di covarianza) σ²_k e peso π_k (la frazione dei punti che le appartengono). Ogni punto ha una probabilità di appartenere a ciascuna componente (assegnazione morbida), e i gruppi possono essere ellittici e di dimensioni diverse. Si addestra con l'algoritmo EM (Expectation-Maximization), che ha la stessa struttura di K-Means a due passi:

E (aspettativa): calcola la responsabilità di ogni componente per ogni punto
    r_ik = π_k · N(x_i | μ_k, σ²_k) / Σ_j π_j · N(x_i | μ_j, σ²_j)

M (massimizzazione): aggiorna i parametri come medie pesate dalle responsabilità
    N_k = Σ_i r_ik        μ_k = Σ_i r_ik x_i / N_k        σ²_k = Σ_i r_ik (x_i − μ_k)² / N_k        π_k = N_k / n

Il passo E è il teorema di Bayes (lo stesso di Naive Bayes) con le componenti come classi nascoste: la «probabilità a posteriori» che il punto i provenga dalla componente k. Il passo M è K-Means con medie pesate invece che «dure». Un risultato teorico garantisce che la verosimiglianza dei dati non decresce mai a ogni iterazione: lo verifichiamo a mano su dati unidimensionali generati da due gaussiane (300 punti da N(−2, 1) e 200 da N(3, 1,5)).

from sklearn.mixture import GaussianMixture

data = np.concatenate([rng2.normal(-2, 1.0, 300), rng2.normal(3, 1.5, 200)])
mu, sd, pi = np.array([-1.0, 1.0]), np.array([1.0, 1.0]), np.array([0.5, 0.5])      # partenza grossolana
normal_pdf = lambda x, m, s: np.exp(-0.5 * ((x - m) / s) ** 2) / (s * np.sqrt(2 * np.pi))

log_likelihood = []
for iteration in range(40):
    weighted = np.stack([pi[k] * normal_pdf(data, mu[k], sd[k]) for k in range(2)], axis=1)     # (n, 2): pi_k N(x_i | mu_k, sd_k)
    log_likelihood.append(np.log(weighted.sum(axis=1)).sum())
    r = weighted / weighted.sum(axis=1, keepdims=True)                              # E: responsabilità
    Nk = r.sum(axis=0)                                                              # M: nuovi parametri
    mu = (r * data[:, None]).sum(axis=0) / Nk
    sd = np.sqrt((r * (data[:, None] - mu) ** 2).sum(axis=0) / Nk)
    pi = Nk / len(data)

print("log-verosimiglianza alle prime 5 iterazioni:", np.round(log_likelihood[:5], 2), "...  finale %.2f" % log_likelihood[-1])
print("non decresce mai: %s" % np.all(np.diff(log_likelihood) > -1e-9))
print("EM a mano:      medie %s   deviazioni %s   pesi %s" % (np.round(mu, 3), np.round(sd, 3), np.round(pi, 3)))
gmm = GaussianMixture(2, random_state=0).fit(data[:, None])
print("scikit-learn:   medie %s   deviazioni %s   pesi %s" % (np.round(gmm.means_.ravel(), 3), np.round(np.sqrt(gmm.covariances_.ravel()), 3), np.round(gmm.weights_, 3)))
print("valori veri:    medie [-2, 3]   deviazioni [1, 1.5]   pesi [0.6, 0.4]")

log-verosimiglianza alle prime 5 iterazioni: [-1619.43 -1078.2  -1075.2  -1074.98 -1074.96] ...  finale -1074.95
non decresce mai: True
EM a mano:      medie [-2.064  2.75 ]   deviazioni [0.888 1.576]   pesi [0.591 0.409]
scikit-learn:   medie [-2.05   2.792]   deviazioni [0.9  1.54]   pesi [0.596 0.404]
valori veri:    medie [-2, 3]   deviazioni [1, 1.5]   pesi [0.6, 0.4]

Il listato parte da stime grossolane (medie −1 e 1, deviazioni 1, pesi 0,5) e per 40 iterazioni ripete: calcolo le probabilità pesate di ciascun punto per ciascuna componente (con la densità gaussiana scritta a mano), le sommo per ottenere la log-verosimiglianza, normalizzo per avere le responsabilità (passo E), e aggiorno medie, deviazioni e pesi come medie pesate (passo M). La log-verosimiglianza passa da −1 619,43 alla prima iterazione a −1 078,2 alla seconda e converge a −1 074,95: sale a ogni passo e non scende mai, come garantito dalla teoria (la verifica np.all(np.diff(...) > -1e-9) dà vero).

I parametri trovati (medie −2,064 e 2,750, deviazioni 0,888 e 1,576, pesi 0,591 e 0,409) sono vicini ai valori che hanno generato i dati (−2 e 3; 1 e 1,5; 0,6 e 0,4), e coincidono con quelli di scikit-learn a meno di piccole differenze di convergenza (−2,050, 2,792; 0,900, 1,540; 0,596, 0,404). Lo scarto dai valori veri è l'errore statistico di un campione di 500 punti, non un difetto dell'algoritmo.

Gruppi ellittici: dove K-Means sbaglia e la miscela no

T = np.array([[2.5, 0.0], [1.8, 0.4]])                                              # matrice che stira e inclina i gruppi
X_ell = np.vstack([rng2.normal(size=(300, 2)) @ T + shift for shift in ([0, 0], [0, 6], [5, 3])])
y_ell = np.repeat([0, 1, 2], 300)

km = KMeans(3, n_init=10, random_state=0).fit(X_ell)
gmm = GaussianMixture(3, covariance_type="full", n_init=5, random_state=0).fit(X_ell)
print("ARI con i gruppi veri:   K-Means %.3f    miscela gaussiana %.3f" % (adjusted_rand_score(y_ell, km.labels_), adjusted_rand_score(y_ell, gmm.predict(X_ell))))

print("\nscelta del numero di componenti con il criterio BIC (più basso = meglio):")
for k in range(1, 7):
    print("   K=%d   BIC %.1f" % (k, GaussianMixture(k, n_init=3, random_state=0).fit(X_ell).bic(X_ell)))

ARI con i gruppi veri:   K-Means 0.555    miscela gaussiana 1.000

scelta del numero di componenti con il criterio BIC (più basso = meglio):
   K=1   BIC 9165.7
   K=2   BIC 8041.7
   K=3   BIC 7196.1
   K=4   BIC 7238.3
   K=5   BIC 7266.5
   K=6   BIC 7305.1

Il listato genera tre gruppi gaussiani allungati e inclinati (applicando una matrice T a punti sferici). K-Means, che assume gruppi sferici, taglia i gruppi con confini rettilinei e ottiene un ARI di 0,555; la miscela con covarianze piene, che modella la forma ellittica di ciascuna componente, li recupera perfettamente (ARI 1,000).

La seconda parte mostra come scegliere K con il criterio BIC (Bayesian Information Criterion): la log-verosimiglianza penalizzata per il numero di parametri, che a differenza della verosimiglianza stessa non cresce sempre con K. Il BIC scende da 9 165,7 (un solo gruppo) a 8 041,7 (due) e raggiunge il minimo, 7 196,1, a K = 3, il valore vero; per K = 4, 5, 6 risale (7 238,3, 7 266,5, 7 305,1). A differenza dell'inerzia di K-Means, che decresce sempre, il BIC ha un minimo vero, ed è per questo un criterio di selezione più principiato (sotto l'ipotesi, da verificare, che i gruppi siano davvero gaussiani).

Stabilità e clustering gerarchico

Un clustering «giusto» dovrebbe essere stabile: se ricampioni i dati e rifai l'analisi, dovresti ottenere gli stessi gruppi. Si misura con il bootstrap: si estraggono campioni con reinserimento, si addestra K-Means su ciascuno, si assegnano i punti originali al clustering così ottenuto e si calcola l'ARI tra questo e il clustering sull'intero dataset. Valori vicini a 1 indicano gruppi robusti. Poi vediamo un'alternativa che non richiede di fissare K: il clustering gerarchico agglomerativo, che parte da ogni punto come gruppo a sé e fonde ripetutamente i due gruppi più vicini, secondo un criterio di collegamento (linkage): single (distanza tra i due punti più vicini dei gruppi), complete (distanza tra i più lontani), average (media), ward (minimizza l'aumento di varianza interna).

from sklearn.cluster import AgglomerativeClustering

X_four, _ = make_blobs(n_samples=600, centers=4, cluster_std=0.9, random_state=3)

def stability(X, k, repetitions=20):
    reference = KMeans(k, n_init=5, random_state=0).fit(X)
    aris = []
    for b in range(repetitions):
        sample = np.random.default_rng(b).integers(0, len(X), len(X))                # campione bootstrap
        km = KMeans(k, n_init=5, random_state=b).fit(X[sample])
        aris.append(adjusted_rand_score(reference.labels_, km.predict(X)))            # confronto sull'intero dataset
    return np.mean(aris), np.std(aris)

print("K    stabilità (ARI medio tra clustering ripetuti)")
for k in (2, 3, 4, 5, 6):
    mean_ari, std_ari = stability(X_four, k)
    print("%-4d  %.3f +- %.3f" % (k, mean_ari, std_ari))

X_moons, y_moons = make_moons(500, noise=0.07, random_state=0)
print("\nclustering gerarchico sulle due lune, ARI con i gruppi veri:")
for linkage in ("ward", "complete", "average", "single"):
    print("   %-9s %.3f" % (linkage, adjusted_rand_score(y_moons, AgglomerativeClustering(2, linkage=linkage).fit_predict(X_moons))))

K    stabilità (ARI medio tra clustering ripetuti)
2     1.000 +- 0.000
3     0.899 +- 0.171
4     1.000 +- 0.000
5     0.834 +- 0.047
6     0.752 +- 0.080

clustering gerarchico sulle due lune, ARI con i gruppi veri:
   ward      0.445
   complete  0.523
   average   0.429
   single    1.000

La funzione stability confronta, con l'ARI, il clustering di riferimento con quelli ottenuti da 20 ricampionamenti. Su quattro gruppi veri, la stabilità è perfetta (1,000) con K = 2 e K = 4, e più bassa con K = 3 (0,899 ± 0,171: a volte due gruppi veri vengono fusi in un modo, a volte nell'altro), K = 5 (0,834) e K = 6 (0,752). Due insegnamenti. Primo, la stabilità segnala i K «naturali» (qui 4), ma anche K = 2 risulta perfettamente stabile, perché due grandi coppie di gruppi sono una struttura robusta: la stabilità non individua un unico K, individua le scale a cui esiste struttura. Secondo, i K troppo grandi sono instabili perché K-Means spezza arbitrariamente gruppi compatti, e questo è un buon segnale di allarme di over-clustering.

Il secondo esperimento è un promemoria sull'importanza del criterio di collegamento: sulle due lune, single linkage ottiene un ARI perfetto (1,000), perché segue le catene di punti vicini (come DBSCAN), mentre Ward (0,445), complete (0,523) e average (0,429) cercano gruppi compatti e tagliano le lune a metà. Anche qui la scelta dell'algoritmo codifica un'ipotesi su che cosa sia un gruppo.

Tre prove da fare

Esercizio 1. Nell'esperimento del rumore, calcola l'ARI dopo aver tolto di proposito la colonna di rumore con un criterio automatico (per esempio misurando, per ogni feature, la silhouette del clustering ottenuto con quella sola feature). Esercizio 2. Con la miscela gaussiana sui gruppi ellittici, confronta covariance_type uguale a «spherical», «diag» e «full»: quale ricostruisce meglio i gruppi, e quale ha il BIC più basso? Esercizio 3. Calcola la stabilità di DBSCAN al variare di ε sulle due lune: quali valori di ε sono stabili?

Per allenarti sul dominio dei profili d'uso delle stazioni di bike sharing c'è uno scheletro di esercizi in Python con tre parti da completare: l'algoritmo di Lloyd scritto da zero con la verifica che la SSE non cresce, l'inizializzazione che finisce in un minimo locale 38 volte peggiore, e DBSCAN etichettato a mano. Finché i punti da completare sono vuoti, stampa quali esercizi mancano.

Che cosa sai fare adesso

Sai scrivere K-Means e DBSCAN da zero e spiegare perché K-Means converge ma a un minimo locale, sai scegliere K con inerzia, silhouette, BIC e stabilità sapendo che nessuno è un oracolo, e sai riconoscere quando la forma dei gruppi o la scala delle feature rende il risultato privo di senso. Sai usare una miscela gaussiana quando i gruppi sono ellittici. Nella lezione successiva, riduzione di dimensionalità, riduciamo il numero di dimensioni conservando l'informazione che conta.

Fonti

Lezione precedente del percorso: Boosting: AdaBoost, Gradient Boosting e cenni su XGBoost. Nella stessa serie rimandano a questa lezione anche: Alberi decisionali: come si sceglie una domanda con Gini ed entropia.