# Ceci requiert le package MASS (de VR).
# 	prcomp utilise var(x)=W/(n-1);
# 	princomp utilise ((n-1)/n)*var(x)=W/n;

par(ask=TRUE);

data(iris3);
#print(iris3);
ir <- rbind(iris3[,,1],iris3[,,2],iris3[,,3]);
lir=log(ir);
n=dim(lir)[1];
lir.species <- factor(c(rep("s",50),rep("c",50),rep("v",50)));
lir.code <- c(rep(2,50),rep(3,50),rep(4,50));

dflir=data.frame(lir);
boxplot(dflir);
pairs(dflir,asp=1,col=lir.code);
pairs(dflir,asp=1,col=lir.code,xlim=c(-2.31, 2.07), ylim=c(-2.31, 2.07));

lir.pca <- princomp(lir,cor=TRUE,scores=TRUE);
print(summary(lir.pca));
#plot(lir.pca);
screeplot(lir.pca,type = "lines");
print(loadings(lir.pca));
lir.pc <- predict(lir.pca);

boxplot(data.frame(lir.pc));
pairs(lir.pc,asp=1,col=lir.code,xlim=c(-3.34, 3.05), ylim=c(-3.34, 2.05));

eqscplot(lir.pc[,1:2],type="n",xlab="1ere composante principale",ylab="2eme composante principale");
text(lir.pc[,1:2],labels=as.character(lir.species),col=lir.code);

eqscplot(lir.pc[,3:4],type="n",xlab="3eme composante principale",ylab="4eme composante principale");
text(lir.pc[,3:4],labels=as.character(lir.species),col=lir.code);

#############################################

# Calcul manuel:

h=3;
if (h==1) { # utilise W/n;
	eig=eigen((n-1)/n*var(lir));
	}
if (h==2) { # utilise W/(n-1);
	eig=eigen(var(lir));
	}
if (h==3) { # utilise CORR;
	eig=eigen(cov2cor(var(lir)));
	}
lambda=eig$values;
O=eig$vectors;
print(lambda);
print(sqrt(lambda));
print(O);

#s=lir.pca$scores;
#xbar=lir.pca$center;
#scale=lir.pca$scale;
#print(s);
#i=2;
#print(s[i,]);
#print(t(O)%*%(diag(1/scale)%*%(lir[i,]-xbar)));