Comments
Description
Transcript
COMPONENTI PRINCIPALI
COMPONENTI PRINCIPALI Mario Romanazzi 1 Introduzione Le componenti principali sono combinazioni lineari delle variabili osservate, ordinate in base ad un criterio di rilevanza informativa. La prima componente principale estrae dai dati la massima quantità d’informazione; le componenti successive (seconda, terza, ecc.) ottimizzano l’informazione residua sotto il vincolo d’incorrelazione con le altre componenti. In generale, avendo osservato p > 1 variabili X 1 , ..., Xp si possono determinare altrettante componenti principali Z1 , ..., Zp , a due a due linearmente indipendenti, tali che le prime 1 ≤ q < p componenti Z1 , ..., Zq spiegano una parte significativa dell’informazione contenuta in X1 , ..., Xp . L’applicazione più importante è la riduzione di dimensionalità dei dati osservati. Infatti la loro proiezione nello spazio delle prime q componenti consente di studiare le proprietà statistiche della distribuzione in uno spazio di dimensione ridotta con una perdita (sperabilmente) limitata d’informazione. Un’altra applicazione è nei modelli di regressione con variabili esplicative fortemente correlate (o addirittura collineari). In tal caso le stime dei parametri prodotte dai metodi standard (minimi quadrati, massima verosimiglianza sotto l’ipotesi di normalità) sono instabili o addirittura non esistono. Un possibile rimedio è la sostituzione delle variabili esplicative X1 , ..., Xp del modello di regressione con le corrispondenti componenti principali Z1 , ..., Zp (o un loro sottoinsieme). Il metodo delle componenti principali è illustrato in tutti i libri di analisi multivariata. Un riferimento autorevole è Mardia et al (1979), Cap. 8. L’argomento presuppone la conoscenza delle nozioni base di algebra delle matrici di cui forniamo un veloce riepilogo nell’Appendice A. L’Appendice B espone la definizione e le principali proprietà di autovalori e autovettori. Per ulteriori approfondimenti il lettore può consultare l’Appendice A del già citato Mardia et al (1979). Una più ampia esposizione particolarmente mirata alle applicazioni di Statistica ed Econometria è Magnus and Neudecker (2001). Alcuni esempi numerici sono sviluppati col programma R (R Development Core Team, 2006). 2 Notazione e risultati preliminari Supponiamo di aver rilevato p variabili numeriche X1 , ..., Xp su un campione di n unità. I dati sono ordinati nella matrice nxp X = (xij ); xij rappresenta il dato della i-esima unità sulla j-esima variabile. Elenchiamo di seguito le statistiche di sintesi più comuni. Pn Il vettore delle medie (centroide della distribuzione) x̄ = (x̄1 , ..., x̄j , ..., xp )T , con x̄j = n−1 i=1 xij . Un’espressione matriciale di x̄ è x̄ = n−1 X T 1n , in cui 1n è il vettore nx1 con tutte le componenti uguali a 1. Il centroide descrive la posizione del campione nello spazio p-dimensionale. La matrice di varianze e covarianze SX = (shk ), quadrata d’ordine p, in cui i termini diagonali Pn shh P = n−1 i=1 (xih − x̄h )2 = var(Xh ) sono le varianze mentre i termini non diagonali shk = n n−1 i=1 (xih − x̄h )(xik − x̄k ) = cov(Xh , Xk ) sono le covarianze della distribuzione. Ogni matrice di covarianza è simmetrica e semi-definita positiva. Essa descrive la dispersione del campione nello spazio p-dimensionale; più precisamente, le deviazioni standard (radici quadrate delle varianze) 1 2 2 NOTAZIONE E RISULTATI PRELIMINARI misurano la dispersione delle variabili attorno alle rispettive medie mentre le covarianze misurano l’interdipendenza lineare delle variabili. Un’espressione matriciale di SX è SX = n−1 n X (xi − x̄)(xi − x̄)T = n−1 XCT XC , i=1 (C) in cui xi indica l’i-esima riga di X (scritta come vettore colonna px1) e XC = X − 1n x̄T = (xij ) (C) indica la matrice dei dati centrati rispetto alle rispettive medie, xij = xij − x̄j . La matrice di correlazione RX = (rhk ), quadrata d’ordine p, in cui i termini diagonali sono identicamente uguali a 1 mentre i termini non diagonali rhk = shk /(shh skk )1/2 sono i coefficienti di correlazione lineare tra le variabili. Un’importante proprietà è −1 ≤ r hk ≤ 1; rhk = 0 nel caso di indipendenza lineare, rhk = ±1 nel caso di massima dipendenza lineare; il segno positivo o negativo di rhk descrive il verso della relazione lineare (concorde o discorde). La matrice di correlazione è la (ST ) matrice di covarianza dei dati standardizzati xij = (xij − x̄j )/(sjj )1/2 ed è pertanto simmetrica e semi-definita positiva. Un’espressione matriciale di RX è T RX = n−1 XST XST , in cui XST = XC (diag SX )−1/2 = (xST ij ) è la matrice dei dati standardizzati; diag SX indica la matrice i cui elementi diagonali sono quelli di SX e quelli non diagonali sono identicamente nulli. Esempio 1 (Indicatori delle economie dell’Unione Europea) Dall’edizione 2009 di Italia in cifre (ISTAT) abbiamo ripreso i seguenti indicatori delle economie dei 27 Paesi dell’Unione Europea: PIL pro-capite a parità di potere d’acquisto (UE27 = 100), deficit pubblico (% sul PIL), debito pubblico (% sul PIL), tasso d’inflazione, tasso d’occupazione. I dati sono riferiti al 2007. Riportiamo di seguito la matrice dati. Paese PIL Deficit Debito Inflazione Occupazione 1 AUSTRIA 123.8 -0.4 59.5 2.2 71.4 2 BELGIO 118.1 -0.3 83.9 1.8 62.0 3 CIPRO 90.6 3.5 59.5 2.2 71.0 4 FINLANDIA 115.8 5.3 35.1 1.6 70.3 5 FRANCIA 109.1 -2.7 63.9 1.6 70.3 6 GERMANIA 114.7 -0.2 65.1 2.3 69.4 7 GRECIA 94.8 -3.5 94.8 3.0 61.4 8 IRLANDA 150.2 0.2 24.8 2.9 69.1 9 ITALIA 101.4 -1.5 103.5 2.0 58.7 10 LUSSEMBURGO 266.2 3.2 7.0 2.7 64.2 11 MALTA 77.7 -1.8 62.2 0.7 54.6 12 PAESI_BASSI 130.9 0.3 45.7 1.6 76.0 13 PORTOGALLO 76.1 -2.6 63.6 2.4 67.8 14 SLOVACCHIA 67.0 -1.9 29.4 1.9 60.7 15 SLOVENIA 89.2 0.5 23.4 3.8 67.8 16 SPAGNA 105.4 2.2 36.2 2.8 65.6 17 BULGARIA 37.3 0.1 18.2 7.6 61.7 18 DANIMARCA 120.0 4.9 26.2 1.7 77.1 19 ESTONIA 67.9 2.7 3.5 6.7 69.4 20 LETTONIA 57.9 0.1 9.5 10.1 68.3 21 LITUANIA 59.5 -1.2 17.0 5.8 64.9 22 POLONIA 53.3 -2.0 44.9 2.6 57.0 3 3 COMBINAZIONI LINEARI 23 REGNO_UNITO 119.1 24 REP_CECA 80.2 25 ROMANIA 42.1 26 SVEZIA 122.2 27 UNGHERIA 62.6 -2.8 -1.0 -2.6 3.6 -5.0 44.2 28.9 12.9 40.4 65.8 2.3 3.0 4.9 1.7 7.9 71.5 66.1 58.8 74.2 57.3 Otteniamo il vettore delle medie, la matrice di varianze e covarianze e la matrice di correlazione con le funzioni R mean, cov e cor. > mean(eu[, 4:8]) PIL Deficit Debito 98.2629630 -0.1074074 43.3000000 Inflazione Occupazione 3.3259259 66.1703704 > round(cov(eu[, 4:8]), 2) PIL Deficit Debito Inflazione Occupazione PIL 1993.66 48.80 3.21 -47.41 95.99 Deficit 48.80 6.98 -28.53 -1.03 9.13 Debito 3.21 -28.53 710.81 -27.39 -35.64 Inflazione -47.41 -1.03 -27.39 5.37 -2.65 Occupazione 95.99 9.13 -35.64 -2.65 35.99 > round(cor(eu[, 4:8]), 2) PIL Deficit Debito Inflazione Occupazione 3 PIL Deficit Debito Inflazione Occupazione 1.00 0.41 0.00 -0.46 0.36 0.41 1.00 -0.41 -0.17 0.58 0.00 -0.41 1.00 -0.44 -0.22 -0.46 -0.17 -0.44 1.00 -0.19 0.36 0.58 -0.22 -0.19 1.00 Combinazioni lineari Una combinazione lineare generalizza il familiare concetto di media aritmetica ponderata. Essa trasforma il vettore delle determinazioni x1 , ..., xj , ..., xp di p variabili numeriche in un valore numerico unico, indicato con z, pari alla somma delle determinazioni numeriche x1 , ..., xj , ..., xp moltiplicate per i coefficienti a1 , ..., aj , ..., ap z = a1 x1 + ... + ap xp . I coefficienti possono assumere valori numerici arbitrari (non tutti nulli). In termini matematici, una Pp combinazione lineare è una funzione t : Rp → R tale che, per ogni vettore x di Rp , t(x) = j=1 aj xj = aT x. La famiglia delle combinazioni lineari è molto ampia e comprende trasformazioni largamente usate nell’analisi dei dati. Somma. La combinazione lineare più comune è la somma, corrispondente ad a = 1p , il vettore con p componenti identicamente uguali a 1. Supponiamo di aver rilevato su un campione di n famiglie le spese mensili divise per capitolo di spesa: abbigliamento, alimentazione, casa, istruzione e cultura, trasporti, tempo libero. La P6spesa totale, pari alla somma delle voci di ogni capitolo, è la combinazione lineare z = 1T6 x = j=1 xj . 4 3 COMBINAZIONI LINEARI Differenza. Il risultato di bilancio di un’impresa, differenza fra il totale delle entrate ed il totale delle uscite, è una combinazione lineare con vettore a = (1, −1)T . Media. Un altro familiare esempio è la media, corrispondente ad a = 1p /p. Quando gli studenti calcolano il voto medio negli esami applicano esattamente questa trasformazione. Esempio 2 (Voto medio degli esami) In un campione di 16 laureati in Economia abbiamo rilevato i voti nelle seguenti discipline: Matematica, Statistica, Economia Politica I e II, Economia Aziendale I e II, Diritto I e II. La tabella sottostante riporta i voti nelle singole discipline insieme col corrispondente voto medio. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 stu mat stat eco1 eco2 eca1 eca2 dir1 dir2 media S1 24 26 30 30 28 27 25 30 27.500 S2 18 21 26 28 30 24 25 19 23.875 S3 26 22 30 28 27 25 22 21 25.125 S4 20 20 21 26 26 23 23 27 23.250 S5 22 23 25 27 30 28 25 28 26.000 S6 21 18 20 19 23 26 21 22 21.250 S7 22 23 28 26 24 28 25 23 24.875 S8 26 24 22 18 24 20 21 20 21.875 S9 18 19 23 25 22 19 20 24 21.250 S10 23 21 24 24 28 23 30 28 25.125 S11 26 25 28 27 29 30 25 28 27.250 S12 24 27 25 23 29 25 25 27 25.625 S13 18 18 21 25 22 19 20 24 20.875 S14 20 20 24 23 26 22 27 23 23.125 S15 26 21 24 21 22 18 19 23 21.750 S16 26 24 22 25 24 23 18 22 23.000 Vale la pena sottolineare la differenza tra la combinazione lineare ”voto medio”, che definisce un’ulteriore variabile (funzione di quelle osservate), e la media campionaria dei voti nelle singole discipline. È ovviamente possibile calcolare la media campionaria anche del ”voto medio”, insieme con le altre usuali statistiche riassuntive. Le medie e le deviazioni standard sono riportate di seguito. [1] Medie mat stat eco1 eco2 eca1 eca2 dir1 dir2 media 22.50 22.00 24.56 24.69 25.88 23.75 23.19 24.31 23.86 [1] Deviazioni Standard mat 3.08 stat 2.73 eco1 3.14 eco2 3.28 eca1 2.92 eca2 3.57 dir1 3.23 dir2 media 3.28 2.15 La distribuzione di una combinazione lineare dipende dalla distribuzione congiunta delle variabili osservate. Ad esempio, se questa è normale, allora ogni combinazione lineare ha una distribuzione normale con media e deviazione standard dipendenti dal vettore delle medie e dalla matrice di varianze e covarianze di partenza. Nel caso campionario, la media e la deviazione standard di ogni combinazione lineare sono semplici funzioni matriciali di x̄ e SX . 5 4 COMPONENTI PRINCIPALI Teorema 3 (Media e varianza delle combinazioni lineari) Per ogni combinazione lineare z = a T x, z̄ = aT x̄, s2Z = aT SX a. Dimostrazione. Sia X la matrice nxp dei dati campionari e sia a il vettore px1 dei coefficienti della combinazione lineare. Il vettore z nx1 dei valori campionari della combinazione lineare è il prodotto (righe per colonne) z = Xa. La media campionaria di z è z̄ = n−1 =n n X zi = n−1 z T 1n i=1 −1 T T a X 1n = aT (n−1 X T 1n ) = aT x̄. La varianza campionaria di z è s2Z = n−1 n X (zi − z̄)2 = n−1 i=1 = aT (n−1 n X aT (xi − x̄)(xi − x̄)T a i=1 n X (xi − x̄)(xi − x̄)T )a = aT SX a. i=1 s2Z T Notiamo che = a SX a ≥ 0 per ogni vettore a il che prova che la matrice di varianze e covarianze S X è semi-definita positiva. Nel caso particolare p = 2, il teorema precedente fornisce le formule, note dai corsi di base z̄ = a1 x̄1 + a2 x̄2 , s2Z 4 = a21 s11 + a22 s22 + 2a1 a2 s12 . Componenti principali Negli esempi della Sezione 4, i coefficienti delle combinazioni lineari sono valori stabiliti in anticipo, in base al problema in esame: di volta in volta si considera la somma, la differenza oppure la media. Un caso completamente diverso è quello in cui viene fissata qualche proprietà della combinazione lineare z ≡ z a e il vettore a dei coefficienti è determinato in modo da soddisfare tali requisiti. L’esempio più importante è proprio la (prima) componente principale, definita come la combinazione lineare normalizzata con massima varianza. Il significato di tale requisito è identificare la combinazione lineare, cioè la proiezione da Rp in R, capace di assorbire la massima quantità d’informazione (misurata dalla varianza) contenuta nei dati osservati. Una combinazione lineare viene definita normalizzata se il vettore dei coefficienti soddisfa il vincolo p X T a a= a2j = 1. j=1 Questa restrizione è necessaria per assicurare l’esistenza e l’unicità della soluzione del problema di massimizzazione della varianza. In base al Teorema 3, dobbiamo determinare il vettore a tale che la funzione var(z a ) = aT SX a ≡ ψ(a) assuma il suo massimo valore, sotto il vincolo aT a = 1. La matrice SX , simmetrica e semi-definita positiva, ha p autovalori non negativi λ1 ≥ λ2 ≥ ... ≥ λp ≥ 0 associati ai corrispondenti autovettori a1 , a2 , ..., ap normalizzati (aTh ah = 1) e a due a due ortogonali (aTh ak = 0, h 6= k). La generica coppia autovalore-autovettore verifica l’equazione matriciale SX a h = λ h a h , h = 1, ..., p. Il numero degli autovalori positivi è uguale al rango di X e di SX . Nel seguente teorema dimostriamo che il vettore a ottimale è l’autovettore a1 di SX associato al massimo autovalore λ1 . 5 PROPRIETÀ STATISTICHE DELLE COMPONENTI PRINCIPALI 6 Teorema 4 (Prima componente principale) Il valore massimo della funzione ψ(a) = a T SX a, sotto il vincolo aT a = 1, è pari al massimo autovalore λ1 di SX e viene raggiunto quando il vettore a coincide col corrispondente autovettore a1 . Dimostrazione. La dimostrazione è basata sulla disuguaglianza di Cauchy-Schwarz secondo la quale, per ogni coppia di vettori u = (u1 , u2 , ..., up )T , v = (v1 , v2 , ..., vp )T di Rp (uT v)2 ≤ (uT u)(v T v). L’uguaglianza vale se e solo se v = λu per qualche scalare λ. Nella funzione ψ(a) poniamo u = a, v = SX a. In base alla disuguaglianza di Cauchy-Schwarz 2 2 a) = (aT SX a) ψ(a)2 ≤ (aT a)(aT SX perchè aT a = 1. L’uguaglianza è soddisfatta se e solo se SX a = λa, cioè se e solo se a è un autovettore di SX . Premoltiplicando l’equazione per aT otteniamo aT SX a = ψ(a) = λaT a = λ, pertanto la funzione ψ(a) viene massimizzata se a è l’autovettore di S X associato al massimo autovalore. Nota 5 Dalla dimostrazione segue che gli autovalori della matrice di covarianza sono le varianze delle componenti principali. In particolare, il massimo autovalore è la varianza della prima, e più importante, componente. Nota 6 Ai p autovettori ah , h = 1, ..., p, della matrice di covarianza corrispondono altrettante componenti principali zh , combinazioni lineari delle variabili originarie, ordinate secondo la varianza (autovalore λ h ) decrescente. Nota 7 Le componenti principali associate agli autovalori più piccoli sono usate per individuare le variabili osservate ridondanti. Un criterio frequentemente usato è di eliminare le variabili osservate aventi massima correlazione assoluta con le ultime s componenti,1 ≤ s ≤ p − 1. Nota 8 Se la matrice dati X ha rango k < p, anche il rango di SX è pari a k e gli ultimi p − k autovalori sono nulli. I corrispondenti autovettori descrivono le relazioni di collinearità tra le variabili osservate. 5 Proprietà statistiche delle componenti principali L’analisi delle componenti principali (ACP) consiste nella determinazione di q < p componenti, basate sui primi q autovettori della matrice di covarianza dei dati osservati. Il valore di q viene scelto in modo da assorbire una quota adeguata d’informazione. Le componenti sono usate per studiare le proprietà del campione in uno spazio di dimensione ridotta rispetto a quello di partenza. Il problema fondamentale è l’interpretazione delle componenti che sono di per sè variabili artificiali, combinazioni lineari di quelle d’origine. L’estrazione delle componenti principali è sempre preceduta da qualche trasformazione dei dati osservati. I dati sono sempre centrati rispetto al vettore delle medie. Inoltre, in molti casi, i dati vengono standardizzati per evitare possibili distorsioni dovute a sperequazioni delle varianze delle variabili d’origine. In questo secondo caso, poichè la matrice di covarianza dei dati standardizzati coincide con la matrice di correlazione dei dati di partenza, le componenti principali sono basate sugli autovalori ed autovettori della matrice di correlazione. Nel seguente teorema ricaviamo medie, varianze e covarianze delle componenti principali. Indichiamo con zh una generica componente. Se i dati sono centrati, zh = XC ah ; se i dati sono standardizzati, zh = XST bh , in cui bh indica l’autovettore della matrice di correlazione associato al suo h-esimo più alto autovalore µh . 5 PROPRIETÀ STATISTICHE DELLE COMPONENTI PRINCIPALI 7 Teorema 9 (Proprietà statistiche) i. Le componenti principali hanno medie nulle z̄h = 0. ii. Se i dati sono centrati, le varianze delle componenti principali sono uguali agli autovalori della matrice di covarianza var(zh ) = λh ; se i dati sono standardizzati, sono uguali agli autovalori della matrice di correlazione var(zh ) = µh . iii. Le componenti principali sono a due a due incorrelate cov(zh , zk ) = 0. Dimostrazione. i. In base al Teorema 3, z̄h = aTh x̄C = 0, perchè x̄C = 0p , il vettore nullo. La dimostrazione nel caso di dati standardizzati è simile. ii. Se i dati sono centrati, la varianza di zh è var(zh ) = aTh SX ah = λh aTh ah = λh . Se i dati sono standardizzati, var(zh ) = bTh RX bh = µh bTh bh = µh . iii. Usando la definizione di covarianza cov(zh , zk ) = n −1 = n−1 n X i=1 n X (zih − z̄h )(zik − z̄k ) zih zik = n−1 zhT zk , i=1 perchè z̄h = z̄k = 0. Se i dati sono centrati, cov(zh , zk ) = n−1 zhT zk = n−1 aTh XCT XC ak = aTh SX ak = λk aTh ak = 0 perchè autovettori distinti sono ortogonali. La dimostrazione nel caso di dati standardizzati è simile. La qualità statistica delle componenti principali si misura mediante la percentuale spiegata della varianza totale. La varianza totale è la traccia della matrice di covarianza (dati centrati) o di correlazione (dati standardizzati) tr(SX ) = tr(RX ) = p X h=1 p X h=1 shh = rhh = p X h=1 p X h=1 λh µh = p 8 6 APPLICAZIONI Supponiamo di considerare tutte le p componenti. In tal caso, la varianza complessivamente assorbita è la varianza totale, tr(SX ) o tr(RX ) = p ce rappresenta pertanto la massima quantità d’informazione Pq estraibile dai dati. Se si considerano le prime q < p componenti, la varianza assorbita è h=1 λh e un indice relativo di qualità è q X λh /tr(SX ). h=1 La formula corrispondente per dati standardizzati è q X µh /tr(RX ) = h=1 q X µh /p. h=1 Per una buona approssimazione il rapporto non dovrebbe essere inferiore al 70%. Per dati Pp standardizzati, un criterio frequentemente usato è di considerare tutti gli autovalori maggiori di 1 = h=1 µh /p, media campionaria degli autovalori. 6 Applicazioni La funzione R che determina le componenti principali di un campione numerico è princomp; l’opzione cor=FALSE centra i dati rispetto alle medie delle variabili mentre cor=TRUE li standardizza . Essa produce una lista contenente tre risultati principali: 1. $loadings, matrice pxp le cui colonne contengono gli autovettori ordinati in modo decrescente, secondo i corrispondenti autovalori, 2. $scores, matrice nxp delle componenti principali ordinate, 3. $sdev, vettore px1 delle deviazioni standard delle componenti principali (uguali alle radici quadrate degli autovalori). La funzione summary, applicata ai risultati di princomp, fornisce la % della varianza spiegata dalle prime q componenti, q = 1, 2, ..., p. Per l’interpretazione delle componenti, sono utili le correlazioni lineari tra le variabili osservate e le componenti principali, ottenibili mediante la funzione cor. 6.1 Indicatori economici dell’Unione Europea Illustriamo di seguito l’ACP degli indicatori economici dell’Unione Europea, considerando prima i dati centrati. Come sappiamo, essa è basata sugli autovalori ed autovettori della matrice di covarianza. > pc <- princomp(eu[, 4:8], cor = FALSE) > summary(pc) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Standard deviation 43.8930764 26.2380377 5.42881267 1.948002663 Proportion of Variance 0.7267847 0.2597027 0.01111792 0.001431505 Cumulative Proportion 0.7267847 0.9864874 0.99760535 0.999036855 Comp.5 Standard deviation 1.5978605134 Proportion of Variance 0.0009631446 Cumulative Proportion 1.0000000000 > round(cor(eu[, 4:8], pc$scores[, 1:5]), 2) 6 APPLICAZIONI 9 Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 PIL 1.00 0.00 0.01 0.00 0.00 Deficit 0.42 0.41 -0.42 -0.60 -0.34 Debito 0.00 -1.00 -0.01 0.00 0.00 Inflazione -0.46 0.44 0.17 0.48 -0.58 Occupazione 0.36 0.23 -0.90 0.07 0.02 La prima componente spiega da sola circa il 73% della varianza totale e le prime due insieme quasi il 99%. La matrice di correlazione delle variabili osservate con le componenti principali mostra che la prima componente ha correlazione pari a 1 con PIL mentre la seconda ha correlazione pari a −1 con Debito. Questo indica che le prime due componenti riproducono pari pari due variabili osservate, PIL e Debito rispettivamente, senza operare una vera sintesi multidimensionale. La ragione risiede nel valore delle varianze di PIL e Debito, molto più alto delle altre variabili. Questo implica, tra l’altro, che il diagramma di dispersione delle prime due componenti è equivalente a quello di PIL e Debito. Il problema può essere superato mediante una standardizzazione preliminare dei dati (le variabili standardizzate hanno tutte varianza unitaria). > pc1 <- princomp(eu[, 4:8], cor = TRUE) > summary(pc1) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 1.4599323 1.2568815 0.7605521 0.64256886 0.54544695 Proportion of Variance 0.4262804 0.3159502 0.1156879 0.08257895 0.05950247 Cumulative Proportion 0.4262804 0.7422307 0.8579186 0.94049753 1.00000000 > round(cor(eu[, 4:8], pc1$scores[, 1:5]), 2) Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 PIL -0.73 -0.33 0.50 -0.31 -0.06 Deficit -0.83 0.29 -0.02 0.36 -0.30 Debito 0.28 -0.86 -0.25 -0.09 -0.33 Inflazione 0.45 0.79 0.04 -0.30 -0.29 Occupazione -0.79 0.15 -0.51 -0.30 0.09 L’ACP dei dati standardizzati produce risultati molto diversi. La % della varianza spiegata dalla prima componente è 43%, dalle prime due 74% e dalle prime tre 86% (valori approssimati). Inoltre, l’esame delle correlazioni con le variabili osservate mostra che 1. La prima componente dipende soprattutto da Deficit, Occupazione e PIL (coefficienti negativi, variabili tra −0.73 e −0.83). Pertanto, Paesi con valori elevati di queste variabili si posizionano sul polo negativo estremo della prima componente e, all’opposto, Paesi con valori bassi si posizionano sul polo positivo estremo. 2. La seconda componente è essenzialmente un contrasto tra Inflazione e Debito (coefficienti rispettivamente uguali a 0.79 e −0.86). Il polo positivo caratterizza Paesi con alta inflazione e basso debito pubblico, quello negativo Paesi con bassa inflazione e alto debito. 3. La terza componente è un contrasto tra PIL e Occupazione (coefficienti rispettivamente uguali a 0.51 e −0.51). 10 6 APPLICAZIONI Alla luce di queste osservazioni, analizziamo il diagramma di dispersione delle prime due componenti. > > > > + > > + n <- dim(eu)[1] col <- rep("black", n) col[eu$Euro == 1] <- "red" plot(pc1$scores[, 1:2], pch = 20, col = col, xlab = "PC1 (43%)", ylab = "PC2 (32%)", main = "Indicatori Economici EU") abline(h = 0, v = 0, col = "red", lty = "dashed") text(pc1$scores[, 1:2], labels = as.character(eu$Sigla), pos = 4, cex = 0.6) Indicatori Economici EU 3 LE ES 2 BU LA 1 SL DE CZ FI LU 0 PC2 (32%) RO IR HU SP SW RS CY PL NE UK −1 AU PO GE FR MA BE −2 GR IT −3 −2 −1 0 1 2 3 PC1 (43%) Ai poli opposti di PC1 troviamo Ungheria (deficit negativo, basso tasso d’occupazione, basso PIL pro capite) e Lussemburgo con Danimarca (alti tassi d’occupazione, specie per Danimarca, e PIL pro capite superiore alla media UE, specie per Lussemburgo). Sul polo negativo di PC2 troviamo Italia, Grecia e Belgio (debito pubblico superiore alla media UE), su quello positivo Lettonia (inflazione molto alta, debito pubblico basso). In generale, i Paesi da più tempo nell’Unione (e quasi tutti quelli aderenti all’euro) occupano la regione triangolare in basso a sinistra, mentre i Paesi di recente adesione occupano la regione in alto a destra. Ulteriore informazione viene fornita dal diagramma di dispersione della prima e della terza componente. 11 6 APPLICAZIONI Indicatori Economici EU 1 PC3 (12%) 2 3 LU MA IR BE SP PL RS RO HU 0 IT SL FI LA CZ BU GR LE ESUK AU GE FR SW NE PO DE −1 CY −3 −2 −1 0 1 2 3 PC1 (43%) La terza componente è dominata da Lussemburgo, in posizione chiaramente periferica, a causa del valore elevatissimo del PIL pro capite (266.2, il valore immediatamente successivo è quello di Irlanda, 150.2; la media UE è 100). Per valutare l’impatto di questo dato anomalo sull’ACP, ripetiamo l’analisi escludendo sia Lussemburgo che Lettonia che ha un dato estremo sul tasso d’inflazione (10.1%). > > > > eu1 <- eu[eu$PIL < 200 & eu$Inflazione < 10, ] col1 <- col[eu$PIL < 200 & eu$Inflazione < 10] pc2 <- princomp(eu1[, 4:8], cor = TRUE) summary(pc2) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 1.5538176 1.2743333 0.67933722 0.57687409 0.40919765 Proportion of Variance 0.4828698 0.3247851 0.09229981 0.06655674 0.03348854 Cumulative Proportion 0.4828698 0.8076549 0.89995471 0.96651146 1.00000000 > round(cor(eu1[, 4:8], pc2$scores[, 1:5]), 2) 6 APPLICAZIONI PIL Deficit Debito Inflazione Occupazione 12 Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 -0.88 -0.28 -0.25 0.00 0.29 -0.68 0.54 0.36 -0.34 0.02 -0.05 -0.92 -0.07 -0.35 -0.12 0.68 0.53 -0.39 -0.30 0.08 -0.84 0.35 -0.34 0.05 -0.25 I risultati mostrano che le prime due componenti arrivano a spiegare circa l’80% della varianza totale. Anche l’interpretazione, basata sulle correlazioni con le variabili osservate, è un po’ diversa. La prima componente è un contrasto tra Inflazione da un lato (coefficiente 0.68) e PIL, Deficit, Occupazione dall’altro (coefficienti −0.88, −0.68, −0.84). La seconda componente è un contrasto tra Deficit e Inflazione da un lato e Debito, dall’altro. Riportiamo di seguito il diagramma di dispersione, con Lussemburgo e Lettonia come punti supplementari. La loro posizione è determinata usando la soluzione ACP trovata per gli altri 25 Paesi, senza che essi abbiano contribuito a determinare autovalori ed autovettori (la logica è quella del campione di prova nell’analisi di regressione). Va sottolineato che i loro dati devono essere preventivamente standardizzati, usando medie e deviazioni standard degli altri 25 Paesi. > ps <- eu[eu$Sigla == "LU" | eu$Sigla == "LE", ] > ps1 <- scale(ps[, 4:8], center = pc2$center, scale = pc2$scale) 13 6 APPLICAZIONI > > + + > > + > > ps2 <- ps1 %*% pc2$loadings[, 1:2] plot(pc2$scores[, 1:2], pch = 20, col = col1, xlim = c(-3.9, 3.5), ylim = c(-3, 3.2), xlab = "PC1 (48%)", ylab = "PC2 (32%)", main = "Indicatori Economici EU", sub = "Lussemburgo, Lettonia punti supplementari") abline(h = 0, v = 0, col = "red", lty = "dashed") text(pc2$scores[, 1:2], labels = as.character(eu1$Sigla), pos = 4, cex = 0.6) points(ps2, pch = 20, col = c("red", "black")) text(ps2, labels = c("LU", "LE"), pos = 4, col = "cyan") Indicatori Economici EU 3 LE ES 2 BU LA SL FI SW 0 SP IR LU RO CZ CY NE RS AU −1 PC2 (32%) 1 DE PL UK GE HU PO FR MA −2 BE GR −3 IT −4 −2 0 2 PC1 (48%) Lussemburgo, Lettonia punti supplementari 6.2 Correlazione lineare e componenti principali L’ACP è un metodo lineare e quindi dipende dalle varianze e dalle correlazioni lineari tra le variabili X1 , ..., Xp . Se X1 , ..., Xp sono incorrelate (le matrici di covarianza e di correlazione sono diagonali), le componenti principali coincidono con le variabili osservate, a meno di un riordinamento secondo la varianza, e l’informazione non è comprimibile per mezzo di combinazioni lineari. Viceversa, se la situazione è di collinearità, le variabili osservate sono ridondanti ed è possibile ridurre la dimensionalità. Il caso limite è quello in cui RX = Jp = 1p 1Tp , la matrice pxp con tutti gli elementi uguali a 1. Tale matrice ha rango 1 e quindi una sola combinazione lineare può assorbire tutta l’informazione. Nel seguente esempio ne calcoliamo autovalori e autovettori. 14 6 APPLICAZIONI Esempio 10 (Autovalori e autovettori della matrice Jp ) Poichè Jp ha rango 1, l’unico autovalore positivo è uguale a p, gli altri sono tutti uguali a 0. L’autovettore u 1 associato al massimo autovalore λ1 = p si ricava usando la definizione. Esso deve verificare l’equazione matriciale 1p 1Tp u1 = (1Tp u1 )1p = pu1 e quindi dev’essere proporzionale al vettore 1p . Imponendo il vincolo di normalizzazione otteniamo u1 = 1p /p1/2 . In termini di componenti principali 1. La prima componente principale spiega il 100% della varianza totale il che significa che l’informazione è essenzialmente uni-dimensionale. 2. La combinazione lineare ottimale (massima varianza) è proporzionale al vettore 1 p , quindi attribuisce pesi uguali alle variabili osservate ed ha il significato di una somma o di una media. 3. Le componenti principali successive alla prima hanno varianza identicamente nulla e sono irrilevanti dal punto di vista statistico. I loro autovettori si possono ricavare usando il vincolo d’incorrelazione. Le situazioni più comuni sono intermedie tra l’indipendenza lineare e la completa collinearità, come nell’esempio seguente. Esempio 11 (Strutture di correlazione) Supponiamo che la 1 1 0 0 1 1 0 0 RX = 0 0 1 −0.5 0 0 −0.5 1 matrice di correlazione di X 1 , ..., X4 sia . Essa indica che le variabili sono divise in due gruppi, {X1 , X2 } e {X3 , X4 }. Le variabili del primo gruppo sono perfettamente correlate, quelle del secondo gruppo sono correlate, ma la correlazione è lontana dal massimo. Inoltre, variabili appartenenti a gruppi diversi sono incorrelate. Ricaviamo autovalori e autovettori con la funzione eigen di R. > r <- matrix(c(1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, -0.5, 0, 0, -0.5, + 1), nrow = 4, ncol = 4, byrow = FALSE) > auto <- eigen(r) > auto$values [1] 2.000000e+00 1.500000e+00 5.000000e-01 9.860761e-32 > auto$vectors [1,] [2,] [3,] [4,] [,1] 0.7071068 0.7071068 0.0000000 0.0000000 [,2] 0.0000000 0.0000000 0.7071068 -0.7071068 [,3] 0.0000000 0.0000000 0.7071068 0.7071068 [,4] 0.7071068 -0.7071068 0.0000000 0.0000000 La matrice ha rango 3 (il quarto autovalore è numericamente uguale a 0). La prima componente principale spiega da sola il 50% della varianza totale ed è essenzialmente la somma delle variabili del primo gruppo. La seconda componente spiega il 37.5% della varianza totale ed è essenzialmente la differenza delle variabili del secondo gruppo. La terza componente spiega il restante 12.5% della varianza totale ed è la somma delle variabili del secondo gruppo. 15 A RICHIAMI DI ALGEBRA DELLE MATRICI A Richiami di algebra delle matrici Una matrice A = (aij ) è una tabella mxn di numeri ordinati in m righe e n colonne. Quando m = n, la matrice è quadrata. Se le matrici A = (aij ) e B = (bij ) hanno lo stesso numero di righe e lo stesso numero di colonne, la loro somma A + B è la matrice C = (cij ), cij = aij + bij . Il prodotto della matrice A = (aij ) per la costante numerica λ è la matrice λA = (λaij ). Se A è mxn e B è nxp, il prodotto (righe per colonne) di A e B è AB = C = (cij ), cij = n X ail blj l=1 Dalla definizione segue che il prodotto righe per colonne di A e B richiede che le matrici siano opportunamente conformate: il numero di colonne di A è uguale al numero di righe di B. L’esistenza di AB non implica in generale l’esistenza di BA e, anche quando esistono entrambi i prodotti, non è detto che risulti AB = BA. Questo prodotto matriciale non gode della proprietà commutativa. L’operatore R del prodotto righe per colonne è %*%. L’esempio ne illustra l’applicazione. > ma <- matrix(c(1, 2, 0, 1, -1, 0), nrow = 2, ncol = 3, byrow = FALSE) > ma [1,] [2,] [,1] [,2] [,3] 1 0 -1 2 1 0 > mb <- matrix(c(1, -1, 0, 0, 1, 1, 0, 2, 0, 1, 1, 3), nrow = 3, + ncol = 4, byrow = FALSE) > mb [1,] [2,] [3,] [,1] [,2] [,3] [,4] 1 0 0 1 -1 1 2 1 0 1 0 3 > ma %*% mb [1,] [2,] [,1] [,2] [,3] [,4] 1 -1 0 -2 1 1 2 3 > md <- matrix(c(5, 0, -1, -3, 1, 1), nrow = 3, ncol = 2, byrow = FALSE) > md [1,] [2,] [3,] [,1] [,2] 5 -3 0 1 -1 1 > ma %*% md 16 A RICHIAMI DI ALGEBRA DELLE MATRICI [1,] [2,] [,1] [,2] 6 -4 10 -5 > md %*% ma [1,] [2,] [3,] [,1] [,2] [,3] -1 -3 -5 2 1 0 1 1 1 La trasposta della matrice A mxn, indicata con AT , è la matrice nxm che si ottiene scrivendo le righe di A come colonne AT = (aji ). Valgono le proprietà (AT )T = A, (A + B)T = AT + B T , (AB)T = B T AT Le matrici con una sola colonna vengono chiamate vettori colonna e le loro trasposte vettori riga. I vettori sono indicati con lettere minuscole. Il vettore (nullo) nx1 con tutti gli elementi uguali a 0 viene indicato con 0n . Il vettore nx1 con tutti gli elementi uguali a 1 viene indicato con 1n . Se x è un vettore nx1, n X xT x = x2i i=1 T 1/2 e la norma di x è k x k= (x x) . La norma è la generalizzazione multidimensionale della familiare nozione di lunghezza di un vettore. Se x e y sono vettori nx1, T x y= n X xi y i i=1 è il prodotto interno di x e y. Il prodotto interno diviso per il prodotto della norma di x e della norma di y, xT y/(k x k k y k), è il coseno dell’angolo formato da x e y. Se xT y = 0, i due vettori sono ortogonali. Se A è quadrata d’ordine m, diag A = B = (bij ) con bii = aii , i = 1, ..., m e bij = 0, i 6= j. Una matrice A quadrata viene detta diagonale se diag A = A. La matrice identità d’ordine m, indicata con Im , è una matrice diagonale con elementi diagonali tutti uguali a 1. Se A e I sono conformate per il prodotto AI = IA = A. Una matrice A quadrata viene detta ortogonale se A T A = AAT = I. Questa proprietà indica che i vettori riga (colonna) di A sono vettori di norma unitaria a due a due ortogonali. La traccia di una matrice A quadrata d’ordine m, indicata con tr A, è la somma dei suoi elementi diagonali tr A = m X aii . i=1 Valgono le seguenti proprietà tr(A + B) = tr A + tr B, tr(λA) = λtr A, tr(AB) = tr(BA). La terza proprietà richiede che AB sia una matrice quadrata. Il rango di riga (colonna) di una matrice è il massimo numero di vettori riga (colonna) linearmente indipendenti. Il rango di una matrice è il minimo valore del rango di riga e del rango di colonna. Per una matrice A mxn il rango r(A) soddisfa la diseguaglianza r(A) ≤ min{m, n}. A RICHIAMI DI ALGEBRA DELLE MATRICI 17 Una matrice A quadrata d’ordine m viene detta non singolare se r(A) = m, in caso contrario viene detta singolare. Se A mxm è non singolare, esiste ed è unica la matrice B non singolare tale che AB = BA = Im . La matrice B, indicata con A−1 , viene detta l’inversa di A. Se tutte le inverse esistono, valgono le proprietà (A−1 )T = (AT )−1 , (AB)−1 = B −1 A−1 . Gli operatori R di trasposizione e inversione di matrici sono t e solve, rispettivamente. Inoltre l’operatore diag, applicato ad una matrice quadrata, ne estrae il vettore degli elementi diagonali; applicato ad un vettore costruisce la matrice diagonale con elementi diagonali uguali agli elementi del vettore. > t(ma) [1,] [2,] [3,] [,1] [,2] 1 2 0 1 -1 0 > mc <- matrix(c(1, -1, 0, 0, 1, 1, 0, 2, 0), nrow = 3, ncol = 3, + byrow = FALSE) > diag(diag(mc), nrow = length(diag(mc))) [1,] [2,] [3,] [,1] [,2] [,3] 1 0 0 0 1 0 0 0 0 > diag(c(1, 1, 2), nrow = 3) [1,] [2,] [3,] [,1] [,2] [,3] 1 0 0 0 1 0 0 0 2 > diag(1, nrow = 2) [1,] [2,] [,1] [,2] 1 0 0 1 > sum(diag(mc)) [1] 2 > solve(mc) [1,] [2,] [3,] [,1] [,2] [,3] 1.0 0.0 0.0 0.0 0.0 1.0 0.5 0.5 -0.5 > mc %*% solve(mc) 18 B AUTOVALORI E AUTOVETTORI [1,] [2,] [3,] [,1] [,2] [,3] 1 0 0 0 1 0 0 0 1 > solve(mc) %*% mc [1,] [2,] [3,] [,1] [,2] [,3] 1 0 0 0 1 0 0 0 1 Una matrice A quadrata viene detta simmetrica se AT = A. Se x è un vettore nx1 e A è quadrata simmetrica d’ordine n, l’espressione n X xT Ax = aij xi xj i,j=1 viene detta forma quadratica in x rispetto alla matrice A. Per n = 2 essa è un polinomio omogeneo di secondo grado a11 x21 + 2a12 x1 x2 + a22 x22 nelle variabili x1 , x2 , con coefficienti dipendenti dalla matrice A. Le forme quadratiche sono classificate in base al loro segno: xT Ax > 0 per ogni x 6= 0n : definita positiva xT Ax ≥ 0 per ogni x : semi definita positiva xT Ax < 0 per ogni x 6= 0n : definita negativa xT Ax ≤ 0 per ogni x : semi definita negativa xT Ax Q 0 indefinita B Autovalori e autovettori Introduciamo questa importante nozione illustrando la motivazione geometrica. Sia x = (x1 , ..., xn )T un vettore nx1. Identifichiamo lo spazio di riferimento di x con Rn , l’insieme delle n-uple di numeri reali. Ogni vettore può essere rappresentato in Rn mediante un segmento orientato uscente dall’origine 0n con l’estremità in corrispondenza del punto Px di coordinate x1 , ..., xn . Esso è caratterizzato dalla sua direzione (la retta orientata uscente dall’origine a cui il vettore appartiene) e dalla sua lunghezza (la norma k x k) (vedi Figura 1). Sia A una fissata matrice quadrata d’ordine n. Si chiamaPtrasformazione lineare l’applicazione ` : Rn → Rn che ad x associa il vettore nx1 y = Ax, con yi = nl=1 ail xl = aTi x, i = 1, ..., n, essendo ai l’i-esimo vettore riga di A. Esempio 12 (Trasformazioni lineari) Consideriamo la trasformazione lineare associata alla matrice 2 1 A= . 0 −1 Il vettore x = (1, 2)T viene trasformato nel vettore (4, −2)T mentre x = (−3, 3)T viene trasformato nel vettore (−3, −3)T . Il vettore x = (1, 0)T viene trasformato nel vettore (2, 0)T = 2x e, a differenza dei casi precedenti, la trasformazione cambia la lunghezza di x, ma non la direzione (vedi Figura 2). 19 3 B AUTOVALORI E AUTOVETTORI 2 u 0 1 −3v −1 v/2 v −3 −2 −u/2 −3 −2 −1 0 1 2 3 Figura 1: Vettori di R2 : u = (2, 2)T sulla retta y = x,v = (2, −1)T sulla retta y = −x/2. Un problema di grande interesse teorico e pratico è l’identificazione delle rette per l’origine invarianti rispetto ad una trasformazione lineare. Il problema ammette una semplice impostazione algebrica. Infatti, per ogni vettore x, il vettore λx che si ottiene moltiplicando x per lo scalare λ rappresenta un vettore uscente dall’origine con la stessa direzione di x, verso uguale o opposto a seconda che sia λ > 0 oppure λ < 0 e lunghezza pari a | λ | volte la lunghezza di x. Pertanto la ricerca delle direzioni invarianti è ricondotta alla ricerca dei vettori x e degli scalari λ che verificano la seguente equazione matriciale Ax = λx. (1) Ogni soluzione λ∗ viene chiamata autovalore di A ed il vettore x∗ viene chiamato autovettore di A associato a λ∗ . Esempio 13 (Autovalori e autovettori) Usando la definizione (1), si verifica facilmente che gli autovalori della matrice A dell’Esempio 12 sono λ1 = 2 associato all’autovettore a1 = (1, 0)T , λ2 = −1 associato all’autovettore a2 = (−1/10−1/2, 3/10−1/2 )T (vedi Figura 2). Notiamo a questo punto che gli autovalori possono essere numeri complessi, anche quando gli elementi di A sono numeri reali, e questo perchè vogliamo trovare ogni possibile vettore x (anche complesso), diverso dal vettore nullo, soddisfacente la (1). Inoltre, gli autovalori possono essere uguali a 0 purchè esista un vettore x∗ , diverso dal vettore nullo, tale che Ax∗ = 0n . Il calcolo degli autovalori e autovettori di una matrice può essere complicato e viene risolto grazie a speciali programmi di calcolo numerico, disponibili anche in R. 20 4 B AUTOVALORI E AUTOVETTORI v 2 u w Az 0 z −2 Aw Au −4 Av −4 −2 0 2 4 Figura 2: La trasformazione x → Ax (tratteggio: vettori trasformati; verde: autovettori di A). Per calcolare gli autovalori ed autovettori di una matrice si può usare la funzione R eigen. Essa produce una lista con componenti 1. $values: vettore degli autovalori, 2. $vectors: matrice avente nelle colonne gli autovettori. Verifichiamo, usando la funzione eigen, che la matrice mc ha un autovalore reale e due autovalori complessi coniugati mentre gli autovettori hanno elementi complessi. > mc <- matrix(c(2, 0, -1, 0, 1, -2, 1, 1, 0), nrow = 3, ncol = 3, + byrow = FALSE) > auto <- eigen(mc) > auto$values [1] 1.770917+0.000000i 0.614542+1.563885i 0.614542-1.563885i > auto$vectors [,1] [,2] [,3] [1,] -0.9362651+0i 0.2497824+0.2819505i 0.2497824-0.2819505i [2,] 0.2782173+0i 0.1169313+0.4744142i 0.1169313-0.4744142i [3,] 0.2144824+0i -0.7870012+0.0000000i -0.7870012+0.0000000i 21 B AUTOVALORI E AUTOVETTORI Riepiloghiamo di seguito alcuni utili risultati riguardanti autovalori e autovettori, supponendo A quadrata d’ordine n. 1. Autovalori e autovettori sono definiti soltanto per matrici quadrate. 2. Gli autovettori sono per definizione diversi dal vettore nullo mentre gli autovalori possono essere uguali a 0. 3. Gli autovettori non sono determinati univocamente; infatti, se x∗ è un autovettore associato all’autovalore λ∗ , lo è anche kx∗ , per ogni scalare k 6= 0. Di qui la consuetudine di ridurli a norma unitaria (quando sono reali). Anche in questo caso, sono comunque definiti a meno di una moltiplicazione per −1. 4. A e AT hanno gli stessi autovalori; gli autovettori possono essere diversi. Pn 5. tr A = i=1 λi . Qn 6. det A = i=1 λi . 7. A è non singolare se e solo se tutti gli autovalori sono diversi da 0. 8. Una matrice singolare ha almeno un autovalore nullo. Le matrici reali simmetriche godono di ulteriori notevoli proprietà. P1. Gli autovalori sono reali. P2. Autovettori associati ad autovalori distinti sono ortogonali. P3. Sia A una matrice simmetrica d’ordine n. Indichiamo con Λ = diag(λ1 , ..., λn ) la matrice diagonale degli autovalori ordinati λ1 ≥ ... ≥ λn e con U la matrice pxp (ortogonale) le cui colonne sono ordinatamente gli autovettori ortonormali corrispondenti agli autovalori. Vale la seguente rappresentazione (teorema degli assi principali) A = U ΛU T . (2) P4 Una matrice simmetrica è definita positiva (semi definita positiva) se e solo se tutti gli autovalori sono positivi (non negativi). Le proprietà P1-P4 sono applicate nell’analisi delle componenti principali, che usa autovalori e autovettori delle matrici di covarianza o di correlazione. Ne riportiamo di seguito un’illustrazione. > b <- matrix(c(9, 4, 2, 4, 4, 2, 2, 2, 1), nrow = 3, ncol = 3, + byrow = FALSE) > print(c("Traccia = ", sum(diag(b))), quote = FALSE) [1] Traccia = 14 > print(c("Determinante = ", det(b)), quote = FALSE) [1] Determinante = 0 > auto <- eigen(b) > auto$values [1] 1.189898e+01 2.101021e+00 -1.561251e-17 22 RIFERIMENTI BIBLIOGRAFICI > print(c("Somma autovalori = Traccia = ", sum(auto$values)), quote = FALSE) [1] Somma autovalori = Traccia = 14 > print(c("Prodotto autovalori = Determinante = ", prod(auto$values)), + quote = FALSE) [1] Prodotto autovalori = Determinante = -3.90312782094782e-16 > auto$vectors [,1] [,2] [,3] [1,] 0.8391211 0.5439447 0.0000000 [2,] 0.4865189 -0.7505327 -0.4472136 [3,] 0.2432595 -0.3752663 0.8944272 > print("Verifica ortogonalità matrice autovettori", quote = FALSE) [1] Verifica ortogonalità matrice autovettori > auto$vectors %*% t(auto$vectors) [,1] [,2] [,3] [1,] 1.000000e+00 6.697659e-17 3.104884e-17 [2,] 6.697659e-17 1.000000e+00 -1.260385e-17 [3,] 3.104884e-17 -1.260385e-17 1.000000e+00 > t(auto$vectors) %*% auto$vectors [,1] [,2] [,3] [1,] 1.000000e+00 6.543838e-17 1.885157e-17 [2,] 6.543838e-17 1.000000e+00 -7.678862e-17 [3,] 1.885157e-17 -7.678862e-17 1.000000e+00 > print("Verifica teorema assi principali", quote = FALSE) [1] Verifica teorema assi principali > auto$vectors %*% diag(auto$values) %*% t(auto$vectors) [1,] [2,] [3,] [,1] [,2] [,3] 9 4 2 4 4 2 2 2 1 Riferimenti bibliografici Magnus JR, Neudecker H (2001) Matrix Differential Calculus with Applications in Statistics and Economics, 2nd edn. Wiley, New York Mardia KV, Kent JT, Bibby JM (1979) Multivariate Analysis, 1st edn. Academic Press, London R Development Core Team (2006) R: A language and environment for statistical computing. URL http://www.R-project.org, ISBN 3-900051-07-0