Statistical Learning Notes PDF

Statistical Learning Enrico Daniele September 2024 - January 2025 Indice 1 La Classificazione 5 1.1 Concetti Introduttivi.....

Statistical Learning Enrico Daniele September 2024 - January 2025 Indice 1 La Classificazione 5 1.1 Concetti Introduttivi........................ 5 1.2 Il problema di Classificazione.................... 5 1.2.1 La variabile casuale Bernoulliana............. 7 Pull out what’s known.................... 7 Légge delle aspettative iterate................ 7 1.3 Classificatore ottimale di Bayes................... 7 1.3.1 Problema 2-class...................... 8 Dimostrazione........................ 8 Dimostrazione........................ 9 1.3.2 Problema multiclass.................... 9 1.4 Modello Bernoulliano-logistico Lineare............... 10 Assunzione 1......................... 11 Assunzione 2......................... 11 Odds............................. 11 Funzione logistca (o anti-logit/sigmoid)........... 11 Funzione soft-max...................... 11 Modelli con costante..................... 12 Assunzione 3......................... 12 1.5 Stima di β̂.............................. 12 1.6 Proprietà di β̂ ed Inferenza..................... 13 1.7 La linearità della Decision Boundary............... 14 1.8 Error rate in-sample......................... 15 1.9 Regressione logistica per i problemi multi-class.......... 15 1  1.9.1 Regressione Logistica One VS All............. 16 1.9.2 Modello Bernoulliano Generalizzato-Logistico (o Multi- nomial Logit)........................ 17 2 Geometria e probabilità nello spazio Multivariato 19 2.1 I concetti di Distanza e di Intorno.................. 19 2.2 Variabili Casuali Multivariate................... 21 2.2.1 Cosa sono le v. c. Multivariate?.............. 21 2.2.2 Momento primo, momento secondo e matrice VCV.... 21 Modello sferico: incorrelazione + omoschedasticità.... 23 Modello diagonale: incorrelazione + eteroschedasticità.. 23 Media campionaria (osservata)............... 23 Matrice di varianze-covarianze non distorta (osservata).. 24 2.3 Dispersione vs direzione....................... 24 2.4 Trasformazioni lineari di variabili casuali inRp......... 27 2.5 Variabile Casuale Normale Multivariata............. 28 2.5.1 Variabile Casuale Normale Multivariata Standard... 28 2.5.2 I density countours..................... 30 2.5.3 Distanza di Mahalanobis................. 30 3 Classificazione: LDA, QDA e kNN 33 3.1 Linear Discriminant Analysis - LDA................ 35 A. LDA:.......................... 35 Dimostrazione........................ 36 3.1.1 Linearità della Decision Boundary............ 36 3.1.2 Stima dei parametri.................... 36 3.2 Quadratic Discriminant Analysis - QDA............. 38 A. QDA:.......................... 38 Dimostrazione........................ 39 3.3 k-Nearest Neigbours (kNN)..................... 40 3.3.1 Classificatore kNN..................... 40 4 Metriche di Performance e Metodi di Thresholding 43 4.1 Alcune misure di performance................... 43 4.2 La curva ROC............................ 44 2 5 Validazione e Selezione del Modello 47 Descrizione locale vs globale................. 48 Parametri e iper-parametri................. 48 Le sorgenti di variazione................... 48 5.1 La questione dell’errore........................ 49 5.2 La Cross Validation......................... 51 Monte Carlo CV - MCCV.................. 53 k-folds CV - kFCV...................... 54 Trade-off nella scelta di k.................. 55 5.3 External Validation......................... 58 Caso 1 - Con EVS...................... 58 Caso 2 - Senza EVS..................... 58 6 Riduzione e Compressione dei dati 61 6.1 Principal Component Analysis - PCA................ 61 6.1.1 Costruzione delle PC.................... 62 6.1.2 Geometria delle PC.................... 65 6.1.3 Proprietà delle PC..................... 69 6.2 PCA................................. 70 7 Il Clustering 75 7.1 K-Means clustering......................... 75 7.1.1 L’iterazione di Lloyd.................... 76 7.1.2 La nozione di dissimilarità................ 78 7.1.3 K-Medoid......................... 80 7.2 Metodi gerarchici.......................... 82 I grafi............................ 82 Single Linkage (Nearest-Neighbour Linkage)....... 84 Complete Linkage...................... 84 Average Linkage....................... 85 7.3 Validazione e scelta di K...................... 86 3  1 La Classificazione 1.1 Concetti Introduttivi Lo Statistical Learning rappresenta la base statistica del Machine Learning. Il ML è una branca dell’Intelligenza Artificiale che raccoglie metodi sviluppati in diverse discipline (come la matematica e l’informatica) per costruire dai dati osservati dei meccanismi di azione, come previsioni e verifiche di ipotesi. Quando lavoriamo con i dati, possiamo distinguere due paradigmi principali: Par. predittivo, dove l’obiettivo è prevedere i casi non osservati sfruttando le relazioni tra le variabili nei dati osservati (tipico del ML); Par. esplicativo, dove l’obiettivo è comprendere le relazioni tra le variabili. Questi approcci si basano su due diverse fonti di conoscenza: Domain-independent Knowledge: strumenti pratici come programmazione, ottimizzazione, EDA (Exploratory Data Analysis), inferenza, ecc.; Domain-dependent Knowledge: conoscenza specifica e visione generale sul problema affrontato. 1.2 Il problema di Classificazione Nel problema di classificazione partiamo da una situazione dove abbiamo un certo numero di oggetti appartenenti a una classe nell’insieme Y. Nei problemi 2-class, queste sono convenzionalmente chiamate 0 e 1, o positivi (+1) e negativi (−1): sarà positiva la classe che è funzione del problema, come i pazienti che presentano una data malattia. Il problema di classificazione consiste nel prevedere, in base alle caratteristiche di un oggetto, la sua classe. Ergo, la quantità da prevedere è una variabile categoriale. In generale, osserviamo un vettore di variabili casuali, X = (X1 ,... , Xp )′ , e vogliamo costruire un classificatore, ossia una funzione Ŷ (·) che prende in input le features osservate e assegna l’oggetto a una certa classe Ŷ (X) ∈ Y. Noi vorremmo che un classificatore assegni le classi correttamente. Di solito si ha una funzione obiettivo da ottimizzare, come una Loss da minimizzare o un’Utilità da massimizzare. Qui vogliamo minimizzare una funzione di perdita. 5  Supponiamo di voler classificare gli individui in maschi e femmine, e di aver osservato una serie di features. Supponiamo di costruire il seguente classificatore: { M → se il nome finisce in ”o” Ŷ (X) := F → altrimenti Date le istanze {Ugo, Maria, Jasmine, Paolo, Gioele}, il classificatore ne classificherà correttamente 4 su 5. Per misurarne la qualità possiamo considerare la Loss 0/1: { 1 → se yi 6= Ŷ (xi ) L(yi , Ŷ (xi )) = I{yi 6= Ŷ (xi )} = 0 → altrimenti La Loss misura la perdita conseguita ogni volta che prendiamo l’istanza i e la as- segniamo a Ŷ (xi ). In altre parole, in base alle caratteristiche di i il classificatore fa una previsione, e la confronteremo con la vera etichetta. A che serve tutto ciò? La classificazione è un problema supervised e questa è solo la prima fase: useremo il classificatore addestrato nel train set per classificare le istanze fuori campione. La qualità di un classificatore dipende dalla distribuzione delle Loss, che a loro volta dipendono dall’incertezza che governa la Y e le X. Possiamo misurare la sua capacità previsiva con l’Expected Prediction Error, ossia il valore atteso della Loss. EP E := E[L(Y, Ŷ (X))] (1) Nel caso della Loss 0/1, L(·) è una v. c. dicotomica; quindi: EP E := E[I{Y 6= Ŷ (X)}] = P (Y 6= Ŷ (X)) (2) Nell’esempio, la Loss media attesa risulta: 1X 1X n n 0+0+0+0+1 L(yi , Ŷ (xi )) = I{yi 6= Ŷ (xi )} = = 0.2 (3) n i=1 n i=1 5 La Loss media, in questo caso, prende il nome di Misclassification Rate, ed è la proporzione di unità misclassificate. Sotto condizioni generali, possiamo determinare una ricetta per costruire un classificatore decente? Analizziamo l’ambiente tipo. Abbiamo: Una Popolazione di v. c., dove le Yi ∈ Y rappresentano le classi della popo- lazione con la loro distribuzione di probabilità; Un vettore di v. c. X (le features) per ogni istanza; Un classificatore Ŷ (X) che è funzione delle features e, quindi, è una v. c.; La Loss, che è la v. c. L(Y, Ŷ (X)). Poiché la Loss assume solo due valori, è una v. c. di Bernoulli. 6 1.2.1 La variabile casuale Bernoulliana La v. c. Bernoulliana Z può assumere solo due valori, fissati per convenzione nell’insieme {0, 1}, con 1 = successo. Se η = P (Z = 1) ⇒ P (Z = 0) = 1 − η. Da ciò si ha che E[Z] = P (Z = 1) = η. La successione ∑n {Z1 ,... , Zn } è un campione casuale con media campionaria Z = n1 i=1 Zi , e corrisponderà alla proporzione di successi. Esso è uno stimatore consistente per η, perché, per la légge dei grandi numeri, lim Z ≈ E[Z] = η. n→∞ Consideriamo una v. c. Reale X e fissiamo l’evento X ≥ 0. Inoltre, consid- eriamo P (X ≥ 0) = η. Si noti che, sotto queste condizioni, imponendole la funzione indicatrice, la si dicotomizza. Infatti, Z = I{X ≥ 0} ∼ Be(η) con: { 0 → se X < 0 Z := 1 → se X ≥ 0 Da ciò si ottiene che E[Z] = E[I{X ≥ 0}] = P (X ≥ 0) = η. Pull out what’s known. Date le v. c. X e Y , e una funzione g(·), allora: E[g(X)Y |X] = g(X)E[Y |X] (4) Légge delle aspettative iterate. Date le v. c. X e Y , allora: E[E[Y |X]] = E[Y ] (5) 1.3 Classificatore ottimale di Bayes C’è un classificatore ottimale che garantisca bassi livelli di rischio? Noi usiamo l’informazione contenuta in X per prevedere Ŷ (X) (facciamo un calcolo con- dizionato). Se formuliamo l’assegnazione sulla base dell’informazione fornita da X = x, il rischio atteso (o errore condizionato) è: E[L(Y, Ŷ (x))| X = x] (6) Si definisce Classificatore Ottimale di Bayes un classificatore Y ∗ (·) per cui la Loss condizionata in un dato punto x è la minore possibile; ossia: E[L(Y, Y ∗ (x))| X = x] ≤ E[L(Y, Ŷ (x))| X = x] ∀ x (7) 7 1.3.1 Problema 2-class Consideriamo il problema seguente: Y ∈ Y , con Y = {0, 1}; P (Y = 1|X) = η(X) e P (Y = 0|X) = 1 − η(X) note (posterior class probabilities - probabilità di appartenere a una classe date le features); Si dimostra che, conoscendo i posterior, si può costruire il Classificatore Bayesiano: { 0 → se η(X) < 12 Y ∗ (X) = 1 → se η(X) ≥ 1 2 Dimostrazione. Siccome la Loss assume due valori precisi, possiamo spezzarla in due componenti che rappresentano i due pezzi della decisione: L(Y, Ŷ ) = I{Y 6= Ŷ (X)} = I{Y = 1}I{Ŷ (X) = 0} + I{Y = 0}I{Ŷ (X) = 1} (8) Il classificatore è ottimale se minimizza l’errore condizionato per ogni valore di X: E[L(Y, Ŷ )|X] = E[I{Y = 1}I{Ŷ (X) = 0}|X] + E[I{Y = 0}I{Ŷ (X) = 1}|X] = = E[I{Y = 1}|X]I{Ŷ (X) = 0}] + E[I{Y = 0}|X]I{Ŷ (X) = 1}] = = P (Y = 1|X)I{Ŷ (X) = 0}] + P (Y = 0|X)I{Ŷ (X) = 1} = = η(X)I{Ŷ (X) = 0}] + [1 − η(X)]I{Ŷ (X) = 1} Si dimostra che il classificatore deve assegnare 0 se η(X) < 1 − η(X), ossia η(X) < 21 , altrimenti assegnerà 1. Ogni classificatore ha una superficie di separazione che divide geometricamente i positivi dai negativi. Tale linea è la decision boundary (frontiera di decisione), e, per p features, ne avrò p − 1. I punti con un colore diverso dalla loro regione, e, quindi, con una classe effettiva diversa dalla prevista, sono i punti misclassificati. 8  Perché questo è il miglior MCR possibile? Per la légge dei valori attesi iterati. Per costruzione, Y ∗ minimizza E[L(Y, Y ∗ (x))|X = x] ∀ x. Poiché ottimizza punto per punto localmente, ottimizzerà anche in generale. Come passiamo dal rischio condizionato locale a quello incondizionato globale (EPE)? Dobbiamo guardare alla nozione di error rate: E[L(Y, Ŷ (·))]. Poiché Y ∗ fornisce la migliore previsione per ogni punto x, allora produce globalmente l’error rate ottimale: E[L(Y, Y ∗ (X))], detto Bayes Error Rate. Si dimostra che il BER rispetto alla Loss 0/1 nel problema 2-class è dato da: E[min{η(X), 1 − η(X)}]. Dimostrazione. Per ottenere il BER dobbiamo ricavare E[L(Y, Y ∗ )]. Usando la légge del valore atteso iterato: E[L(·)] = E[E[L(·)|X]] (9) Calcoliamo innanzitutto E[L(Y, Y ∗ )|X]. Prima abbiamo ricavato: E[L(Y, Ŷ )|X] = η(X)I{Ŷ (X) = 0}] + [1 − η(X)]I{Ŷ (X) = 1} (10) Da cui otteniamo: E[L(Y, Y ∗ )|X] = η(X)I{Y ∗ (X) = 0}] + [1 − η(X)]I{Y ∗ (X) = 1} (11) Quando η(X) < 12 , Y ∗ = 0, quindi E[L(Y, Y ∗ )|X] = η(X); Quando η(X) ≥ 12 , Y ∗ = 1, quindi E[L(Y, Y ∗ )|X] = 1 − η(X). Possiamo concludere che il rischio condizionato di Y ∗ è min{η(X), 1 − η(X)}. Ora, poiché il rischio incondizionato è uguale al valore atteso del rischio condizion- ato, per la legge delle aspettative iterate: E[L(Y, Y ∗ (X))] = E[E[L(Y, Y ∗ (X))|X]] = E[min{η(X), 1 − η(X)}] (12) Dato l’input x, il rischio condizionato sarà: E[L(Y, Y ∗ (x))|X = x] = min{η(x), 1 − η(x)} (13) 1.3.2 Problema multiclass Analizziamo il problema seguente. Y ∈ Y , con Y = {Y1 ,... , YK }. Per brevità Y = {1,... , K}; Supponiamo di conoscere i posterior: P (Y = k|X) = ηk (X); Consideriamo la Loss 0/1 : L(Y, Ŷ (X)) = I{Y 6= Ŷ (X)}. 9 Il classificatore ottimale di Bayes è dato da Y ∗ (X) = argmax ηk (X), ossia k=1,..., K Y ∗ (x) = k se la P (Y = k|X = x) è massima sull’insieme di tutti i posterior calcolati in x. Supponiamo che Y = {A, B, C}, X = {1, 2, 3, 4} e che i posterior siano: x=1 x=2 x=3 x=4 ηA (x) = P (Y = A|X = x) 0.0000 0.0090 0.2531 0.7895 ηB (x) = P (Y = B|X = x) 0.3226 0.8442 0.7420 0.2105 ηC (x) = P (Y = C|X = x) 0.6774 0.1467 0.0049 0.0000 Si calcola facilmente che: x=1 x=2 x=3 x=4 Y ∗ (x) C B B A Si dimostra che, in corrispondenza del punto x, il rischio condizionato del classificatore Bayesiano è: E[L(Y, Y ∗ )|X = x] = 1 − max η(x) = 1 − max P (Y = k|X = x) (14) k k Da ciò ne consegue che il BER sotto la Loss 0/1 nel problema multiclass è dato da: E[L(Y, Y ∗ (X))] = 1 − E[max ηk (X)] = 1 − E[max P (Y = k|X)] (15) k k Come si costruisce il vero classificatore ottimale? Dovrei conoscere la distribuzione congiunta di X e Y (ossia i posterior). La determinazione dell’EPE è importante non solo per quantificare il livello di errore che ci dobbiamo aspettare, ma anche per scegliere tra metodi e algoritmi alternativi. Tuttavia, il calcolo analitico del ris- chio atteso è spesso impossibile, anche quando si conoscono i posterior. Ci sono due meta-modi per costruire un classificatore: 1. Osserviamo sul train set le Y e le X e stimiamo la distribuzione condizionata di Y |X (i posterior) attraverso vari metodi, come i modelli di regressione. A questo meta-modo corrispondono i classificatori LDA, QDA e kNN; 2. Costruiamo direttamente un modello di regressione che approssima η(X). A tale meta-modo corrispondono i classificatori di tipo logistico. 1.4 Modello Bernoulliano-logistico Lineare ′ ′ Tipicamente, in un problema 2-class, Xn = {(y1 , x1 ),... , (yn , xn )} è il nostro ′ dataset , dove yi ∈ {0, 1}, e xi contiene features numeriche. 10 Assunzione 1. La Popolazione obbedisce al modello Y |X ∼ Be(η(X)). Poiché η(X) = E[Y |X], η(X) è una funzione di regressione. Inoltre, i posterior sono: η(X) = P (Y = 1|X) e 1 − η(X) = P (Y = 0|X) (16) Assunzione 2. Il campione {(Yi , Xi )} è una successione di v. c. iid, dove Yi |Xi ∼ Be(η(Xi )) ∀ i = 1,... , n. Quindi: ( η(Xi ) → se yi = 1 P (Yi = yi |Xi ) = (17) 1 − η(Xi ) → se yi = 0 In forma compatta, la funzione di probabilità condizionata di Yi sarà: P (Yi = yi |Xi = xi ) = η(xi )yi [1 − η(xi )]1−yi (18) Per uno dei principî base della probabilità, possiamo fattorizzare la probabilità con- giunta di Xi e Yi nella condizionata per la condizionante: P (Yi , Xi ) = P (Yi |Xi )P (Xi ) (19) Solo la parte condizionale del modello dipende da η. Quindi, quando stimeremo η, l’altra la potremo ignorare. In questo tipo di problemi si tende a considerare solo la parte condizionata, perché, generalmente, è quella che dipende dal parametro incognito. Facciamo una parentesi sulle funzioni logit e logistic. Odds. L’odds di un evento è il rapporto fra la sua probabilità di successo e quella di insuccesso. Quindi, sia A un evento qualsiasi con P (A) = p, allora: p odds(A) = (20) 1−p Quando il denominatore è piccolo, gli odds perdono di significato, perché il rap- porto assume variazioni non più interpretabili. Ecco perché si passa al logaritmo (trasformazione monotòna crescente: cambia la scala, non l’interpretazione):   p logit(p) = ln ∈ R (21) 1−p Funzione logistca (o anti-logit/sigmoid). È l’inversa del logit: eδ 1 Se δ = logit(p) → logistic(δ) = p = = (22) 1 + eδ 1 + e−δ Funzione soft-max. È un metodo di ponderazione. Dati i numeri {z1 ,... , zM }: e zj sof tmax(zj ) = PM ∈ (0, 1) (23) m=1 e zm 11 Modelli con costante. Dati i due modelli: yi = β1 + β2 xi2 +... + βp xip e yi = β1 xi1 + β2 xi2 +... + βp xip (24) Il primo equivale al secondo, dove si pone che x1 = 1 ∀ i. Possiamo considerare l’intercetta come il coefficiente di un regressore che vale 1 per ogni osservazione. Assunzione 3. Il log-odds dell’evento {Yi = 1|Xi } è una funzione lineare dei predittori: η(Xi ) Xp log[ ] = β1 Xi1 +... + βp Xip = βj Xij = g(Xi ; β) (25) 1 − η(Xi ) j=1 In tale espressione, β è il vettore dei parametri ed η(X) è la funzione delle features che controlla la probabilità di Y = 1. L’obiettivo è stimare β per poter stimare i logaritmi degli odds, e, quindi, gli η(·). Infine, mettendoli nella Y , ottengo i posterior delle due classi. Dall’equazione precedente ricavo: ∑p β X eg(Xi ; β) e j=1 j ij η(Xi ; β) = logistic[g(Xi ; β)] = g(X ; β) = ∑p (26) 1+e i β X 1 + e j=1 j ij Come costruiamo il classificatore Bayesiano? Supponiamo di calcolare lo stimatore β̂ per β e la sua stima b usando un training set. Allora, la stima di η(·) sull’i-mo punto osservato è: ∑ p bj xij e j=1 η(xi ; b) = ∑p (27) bj xij 1+e j=1 Quindi, costruisco il seguente classificatore: ( 1 → se η(xi ; b) ≥ 1 ŷ(xi ; b) := 2 (28) 0 → se η(xi ; b) < 1 2 Si noti che, quando η(xi ; b) ≥ 12 , allora logit[η(xi ; b)] ≥ 0. Quindi, lo possi- amo riscrivere anche così: ( 1 → se logit[η(xi ; b)] ≥ 0 ŷ(xi ; b) := (29) 0 → se logit[η(xi ; b)] < 0 1.5 Stima di β̂ Il metodo di stima più diffuso è quello ML, perché assicura asintoticamente stime non distorte, consistenti, normali ed efficienti. P [(Yi , Xi ); β] = P (Yi |Xi ) ∗ P (Xi ) = η(Xi ; β)Yi [1 − η(Xi ; β)]1−Yi ∗ P (Xi ) 12 Lo stimatore ML β̂ è definito come: Y n β̂ := argmax Ln (β) dove Ln (β) := P [(Yi , Xi ); β] (30) β i=1 Il termine P (Xi ) non coinvolge β; quindi, l’espressione da risolvere sarà: Y n β̂ := argmax { η(Xi ; β)Yi [1 − η(Xi ; β)]1−Yi } (31) β i=1 Procediamo con la log-Verosimiglianza: Y n X n ln (β) = log{ η(·)Yi [1 − η(·)]1−Yi } = log{η(·)Yi [1 − η(·)]1−Yi } = i=1 i=1 X n = {Yi log[η(·)] + (1 − Yi )log[1 − η(·)]} = i=1 X n X p ∑p βj Xij = [Yi βj Xij − log(1 + e j=1 )] i=1 j=1 Dato il training set X, la stima di β̂ viene calcolata risolvendo: X n X p ∑p βj xij b := argmax [yi βj xij − log(1 + e j=1 )] (32) β i=1 j=1 Tale espressione non ha una forma chiusa. Si procede o tramite algoritmi: IRLS (Iteratively Reweighted Least Squares): a ogni step si risolve un prob- lema di minimi quadrati pesati (funziona, ma senza garanzie di convergenza); Coordinate-descent (dà garanzie di convergenza). Oppure tramite software (qui R): stats::glm consente di calcolare b perché il modello di regressione logistico è un caso particolare di ”modello lineare generalizzato”; glmnet::glmnet si basa su un’implementazione dell’algoritmo coordinate-descent, ed è molto indicato nei problemi di Alta Dimensionalità. 1.6 Proprietà di β̂ ed Inferenza Se i dati sono generati da un modello che approssima quello Bernoulliano-Logistico, i parametri stimati convergono in probabilità ai veri parametri che generano i dati, 13 da cui possiamo aspettarci che le stime dei posterior convergano ai veri valori della Popolazione per campioni molto grandi. Quindi, asintoticamente, mi devo as- pettare che la classificazione raggiunge il BER (Consistenza). Inoltre, si dimostra che, sotto certe condizioni, la distribuzione di β̂ è una Normale multivariata. I coefficienti del modello misurano l’impatto marginale delle covariate sui log- odds in favore di Y = 1: ∂ logit[η(Xi ; β)] = βj (33) ∂Xj Da ciò si dimostra che 100 · (eβj − 1) misura il tasso di variazione percentuale dell’odds in favore di Y = 1 quando il predittore Xj varia di una unità. 1.7 La linearità della Decision Boundary Per definizione, i modelli sono falsi, perché sono solo un’astrazione mentale per descrivere una realtà. Quindi non esistono modelli veri, ma solo più o meno buoni. Molto è legato al tipo di separazione che le classi richiedono. Ogni modello ha una sua struttura matematico-geometrica e probabilistica che adegua a separare certe geometrie di dati e non altre. Perché il modello Bernoulliano-Logistico è classificato come lineare? La lin- earità non sta solo nella funzione g(·), ma anche nel modo in cui separa le classi. Prendiamo il dataset baby.dat dove Y = {0, 1} = {not-Survive, Survive} e sup- poniamo di aver classificato usando solo due variabili: Weight ed Age; quindi, con X = (1, X2 , X3 ). Assegniamo xi a Y = 1 se logit[η(xi ; β)] ≥ 0 o se η(xi ; β) ≥ 12. La funzione g(·) è data da: logit[η(xi ; β)] = β1 + β2 W eight + β3 Age = β1 + β2 X2 + β3 X3 (34) Lo spazio delle features, ossia l’insieme di tutte le combinazioni dei loro possibili valori, in linea generale, è R2 ; quindi, la Decision Boundary è data da: B = {x ∈ R2 : β1 + β2 x2 + β3 x3 = 0} (35) Da questa, possiamo ottenere l’equazione della linea: x3 = β1 β3 − β3 x2. β2 General- izzando, la DB è quell’iperpiano che separa le classi: η(x; β) B = {x ∈ Rp : log[ ] = 0} (36) 1 − η(x; β) Sostituendo il log-odds col suo modello lineare, otteniamo: B = {x ∈ Rp : β1 x1 +... + βp xp = 0} (37) 14 1.8 Error rate in-sample Analizziamo il problema della stima del tasso di errore. La Loss è una v. c. Bernoul- liana che assumerà i valori di 1 o 0 a seconda che il classificatore misclassifichi o meno. Poiché la Loss è un indicatore, la posso riscrivere come E[I{Y 6= Ŷ }]. Se vogliamo ottenere una stima di questa quantità astraendola alla Popolazione, ma usando le informazioni contenute nel train set, cosa possiamo fare? Applichiamo il principio analogico dell’Inferenza: non sappiamo che fare? Facciamo sul campione l’analogo campionario di quello che faremmo sulla Popo- lazione. Fare l’analogo campionario significa operare con la distribuzione empirica. Quindi: Sostituiamo Y e Ŷ con le loro versioni campionarie, Yi e Ŷi ; ∑n Sostituiamo E[·] con n1 i=1 [·] Stimiamo la media campionaria delle loss: 1X 1X n n L¯n = L(Yi , Ŷi ) = I{Yi 6= Ŷi } (38) n i=1 n i=1 ∑n La versione osservata di tale oggetto sarà n1 i=1 I{yi 6= ŷi }, e corrisponde alla proporzione dei punti misclassificati nel train set. Si dimostra, però, che tale stima- tore dell’Error Rate avrebbe proprietà statistiche scadenti: ha un sacco di bias e una varianza enorme. 1.9 Regressione logistica per i problemi multi-class Nel contesto multi-class, la Y diventa una variabile categoriale multipla in K livelli. Una cosa che possiamo sempre fare è dummificarla, ossia trasformarla in K − 1 variabili dicotomiche. Supponiamo di avere n oggetti: l’oggetto i-mo apparterrà a una delle K classi. Definiremo K − 1 dummy tali che la j-ma dummy sarà 1 se la classe dell’oggetto i-mo è j, 0 altrimenti: ( (j) 1 → se Yi = j Yi = (39) 0 → altrimenti Quindi, definiamo le variabili binarie Y (1) ,... , Y (K−1). Ce ne servono K − 1 perché, se valessero tutte 0, l’oggetto automaticamente apparterrà alla classe K. Sulla base di tale meccanismo, il primo metodo di classificazione logistica multi- class si chiama Regressione Logistica One VS All. 15 1.9.1 Regressione Logistica One VS All Consideriamo la k-ma componente di Y , ossia Y k , e denotiamo: I posterior della classe k-ma date le features: η k (X) = P (Y k = 1|X); k ∑p I log-odds: log( 1−η η (X) k (X) ) = j=1 βjk Xj ; Il vettore di coefficienti β (k) che modella i log-odds in favore della classe k. Dobbiamo stimare K · p coefficienti, corrispondenti a {β (1) ,... , β (K) }. Avrò un vettore di parametri per ogni dicotomizzazione. La procedura è questa: abbiamo in input K dummy e la matrice delle features X. Per ogni classe stimiamo i parametri del modello logistico e andiamo a prevedere i posterior. Otterremo k vettori di posterior; poiché, però, questi son fatte con re- gressioni separate, può capitare che alcune somme non facciano 1. Quindi, bisogn- erà normalizzarle. Alla fine otterremo una tabella dove sulle righe abbiamo le os- servazioni e sulle colonne i posterior per ogni classe. Riga per riga guarderemo al posterior maggiore, e assegneremo l’oggetto alla relativa classe. Tale metodo ha una enorme inefficienza, dovuta al fatto di dover stimare una quantità smisurata di parametri per K classi separandola in K − 1 problemi dis- tinti, producendo una varianza più alta nelle stime e, quindi, nella classificazione. Inoltre, per funzionare bene, ci deve essere una certa indipendenza fra i posterior di ogni unità, ma tale ipotesi non sempre è plausibile. 16 1.9.2 Modello Bernoulliano Generalizzato-Logistico (o Multinomial Logit) La Multinomial Logit è una distribuzione semplice su K categorie: la Y può as- sumere K valori diversi, ognuno con probabilità η 1 (X),... , η K−1 (X). L’assunzione 1 del modello Bernoulliano-Logistico viene sostituita con la sua generalizzazione: Y |X ∼ GenBer[η 1 (X),... , η K−1 (X)] dove P (Y = k|X) = η k (X) (40) ∑K−1 La probabilità di appartenere alla k-ma classe sarà η (X) = 1 − K k=1 k η (X), ossia il complemento a uno degli altri posterior. Come nell’ass. 3 del modello Bernoulliano-Logistico, assumiamo che: η k (Xi ) Xp ∀ k = 1,... , K − 1 : log[ K ] = g(Xi ; β k ) = β1k Xi1 +... + βpk Xip = βjk Xij η (Xi ) j=1 (41) Ossia, consideriamo il logaritmo del rapporto tra il posterior della classe k-ma e quello della baseline. Avremo, quindi, K − 1 funzioni di regressione con cui, una volta note, calcoleremo i posterior. Per tutti i k 6= K avremo: k eg(Xi ; β ) η k (Xi ) = PK (42) g(Xi ; β k ) k=1 e Invece, per la reference class (k = K), avrò: 1 η K (Xi ) = PK−1 k) (43) 1+ k=1 eg(Xi ; β Per ottenere il classificatore Bayesiano, otteniamo le stime b1 ,... , bK , e calcol- iamo l’assegnazione: yˆi = argmax η(xi ; bk ) (44) k=1,..., K Si può dimostrare che, in questo caso, i posterior si calcolano come il softmax delle funzioni di regressione. I parametri da stimare sono di meno, p · (K − 1), perché abbiamo K − 1 funzioni di regressione, ognuna con p features (nel One VS All sono K · p). 17 18  2 Geometria e probabilità nello spazio Multivariato 2.1 I concetti di Distanza e di Intorno Quando parliamo di probabilità multivariata non è che distribuzione hanno i val- ori sulla X, ma che distribuzione hanno i punti in uno spazio più complicato. Che probabilità abbiamo che i punti stiano in una zona piuttosto che in un’altra? Cominciamo a vedere come i concetti geometrici vengono generalizzati ad am- bienti di lavoro che non possiamo nemmeno vedere, come la quarta dimensione. Il concetto fondamentale di ogni geometria è la distanza. Qualsiasi distanza deve obbedire ai seguenti assiomi (per ogni {a, b, c}): Identità: a = b ⇐⇒ dist(a, b) = 0; Simmetria: dist(a, b) = dist(b, a); Sub-additività: dist(a, b) ≤ dist(a, c) + dist(c, b). Tale assioma af- ferma che, dati i tre oggetti a, b e c, la distanza tra a e b non può essere maggiore della somma delle distanze tra a e c, e tra c e b. Supponiamo di avere 2 punti posti su una retta: x1 = 10 ∧ x2 = 6. Quindi, abbiamo che x1 , x2 ∈ R. Quant’è la distanza? Diremmo banalmente 4. La nozione di distanza non presuppone una direzione: in mente, quando pensiamo alle distanze, abbiamo il concetto di simmetria (le distanze non simmetriche non sono distanze). Infatti, la distanza tra x1 a x2 è la stessa che tra x2 e x1. La misura che abbiamo intuitivamente usato è: dist(x1 , x2 ) = |x1 − x2 | = |10 − 6| = 4 (45) Consideriamo i punti x1 , x2 ∈ R2 , con x1 = (3, 2) e x2 = (8, 6). Possiamo calcolare la distanza come lunghezza del segmento x1 x2 col Teorema di Pitagora: q dist(x1 , x2 ) = (8 − 3)2 + (6 − 2)2 = 6.4 (46) Ma se siamo in dimensioni maggiori? La generalizzazione del teorema di Pitagora nel caso q-dimensionale si chiama Distanza Euclidea. Dati i punti a, b ∈ Rq , dove a = (a1 ,... , aq ) e b = (b1 ,... , bq ), la Distanza Euclidea tra a e b è: v q uX u q dist2 (a, b) = (a1 − b1 )2 +... (aq − bq )2 =t (aj − bj )2 (47) j=1 19 In generale, possiamo definire una classe di distanze che include anche quella eu- clidea: la classe delle distanze di Minkowski (dove p è l’iperparametro): v q uX 1 u q = t p distp (a, b) = p |a1 − b1 |p +... + |aq − bq |p ) p |aj − bj |p (48) j=1 Quando p = 1 abbiamo la distanza di Manhattan (o distanza Taxicab): X q dist1 (a, b) = |aj − bj | (49) j=1 Quando p → ∞, abbiamo la distanza di Chebyshev: dist∞ (a, b) = lim distp (a, b) = max{|aj − bj |} (50) p→∞ j Se applichiamo queste diverse forme di distanza allo stesso spazio geometrico, come lo spazio Euclideo bidimensionale, che cambia? Cambia il concetto di Intorno. L’intorno è l’insieme dei punti che stanno nello spazio Rq che, in base a una certa nozione di distanza, stanno a una distanza minore di ε; ossia: B(a; ε) = {x ∈ Rq |distp (x, a) < ε} (51) Si dice ”palla centrata su a di raggio ε”. Si dice ”palla aperta” se non se ne definisce un bordo, quindi un limite (non c’è l’uguale nella disuguaglianza). Supponiamo di essere in uno spazio R2. Mettiamo di aver preso tanti punti (quelli in grigio) e che il nostro a sia fissato in (0, 0). Abbiamo stabilito che ε = 1. Per tutti i punti del grafico abbiamo calcolato la distanza di Taxicab, quella Euclidea e quella di Chebyshev. Abbiamo colorato di rosso solo quei punti la cui distanza era minore di ε. Lo spazio rosso risultante è la visione grafica dell’Intorno. Può risultare controintuitivo pensare che i punti sul bordo dell’Intorno calcolato con Taxicab o con Chebyshev siano a distanza 1, ma è così. 20 2.2 Variabili Casuali Multivariate 2.2.1 Cosa sono le v. c. Multivariate? Una v. c. è una funzione che va dallo spazio degli eventi a un insieme matem- aticamente ben compreso. Nel corso di studi, abbiamo visto come si calcolano le probabilità, come si rappresenta la loro legge di probabilità e alcune caratteris- tiche annesse come i momenti primo (media) e secondo (variabilità, ossia quanto è dispersa la distribuzione di probabilità intorno al centro di massa). La varianza dipende dai primi due momenti: V ar[X] = E2 [X] − E[X]2. Per una v. c. a media 0, la varianza coincide col momento secondo. Trasportiamo queste nozioni al caso multivariato, iniziando dalle basi con R2. Adesso, la variabile è una funzione che trasforma gli eventi in due numeri. An- che nel caso multivariato, sotto certe condizioni, i concetti di momento primo e secondo caratterizzano la stessa questione del caso univariato. Anche nel caso multivariato il momento primo non sempre caratterizza il centro di massa della distribuzione (lo fa quando abbiamo distribuzioni unimodali e simmetriche). In- oltre, il momento primo ci dice più o meno dove si ammassa la maggior parte dei punti che hanno la probabilità più alta. 2.2.2 Momento primo, momento secondo e matrice VCV Come si definisce il momento primo nel caso multivariato? Il momento primo di un vettore corrisponde al vettore dei momenti primi. In altre parole, non è al- tro che il vettore dei momenti marginali: prendiamo la distribuzione congiunta e ignoriamo completamente una delle due marginali:     E[X1 ] µ1 E[X] := = =µ (52) E[X2 ] µ2 Il vettore µ, nel nostro caso (µ1 , µ2 ), definisce il punto intorno al quale si con- centra parecchia massa centrale di probabilità. E per il momento secondo? Nel caso univariato, esso ha a che fare con la nozione di quadrato. Se X è un vettore, il quadrato di un vettore non esiste. Dob- biamo trovare una nozione di momento secondo che ricostruisca un concetto di variabilità simile a quello che la varianza svolge nel caso univariato. Il momento secondo di un vettore di v. c. è il valore atteso del suo prodotto esterno, ossia:     ′ X12 X1 X2 E[X12 ] E[X1 X2 ] E[XX ] = E 2 = (53) X2 X1 X2 E[X2 X1 ] E[X22 ] 21 Se volessimo calcolare la varianza marginale di X1 , la calcoleremmo in questo modo: σ12 = V [X1 ] = E[(X1 − µ)2 ] (54) Inoltre, notiamo che, quando definiamo la Covarianza fra due variabili, questa è uguale al valore atteso del prodotto delle due variabili centrate: Cov(X1 , X2 ) = σ12 = E[(X1 − µ1 )(X2 − µ2 )] (55) Sapendo tutto ciò, svolgiamo i calcoli:     ′ E[(X1 − µ1 )2 ] E[(X1 − µ1 )(X2 − µ2 )] σ12 σ12 E[(X−µ)(X−µ) ] = = E[(X2 − µ2 )(X1 µ1 )] E[(X2 − µ2 )2 ] σ21 σ22 (56) Tale espressione indica la matrice di varianze e covarianze Σ del vettore X, dove sulla Diagonale principale abbiamo le varianze marginali e altrove le covarianze. Essa ci informa sulla dispersione che la distribuzione produce intorno al vettore medio. Estrapolando al caso generale con p variabili, avremo:     E[X1 ] µ1 E[X] := ...  = ... = µ (57) E[Xp ] µp Il momento secondo sarà una matrice p · p:   ′ E[X12 ]... E[X1 Xp ] E[XX ] = E[Xr Xq ] .........  (58) E[Xp X1 ]... E[Xp2 ] Il vettore scarto sarà:   X 1 − µ1 (X − µ) = ...  (59) X p − µp La matrice di varianze e covarianze p · p, invece, è:   σ12... σ1p Σ = σrq = ......... (60) σp1... σp2 Quindi, nel caso multivariato la dispersione dei punti intorno alla media è caratter- izzata essenzialmente da due cose: le varianze marginali e i legami lineari a coppie delle variabili. Vediamo le proprietà di Σ: È simmetrica, perché la covarianza a coppie è simmetrica. Questo implica che i suoi elementi sono meno di p2. infatti, al massimo vi saranno: p·p−p p(p + 1) p + = (61) # varianze marginali 2 2 # covarianze off-diagonal # elementi totali È semidefinita positiva, ossia i suoi autovalori sono tutti reali e non negativi. Vediamo alcune situazioni particolari. 22 Modello sferico: incorrelazione + omoschedasticità.     σ2 0... 0 1 0... 0 ... σ2 0  0 V [X] = ...  = σ 2 ... 1...  = σ 2 Ip (62) ............ ............ 0...... σ2 0...... 1 Tale matrice è detta ”sparsa”, nel senso che la maggior parte dei suoi elementi sono nulli. Il numero di parametri che la caratterizza è solo 1, ossia σ 2. Modello diagonale: incorrelazione + eteroschedasticità.   σ12 0... 0 ... σ22... 0 V [X] =  ...  (63)......... 0...... σp 2 Anche in questo caso si tratta di una matrice sparsa, con # parametri = p. Ce ne sta un altro di caso, in cui tutto può essere o meno correlato. La corre- lazione, nel caso univariato, fra due variabili è data da: cov[Xr , Xq ] σrq ρrq = Cor[Xr , Xq ] = p = (64) V [Xr ]V [Xq ] σr σq La correlazione di un vettore è una matrice di correlazione p · p (simmetrica):   1 ρ12... ρ1p ρ21 1... ρ2p   cor[X] = R =   (65)............ ρp1...... 1 Poiché vogliamo analizzare le questioni anche dal punto di vista pratico, estraendo campioni da popolazioni multivariate e vedere da cosa si caratterizza la loro disper- sione, tutte queste quantità le dobbiamo stimare, e procederemo con l’equivalente campionario delle quantità studiate. Media campionaria (osservata). Ricordiamo che il momento primo di un vet- tore di variabili è uguale al vettore dei loro momenti primi:    1 Pn  1X n x̄1 n i=1 xi1 x̄n =   xi =... = ...  (66) n i=1 x̄ 1 Pn p x n i=1 ip 23 Matrice di varianze-covarianze non distorta (osservata). 1 X n ′ Sn = (xi − x̄n )(xi − x̄n ) (67) n − 1 i=1 Per ricavare da questa la matrice di correlazione campionaria, dobbiamo costruirci una matrice diagonale Q (p · p) con i reciproci delle radici delle varianze. √1  S11 0... 0  0 √1......    Q= S22  (68) ............  0...... √1 Spp La matrice di correlazione campionaria p · p, con gli 1 e le correlazioni di Pearson a coppie, sarà il prodotto ”sandwich”: Rn = QSn Q. Tutti questi sono stimatori che hanno analoghe proprietà alle corrispondenti univariate: sono non distorti (in media stimano la quantità giusta da stimare) e hanno lo stesso comportamento asintotico. Se prendiamo un numero sufficiente- mente grande di osservazioni, le stime campionarie convergono ai veri valori della p p popolazione: X̄n → µ; Sn → Σ; Rn → R. p 2.3 Dispersione vs direzione Quando parliamo di v. c., il modo più semplice per capire la casualità è speri- mentare con la casualità. Prendiamo un vettore di 2 v. c. da una distribuzione che ha come valore medio il punto (0, 0) e, in particolare, con una matrice VCV sferica e con varianza unitaria. Quindi: Σ = σ 2 I2 , con σ 2 = 1. Per vedere i campioni, abbiamo generato una matrice di dati, dove ogni punto è una riga e ha coordinate (X1 , X2 ). Calcoliamo le 2 medie campionarie e costru- iamo il vettore medio. Notiamo che, geometricamente, esso sta al centro della dis- tribuzione, e che i punti stanno più azzeccati verso il centro e poi, allontanandosi, si diradano sempre di più (questo effetto è catturato dalla matrice VCV). Come dobbiamo pensare quando pensiamo alle marginali? Quando marginal- izziamo rispetto a X1 , stiamo fregandocene dell’altra, e stiamo vedendo che valore hanno i punti sull’asse delle ascisse, come se si proiettassero tutti sull’asse orizzon- tale. Volendo, Possiamo proiettare i punti anche su una direzione. Una direzione è qualsiasi retta passante per un punto. Per un punto passano infinite direzioni, per due punti ne passa solo una. Nel campo univariato possiamo guardare ai dati 24 in una sola direzione, ma nel multivariato, possiamo guardare alla nuvola di punti attraverso infinite direzioni. Prendiamo due direzioni passanti per il punto centrale, e proiettiamo lo scatter lungo le due direzioni. Otteniamo, così, due stripchart statisticamente uguali: nel caso sferico, in qualunque direzione guardiamo, passando dal centro, la dispersione dei dati rispetto al centro è omogenea (non cambia). Ora prendiamo un caso dove la matrice VCV è densa. Facciamo un esempio dove abbiamo due variabili estremamente negativamente correlate e con varianze diverse. Avremo uno scatter dei punti approssimativamente ellissoidale.     3 −4.7 1 −0.9 Σ= ⇒R= (69) −4.7 9 −0.9 1 Se c’è una forte correlazione fra le marginali, lo scatter non avrà più una forma sferica, ma diventerà pesantemente direzionato, e la dispersione non sarà uniforme in tutte le direzioni. Se la correlazione è forte ci sono direzioni in cui c’è una forte dispersione rispetto al centro, e direzioni in cui è bassissima. Le direzioni che hanno la minore e maggiore varianza sono nello spettro della matrice VCV, ossia nei suoi autovalori e autovettori. La matrice VCV contiene tutta l’informazione sulla geometria dei dati che ha a che fare con la nozione di dispersione intorno al momento primo della dis- 25 tribuzione. La matrice VCV la definiamo come il momento secondo del vettore scarto, dove il vettore scarto è il vettore centrato X − µ. Quando abbiamo introdotto i vettori, abbiamo parlato di versore. I versori sono i vettori che hanno le direzioni delle coordinate del sistema. Ovviamente non esistono solo le direzioni determinate dagli assi di riferimento: abbiamo variabilità e dispersione in tutte le varie direzioni. Ciò che determina il fatto che ci sono alcune direzioni in cui c’è più o meno varianza è legato alla struttura di covarianza dei dati. Se abbiamo più struttura di correlazione, avremo direzioni in cui c’è un’estrema varianza ed altre in cui ce n’è molto meno. Se, invece, c’è incorrelazione, le varianze sono dominate essenzialmente da quello che succede nelle direzioni principali, e in tutte le altre direzioni la dispersione sarà simile. Non possiamo giudicare quanta variabilità ci sia in un dataset multivariato guardando allo scatterplot, perché sarebbe infattibile per un essere umano oltre la terza dimensione. Quindi, come facciamo a scrivere una funzione della matrice VCV che ci dia un’indicazione della variabilità reale? Ci sono due metodi. Il primo è considerare la Varianza Totale, ossia la somma delle varianze marginali: X p X p T V (X) = traccia(Σ) = σj2 = V [Xj ] (70) j=1 j=1 Questa corrisponde alla traccia della matrice (la traccia di una matrice quadrata corrisponde alla somma degli elementi sulla diagonale principale). Logicamente, la TV cattura completamente la nozione di dispersione multivariata solo nel caso sferico e nel caso diagonale. Infatti, tale misura trascura sistematicamente tutto ciò che c’è al di fuori della diagonale principale. In tutti gli altri casi, quando c’è correlazione fra i dati, ci saranno alcune di- rezioni in cui c’è più o meno dispersione, e questa discrepanza è tanto più forte quanto più forte è la struttura di correlazione. La misura più generale di disper- sione in questi casi è la Varianza Generalizzata: GV (X) = det(Σ) (71) La GV tiene conto di tutte le varianze e le covarianze. La GV sintetizza, quindi, la dispersione in tutte le direzioni. Se il determinate è grande, allora avremo una forte dispersione intorno al vettore medio; altrimenti, se è prossimo allo zero, tutti i punti stanno ammassati intorno al vettore medio. Si può dimostrare che, se prendiamo uno scatter multidimensionale e ne cal- coliamo la matrice VCV, la radice quadrata del determinante è una quantità pro- 26 porzionale al volume della struttura geometrica che contiene tutti i punti osser- vati. Quindi, quanto maggiore è il determinante, tanto più i punti hanno bisogno di una superficie maggiori. Quando il determinante, invece, è prossimo allo zero, stiamo operando con uno scatter di punti che può essere incapsulato in un’urna convessa, ossia in una struttura geometrica con un volume molto piccolo. Un caso limite è quando abbiamo distribuzioni con matrice VCV con una GV pari a 0. Ciò significa che il determinate sarà 0: si parlerà di distribuzione a varianza singolare o di distribuzione singolare. Tutto starà ammassato sul vettore medio, ossia, che per qualunque direzione passante per il vettore medio non c’è variabilità. 2.4 Trasformazioni lineari di variabili casuali inRp Spesso non lavorariamo coi dati originali, ma con quelli trasformati. Una classe di trasformazioni essenziale è quella delle trasformazioni lineari. Che significa una trasformazione lineare nel caso multivariato? Significa un sistema lineare. Supponiamo di avere un vettore di v. c. X ∈ Rp , con E[X] = µ e V [X] = Σ. Supponiamo anche di prendere un vettore di costanti t ∈ Rp e una matrice A costante (di coefficienti) di dimensione p · p. Una trasformazione lineare di X in dimensione p non è altro che un sistema lineare: Y = t +A·X (72) p·1 p·1 p·p p·1 Come succede nel caso univariato, c’è un modo semplice per andare dai momenti originali a quelli delle trasformate e viceversa. Il vettore medio è un operatore lin- eare affine, ossia trasforma il momento così come si trasforma la variabile. La vari- anza, invece, è più complessa: così come la varianza di una trasformazione lineare nel caso bivariato è fatta come coefficienti di trasformazione angolare al quadrato per la varianza originale, nel caso multivariato è simile: ′ E[Y ] = t + Aµ e V [Y ] = AΣA (73) Si può notare come, fissati i due coefficienti t ed A, i primi due momenti di Y siano funzioni dei primi due momenti di X. Nel caso univariato, le trasformazioni lineari fanno due operazioni: shift (spostano il centro della distribuzione) e scale (allargano o restringono il range di misurazione dei dati). Nel caso multivariato ne fanno molte di più (tra le più importanti): il mirroring (rotazione) e lo stretching (la deformazione). 27 2.5 Variabile Casuale Normale Multivariata Come nel caso univariato, esiste una distribuzione Normale anche nel caso multi- variato: prende il nome di Normale Multivariata. Tale distribuzione è controllata da due parametri, che coincidono col vettore medio e con la matrice VCV. Scriver- emo: X ∼ Np (µ, Σ). La sua funzione di densità è: 1 ′ 1 Σ−1 (x−µ) ϕ(x; µ, Σ) = p e− 2 (x−µ) (74) det(Σ) · (2π)p Tale espressione è la generalizzazione della Normale Univariata; infatti, sostituendo 1 a p, si ottiene la sua forma più semplice. La forma quadratica posta all’esponenziale dopo − 12 prende il nome di Distanza di Mahalanobis. Tale quantità misura la dis- tanza dei punti rispetto al centro. La distribuzione Normale appartiene alla classe delle distribuzioni ellittico-simmetriche, che sono caratterizzate da due parametri: uno di centralità e uno di scatter. Nel caso sferico i density contours set della Normale (quando tagliamo la testa della Normale e proiettiamo la forma che ci rimane in giù) hanno forma circolare in due dimensioni, sferica in tre dimensioni e ipersferica in più dimensioni. Man mano che i cerchi si avvicinano al centro, staremo in zone di densità sempre più alte (e viceversa). Ovviamente, man mano che ci allontaniamo, i cerchi catturano una massa di probabilità di punti sempre maggiore. Quando i contours diventano di dimensioni infinite, stiamo calcolando una massa di probabilità pari a 1. quanto più è forte l’eteroschedasticità e la struttura di correlazione, i density contours diventano sempre più ellittici. Inoltre, quando abbiamo incorrelazione ed omoschedasticità, la matrice VCV ha una struttura diagonale, e le ellissi saranno parallele a uno o più assi. 2.5.1 Variabile Casuale Normale Multivariata Standard Anche nel caso Multivariato esiste la Normale Standard: è quella centrata sull’origine dello spazio euclideo, (0,... , 0), e avente matrice VCV pari alla matrice di iden- tità I in p dimensioni: Z ∼ N (0, Ip ). Quindi, la NMS è una distribuzione sferica e produce contours set circolari, sferici ed ipersferici. Molte distribuzioni multivariate si basano sulla perturbazione della Normale Standard attraverso delle trasformazioni lineari (nel caso multivariato si chiama trasformazione affine). Inoltre, ve ne sono altre (come la T-Student) che si otten- gono perturbando, per esempio, la matrice VCV di una Normale. 28  Modificando appena la Normale, la si può trasformare in quella che si chiama Normale Generalizzata, la quale comprende tutta una serie di Modelli Multivariati. Tantissimi teoremi limite si basano sulla Normale. Infatti, esistono le leggi cen- trali del limite anche in Rp , e molte statistiche multivariate, in qualche modo asin- toticamente, hanno distribuzioni legate alla Normale. La Normale Multivariata ha una chiara nozione di centralità, grazie alla pre- senza di un punto di simmetria (la media), intorno al quale, muovendoci in deter- minate direzioni, possiamo prevedere con precisione il comportamento della dis- tribuzione. Questo punto è, quindi, definibile come il suo centro di massa. Perché il centro di massa è una misura di posizione così chiara? Una distribuzione non simmetrica non presenta un centro di massa ben definito. Inoltre, una dis- tribuzione con code pesanti, anche se simmetrica, può creare difficoltà. In questi casi, gran parte della probabilità è concentrata agli estremi, e il campionamento da tale distribuzione può portare a una media instabile o poco rappresentativa. La Normale Multivariata possiede una proprietà unica: le sue code decadono rapidamente verso 0, seguendo una velocità esponenziale rispetto al quadrato della distanza dal centro. Questo comportamento garantisce una stabilità intrinseca alla media, rendendola una misura di posizione robusta e affidabile. Non ci serve caratterizzare tutte le possibili Normali, perché conoscendo quella Standard, le possiamo produrre tutte. Se vogliamo rappresentare una Normale con vettore medio µ e matrice VCV Σ, possiamo rappresentarla come l’equivalente stocastico di questa trasformazione lineare. Cioè, se prendiamo la v. c. Multi- variata X dove X = µ + QZ, dove Z è una Normale Sferica, per le proprietà precedenti dei momenti delle trasformazioni lineari il valore atteso di X è uguale a ′ ′ µ + Q · 0 = µ, e la varianza è il prodotto sandwich QIp Q = QQ = Σ. ′ Come facciamo a trovare una Matrice Q tale che QQ = Σ?Questo si chiama problema di fattorizzazione di matrici, e vi sono tre strade: Fare la decomposizione spettrale di Σ (autovettori e autovalori), e moltipli- care la matrice degli autovettori per la radice della matrice degli autovalori; Fare la decomposizione di Cholesky e prendere il triangolo sinistro; Fare la decomposizione a valori singolari (SVD). 29 2.5.2 I density countours I density contours della Normale saranno di forma sferica o ellissoidale. Il centro dell’ellisse è determinato da µ e il volume e la direzione da Σ (il volume dipende dal determinante). Ogni density contour rappresenta una regione di punti intorno alla media e, quindi, conterrà una massa dei punti con una certa probabilità. Chiameremo Eα l’ellissoide di centro µ tale che, sotto una Normale, la proba- bilità che X appartenga quest’ellissoide sia uguale ad α: P (X ∈ Eα ) = α. Calco- lare questa probabilità significa, ad esempio in due dimensioni, calcolare la superfi- cie dell’ellissoide fino ad arrivarne a uno che misura α, che sarà pari alla probabilità che un punto estratto da una Normale stia là dentro. Possiamo dire che gli ellissoidi sono i corrispondenti dei quantili del caso univariato. 2.5.3 Distanza di Mahalanobis Possiamo costruire una nozione di distanza simile a quella euclidea, cioè, che alla stessa distanza, troviamo lo stesso livello di densità qualunque sia la direzione in cui ci muoviamo? Sì, e questa distanza si chiama Distanza di Mahalanobis. Supponiamo di avere due vettori in dimensione p, xr e xq. Si definisce dis- tanza di Mahalanobis tra xr e xq , in un sistema descritto dalla matrice Σ: q distM (xr , xq ) := (xr − xq )′ Σ−1 (xr − xq ) (75) Noi siamo interessati alla distanza tra un punto generico xi e il centro: q distM (xi , µ) := (xi − µ)′ Σ−1 (xi − µ) (76) Spesso ci riferiamo alla sua forma quadratica (Squared Mahalanobis Distance): ′ smd(xi , µ) = (xi − µ) Σ−1 (xi − µ) (77) La Distanza di Mahalanobis è una generalizzazione della distanza Euclidea. Vedi- amo perché essa si riconduce alla distanza Euclidea nel caso sferico. q r 1 distM (xi , xj ) : = (xi − xj )′ (σ 2 Ip )−1 (xi − xj ) = (xi − xj )′ ( 2 Ip )(xi − xj ) = σ v q uX p 1 1u 1 = (xi − xj ) (xi − xj ) = t ′ (xrj − xqj )2 = dist2 (xr , xq ) σ σ j=1 σ Quando calcoliamo la distanza di Mahalanobis con una matrice VCV sferica, sti- amo praticamente calcolando una quantità proporzionale alla distanza euclidea, e la proporzionalità è data dal fattore σ1 Ip. 30  Quando ci allontaniamo dal caso sferico, non la dobbiamo più pensare come a una distanza euclidea. È sempre una distanza nel senso che più e grande, tanti più punti stanno in quell’intorno, però non sarà simmetrica: dipenderà dalla direzione dove ci si muove, e sarà Σ a fare sì che i valori delle distanze di Mahalanobis vengano deformati a seconda della direzione. Qual è la relazione importante tra le Distanze di Mahalanobis e le dimensioni degli ellissoidi? Supponiamo di osservare un insieme di punti in R2 , e di calcolare la distanza di Mahalanobis con un’appropriata matrice VCV (stimando la matrice VCV osservata), e vogliamo calcolare l’SMD tra un punto osservato e il vettore medio stimato: smd(xi , x̄). Coloriamo in blu i punti che stanno a una distanza di Mahalanobis al quadrato non più di 2.77 (Perché questo numero? Perché l’smd segue una distribuzione Chi-quadrato con p gradi di libertà, e il 75-mo percentile è 2.77 quando p = 2). In sintesi, E0.75 è la regione di confidenza dove si troverà il 75% dei punti osservati. Come facciamo a sapere che è l’ellissoide di livello 0.75? Perché conosciamo la distribuzione, e anche a quale distanza troveremo l’ellissoide di dimensione 0.75. La distanza di Mahalanobis, anche in uno spazio non sferico, si comporta uguale in tutte le direzioni, anche se lo spazio è stato contratto, perché Σ corregge tutta la deformazione: rende sferico uno spazio non sferico. 31 32  3 Classificazione: LDA, QDA e kNN Nel contesto LDA, QDA e kNN non partiamo dalla costruzione di un modello per i posterior, ma da un modello per la distribuzione congiunta delle features, e solo da lì la marginalizziamo per ottenere i posterior. Tali metodi non funzionano nel caso di features non numeriche (come Maschio/Femmina). In realtà, dovrebbero essere anche continue, ma riescono a funzionare bene anche nel caso di features a interi, basta che abbiano un range di valori plausibile (Se abbiamo una dummy con 100 mila osservazioni, combinereremmo solo disastri). Supponiamo di conoscere, o di poter derivare, le prior probability (o prevalenze), ossia la probabilità incondizionata di appartenere a una classe: πk = P (Y = k). Supponiamo di conoscere le distribuzioni condizionate di classe (le Class Don- ditional Censity), ossia la distribuzione di X|Y = k. Non sappiamo quale sia davvero la distribuzione, ma sappiamo che essa è rappresentata da una certa fun- zione densità fk (x), dove x ∈ Rp e k = 1,... , K. Supponiamo di avere un intorno multivariato B intorno a x, B(x). Per cal- colare la probabilità che X appartenga a questo intorno, dovrei conoscere la dis- tribuzione delle X e vedere quanta massa di X sta dentro quest’Intorno. Per la legge delle probabilità totali, poiché le classi sono esaustive (un oggetto deve stare per forza in una classe), questa la possiamo scomporre come segue: X K P {X ∈ B(x)} = P {X ∈ B(x)|Y = k)}P (Y = k) (78) k=1 Se schiacciamo un intorno fino a farlo avere un’ampiezza infinitesima, il concetto di probabilità si trasforma in quello di densità. Infatti, la densità in un punto x è il limite della probabilità di appartenenza di X ad un intorno di x (se la distribuzione di X è continua). Infatti, la densità è una derivata (derivata di Radon-Nikodym). Quindi, la probabilità di stare in un intorno infinitesimale di x è proporzionale alla densità in x, la quale la possiamo rappresentare come una densità incondizion- ata delle X, m(x), data dalla somma delle probabilità incondizionate di X, i πk , per le densità condizionate calcolate in x: X K m(x) = πk fk (x) (79) k=1 Questa si chiama Mistura di Probabilità, ed è una combinazione convessa di dis- tribuzioni di probabilità. Una combinazione convessa è una combinazione lineare che ha tutti i coefficienti tra 0 e 1. 33  Da questa noi possiamo ricavare i posterior, ossia la probabilità che un punto appartenga a una classe dato che X stia in un intorno dato. Applicando l’identità delle probabilità condizionate, possiamo scrivere: P {X ∈ B(x)|Y = k}P (Y = k) P {Y = k|X ∈ B(x)} = (80) P {X ∈ B(x)} Restringendo al limite l’intorno di x, la probabilità condizionata che un oggetto appartenga alla classe k dato che le sue features stanno in una certa regione sarà: πk fk (x) πk fk (x) τk (x) = = PK (81) m(x) k=1 πk fk (x) Tali quantità si chiamano Posterior Weights (o Posterior Density) e sono il rap- porto tra il contributo della classe k-ma nella distribuzione di x fratto l’intera dis- tribuzione di x. Tau non misura le probabilità, ma le densità dei posterior, ossia quanto è densa di punti la probabilità condizionata di Y |X in un intorno infinites- imale di x. Il classificatore ottimale Bayesiano assegnerà le istanze alla classe che massimizza i posterior: Y ∗ (x) = argmax P {Y = k|X = x} = argmax τk (x) (82) k k Supponiamo di avere due classi e di confrontarne i posterior. Assegneremo l’unità i alla classe 1 se τ1 ≥ τ2. Notiamo, però, che τ1 e τ2 hanno lo stesso denominatore, che non dipende dalle classi, ma dalla distribuzione incondizionata delle features. Quindi, assegneremo a k quando il prodotto πk fk (x) è il maggiore possibile: π1 f1 (x) π2 f2 (x) τ1 (x) ≥ τ2 (x) =⇒ ≥ =⇒ π1 f1 (x) ≥ π2 f2 (x) (83) m(x) m(x) Applicandovi una funzione monotòna crescente, non cambia l’ordinamento: Y ∗ (x) = argmax πk fk (x) = argmax {log(πk ) + log[fk (x)]} (84) k k Nella pratica, però, osserviamo solo dei dati sulle Y e sulle X; quindi, i πk le fk (x) li dobbiamo stimare usando le informazioni nel train set. Ci sono 3 approcci: 1. LDA e QDA. Non conoscendo le funzioni di densità fk (·), assumiamo la loro Gaussianità, così traduciamo il problema della stima delle funzioni in un problema di stima dei loro parametri. Quindi, abbiamo K distribuzioni Normali, ognuna col suo vettore medio e la sua matrice VCV. A seconda di come parametrizziamo la matrice VCV e di quanto sia grande p, possi- amo avere una enormità di parametri. In generale, non è buona prassi usare i metodi classici nei casi ove il numero di features ecceda quello delle os- servazioni (Alta Dimensionalità). Questi approcci non scalano bene con p, 34  quando abbiamo dataset piccolissimi, o dataset grandi ma con un p enorme, o ancora quando abbiamo un numero enorme di classi, perché all’aumentare di K abbiamo più matrici VCV da stimare e più vettori di medie. 2. Metodi Kernel Density (kNN: k-Near Neighbours). Anche qui, però, la scalabilità rispetto a p non è grandiosa, perché anche i metodi kernel, quando aumentiamo la dimensionalità, iniziano a diventare instabili e computazional- mente complicati; 3. Naive Bayes Classifier. Tale metodo è simile al kNN, ma risolve il problema della scalabilità rispetto a p; tuttavia, si basa su una forte ipotesi di indipen- denza condizionale tra le features. 3.1 Linear Discriminant Analysis - LDA L’LDA è un mondo in cui assumiamo la Gaussianità delle f1 (·),... , fK (·). Il problema di apprendimento, quindi, diventa un problema di stima dei parametri A. LDA: in ogni classe k la distribuzione congiunta delle features X è una Nor- male Multivariata, centrata sul vettore medio µk e avente matrice VCV Σ, comune a tutte le classi (omoschedasticità); ossia: X|Y = k ∼ N (µk ; Σ) ∀ k = 1,... , K (85) Tale assunzione ha delle implicazioni: Le Class Conditional Densities sono descritte dalla funzione fk (x) = ϕ(x; µk , Σ); Poiché la Normale è una distribuzione ellittico-simmetrica, se le features sono poco rappresentabili da scatter simili occorre cambiare classificatore; Poiché c’è omoschedasticità, la dispersione delle features ha caratteristiche geometriche omogenee in tutte le classi. Sotto l’A. LDA il classificatore Bayesiano sarà dato da: Y ∗ (x) = argmax δk (x) (86) k ′ 1 ′ −1 con δk (x) = x Σ−1 µk − µ Σ µk + log(πk ) (87) 2 k La funzione δk (x) si chiama Linear Discriminant Function, e determina com’è fatta la decision boundary. Essa è una funzione lineare affine delle x, perché l’unica parte che ne dipende è il primo termine, a cui ne viene sottratta una indipendente. 35 Dimostrazione. Ricordiamo che: 1 ′ 1 Σ−1 (x−µk ) ϕ(x; µk , Σ) = p e− 2 (x−µk ) (88) det(Σ) · (2π)p Vogliamo risolvere il problema di ottimo. Scriviamo: log(πk ) + log[fk (x)] = log(πk ) + log[ϕ(x; µk , Σ)] = " # 1 ′ −1 −1 = log(πk ) + log p e 2 (x−µ k ) Σ (x−µ k ) = det(Σ) · (2π)p p 1 1 ′ = log(πk ) − log(2π) − log[det(Σ)] − (x − µk ) Σ−1 (x − µk ) 2 2 2 Fissiamo con a i termini indipendenti da da k. Inoltre, considerando la propri- età delle matrici trasposte, per cui la trasposta di una differenza è uguale alla dif- ′ ′ ′ ferenza delle trasposte, ossia (A−B) = A −B , possiamo espandere il quadrato dell’ultimo termine: 1 ′ −1 ′ ′ =... = log(πk ) + a − (x Σ x + µk Σ−1 µk − 2x Σ−1 µk ) = 2 1 ′ ′ = log(πk ) + a + b − (µk Σ−1 µk − 2x Σ−1 µk ) 2 Poiché i termini costanti sono ininfluenti nei problemi di ottimizzazione, elimi- nandoli e semplificando, il problema di ottimizzazione diventa: 1 ′ −1 ′ Y ∗ (x) := argmax [log(πk ) − µ Σ µk + x Σ−1 µk ] = argmax δk (x) (89) k 2 k k 3.1.1 Linearità della Decision Boundary Poiché l’LDA produce Decision Boundaries lineari, è un metodo di classificazione lineare, e si può dimostrare che produce decision boundaries che consistono in rette (p = 2), piani (p = 3) e iper-piani (p > 3). Fissate due classi {l, m}, il classificatore assegnerà un punto x alla classe l se δl (x) ≥ δm (x). La Decision Boundary sarà data dall’insieme: Bl, m = {x| δl (x) = δm (x)} (90) Ovviamente, la probabilità di avere due punti osservati esattamente lungo la Deci- sion Boundary è 0, quindi nella pratica il problema non si pone. Se, per problemi di arrotondamento, dovesse succedere, lo assegneremo a caso. 3.1.2 Stima dei parametri Quanti parametri dobbiamo stimare? 36  K − 1 prevalenze, π1 ,... , πK−1 ; K vettori medi di dimensione p, ossia µ1 ,... , µK ; Una matrice VCV composta da p(p+1) 2 parametri. In totale, avremo {K(p + 1) + p(p+1) 2 − 1} parametri; quindi, tale metodo scala linearmente rispetto a K e quadraticamente rispetto a p. Al solito abbiamo un train set, Xn = {(y1 , x1 ),... , (yn , xn )}, dove osservi- amo n vettori di features xi , e le yi avranno valori nell’insieme Y = {1,... , K}. Nella trattazione useremo wik = I{yi = k} per indicare il peso assegnato all’istanza i di appartenere alla classe k (stiamo bernoullianizzando il problema): ( 1 → se yi = k wik = I{yi = k} = (91) 0 → altrimenti Fissiamo un po’ di notazione. Per ogni classe possiamo riscrivere le size: X n X n nk = I{yi = k} = wik (92) i=1 i=1 ∑K Quindi, possiamo stimare i prior estraendo la media (si noti che k=1 π̂k = 1): 1X n nk π̂k = I{yi = k} = (93) n i=1 n Poiché dobbiamo stimare le medie dei gruppi condizionate ai gruppi, possiamo riscrivere i centri di classe µk come una media pesata dei vettori osservati: 1 X Xn 1 µ̂k = xi = P n wik xi (94) nk i:y =k i=1 wik i=1 i Ora possiamo stimare le Within Class Covariance Matrix. Innanzitutto, defini- amo lo stimatore della matrice VCV non distorto per ogni classe: 1 X 1 Xn Σ̂k = (xi − µ̂k )(xi − µ̂k )′ = Pn wik (xi − µ̂k )(xi − µ̂k )′ nk − 1 i:y =k ( i=1 wik ) − 1 i=1 i (95) Anche Σ̂k può essere visto come l’analogo pesato della covarianza campionaria (stats::cov.wt()). Ora, però, dobbiamo introdurre uno stimatore pooling di Σ: 1 XK Σ̂ = (nk − 1)Σ̂k (96) n − K k=1 37 Questa è la media pesata delle matrici VCV classe per classe, che darà più peso ai gruppi più numerosi (si presuppone che, all’aumentare della numerosità, si fanno stime migliori. Sostituendo i parametri della popolazione con le stime, avremo: ′ −1 1 ′ −1 δk (x) = x Σ̂ µ̂k − µ̂ Σ̂ µ̂k + log(π̂k ) (97) 2 k Massimizzando i δk (x), massimizzeremo i posterior punto per punto. Che ci as- pettiamo da questo classificatore? Quali sono le garanzie quando l’A. LDA è sod- disfatta? Cosa accade quando viene meno? E, nella pratica, quando la popolazione devia dall’A. LDA, cosa ci dobbiamo aspettare? Quando applichiamo tale metodologia, nella funzione discriminante ci stiamo mettendo le stime ML, che, per grandi campioni, convergono ai veri parametri che hanno generato la popolazione. Stiamo facendo un ragionamento ipotetico: se la popolazione fosse coerente col modello LDA, le stime sarebbero stime ML, le quali sono consistenti, asintoticamente non distorte e asintoticamente Normali. Poiché asintoticamente stimiamo i parametri ML, quando li mettiamo nella funzione discriminante (continua per questi parametri), per il teorema di continu- ità e di convergenza massimizziamo i posterior, e, quindi, otteniamo il classificatore ottimale Bayesiano e raggiungiamo il BER. Anche quando l’A. LDA non è soddisfatta ci sono risultati teorici e sperimen- tali che mostrano che l’LDA assicura risultati estremamente competitivi e ”stabili” in moltissime situazioni di non-Gaussianità. 3.2 Quadratic Discriminant Analysis - QDA La QDA non è altro che una generalizzazione dell’LDA: l’unica cosa che cambia è, che mentre, nell’LDA assumevamo l’omoschedasticità delle classi, nella QDA assumiamo che lo scatter si muove liberamente (ogni classe può avere il suo scat- ter). Questo introduce una flessibilità maggiore, il che non sempre si traduce in un vantaggio nella pratica. Ciò che cambia, quindi, è la regola di allocazione delle osservazioni, ossia la frontiera di decisione. A. QDA: in ogni classe k la distribuzione congiunta delle features X è una Nor- male Multivariata centrata sul vettore medio µk e con matrice VCV Σk ; ossia: X|Y = k ∼ N (µk ; Σk ) ∀ k = 1,... , K (98) L’implementazione del classificatore ottimale Bayesiano è dato da: Y ∗ (x) := argmax δk (x) (99) k 38 Qui δk (x) si chiama Funzione Discriminante Quadratica, e avrà questa forma: 1 1 δk (x) = log(πk ) − log[det(Σk )] − (x − µk )′ Σ−1 k (x − µk ) (100) 2 2 Essa descriverà delle forme quadriche, paraboloidi (o iper-paraboloidi se p > 3). Dimostrazione. Riscriviamo log(πk ) + log[fk (x)] usando l’A. QDA: =... = log(πk ) + log[ϕ(x; µk , Σk )] = p 1 1 = log(πk ) − log(2π) − log[det(Σk )] − (x − µk )′ Σ−1 k (x − µk ) 2 2 2 Trascurando i termini che non dipendono da k: 1 1 δk (x) = log(πk ) − log[det(Σk )] − (x − µk )′ Σ−1 k (x − µk ) (101) 2 2 Quindi, il problema argmax {log(πk )+log[ϕ(x; µk , Σk )]} equivale a argmax δk (x). k k Analogamente all’LDA, la Decision Boundary tra le due classi l ed m è l’insieme: Bl, m = {x | δl (x) = δm (x)} (102) Quanti parametri abbiamo da stimare? K − 1 prevalenze, π1 ,... , πK−1 ; K vettori medi di dimensione p, ossia µ1 ,... , µK ; K matrici VCV, composte da p(p+1) 2 parametri 2 In totale, avremo {K p +3p+2 2 − 1} parametri. Anche tale metodo scala linear- mente rispetto a K e quadraticamente rispetto a p. La stima dei πk e µk saranno, anche qui, π̂k e µ̂k. La stima delle Σk è data dalle Within-Classs Covariance Matrix Σ̂k. Sostituendo i parametri della popolazione con le stime, otteniamo: 1 1 −1 δk (x) = log(π̂k ) − log[det(Σ̂k )] − (x − µ̂k )′ Σ̂k (x − µ̂k ) (103) 2 2 Massimizzando la Funzione Discriminante Quadratica per ogni punto del train set, otteniamo la classificazione, con garanzie simili a quelle dell’LDA. 39 3.3 k-Nearest Neigbours (kNN) Si parla di metodi NN perché, in base alla stessa logica, si possono costruire metodi per stimare densità multivariate, metodi di regressione basati sul kNN, e così via. I metodi kNN non fanno alcuna assunzione di tipo parametrico sulla Popolazione; tuttavia, anche essi partono dalla stima delle densità di classe condizionate per ot- tenere una stima dei posterior, ma lo fanno in modo non parametrico (senza as- sumere una precisa forma distribuzionale). Attenzione: qui k non denota una classe, ma una cosa che si chiama ordine del neighbour. Consideriamo dei vettori osservati di features {x1 ,... , xn }, e consideriamo un qualunque punto x nello spazio delle features. Scelta la distanza desiderata (eu- clidea, di Manhattan e qualsiasi altra), chiameremo di (x) = dist(xi , x), ossia la distanza fra xi e x. Indichiamo con {(1),... , (n)} i ranks delle distanze ordi- nante in senso non-decrescente in modo che di (x) ≤ d(i+1) (x) ∀ i = 1, 2,... Il kNN è l’insieme delle k istanze del training set più vicine a x, ovverosia: Nk (x) = {xi |di (x) ≤ d(k) (x)} (104) La definizione del kNN dipende da due tunings (iper-parametri): k determina quanto localmente guardiamo intorno ad x (anche bandwidth, perché ha lo stesso ruolo della bandwidth nel kernel density); dist(·) determina la geometria dell’intorno. 3.3.1 Classificatore kNN Il kNN si basa su un modo di metrizzare lo spazio indipendente dallo scaling delle distanze tra i punti. Il kNN determina una nozione di neighbouring depurata dalla scala numerica effettiva che governa le distanze degli oggetti considerati. Ci sono 2 iper-parametri che governano il kNN, ossia k e il tipo di distanza considerata. Per evitare confusioni, consideriamo L classi, Y = {1,... , L}. Sia Xn = {(y1 , x1 ),... , (yn , xn )} il train set, e x è un punto nello spazio delle features da assegnare ad una delle L classi. Il kNN usa la distribuzione congiunta delle features delle classi osservate per stimare i posterior, ma lo fa in modo non parametrico. Chiamiamo ηl (x) il poste- rior della classe l quando l’input è x: ηl (x) = P (Y = l|X = x) (105) 40 Si può dimostrare che una stima statisticamente valida di questo posterior può es- sere calcolata con uno stimatore di classe kernel: 1X n η̂l (x) = I{xi ∈ Nk (x)}I{yi = l} (106) k i=1 Il fattore I{yi = l} varrà 1 se l’istanza i è di classe l. Invece, il termine I{xi ∈ Nk (x)} varrà 1 se xi sta nel neighbour di x di dimensione k. Quindi, abbiamo una somma di 0 e 1. Quindi, η̂l (x) calcola la proporzione di osservazioni nel train set di classe l che sta nel k-mo neighbour e, quindi, la stima di questa probabilità. Nella pratica, il kNN implementa il classificatore ottimale Bayesiano, perché, se prendiamo un input x, guardiamo al suo neighbour e stimiamo la probabilità a posteriori di trovare nel neighbour osservazioni di classe l, attribuiremo x alla classe che massimizza il posterior, quindi Ŷ (x) = argmax η̂l (x). l Il kNN è un classificatore di majority rule: stima il posterior classe per classe e attribuisce il punto x alla classe che lo massimizza. Se quella quantità la molti- plichiamo per una costante indipendente da l, il problema resterebbe invariato; quindi, l’equazione equivale a Ŷ (x) = argmax k η̂l (x). l Se guardiamo a cosa corrisponde k η̂l (x), essa non è altro che una somma di indicatori (sono conteggi). Quella sommatoria, quindi, per ogni classe k conta quanti elementi di classe l stanno in un neighbour di ordine k del punto x. In termini di geometria, il kNN separa le classi in un modo che consente di riprodurre ogni forma di frontiera di decisione. Tutto dipende da k: se è troppo piccolo, si ha undersmoothing (o overfitting - si ha una stima di densità dei dati troppo locale); viceversa, si omogenizzano regioni troppo diverse rischiando di mischiare le classi. Il valore k1 esprime la complessità con cui rappresentiamo i dati. La performance del kNN dipende molto da k, e non ci sono molti casi dove abbiamo garanzie analitiche. Per k = 1, nel problema 2-class si dimostra che asin- toticamente l’error rate del kNN al massimo equivale al doppio del Bayes Error. 41 42  4 Metriche di Performance e Metodi di Thresholding 4.1 Alcune misure di performance Spesso l’MCR non è una misura di nostro interesse. Facciamo un esempio. Classi Previste Classi Vere Negativi Positivi Totale Negativi 88 6 94 Positivi 4 2 6 Tot 92 8 100 Tale classificatore classifica correttamente 88 negativi su 94 e 2 positivi su 6. Si noti che c’è un fortissimo sbilanciamento tra le classi (i negativi sono il 94%). Sebbene il classificatore abbia un’accuratezza del 90%, possiamo assicurarcene una del 94% semplicemente assegnando ogni istanza alla classe prevalente. Questo è il paradosso dell’accuratezza: il classificatore triviale assicura un MCR solo del 6%. C’è un altro motivo per cui l’MCR potrebbe non essere un buon indicatore delle performance di un classificatore. Ad esempio, quando operiamo nel caso diagnostico, un falso positivo e un falso negativo hanno due importanze molto differenti. Si ricordi che l’MCR usa una funzione di perdita 0/1, che assegna 1 in- dipendentemente dall’errore commesso (non vorremmo somministrare farmaci ai sani e trascurare gli ammalati! Se dobbiamo scegliere, è più grave il falso negativo). Ci sono molte misure che cercano di correggere il difetto dell’MCR: F P R = FNP. È la frazione dei positivi misclassificati, e corrisponde alla probabilità di