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

par(ask=TRUE);

data(crabs);
lcrabs <- log(crabs[,4:8]);
#print(lcrabs);
n=dim(lcrabs)[1];

crabs.grp <- factor(c("B","b","O","o")[rep(1:4,each=50)]);
crabs.codegrp <- c("Blue","lightblue","orangered1","orange")[rep(1:4,each=50)];

boxplot(lcrabs);
pairs(lcrabs,asp=1,col=crabs.codegrp);
pairs(lcrabs,asp=1,col=crabs.codegrp,xlim=c(1.80, 4.01), ylim=c(1.80, 4.01));

lcrabs.pca <- princomp(lcrabs,scores=TRUE);
print(summary(lcrabs.pca));
#plot(lcrabs.pca);
screeplot(lcrabs.pca,type = "lines");
print(loadings(lcrabs.pca));
lcrabs.pc <- predict(lcrabs.pca);
#dimnames(lcrabs.pc) <- list(NULL,paste("PC",1:5,sep=""));

boxplot(data.frame(lcrabs.pc));
pairs(lcrabs.pc,asp=1,col=crabs.codegrp,xlim=c(-.95, 1.66), ylim=c(-.95, 1.66));

hp=3;
if (hp==1) { # avec labels, orthonormé
	eqscplot(lcrabs.pc[,1:2],type="n",xlab="1ere composante 	principale",ylab="2eme composante principale");
	text(lcrabs.pc[,1:2],labels=as.character(crabs.grp),col=crabs.codegrp);
	legend("topright", c("B: B-M","b: B-F","O: O-M","o: O-F"), 	text.col=c("Blue","lightblue","orangered1","orange"));

	eqscplot(lcrabs.pc[,4:5],type="n",xlab="4eme composante 	principale",ylab="5eme composante principale");
	text(lcrabs.pc[,4:5],labels=as.character(crabs.grp),col=crabs.codegrp);
	legend("topright", c("B: B-M","b: B-F","O: O-M","o: O-F"), 	text.col=c("Blue","lightblue","orangered1","orange"));
	}
if (hp==2) { # avec labels, non orthonormé
	plot(lcrabs.pc[,1:2],type="n",xlab="1ere composante 	principale",ylab="2eme composante principale");
	text(lcrabs.pc[,1:2],labels=as.character(crabs.grp),col=crabs.codegrp);
	legend("topright", c("B: B-M","b: B-F","O: O-M","o: O-F"), 	text.col=c("Blue","lightblue","orangered1","orange"));

	plot(lcrabs.pc[,4:5],type="n",xlab="4eme composante 	principale",ylab="5eme composante principale");
	text(lcrabs.pc[,4:5],labels=as.character(crabs.grp),col=crabs.codegrp);
	legend("topright", c("B: B-M","b: B-F","O: O-M","o: O-F"), 	text.col=c("Blue","lightblue","orangered1","orange"));
	}
if (hp==3) { # sans label, orthonormé
	eqscplot(lcrabs.pc[,1:2],xlab="1ere composante principale",
	ylab="2eme composante principale",col=crabs.codegrp);
	legend("topright", c("B-M","B-F","O-M","O-F"), 	pch=1,col=c("Blue","lightblue","orangered1","orange"));

	eqscplot(lcrabs.pc[,4:5],xlab="4eme composante principale",
	ylab="5eme composante principale",col=crabs.codegrp);
	legend("topright", c("B-M","B-F","O-M","O-F"), 	pch=1,col=c("Blue","lightblue","orangered1","orange"));
	}
if (hp==4) { # sans labels, non orthonormé
	plot(lcrabs.pc[,1:2],xlab="1ere composante principale",
	ylab="2eme composante principale",col=crabs.codegrp);
	legend("topright", c("B-M","B-F","O-M","O-F"), 	pch=1,col=c("Blue","lightblue","orangered1","orange"));

	plot(lcrabs.pc[,4:5],xlab="4eme composante principale",
	ylab="5eme composante principale",col=crabs.codegrp);
	legend("topright", c("B-M","B-F","O-M","O-F"), 	pch=1,col=c("Blue","lightblue","orangered1","orange"));
	}

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

#print(splom(~ lcrabs.pc[,1:3], groups=crabs.grp, panel=panel.superpose,col=c("Blue","lightblue","orangered1","orange"), key=list(text=list(c("B-M","B-F","O-M","O-F")), points=Rows(trellis.par.get("superpose.symbol"),1:4),columns=4,col=c("Blue","lightblue","orangered1","orange"))));

sex <- crabs$sex;levels(sex)<-c("Female","Male");
sp <- crabs$sp;levels(sp)<-c("Blue","Orange");
#print(splom(~lcrabs.pc[,1:3] | sp*sex,cex=0.5,pscales=0));

#xgobi(lcrabs, colors=c("SkyBlue","SlateBlue","Orange","Red")[rep(1:4,each=50)])

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

# Calcul manuel:

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

#s=lcrabs.pca$scores;
#xbar=lcrabs.pca$center;
#scale=lcrabs.pca$scale;
##print(s);
#i=1;
#v=c(lcrabs[i,1],lcrabs[i,2],lcrabs[i,3],lcrabs[i,4],lcrabs[i,5]);
#print(s[i,]);
#print(t(O)%*%(diag(1/scale)%*%(v-xbar)));