// version du 10/11/12

// Display warning for floating point exception
ieee(1)
// calcul l'incertitude sur les coefficients d'une regression lineaire
// de facon statistique
// les donnees sont fournies par le tableau resp1.dat
// la table de student est fournie dans le tableau TableStudent.dat

clear;
clf;

// lecture du fichier de donnees
resp=read('resp1.dat',-1,2);
// lecture de la table de Student
// student=read('student.dat',-1,1);

// affichage des points experimentaux
scf(1);
clf;
a=get("current_axes");
a.axes_visible="on";
a.x_location="bottom";
a.x_label.text="x";
a.y_label.text="y";
plot2d(resp(:,1),resp(:,2),-1);


// selection des points
taille=size(resp);
nbtest = taille(1)-1;

for test = 1:nbtest-1;

  // np1 np2 : les bornes des points consideres
  np1 = 1;
  np2 = 2+test;

  val1 = resp(np1:np2,1);
  val2 = resp(np1:np2,2);
  nb = size(mtlb_double(val1),1);
  // calcul des valeurs de la droite des moindres carres
  for n = 1:nb;
    v11(1,n) = mtlb_double(val1(n,1))^2;  v22(1,n) = mtlb_double(val2(n,1))^2;  v12(1,n) = mtlb_double(val1(n,1))*mtlb_double(val2(n,1));
  end;
  denomi = mtlb_s(nb*mtlb_sum(v11),sum(mtlb_double(val1),firstnonsingleton(mtlb_double(val1)))^2);
  // pente
  p = mtlb_s(nb*mtlb_sum(v12),sum(mtlb_double(val1),firstnonsingleton(mtlb_double(val1)))*sum(mtlb_double(val2),firstnonsingleton(mtlb_double(val2))))/denomi;
  // ordonnee a  l'origine
  c = mtlb_s(sum(mtlb_double(val2),firstnonsingleton(mtlb_double(val2)))*mtlb_sum(v11),mtlb_sum(v12)*sum(mtlb_double(val1),firstnonsingleton(mtlb_double(val1))))/denomi;

  //incertitudes sur ces coefficients
  // residus
  for n = 1:nb;
    res(1,n) = mtlb_s(mtlb_double(mtlb_e(val2,n)),mtlb_a(p*mtlb_double(mtlb_e(val1,n)),c));
    res2(1,n) = res(n)^2;
  end;
  // ecart type estime ^2
  sige2 = mtlb_sum(res2)/(nb-2);
  // variances sur p et c ^2
  sigp2 = (sige2*nb)/denomi;
  sigc2 = (sige2*mtlb_sum(v11))/denomi;
  // incertitudes
  if n>2 then
      deltap = cdft("T",n-2,0.025,0.975)*sqrt(sigp2);
      deltac = cdft("T",n-2,0.025,0.975)*sqrt(sigc2);
  else
    deltap = 0;  deltac = 0;
  end;
  // affichage des resultats
  aff0(test)=np2;
  aff(test,1) = p;
  aff(test,2) = p-deltap;
  aff(test,3) = p+deltap;
  aff1(test,1) = c;
  aff1(test,2) = c-deltac;
  aff1(test,3) = c+deltac;
end;
'evolution de la pente avec le nombre de points'
aff
'evolution de l ordonnee a l origine avec le nombre de points'
aff1

// affichage de l'evolution de la pente 
// avec l'augmentation du nombre de points consideres
scf(2);
clf;
a=get("current_axes");
a.axes_visible="on";
a.x_location="bottom";
a.x_label.text="nb points";
a.y_label.text="pentes";
plot2d(aff0,aff(:,1),style=1);
plot2d(aff0,aff(:,2),style=2);
plot2d(aff0,aff(:,3),style=2);
//
// affichage de l'evolution de l'ordonnee a l'origine 
// avec l'augmentation du nombre de points consideres
scf(3);
clf;
a=get("current_axes");
a.axes_visible="on";
a.x_location="bottom";
a.x_label.text="nb points";
a.y_label.text="ordonnee a l origine";
plot2d(aff0,aff1(:,1),style=1);
plot2d(aff0,aff1(:,2),style=2);
plot2d(aff0,aff1(:,3),style=2);