""" Corso Machine Learning Engineer - Modulo 03, Lezione 8 "Clustering: K-Means e DBSCAN senza etichette" OBIETTIVO DI QUESTO FILE ------------------------ Tre esercizi su un dominio nuovo rispetto all'articolo: raggruppare le stazioni di un servizio di bike sharing in base al loro profilo d'uso. Ogni stazione è descritta da due numeri: prelievi nella fascia 7-10 e prelievi nella fascia 17-20, in centinaia. 1. l'assegnazione ai centroidi, scritta da zero; 2. l'aggiornamento dei centroidi, e la dimostrazione numerica che la SSE scende a ogni mezzo passo - più un'inizializzazione che la fa finire in un MINIMO LOCALE 38 volte peggiore del minimo globale. È il motivo per cui `n_init` esiste; 3. DBSCAN a mano: etichettare nucleo / bordo / rumore, e vedere che cambiare min_pts da 3 a 2 trasforma due punti di rumore in un gruppo. COME SI USA ----------- python3 esercizi.py Il file gira anche così com'è: stampa quali esercizi mancano ancora. AMBIENTE DICHIARATO ------------------- Python 3.12+ e NumPy 2.3. Nessun altro pacchetto. python3 -m venv .venv && source .venv/bin/activate pip install "numpy>=2.0" Tutti e tre gli esercizi sono DETERMINISTICI: nessun seed, nessuna inizializzazione casuale. I centroidi di partenza sono scelti a mano apposta, perché il punto è proprio vedere che cosa cambiano. """ import numpy as np # ---------------------------------------------------------------------- # Nove stazioni, tre quartieri. # # Quartiere centro : (0,0) (1,0) (0,1) # Quartiere università : (10,0) (11,0) (10,1) # Quartiere stazione FS : (20,0) (21,0) (20,1) # # Tre gruppi evidenti a occhio. Il punto dell'esercizio è che K-Means può # non trovarli, e non lo dice. # ---------------------------------------------------------------------- STAZIONI = np.array([ [0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [10.0, 0.0], [11.0, 0.0], [10.0, 1.0], [20.0, 0.0], [21.0, 0.0], [20.0, 1.0], ]) INIT_FORTUNATA = np.array([[0.0, 0.0], [10.0, 0.0], [20.0, 0.0]]) INIT_SFORTUNATA = np.array([[0.0, 0.0], [1.0, 0.0], [15.0, 0.0]]) # ---------------------------------------------------------------------- # ESERCIZIO 1 - Assegnare ogni stazione al centroide più vicino # # Prima metà del passo di Lloyd. Per ogni punto, trova l'indice del centroide # a distanza euclidea minima. # # Domanda: perché non serve la radice quadrata? # ---------------------------------------------------------------------- def assegna(X, centroidi): """Indice del centroide più vicino, per ogni riga di X. Parametri --------- X : np.ndarray (n, d) centroidi : np.ndarray (k, d) Ritorna ------- np.ndarray (n,) di interi tra 0 e k-1, oppure None se non svolto. Suggerimento senza cicli: `X[:, None, :] - centroidi[None, :, :]` ha forma (n, k, d); sommando i quadrati su `axis=2` ottieni le distanze al quadrato (n, k), e `argmin(axis=1)` fa il resto. La radice quadrata NON serve: è una funzione crescente, quindi non cambia quale sia il minimo. Toglierla risparmia n*k radici a ogni iterazione, e in un ciclo che gira cento volte su un milione di punti non è poco. """ # TODO (esercizio 1): sostituisci `return None` con l'assegnazione. return None # ---------------------------------------------------------------------- # ESERCIZIO 2 - Aggiornare i centroidi, e guardare la SSE # # Seconda metà del passo di Lloyd: ogni centroide diventa la media dei punti # che gli sono stati assegnati. # # La media NON è una scelta estetica: è il punto che minimizza la somma dei # quadrati delle distanze. Derivando rispetto a c la funzione # somma_i ||x_i - c||^2 si ottiene -2 * somma_i (x_i - c) = 0, cioè # c = media degli x_i. È per questo che K-Means e la SSE vanno insieme, e che # con una metrica diversa (Manhattan) il centro giusto sarebbe la mediana. # # Domanda: che cosa fai con un cluster rimasto senza punti? # ---------------------------------------------------------------------- def aggiorna(X, etichette, k, centroidi_precedenti): """Nuovi centroidi come media dei punti assegnati. Parametri --------- X : np.ndarray (n, d) etichette : np.ndarray (n,) di interi k : int centroidi_precedenti : np.ndarray (k, d) - da riusare per i cluster vuoti. Ritorna ------- np.ndarray (k, d), oppure None se l'esercizio non è ancora svolto. Il cluster VUOTO è il caso che fa esplodere le implementazioni ingenue: la media di zero punti è `nan`, e da quel momento ogni distanza diventa `nan` e `argmin` restituisce 0 per tutti. Qui la regola è semplice e dichiarata: un cluster vuoto tiene il centroide precedente. Le librerie serie fanno di meglio (riposizionano il centroide sul punto più lontano), ma la cosa importante è non lasciare un `nan` in giro. """ # TODO (esercizio 2): sostituisci `return None` con l'aggiornamento. return None def sse(X, etichette, centroidi): """Somma dei quadrati degli scarti dai rispettivi centroidi. Già pronta.""" scarti = X - centroidi[etichette] return float((scarti ** 2).sum()) def lloyd(X, centroidi_iniziali, max_iterazioni=20): """Il ciclo completo. Già pronto: usa le tue due funzioni. Ritorna (etichette, centroidi, storia_sse) oppure None se mancano i TODO. """ centroidi = centroidi_iniziali.copy() k = len(centroidi) storia = [] etichette = None for _ in range(max_iterazioni): nuove = assegna(X, centroidi) if nuove is None: return None storia.append(sse(X, nuove, centroidi)) # dopo l'assegnazione nuovi_centroidi = aggiorna(X, nuove, k, centroidi) if nuovi_centroidi is None: return None storia.append(sse(X, nuove, nuovi_centroidi)) # dopo l'aggiornamento fermo = etichette is not None and np.array_equal(nuove, etichette) etichette, centroidi = nuove, nuovi_centroidi if fermo: break return etichette, centroidi, storia # ---------------------------------------------------------------------- # ESERCIZIO 3 - DBSCAN a mano # # Sette fermate su una pista ciclabile rettilinea, alle distanze in chilometri # # 0 , 1 , 2 , 3 , 8 , 9 , 20 # # Con eps = 1.5 e min_pts = 3 (il punto conta se stesso nel vicinato), # etichetta ciascun punto come: # # NUCLEO : ha almeno min_pts punti entro eps (se stesso compreso) # BORDO : non è nucleo, ma ha almeno un NUCLEO entro eps # RUMORE : nessuna delle due # # Domanda: che cosa succede alle fermate 8 e 9 se min_pts scende a 2? # ---------------------------------------------------------------------- FERMATE = np.array([0.0, 1.0, 2.0, 3.0, 8.0, 9.0, 20.0]) def classifica_dbscan(punti, eps, min_pts): """Etichetta ogni punto come "nucleo", "bordo" o "rumore". Parametri --------- punti : np.ndarray (n,) - coordinate su una retta. eps : float - raggio del vicinato, estremi COMPRESI (<=, non <). min_pts : int - quanti punti servono nel vicinato, se stesso incluso. Ritorna ------- lista di n stringhe, oppure None se l'esercizio non è ancora svolto. Procedi in due passate, e l'ordine conta: prima decidi chi è nucleo GUARDANDO SOLO I CONTEGGI, poi decidi chi è bordo guardando i nuclei già trovati. Se mescoli le due passate rischi di dichiarare bordo un punto vicino a un altro bordo, e DBSCAN non lo permette. """ # TODO (esercizio 3): sostituisci `return None` con le due passate. return None # ---------------------------------------------------------------------- # Runner # ---------------------------------------------------------------------- def _sezione(titolo): print() print(titolo) print("-" * len(titolo)) def _da_fare(numero): print(f" [ ] Esercizio {numero} non ancora svolto: completa il TODO qui sopra.") def esegui_esercizi_1_e_2(): _sezione("Esercizi 1 e 2 - K-Means e il minimo locale") if assegna(STAZIONI, INIT_FORTUNATA) is None: _da_fare(1) return 0 if aggiorna(STAZIONI, np.zeros(len(STAZIONI), dtype=int), 3, INIT_FORTUNATA) is None: _da_fare(2) return 1 for nome, init in (("inizializzazione fortunata ", INIT_FORTUNATA), ("inizializzazione sfortunata", INIT_SFORTUNATA)): etichette, centroidi, storia = lloyd(STAZIONI, init) print(f" {nome}: centroidi di partenza {init[:, 0].tolist()}") print(f" SSE lungo il ciclo: {[round(s, 4) for s in storia]}") print(f" centroidi finali:" f" {[tuple(np.round(c, 4)) for c in centroidi]}") print(f" etichette finali: {etichette.tolist()}") print(f" SSE finale: {sse(STAZIONI, etichette, centroidi):.4f}") print(" Stessi dati, stesso algoritmo, stesso numero di cluster: due") print(" risultati diversi. K-Means converge sempre, ma a un minimo") print(" LOCALE. L'unica difesa è ripartire da inizializzazioni diverse") print(" e tenere la SSE più bassa: è esattamente cio' che fa n_init,") print(" ed è il motivo per cui n_init=1 è quasi sempre un errore.") return 2 def esegui_esercizio_3(): _sezione("Esercizio 3 - DBSCAN a mano") for min_pts in (3, 2): etichette = classifica_dbscan(FERMATE, eps=1.5, min_pts=min_pts) if etichette is None: _da_fare(3) return False print(f" eps = 1.5 , min_pts = {min_pts}") for x, e in zip(FERMATE, etichette): print(f" fermata al km {x:5.1f} -> {e}") print(" Due punti (8 e 9) passano da RUMORE a NUCLEO solo perché") print(" min_pts è sceso di uno. DBSCAN non ha bisogno di sapere quanti") print(" gruppi ci sono, ma ha bisogno di eps e min_pts, e quei due") print(" numeri decidono tutto. Non sono meno arbitrari di K: sono solo") print(" arbitrari in un'altra unità di misura.") return True def main(): print("Corso MLE - Modulo 03, Lezione 8 - esercizi di codice") fatti = esegui_esercizi_1_e_2() if esegui_esercizio_3(): fatti += 1 print() print(f"Esercizi completati: {fatti} su 3.") if fatti < 3: print("Completa i TODO e rilancia: python3 esercizi.py") if __name__ == "__main__": main() # ====================================================================== # SOLUZIONI DI RIFERIMENTO # ====================================================================== # # --- Esercizio 1 ------------------------------------------------------- # def assegna(X, centroidi): # differenze = X[:, None, :] - centroidi[None, :, :] # distanze2 = (differenze ** 2).sum(axis=2) # return distanze2.argmin(axis=1) # # --- Esercizio 2 ------------------------------------------------------- # def aggiorna(X, etichette, k, centroidi_precedenti): # nuovi = centroidi_precedenti.copy().astype(float) # for j in range(k): # membri = X[etichette == j] # if len(membri) > 0: # nuovi[j] = membri.mean(axis=0) # return nuovi # # Atteso con l'INIZIALIZZAZIONE FORTUNATA, centroidi (0,0), (10,0), (20,0): # l'assegnazione è subito quella giusta, tre gruppi da tre. # I centroidi diventano (1/3, 1/3), (31/3, 1/3), (61/3, 1/3). # SSE di ciascun gruppo, a mano: per il gruppo (0,0),(1,0),(0,1) con # centroide (1/3,1/3) -> 2/9 + 5/9 + 5/9 = 12/9 = 1.3333. # Gli altri due gruppi sono traslazioni dello stesso triangolo, quindi # hanno la stessa SSE. # SSE FINALE = 3 * 4/3 = 4.0000 # # Atteso con l'INIZIALIZZAZIONE SFORTUNATA, centroidi (0,0), (1,0), (15,0): # due centroidi si dividono il quartiere centro, il terzo copre università # E stazione FS insieme. Al secondo giro nulla si sposta più: è un punto # fisso. # Centroidi finali: (0, 0.5), (1, 0), (46/3, 1/3) = (15.3333, 0.3333). # SSE a mano: # cluster {(0,0),(0,1)} intorno a (0,0.5) -> 0.25 + 0.25 = 0.5 # cluster {(1,0)} intorno a (1,0) -> 0 # cluster dei sei punti intorno a (15.3333, 0.3333): # (10,0) -> 28.4444 + 0.1111 = 28.5556 # (11,0) -> 18.7778 + 0.1111 = 18.8889 # (10,1) -> 28.4444 + 0.4444 = 28.8889 # (20,0) -> 21.7778 + 0.1111 = 21.8889 # (21,0) -> 32.1111 + 0.1111 = 32.2222 # (20,1) -> 21.7778 + 0.4444 = 22.2222 # somma = 152.6667 # SSE FINALE = 153.1667, cioè 38 volte il minimo globale. # # Due cose da portare via. (1) La SSE scende a ogni mezzo passo, sempre: # l'assegnazione non può aumentarla (ogni punto va al centroide più vicino) # e nemmeno l'aggiornamento (la media minimizza la somma dei quadrati). Da # qui la convergenza in un numero finito di passi: le assegnazioni possibili # sono finite e la SSE non risale mai. (2) Convergere non significa trovare # la soluzione giusta: la SSE finale di 153.17 è un minimo locale perfetto. # Nessun criterio interno avvisa. L'unico segnale è confrontare più avvii, # ed è quello che fa `n_init`. # # --- Esercizio 3 ------------------------------------------------------- # def classifica_dbscan(punti, eps, min_pts): # distanze = np.abs(punti[:, None] - punti[None, :]) # vicini = distanze <= eps # nucleo = vicini.sum(axis=1) >= min_pts # etichette = [] # for i in range(len(punti)): # if nucleo[i]: # etichette.append("nucleo") # elif np.any(vicini[i] & nucleo): # etichette.append("bordo") # else: # etichette.append("rumore") # return etichette # # Atteso con eps = 1.5 e min_pts = 3 (vicinati, estremi compresi): # km 0 : vicini {0, 1} -> 2 punti -> non nucleo; ha il nucleo 1 # entro 1.5 -> BORDO # km 1 : vicini {0, 1, 2} -> 3 punti -> NUCLEO # km 2 : vicini {1, 2, 3} -> 3 punti -> NUCLEO # km 3 : vicini {2, 3} -> 2 punti -> non nucleo; ha il nucleo 2 # entro 1.5 -> BORDO # km 8 : vicini {8, 9} -> 2 punti -> non nucleo; il suo unico # vicino (9) non è nucleo -> RUMORE # km 9 : vicini {8, 9} -> idem -> RUMORE # km 20 : vicini {20} -> 1 punto -> RUMORE # # Atteso con eps = 1.5 e min_pts = 2: # km 0, 1, 2, 3 : tutti NUCLEO (ognuno ha almeno due punti nel vicinato) # km 8, 9 : NUCLEO - ed ecco il gruppo che prima non esisteva # km 20 : RUMORE (resta solo con se stesso) # # Un solo punto di differenza in min_pts e il numero di gruppi passa da uno a # due. La morale è la stessa di K-Means con K, solo travestita: DBSCAN non # chiede quanti gruppi vuoi, ma chiede che densità chiami "gruppo". Il modo # onesto di scegliere eps è il grafico delle distanze al k-esimo vicino # ordinate, cercando il gomito; il modo onesto di scegliere min_pts è partire # da 2*d (con d = numero di feature) e giustificare gli scostamenti. # Nota sull'ordine delle due passate: il punto al km 3 è bordo perché il suo # vicino al km 2 è nucleo. Se avessi deciso "bordo" guardando i vicini già # etichettati invece dei nuclei, il km 0 sarebbe diventato bordo di un bordo, # cosa che DBSCAN non ammette: la raggiungibilità parte sempre da un nucleo.