Laboratori + Soluzioni esami di Analisi Esplorativa

Laboratorio 1

Costruzione matrice X

X=matrix(c(3,4,2,6,8,2,5,5,5.5,4,7,10,5,7.5), ncol=2, byrow=FALSE)

Matrice trasposta

t(X)

Statistiche di sintesi

summary(X)

Calcolare la media e la varianza per la prima variabile

mean(X[,1])

n=nrow(X)

var(X[,1])*(n-1)/n

Calcolare il vettore delle medie

apply(X, MARGIN=2, FUN="mean")

Calcolare la matrice di varianze/covarianze

S=var(X)*(n-1)/n

Calcolare la matrice di correlazione

R=cor(X) è capace solo di percepire la relazione lineare

Costruire il diagramma di dispersione

plot(X)

Statistiche di sintesi - Dati del cervello

library("MASS")

Y=Animals

boxplot.stats(Animals$brain) # baffo sx e dx sono rispettivamente primo e ultimo valore assunto da $stats. n=numerosità, conf=IC, out=outliers

Identificare il nome degli outliers

rownames(Y)[which(Y$brain %in% boxplot.stats(Y$brain)$out)]

Costruire il diagramma di dispersione della relazione brain e body e valutare la correlazione

plot(brain~body, Y)

cor(Y)

Trasformare in log le variabili per aumentare la correlazione

plot(log(brain)~log(body), Y)

cor(log(Y))

Notiamo 3 osservazioni che non si comportano come le altre. Potrebbe trattarsi di outliers bivariati, siccome non riusciamo ad identificare queste osservazioni anomale con l’utilizzo del boxplot, risolviamo con il bagplot

library(aplpack)

bagplot(log(Y))

Ristimo il plot escludendo i 3 outlier. La correlazione aumenta!

plot(log(brain)~log(body),Y[-which(rownames(Y) %in% c("Brachiosaurus","Triceratops","Dipliodocus")),])

cor(log(Y[-which(rownames(Y) %in% c("Brachiosaurus","Triceratops","Dipliodocus")),]))

Diagramma tridimensionale colorando le unità a seconda del sesso

measure=read.csv("https://raw.githubusercontent.com/aldosolari/AE/master/dati/measure.csv")

library(lattice)

cloud(chest ~ waist + hips, group=gender, data=measure)

Matrice dei diagrammi di dispersione

plot(measure[,-4], col=(measure$gender=="male")+1) # 4 è il numero della colonna gender

Rappresentare le unità tramite la tecnica delle facce di Chernoff

library("TeachingDemos")

faces(measure[,-column.gender], scale=T) # scale=T standardizza i dati

Rappresentare le unità tramite la tecnica delle stelle

stars(measure[,-column.gender], scale=T, key.loc = c(12, 10)) # per posizionare la legenda alle coordinate (x,y)

Rappresentare il diagramma di dispersione per latitudine e longitudine dei terremoti

library(datasets)

X=quakes

plot(lat~long, X, pch='.')

Costruire il diagramma di dispersione lat/long condizionato a 5

coplot(lat~long | depth, data=quakes, given.v=co.intervals(quakes$depth, number=5, overlap=0), rows=1)

Descrizione di Tizio e Caio

Tizio è alto 180 cm e pesa 70 Kg, Caio è alto 160 cm e pesa 50 Kg

Tizio=c(180,70)

Caio=c(160,50)

X=rbind(Tizio, Caio)

colnames(X)=c("Altezza","Peso")

Spazio delle variabili

plot(X, xlim=c(0,200), ylim=c(0,200))

text(x=X[,"Altezza"], y=X[,"Peso"], labels=row.names(X), pos=3)

barx=matrix(colMeans(X), ncol=1) # vettore delle medie

baricentro=t(barx) # baricentro (vettore delle medie trasposto)

points(baricentro, pch=19)

Spazio delle osservazioni

tX=t(X)

plot(tX, xlim=c(0,200), ylim=c(0,200), pch=".")

text(x=tX[,"Tizio"], y=tX[,"Caio"], labels=row.names(tX), pos=4)

arrows(x0=0,y0=0,x1=tX[,"Tizio"], y1=tX[,"Caio"])

Vettori scarto dalla media (colonna 1 e 2)

X=matrix(c(4,1,-1,3,3,5), byrow=T, ncol=2, nrow=3)

