Laboratori + Soluzioni esami di Analisi Esplorativa

DOPPIA ENTRATA CHE INCROCIA I GRUPPI OTTENUTI CON LA VARIABILE BLOCCO

paesi=read.table("https://raw.githubusercontent.com/aldosolari/AE/master/dati/paesi.txt", header=T)

X=paesi[,-c(1,11)]

n=nrow(X)

p=ncol(X)

km=kmeans(X, centers=X[c(1,25,26),], algorithm="Lloyd") # K-medie

table(km$cluster) # numerosità dei gruppi

table(km$cluster, paesi$blocco) # tabella gruppi e blocco

#DETERMINARE I CENTROIDI, LE SOMME DEI QUADRATI WITHIN W E BETWEEN B, IL VALORE DELL'INDICE CH

km$centers # centroidi

W=km$tot.withinss #somme dei quadrati within

B=km$betweenss #somme dei quadrati between

K=3

CH=(B/(K-1))/(W/(n-K)) # indice CH

#COSTRUIRE IL GRAFICO SILHOUETTE BASATO SULLA DISTANZA EUCLIDEA

library(cluster)

D = dist(X, method="euclidean") # matrice delle distanze euclidee

sil=silhouette(x=km$cluster, dist=D)

row.names(sil)=paesi$Country

plot(sil)

#STANDARDIZZARE I DATI SWISS E UTILIZZARE L'ALGORITMO DELLE K-MEDIE PER K=2,4,6 INIZIALIZZANDO I CENTROIDI CON LE OSSERVAZIONI

DI RIGA 1,2,...,K. RIPORTARE PER CIASCUN VALORE DI K IL RISPETTIVO VALORE DELL’INDICE CH

<code>
library(datasets)
X=swiss
n=nrow(X)
p=ncol(X)
S=var(X)*((n-1)/n)
Z=scale(X,center=T,scale=diag(S)^(1/2))
K=c(2,4,6)
CH=vector()
library(cluster)
for (k in 1:length(K)){
  km=kmeans(Z, centers=Z[1:K[k],])
  CH[k]=(km$betweenss/(K[k]-1))/(km$tot.withinss/(n-K[k]))
}
rbind(K,CH)
</code>

#CALCOLARE LA MATRICE DI DISTANZE TRA LE N UNITÀ STATISTICHE UTILIZZANDO LA DISTANZA EUCLIDEA

<code>
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)
p=ncol(X)
D2=dist(X,method="euclidean")
</code>

#la matrice di distanze di Minkowski m≥1 è invariante rispetto alle traslazioni
#la matrice di distanze Euclidee è invariante rispetto alle permutazioni
#la distanza Euclidea è invariante rispetto alle rotazioni, mentre la distanza di Manhattan no

#DISTANZA DI MAHALANOBIS

