# Datos de atletismo:

records.pista <- data.frame(read.table("Documents/CIMATOLD/CLASSES/STAT-MULTIVAR/Data sets/T1-9.txt",header=F));
names(records.pista) <- c("100m(s)","200m(s)","400m(s)","800m(min)", "1500m(min)", "3000m(min)","Maraton","Pais");
print(records.pista);

attach(records.pista);
Xnp <- as.matrix(records.pista[,1:7]);
Paises <- records.pista[,8];
detach("records.pista");

dimX <- dim(Xnp);
n    <- dimX[1];
p    <- dimX[2];
S    <- var(Xnp);
Dm12       <- matrix(0,nrow=p,ncol=p);
diag(Dm12) <- 1/sqrt(diag(S));
R          <- Dm12%*%S%*%Dm12;

x.bar      <- apply(Xnp,2,mean);
Znp        <- matrix(0,nrow=n,ncol=p);
for(i in 1:n){
	Znp[i,] <- Dm12%*%(Xnp[i,]-x.bar);	
}

eigs.R <- eigen(R);
eigs.S <- eigen(S);



princomp.scoresR <- matrix(0,nrow=n,ncol=p); # Componentes principales ('scores')
corr.pcompi.varkR <- matrix(0,nrow=p,ncol=p);
prop.varsR        <- rep(0,p);

princomp.scoresS <- matrix(0,nrow=n,ncol=p); # Componentes principales ('scores')
corr.pcompi.varkS <- matrix(0,nrow=p,ncol=p);
prop.varsS        <- rep(sum(eigs.S$values),p);
diagS             <- diag(S);

for(i in 1:p){

		princomp.scoresR[,i]  <- Znp%*%eigs.R$vectors[,i];
		corr.pcompi.varkR[,i] <- sqrt(eigs.R$values[i])*eigs.R$vectors[,i];
		prop.varsR[i]         <- eigs.R$values[i]/p;

		princomp.scoresS[,i]  <- Xnp%*%eigs.S$vectors[,i];
		corr.pcompi.varkS[,i] <- (sqrt(eigs.S$values[i])*eigs.S$vectors[,i])/sqrt(diagS);
		prop.varsS[i]         <- eigs.S$values[i]/prop.varsS[i];

}


# Colocando nombres a filas y columnas de las matrices resultantes:
names.corrs <- list(names(records.pista)[1:p],c("PrinComp 1","PrinComp 2", "PrinComp 3", "PrinComp 4","PrinComp 5","PrinComp 6","PrinComp 7"));
names.scores <- list(Paises,c("PrinComp 1","PrinComp 2", "PrinComp 3", "PrinComp 4","PrinComp 5","PrinComp 6","PrinComp 7"));

dimnames(corr.pcompi.varkR) <- names.corrs;
dimnames(corr.pcompi.varkS) <- names.corrs;

dimnames(princomp.scoresR) <- names.scores;
dimnames(princomp.scoresS) <- names.scores;

pca.atletismo <- princomp(covmat=R,scores=TRUE);
loadings.pca(pca.atletismo) # Deberia ser igual a matriz 'corr.pcompi.varkR', pero no lo es porque el paquete re-escala los
								# vectores en cada columna para que estos tengan norma =1 

# Para obtener los resultados del paquete, hagan:

l1 <- matrix(corr.pcompi.varkR[,1],nrow=7,ncol=1);
corr.pcompi.varkR[,1]/sqrt(t(l1)%*%l1)