tildex1=matrix(X[,1]-mean(X[,1]), ncol=1)

tildex2=matrix(X[,2]-mean(X[,2]), ncol=1)

Devianza

dev1=t(tildex1) %*% tildex1 # = (nrow(X)-1)*var(X[,1])

dev2=t(tildex2) %*% tildex2

Codevianza

cod=t(tildex1) %*% tildex2 # = (nrow(X)-1)*var(X)[1,2]

Correlazione

cor=cod/sqrt(dev1*dev2) # = cor(X)[1,2]

Angolo tra tildex1 e tildex2

acos(cor) # in radianti

acos(cor)*(180/pi) # in gradi

Laboratorio 2

Diagramma dispersione X

X=matrix(c(2,3,3,4,4,5,6,6,7,8,7,8,10,6,8,10,12,13,11,12), ncol=2)

n=nrow(X)

plot(X, xlim=c(-4,13), ylim=c(-4,13), asp=1, bg=heat.colors(n), pch=21, cex=2)

abline(h=0)

abline(v=0)

Vettore delle medie

uno=matrix(rep(1,n), ncol=1)

xbar=(1/n) * t(X) %*% uno

Matrice identità

I=diag(rep(1,n))

Matrice di centramento

H=I-(1/n)*uno %*% t(uno)

Simmetria

t(H) = H

sum(t(H)-H)

Idempotenza

HH = H

sum(H%*%H-H)

Matrice dei dati centrati

C=H%*%X

centrare una matrice già centrata, non produce alcun effetto:

sum(H%*%C-C)

Matrice di varianze/covarianze

S=(1/n)*t(C) %*% C

Matrice di correlazione X

R=D^(−1/2)SD^(−1/2)

D=diag(diag(S)^(-.5))

R=D %*% S %*% D

Costruire la matrice di varianze/covarianze X

S=D^(1/2)RD^(1/2)

D2=diag(diag(S)^(.5))

S=D2 %*% R %*% D2

Matrice dei dati standardizzati

Z=C %*% D

Matrice di var/covar dei dati centrati

S.C=(1/n)*t(H %*% C) %*% (H %*% C)

Matrice di var/covar dei dati standardizzati

S.Z=(1/n)*t(H %*% Z) %*% (H %*% Z)

Matrice di correlazione dei dati centrati

R.C=diag(diag(S.C)^(-.5)) %*% S.C %*% diag(diag(S.C)^(-.5))

Matrice di correlazione dei dati standardizzati

S.Z=diag(diag(S.Z)^(-.5)) %*% S.Z %*% diag(diag(S.Z)^(-.5))

Varianza totale

S=matrix(c(2.2,0.4,0.4,2.8),ncol=2)

sum(diag(S))

Varianza generalizzata

det(S)

Indice relativo di variabilità

det(S)/prod(diag(S))

Autovalori di S

Lambda1=eigen(S)$value[1]

Lambda2=eigen(S)$value[2]

Lambda=diag(eigen(S)$values) # matrice degli autovalori

Autovettori normalizzati

v1=eigen(S)$vectors[,1, drop=F]

v2=eigen(S)$vectors[,2, drop=F]

V=eigen(S)$vectors # matrice degli autovettori

t(v1) %*% v1 # gli autovettori sono di lunghezza unitaria

t(v1) %*% v2 # gli autovettori sono ortogonali

Verifico la decomposizione spettrale

S-V %*% Lambda %*% t(V)

Verifico l'ortogonalità di V

V %*% t(V)-diag(c(1,1))

Calcolare S^(1/2)

SqrtS=V %*% Lambda^(1/2) %*% t(V) # S^(1/2)

InvS=solve(S) # S^(-1)

InvS %*% S # verifico che è l'inversa di S

la varianza totale di S è uguale alla somma degli autovalori di S

sum(diag(S)) == sum(diag(Lambda))

la varianza generalizzata di S è uguale al prodotto degli autovalori di S

det(S) == prod(diag(Lambda))

Matrice dei dati ortogonalizzati Ztilde

X=matrix(c(2,3,3,4,4,5,6,6,7,8,7,8,10,6,8,10,12,13,11,12),ncol=2)

n=nrow(X)

Ztilde = Xtilde S^(−1/2)

uno=matrix(rep(1,n),ncol=1)

I=diag(rep(1,n))

