function [theta,MSE,pval]=regression(M,Y,alpha)
  // regression lineaire avec table d'anova
  // M = matrice (n,p) des regresseurs, Y = observations (n,1)
  // theta = estimateur des parametres
  // MSE   = Mean Square Error estimateur de la variance
  // pval    = p-valeur du test de Ho : "pas de regresseurs"
  [n,p]=size(M);
  theta = inv(M'*M)*M'*Y; // estimateur des parametres
  // calcul des sommes de carres
  YE = M*theta;                 // projete sur E
  SSE = sum((Y - YE).^2);       // ||Y-YE||^2
  MSE = SSE/(n-p);
  printf("\n");
  printf("Estimateur de theta:\n");
  printf('theta(%d)=%f\n', [1:p]',theta);
  printf("Estimateur de sigma2: %f\n", MSE);
  SSM = sum((YE - mean(Y)).^2); // ||YE-YH||^2
  MSM = SSM/(p-1);
  F = MSM/MSE;                      // Fisher
  q = cdff("F",p-1,n-p,1-alpha,alpha); // quantile d'ordre 1-alpha
  [pp,pval] = cdff("PQ",F,p-1,n-p);  // p-valeur
  // Affichage de la table d'anova
  printf("\n");
  printf("TABLE D''ANALYSE DE LA VARIANCE\n");
  printf("Source   SS      DF    MS      Fisher  p-valeur\n");
  printf("-----------------------------------------------\n");
  printf("Model  %6.1f   %3d  %6.1f    %4.1f    %f\n",...
		 SSM,p-1,MSM,F,pval);
  printf("Error  %6.1f   %3d  %6.1f\n",SSE,n-p,MSE);
  printf("-----------------------------------------------\n");
  printf("Quantile de F(%d,%d) d''ordre %f : %f\n",p-1,n-p,1-alpha,q);
endfunction;