<code>
library("MASS")
X=as.matrix(log(Animals[-c(6,16,26),]))
colnames(X)=c("logbody","logbrain")
n=nrow(X)
p=ncol(X)
xbar=matrix(colMeans(X),
</code>

nrow=p, ncol=1) # vettore medie
S=var(X)*((n-1)/n) # matrice di var/cov
#Verificare numericamente la presenza di valori anomali calcolando il quadrato della distanza di mahalanobis dal baricentro
confrontandola con il quantile 95% di una chi-quadrato
(X[1,]-xbar) %*% solve(S) %*% (X[1,]-xbar) #distanza di Mahalanobis al quadrato per la prima osservazione
dM2=apply(X,MARGIN=1, function(u) t(u-xbar) %*% solve(S) %*% (u-xbar)) #distanza di Mahalanobis al quadrato per le osservazioni
which(dM2 > qchisq(0.95,df=2)) #il quadrato della distanza di Mahalanobis dal baricentro calcolata sui dati originali è uguale al quadrato della distanza Euclidea dall'origine calcolata sui dati ortogonalizzati
#il quadrato della distanza di Mahalanobis dal baricentro è invariante per trasformazioni di scala (es: kg -> gr)
#CLUSTER ANALYSIS
X=matrix(c(1,3,2,4,1,5,5,5,5,7,4,9,2,8,3,10), ncol=2, nrow=8, byrow=T)
n=nrow(X)
colnames(X)=c("x1","x2")
rownames(X)=1:n
D=dist(X,method="euclidean")
hc.single=hclust(D,

method="single") hc.complete=hclust(D, method="complete") hc.average=hclust(D, method="average") plot(hc.single, hang=-1) plot(hc.complete, hang=-1) plot(hc.average, hang=-1) #MATRICE DI DISSIMILARITÀ CON L'INDICE DI GOWER library(cluster) flowerrow.names(flower)=c("begonia","broom","camellia","dahlia","forget-me-not","fuchsia","geranium","gladiolus","heather","hydrangea","iris","lily","lily-of-the-valley","peony","pink carnation","redrose","scotch rose", "tulip") D=daisy(flower, metric="gower", type=list(symm=c(1,2),asymm=3)) #MATRICE DI DISTANZA EUCLIDEA C1=read.table("http://azzalini.stat.unipd.it/Libro-DM/C1.dat", sep="", header=T) D = dist(C1[,c("x1","x2")], method = "euclidean") #LABORATORIO 4 #Si consideri ilmodello fattoriale con k=1 fattore e p=3 variabili. Si generino n=20 osservazioni iid da questo modello (conoscendo psi, Lambda, f, u), salvando ciascuna osservazione nella riga i-sima della matrice Xn=20p=3 Lambda=matrix(c(0.5,0.2,-0.1), ncol=1) Psi=diag(c(0.0225, 0.0025, 0.008)) X=matrix(NA,nrow=n,ncol=p) for (i in 1:n){ f=rnorm(1) u=rnorm(p, mean=0, sd=sqrt(diag(Psi)) ) X[i,]=Lambda %*% f+u } MATRICE DI VAR/COVAR Σ=ΛΛ′+Ψ Sigma=Lambda %*% t(Lambda) + Psi Per ciascuna osservazione x, si considerari la trasformazione di scala z=Ax per A=D^(−1/2) dove D=diag(Σ), ottenendo la matricedei dati trasformati Y=XA. Calcolare i pesi fattoriali Λy=AΛ e le varianze specifiche Ψy=AΨA′ A=diag(diag(Sigma)^(-1/2)) Y=X %*% A Lambday=A %*% Lambda Psiy=A %*% Psi %*% t(A) pesi fattoriali diag(Psiy) varianze specifiche di y Corrx=A %*% Sigma %*% t(A) Stimare con il metodo della massima verosimiglianza il

<!-- modello fattoriale con k=1 fattore (funzione factanal(), argomentifactors=1 e rotation="none"), ed ottenere -->
<!-- #1) la stima della matrice dei pesi fattoriali Λ^af=factanal(x=X, factors=1, rotation="none")hatLambda=af$loadings[,] -->
<!-- #2) la stima delle comunalità h^2ihatH2=apply(af$loadings^2,1,sum) -->
<!-- #3) la stima delle varianze specifiche ψ^ihatPsi=diag(af$uniquenesses) -->
<!-- STIMARE CON IL METODO ML IL MODELLO FATTORIALE OTTENENDO LA STIMA DELLA MATRICE DEI PESI FATTORIALIR = matrix(c(1.000, 0.439, 0.410, 0.288, 0.329, 0.248,0.439, 1.000, 0.351, 0.354, 0.320, 0.329,0.410, 0.351, 1.000, 0.164, 0.190,0.181,0.288, 0.354, 0.164, 1.000, 0.595, 0.470,0.329, 0.320, 0.190, 0.595, 1.000, 0.464,0.248, 0.329, 0.181, 0.470, 0.464, 1.000),byrow=T, ncol=6)colnames(R)=c("Gaelic", "English", "History", "Arithmetic", "Algebra", "Geometry")n = 220 -->
<!-- MODELLO FATTORIALE CON 2 FATTORI SENZA ROTAZIONEaf=factanal(covmat=R, factors=2, -->

rotation="none", n.obs=n)
Lambda=af$loadings[,] # stima pesi fattoriali

# MODELLO FATTORIALE CON 2 FATTORI CON ROTAZIONE VERIMAX
af2=factanal(covmat=R, factors=2, rotation="varimax", n.obs=n)
af2$rotmat # matrice di rotazione
af2$loadings # matrice dei pesi fattoriali

# Si riporti il p-value del test rapporto di verosimiglianza per verificare H0:k=1 fattore è sufficiente, 
# e si concluda ad un livello di significatività del 5% se è opportuno utilizzare il modello fattoriale ad 1 fattore
factanal(covmat=R, factors=1, n.obs = n)$PVAL
p=6
k=1
af=factanal(covmat=R, factors=k, n.obs=n)
hatLambda = af$loadings[,]
hatPsi=diag(af$uniqueness)
fit=hatLambda %*% t(hatLambda) + hatPsi
t=n*log(det(fit)/det(R))
gdl=((p-1)^2-p-k)/2
pchisq(t,lower.tail=FALSE, df=gdl)

library(bootstrap)
X=s
corn=nrow(X)
p=ncol(X)
R=cor(X) #matrice di correlazione

