function c=covariance(echant)
//      retourne la matrice de covariance
//      de l'echantillon echant. L'entree est une matrice 
//      a n lignes (individus) et d colonnes (variables)

moyenne = mean(echant,"r");
echcent = echant - ones(size(echant,"r"),1)*moyenne;
c=(echcent'*echcent)/size(echant,"r");

function c=correlation(echant)
//      retourne la matrice de correlation
//      de l'echantillon echant. L'entree est une matrice 
//      a n lignes (individus) et d colonnes (variables)

cov=covariance(echant);
ectinv=(1)./sqrt(diag(cov));
c=diag(ectinv)*cov*diag(ectinv);

function c=ACP_Reduite(echant);
//      Affiche les projections des individus et des variables 
//      sur les 3 premiers plans principaux de l'ACP reduite.
//      Retourne les pourcentages d'inertie expliquee.
//      L'entree est une matrice 
//      a n lignes (individus) et d colonnes (variables)

//
// Calculs de moments
//
moyenne = mean(echant,"r");            // moyennes par variable
ni = size(echant,"r");                 // nombre d'individus
nv = size(echant,"c");                 // nombre de variables
echred = echant - ones(ni,1)*moyenne;  // echantillon centre
mat_cov =(echred'*echred)/ni;          // matrice de covariance
ecti = diag((1)./sqrt(diag(mat_cov))); // inverses des ecart-types  
echred = echred*ecti;                  // echantillon reduit
mat_cor = ecti*mat_cov*ecti;           // matrice de correlations
//
// Calculs des projections sur les axes principaux
//
[D,U]=bdiag(mat_cor);                  // diagonalisation
[c,k]=sort(diag(D));                   // tri des valeurs propres
c=round((c/sum(c,"r"))*100);           // pourcentages d'inertie expliquee
axes = U(:,k);                         // reordonner les vecteurs propres
indiv = echred*axes;                   // nouvelles coordonnees  
//
// Representation graphique plan 1-2
//
e= max(abs(indiv));                    // echelle de representation
xset("window",0); xbasc();             // preparation du graphique  
xset("font",2,3);
plotframe([-e,-e,e,e],[2,10,2,10],[%f,%f],..
   ["Projection sur le plan 1-2","1er Axe","2eme Axe"]);
plot2d(indiv(:,1),indiv(:,2),0,"000"); // afficher les individus
cvar = axes([1,2],:);                  // coordonnees des variables
cvar = cvar.*.[1,0];                   // insertion de l'origine
cvar = cvar';
                                       // Afficher les variables
plot2d(cvar(:,1),cvar(:,2),5*ones(1,nv),"000")


function c=ACP_Reduite2(echant);
//      Affiche les projections des individus et des variables 
//      sur les 3 premiers plans principaux de l'ACP reduite.
//      Retourne les pourcentages d'inertie expliquee.
//      L'entree est une matrice 
//      a n lignes (individus) et d colonnes (variables)

//
// Calculs de moments
//
moyenne = mean(echant,"r");                // Moyennes par variable
n_indiv = size(echant,"r");                // Nombre d'individus
n_var = size(echant,"c");                  // Nombre de variables
echred = echant - ones(n_indiv,1)*moyenne; // Echantillon centre
mat_cov =(echred'*echred)/n_indiv;         // Matrice de covariance
ect_inv = diag((1)./sqrt(diag(mat_cov)));  // Inverses des ecart-types  
echred = echred*ect_inv;                   // Echantillon reduit
mat_cor = ect_inv*mat_cov*ect_inv;         // Matrice de correlations
//
// Calculs des projections sur les axes principaux
//
[D,U]=bdiag(mat_cor);                      // Diagonalisation
[c,k]=sort(diag(D));                       // Tri des valeurs propres
c=round((c/sum(c,"r"))*100);               // Pourcentages d'inertie expliquee
axes = U(:,k);                             // Reordonner les vecteurs propres
proj_indiv = echred*axes;                  // Nouvelles coordonnees des individus 
x_indiv = proj_indiv(:,1);                 // Coordonnees des individus 1er axe
y_indiv = proj_indiv(:,2);                 // Coordonnees des individus 2eme axe
z_indiv = proj_indiv(:,3);                 // Coordonnees des individus 3eme axe
x_var = axes(1,:);                         // Coordonnees des variables 1er axe
y_var = axes(2,:);                         // Coordonnees des variables 2eme axe
z_var = axes(3,:);                         // Coordonnees des variables 3eme axe
//
// Representation graphique plan 1-2
//
e= 5;                                      // Echelle de representation
xset("window",0);
xset("default");
xbasc();
xset("font",2,3);
plotframe([-e,-e,e,e],[2,10,2,10],[%f,%f],["Projection sur le plan 1-2","1er Axe","2eme Axe"]);
plot2d(x_indiv,y_indiv,0,"000");           // Afficher les individus
for i=1:n_var,                             // Afficher les variables par un trait de couleur
    plot2d([0,x_var(i)],[0,y_var(i)],i,"000");
end;
//
// Representation graphique plan 1-3
//
xset("window",1);
xset("default");
xbasc();
xset("font",2,3);
plotframe([-e,-e,e,e],[2,10,2,10],[%f,%f],["Projection sur le plan 1-3","1er Axe","3eme Axe"]);
plot2d(x_indiv,z_indiv,0,"000");           // Afficher les individus
for i=1:n_var,                             // Afficher les variables par un trait de couleur
    plot2d([0,x_var(i)],[0,z_var(i)],i,"000");
end;
//
// Representation graphique plan 2-3
//
xset("window",2);
xset("default");
xbasc();
xset("font",2,3);
plotframe([-e,-e,e,e],[2,10,2,10],[%f,%f],["Projection sur le plan 2-3","2eme Axe","3eme Axe"]);
plot2d(y_indiv,z_indiv,0,"000");           // Afficher les individus
for i=1:n_var,                             // Afficher les variables par un trait de couleur
    plot2d([0,y_var(i)],[0,z_var(i)],i,"000");
end;






