Analisi delle Componenti Principali
Dati Marks
- Caricare i dati
marks
marks <- read.table("http://www.maths.leeds.ac.uk/~charles/mva-data/openclosedbook.dat",header = TRUE)
X = as.matrix(marks)
# assegno i nomi alle variabili
colnames(X) <- c("Mechanics", "Vectors", "Algebra", "Analysis", "Statistics")
# guardo le prime righe
head(X) Mechanics Vectors Algebra Analysis Statistics
[1,] 77 82 67 67 81
[2,] 63 78 80 70 81
[3,] 75 73 71 66 81
[4,] 55 72 63 70 68
[5,] 63 63 65 70 63
[6,] 53 61 72 64 73n = nrow(X)
p = ncol(X)- Calcolare la matrice dei dati centrati \(\tilde{X}\) come trasformazione lineare \(\underset{n\times p}{X}\underset{p\times q}{A'} + \underset{n\times 1}{1}\underset{1\times q}{b'}\) con \(q=p\), \(A = \underset{p\times p}{I}\) e \(b = -\underset{p\times 1}{\bar{x}}\)
A = diag(rep(1,p))
one.n = matrix(rep(1,n))
b = (1/n) * t(X) %*% one.n
Xtilde = X %*% t(A) + one.n %*% (-t(b))- Calcolare il voto medio di ciascun studente come combinazione lineare \(\underset{n\times 1}{y} = \underset{n\times p}{X}\underset{p\times 1}{a}\) con \(a_j = 1/p\), \(j=1,\ldots,p\)
a = matrix(rep(1/p, p), ncol=1)
y = X %*% a- Calcolare la prima componente principale di \(\tilde{X}\) come \(\underset{n\times 1}{y_1} = \underset{n\times p}{\tilde{X}} \underset{p\times 1}{v_1}\) dove \(v_1\) è il primo autovettore di \(S=\frac{1}{n} \tilde{X}' \tilde{X}\) associato all’autovalore più grande \(\lambda_1\). Verificare che la varianza di \(y_1\) è pari a \(\lambda_1\) e che è maggiore della varianza della combinazione lineare normalizzata \(\underset{n\times 1}{y} = \underset{n\times p}{\tilde{X}}\underset{p\times 1}{a}\) con \(a_j = 1/\sqrt{p}\), \(j=1,\ldots,p\).
# decomposizione spettrale di S
S = (1/n) * t(Xtilde) %*% Xtilde
eigen = eigen(S)
Lambda = diag(eigen$values)
V = eigen$vectors
# pesi (loadings) della 1ma componente principale
v1 = V[,1, drop=FALSE]
v1 [,1]
[1,] -0.5054457
[2,] -0.3683486
[3,] -0.3456612
[4,] -0.4511226
[5,] -0.5346501# punteggi (scores) della 1ma componente principale
y1 = Xtilde %*% v1
# varianza di y1
var(y1) * (n-1)/n # coincide con Lambda[1,1] [,1]
[1,] 679.1831# confronto con altra combinazione lineare normalizzata
a = matrix(rep(1/sqrt(p), p), ncol=1)
y = Xtilde %*% a
var(y) * (n-1)/n [,1]
[1,] 662.6463- Calcolare le \(p\) componenti principali \(\underset{n\times p}{Y} = \underset{n\times p}{\tilde{X}} \underset{p\times p}{V}\). Verificare che il vettore medio di \(Y\) è nullo, la matrice di varianze/covarianze \(S^Y\) di \(Y\) è pari a \(\Lambda\), che la varianza totale e generalizzata di \(S^Y\) è pari a quella di \(S\)
# p componenti principali Y
Y = Xtilde %*% V
# vettore medio di Y
round(
(1/n) * t(Y) %*% one.n
, 8) [,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0# matrice di varianze covarianze di Y
S_Y = (1/n) * t(Y) %*% Y # coincide con Lambda
round(
S_Y
, 8) [,1] [,2] [,3] [,4] [,5]
[1,] 679.1831 0.0000 0.0000 0.00000 0.00000
[2,] 0.0000 199.8144 0.0000 0.00000 0.00000
[3,] 0.0000 0.0000 102.5684 0.00000 0.00000
[4,] 0.0000 0.0000 0.0000 83.66873 0.00000
[5,] 0.0000 0.0000 0.0000 0.00000 31.78791# varianza totale di S_Y
sum(diag(S_Y)) # coincide con sum(diag(S_Y))[1] 1097.022# varianza generalizzata di S_Y
det(S_Y) # coincide con det(S)[1] 37021339491- Calcolare le componenti principali con il comando
princomp()eprcomp()
pca = princomp(X)
summary(pca) # Standard deviation coincide con sqrt(diag(Lambda))Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
Standard deviation 26.061142 14.1355705 10.12760414 9.14706148 5.63807655
Proportion of Variance 0.619115 0.1821424 0.09349705 0.07626893 0.02897653
Cumulative Proportion 0.619115 0.8012575 0.89475453 0.97102347 1.00000000# pesi
pca$loadings[,] # coincide con V Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
Mechanics 0.5054457 0.74874751 0.2997888 0.296184264 0.07939388
Vectors 0.3683486 0.20740314 -0.4155900 -0.782888173 0.18887639
Algebra 0.3456612 -0.07590813 -0.1453182 -0.003236339 -0.92392015
Analysis 0.4511226 -0.30088849 -0.5966265 0.518139724 0.28552169
Statistics 0.5346501 -0.54778205 0.6002758 -0.175732020 0.15123239# punteggi
head( pca$scores ) # coincide con head(Y) Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
[1,] 66.32077 6.447125 7.0736275 -9.6463833 5.4557651
[2,] 63.61810 -6.754424 0.8599283 -9.1490636 -7.5656517
[3,] 62.92626 3.080258 10.2297139 -3.7238434 -0.3841125
[4,] 44.53775 -5.577218 -4.3780192 -4.4816746 4.4065605
[5,] 43.28425 1.133228 -1.5314139 5.8059805 0.7378218
[6,] 42.55249 -10.972900 4.8671678 -0.4788987 -7.1021171pca2 = prcomp(X, center = TRUE)
summary(pca2) # Standard deviation coincide con sqrt(diag(Lambda)*(n/n-1))Importance of components:
PC1 PC2 PC3 PC4 PC5
Standard deviation 26.2105 14.2166 10.1856 9.19948 5.67039
Proportion of Variance 0.6191 0.1821 0.0935 0.07627 0.02898
Cumulative Proportion 0.6191 0.8013 0.8948 0.97102 1.00000# pesi
pca2$rotation[,] # coincide con V PC1 PC2 PC3 PC4 PC5
Mechanics -0.5054457 -0.74874751 0.2997888 -0.296184264 -0.07939388
Vectors -0.3683486 -0.20740314 -0.4155900 0.782888173 -0.18887639
Algebra -0.3456612 0.07590813 -0.1453182 0.003236339 0.92392015
Analysis -0.4511226 0.30088849 -0.5966265 -0.518139724 -0.28552169
Statistics -0.5346501 0.54778205 0.6002758 0.175732020 -0.15123239# punteggi
head( pca2$x ) # coincide con head(Y) PC1 PC2 PC3 PC4 PC5
[1,] -66.32077 -6.447125 7.0736275 9.6463833 -5.4557651
[2,] -63.61810 6.754424 0.8599283 9.1490636 7.5656517
[3,] -62.92626 -3.080258 10.2297139 3.7238434 0.3841125
[4,] -44.53775 5.577218 -4.3780192 4.4816746 -4.4065605
[5,] -43.28425 -1.133228 -1.5314139 -5.8059805 -0.7378218
[6,] -42.55249 10.972900 4.8671678 0.4788987 7.1021171- Calcolare la correlazione tra il voto centrato sullo 0 dell’esame in
Statistics(quinta colonna della matrice \(\tilde{X}\)) e i punteggi della prima componente principale (prima colonna della matrice \(Y\)). Si ricordi la formula \[\frac{v_{jk}\sqrt{\lambda_k}}{\sqrt{s_{jj}}}\] per la \(j\)-sima colonna di \(\tilde{X}\) e la \(k\)-sima colonna di \(Y\).
V[5,1]*pca$sdev[1]/sqrt(diag(S)[5]) Comp.1
-0.8121103 # corrisponde alla correlazione calcolata sulla base dei dati
cor(Xtilde[,"Statistics"],Y[,1])[1] -0.8121103- Scegliere il numero \(q\) di componenti principali utilizzando i criteri (i) proporzione di varianza spiegata dalle prime \(q\) componenti superiore all’80%; (ii) varianza spiegata da ciascuna componente maggiore di \(c=\frac{1}{p}\sum_{j=1}^{p}\lambda_j\) (iii) utilizzando lo scree plot.
# scelta di q: proporzione di varianza spiegata > 80%
summary(pca)Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
Standard deviation 26.061142 14.1355705 10.12760414 9.14706148 5.63807655
Proportion of Variance 0.619115 0.1821424 0.09349705 0.07626893 0.02897653
Cumulative Proportion 0.619115 0.8012575 0.89475453 0.97102347 1.00000000# scelta di q: varianza spiegata > c
c = mean(pca$sdev^2)
pca$sdev^2 > cComp.1 Comp.2 Comp.3 Comp.4 Comp.5
TRUE FALSE FALSE FALSE FALSE # scelta di q: scree plot
plot(pca, type="line")
Biplot
Con i dati utilizzati nell’esempio precedente, costruire il biplot con il comando biplot().
# biplot
biplot(pca)
La funzione biplot() rappresenta
- gli AUTOVETTORI NON STANDARDIZZATI = AUTOVETTORI * SQRT(AUTOVALORI), moltiplicati poi per un fattore 0.8 * sqrt(n-1)
- i PUNTEGGI RISCALATI (in modo che la somma dei valori al quadrato sia pari a 1)
In R:
# AUTOVETTORI (PESI)
V = eigen(S)$vectors
V[,1] = -V[,1]
# AUTOVALORI (VARIANZE)
lambda = eigen(S)$values
# AUTOVETTORI NON STANDARDIZZATI
VL = V %*% diag(sqrt(lambda))
# PUNTEGGI
Y = Xtilde %*% V
# BIPLOT
PCA = princomp(Xtilde)
biplot(PCA)
# sovrascrivo le frecce in blu per le prime due variabili
arrows(0, 0,
VL[1,1] * 0.8 * sqrt(n - 1),
VL[1,2] * 0.8 * sqrt(n - 1),
lwd = 1, angle = 30, length = 0.1, col = 4)
arrows(0, 0,
VL[2,1] * 0.8 * sqrt(n - 1),
VL[2,2] * 0.8 * sqrt(n - 1),
lwd = 1, angle = 30, length = 0.1, col = 4)
# PUNTEGGI RISCALATI (la somma dei valori al quadrato di ciascuna colonna è 1)
Ystd = Y %*% diag(1/(sqrt(lambda)*sqrt(n)))
colSums(Ystd^2)[1] 1 1 1 1 1labs = as.character(1:n)
plot(Ystd[,1:2], asp=1, type = "n")
text(Ystd[,1:2], labs, col=4)
Dati Wine
- Importare i dati ed effettuare l’analisi delle componenti principali dei dati standardizzati \(\underset{n\times p}{Z}\) (escludendo la variabile
type), calcolando la matrice dei punteggi \(\underset{n\times p}{Y} = \underset{n\times p}{X}\underset{p\times p}{Z}\) e la matrice dei pesi \(\underset{p\times p}{V}\). Calcolare la percentuale di varianza spiegata dalle prime \(q\) componenti principali, per \(q=1,\ldots,10\).
rm(list=ls())
wine <- read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", header = FALSE)
colnames(wine)<-c("Type","Alcohol",
"Malic","Ash",
"Alcalinity","Magnesium",
"Phenols","Flavanoids",
"Nonflavanoids","Proanthocyanins",
"Color_int","Hue","Dilution","Proline")
wine$Type = factor(wine$Type)
X = as.matrix(wine[,-1])
n = nrow(X)
p = ncol(X)
# PCA
pca = princomp(X, cor=T)
V = pca$loadings
Y = pca$scores
# prop. di varianza spiegata dalle prime q cp
q <- 10
cumsum(pca$sdev^2/sum(pca$sdev^2))[1:q] Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
0.3619885 0.5540634 0.6652997 0.7359900 0.8016229 0.8509812 0.8933680 0.9201754 0.9423970 0.9616972 - Costruire il diagramma di dispersione dei punteggi delle prime due componenti principali, colorando i punti con colori diversi a seconda della tipologia di vino (variabile
type).
plot(Y[,1:2], col=as.numeric(wine$Type), asp=1)
- Costruire il biplot con la funzione
autoplotdella libreriaggfortify.
library(ggfortify)
autoplot(pca, data = wine, colour = "Type",
loadings=TRUE, loadings.label = TRUE) + theme_bw()
Dati Face
- Importare i dati e visualizzare l’immagine con il comando
image. Effettuare l’analisi delle componenti principali dei dati centrati \(\underset{n\times p}{\tilde{X}}\), calcolando la matrice dei punteggi \(\underset{n\times p}{Y} = \underset{n\times p}{\tilde{X}}\underset{p\times p}{V}\) e la matrice dei pesi \(\underset{p\times p}{V}\).
rm(list=ls())
face <- read.table("https://raw.githubusercontent.com/aldosolari/AE/master/docs/dati/face.txt", header=FALSE)
X = as.matrix(face)
n = nrow(face)
p = ncol(face)
# visualizza immagine
image(X, col=gray(0:255/255), asp=p/n)
# PCA
pca = princomp(X, cor=F)
V = pca$loadings
Y = pca$scores
xbar = matrix(pca$center, ncol=1)- Ottenere la migliore approssimazione per \(\underset{n\times p}{\tilde{X}}\) di rango \(q\), \(\underset{n\times p}{A_q} = \underset{n\times q}{Y_q} \underset{q\times p}{V_q'}\), con \(q=10\). Costruire l’immagine compressa \(\underset{n\times p}{C} = \underset{n\times p}{A_q} + \underset{n\times 1}{1}\underset{1\times p}{\bar{x}}\), assicurandosi che tutti gli elementi di \(\underset{n\times p}{C}\) siano compresi tra 0 e 1.
q = 10
Yq = Y[,1:q]
Vq = V[,1:q]
# migliore approssimazione di rango q
Aq = Yq %*% t(Vq)
# compressione immagine
one.n = matrix(rep(1,n), ncol=1)
face2 = Aq + one.n %*% t(xbar)
# forzo i valori tra 0 e 1
face2 <- pmax(pmin(face2, 1), 0)- Visualizzare l’immagine compressa e confrontarla con l’immagine originale calcolando il fattore di riduzione in termini di pixels e bytes (utilizzando il comando
object.size)
# visualizza immagine compressa
image(face2, col=gray(0:255/255), asp=p/n)
# salve immagine compressa
# library(png)
# writePNG(face,"face.png")
# confronta pixels utilizzati
pixels = prod(dim(face))
pixels2 = prod(dim(Yq)) + prod(dim(Vq)) + prod(dim(xbar))
round(pixels/pixels2, 2) # fattore di riduzione[1] 11.02