# Utilizzando la funzione factanal(), stimare con il metodo ML il modello fattoriale con k=2 fattori 
# eseguendo la rotazione varimax
af2=factanal(covmat=R, factors=2, rotation="varimax", n.obs=n)
af2$loadings # matrice dei pesi fattoriali
af2$rotmat # matrice di rotazione

<pre>matrice dei pesi fattorialiaf=factanal(X, factors=2, rotation = "varimax", scores="regression")af$loadings</pre> <pre>#Stimare i punteggi fattoriali con il metodo di Thompsonpuntfat=af$scores</pre> <pre>#Calcolo la stima iniziale della comunalità come hi^2=max(rij)R0=R-diag(rep(1,p)) # sostituisco 1 sulla diagonale di R con 0h2=apply(abs(R0),2,max)Psi=diag(1-h2) # calcolo la stima di PsiRstar=R-Psi # calcolo la matrice di correlazione ridotta R*"life"=structure(.Data = list(c(63., 34., 38., 59., 56., 62., 50., 65., 56., 69., 65., 64., 56., 60., 61., 49., 59., 63., 59., 65., 65., 64.,64., 67.,61., 68., 67., 65., 59., 58., 57.), c(51., 29., 30., 42., 38., 44., 39., 44., 46., 47., 48., 50., 44., 44., 45., 40., 42., 44., 44., 48., 48., 63.,43.,45., 40., 46., 45., 46., 43., 44., 46.), c(30., 13., 17., 20., 18., 24., 20., 22., 24., 24., 26., 28., 25., 22., 22., 22., 22., 23., 24., 28., 26.,21.,21., 23., 21., 23., 23., 24., 23., 24., 28.), c(13., 5., 7., 6., 7., 7., 7., 7., 11., 8.,</pre>

9., 11., 10., 6., 8., 9., 6., 8., 8., 14., 9., 7., 6., 8., 10.,8.,8., 9., 10., 9., 9.), 
c(67., 38., 38., 64., 62., 69., 55., 72., 63., 75., 68., 66., 61., 65., 65., 51., 61., 67., 63., 68., 67., 68.,68., 74., 67.,75., 74., 71., 66., 62., 60.), 
c(54., 32., 34., 46., 46., 50., 43., 50., 54., 53., 50., 51., 48., 45., 49., 41., 43., 48., 46., 51., 49., 47.,47., 51.,46., 52., 51., 51., 49., 47., 49.), 
c(34., 17., 20., 25., 25., 28., 23., 27., 33., 29., 27., 29., 27., 25., 27., 23., 22., 26., 25., 29., 27., 25.,24.,28., 25., 29., 28., 28., 27., 25., 28.), 
c(15., 6., 7., 8., 10., 14., 8., 9., 19., 10., 10., 11., 12., 9., 10., 8., 7., 9., 8., 13., 10., 9., 8., 10.,11.,10., 10., 10., 12., 10., 11.)), 
class = "data.frame", 
names = c("m0", "m25", "m50", "m75", "w0", "w25", "w50", "w75"),
row.names = c("Algeria", "Cameroon", "Madagascar", "Mauritius", "Reunion",

"Seychelles", "South Africa(C)", "SouthAfrica(W)","Tunisia", "Canada", "Costa Rica", "Dominican Rep", "El Salvador", "Greenland", "Grenada", "Guatemala","Honduras","Jamaica", "Mexico", "Nicaragua", "Panama", "Trinidad(62)", "Trinidad (67)","United States (66)", "United States (NW66)", "UnitedStates (W66)", "United States (67)", "Argentina","Chile", "Columbia", "Ecuador"))#test per il numero di fattori incorporati nell'approccio alla massima verosimiglianzasapply(1:3, function(f) factanal(life, factors = f, method="mle")$PVAL) #sono suggeriti 3 fattorifactanal(life, factors = 3, method="mle")#Le stime dei punteggi fattoriali secondo il metodo di Thomson

sono:scores=factanal(life,factors=3,method="mle",scores="regression")$scores #ESAME 6/2/19 library(datasets) state.centerstate.x77 summary(prcomp(state.x77)) #analisi in componenti principali n=nrow(state.x77) p=ncol(state.x77) round(apply(state.x77,2,var)*((n-1)/n),1) #varianze delle variabili presenti arrotondate al 1° decimale c=mean(eigen(cor(state.x77))$values/p) #media delle % di varianza spiegata da ciascuna componente principale sum(eigen(cor(state.x77))$values/p > c) #numero di componenti principali con varianza spiegata superiore a c #Sia Q la matrice dei dati standardizzati ottenuta a partire da Y . Si svolga l'analisi delle c