d=10;
rho=0.5;

// Fabrique une matrice avec 1 sur la diagonale
// et rho en dehors de la diagonale

Rho = (1 - rho) *eye(d,d) + rho * ones(d,d)

// spec calcule les valeurs propre de Gamma (son spectre)
// Toutes les valeurs propres sont positives lorsque rho=0.50

spec(Rho)

// ce n'est pas le cas lorsque  rho=-0.20 ! Pourquoi ?
rho=-0.2;
Rho1 = (1 - rho) *eye(d,d) + rho * ones(d,d);
spec(Rho1)

// fabrique la matrice de terme general 
//      sigma_i * Rho_ij * sigma_j 
sigma = 0.3 * ones(1,d); 
Gamma =    ********* ... ********** ;

// Calcul d'une racine carrée de cette matrice
Sigma = sqroot(Gamma);

// On verifie, on sait jamais, on doit obtenir un truc petit ...
norm(Sigma*Sigma' - Gamma)

// *Attention* Sigma = sqroot(Gamma)
// n'est carrée que si Gamma est inversible.
// Lorsque la matrice Gamma est non inversible (ex: rho = 1).
// le "sqroot" de scilab renvoie une matrice dxn (n<d) au lieu de dxd
// ou n est le rang de Gamma

rho=1;
Rho2 = (1 - rho) *eye(d,d) + rho * ones(d,d);
spec(Rho2)
Gamma2 = diag(sigma) * Rho2 * diag(sigma);

// size(Sigma) donne la taille de Sigma2
Sigma2=sqroot(Gamma2);
size(Sigma2) // Sigma2 n'est pas carrée

// Mais, on a toujours Sigma2*Sigma2' = Gamma2
norm(Sigma2*Sigma2' - Gamma2)