H=I-(1/n)*uno %*% t(uno)

Xtilde=H %*% X

S=(1/n)*t(Xtilde) %*% Xtilde

Lambda=diag(eigen(S)$values)

V=eigen(S)$vectors

InvSqrtS=V %*% diag(diag(Lambda)^(-.5)) %*% t(V) # S^(1/2)

Ztilde= Xtilde %*% InvSqrtS

(1/n)*t(Ztilde) %*% uno # vettore delle medie

(1/n)*t(H %*% Ztilde) %*% (H %*% Ztilde) # matrice di varianze e covarianze

Decomposizione in valori singolari

SVD=svd(Xtilde)

Rango di Xtilde

r=qr(Xtilde)$rank

Verifico SVD

Xtilde-SVD$u %*% diag(SVD$d) %*% t(SVD$v)

Calcolare la 1° componente principale di Xtilde

marks=read.table("http://www.maths.leeds.ac.uk/~charles/mva-data/openclosedbook.dat",header=T)

X=as.matrix(marks)

colnames(X)=c("Mechanics", "Vectors", "Algebra", "Analysis", "Statistics")

n=nrow(X)

p=ncol(X)

Xtilde = X %*% t(diag(rep(1,p))) + matrix(rep(1,n)) %*% (-t((1/n) * t(X) %*% matrix(rep(1,n))))

S=(1/n)*t(Xtilde) %*% Xtilde

Lambda=diag(eigen(S)$values) # l'autovalore più grande è il primo

V=eigen(S)$vectors

v1=V[,1, drop=F] # pesi (loadings) della 1° componente principale

y1=Xtilde %*% v1 # punteggi (scores) della 1° componente principale

var(y1)*(n-1)/n # La varianza di y1 coincide con Lambda[1,1]

Calcolare le p componenti principali Y=Xtilde*V

Y = Xtilde %*% V

Verificare che il vettore medio di Y è nullo

(1/n) * t(Y) %*% uno

La matrice di var/covar di Y è pari a Lambda

SY=(1/n) * t(Y) %*% Y

round(SY,3)==round(Lambda,3)

Calcolo le componenti principali

pca=princomp(X) # Standard deviation = sqrt(diag(Lambda))

summary(pca)

Calcolo pesi e punteggi

pca$loadings # pesi: coincide con V

pca$scores # punteggi: coincide con head(Y)

pca=prcomp(X) # uguale all'altro comando ma con varianza corretta

summary(pca) # Standard deviation coincide con sqrt(diag(Lambda)*(n/n-1))

pca$rotation[,] # pesi: coincide con V

head(pca$x) # punteggi: coincide con head(Y)

Biplot

biplot(pca)

Calcolare la correlazione tra il voto centrato sullo 0 dell’esame in Statistics e i punteggi della 1° componente principale

cor(Xtilde[,"Statistics"],Y[,1])

Scegliere il numero q di componenti principali utilizzando i criteri:

• Proporzione di varianza spiegata dalle prime q componenti superiore all’80%

summary(pca) # scelgo le componenti 1 e 2

• Varianza spiegata da ciascuna componente maggiore del λ medio

pca$sdev^2 > mean(pca$sdev^2) # scelgo solo la componente 1

• Utilizzando lo scree plot

plot(pca, type="line")

Ottenere la migliore approssimazione per Xtilde di rango 10

face=read.table("https://raw.githubusercontent.com/aldosolari/AE/master/dati/face.txt")

X=as.matrix(face)

n=nrow(face)

p=ncol(face)

pca=princomp(X, cor=F)

V=pca$loadings

Y=pca$scores

xbar=matrix(pca$center, ncol=1)

q=10

Yq=Y[,1:q]

Vq=V[,1:q]

Aq=Yq %*% t(Vq) # migliore approssimazione di rango q

one.n = matrix(rep(1,n), ncol=1)

face2 = Aq + one.n %*% t(xbar) # compressione immagine

face2=pmax(pmin(face2, 1), 0) # forzo i valori tra 0 e 1

pixels=prod(dim(face))

pixels2=prod(dim(Yq)) + prod(dim(Vq)) + prod(dim(xbar))

round(pixels/pixels2,2) # fattore di riduzione pixels

size=object.size(face)

size2=object.size(Yq) + object.size(Vq) + object.size(xbar)

round(size/size2)