// e donne l'effort dans la promotion du p lanning familial et c le
// changement (la baisse) )du taux de natalite dans 20 pays sud
// americains dans les annees 1965-75
function [e,c]=donnees()
  e=[0 0 16 16 21 15 14 6 13 9 3 7 23 4 0 19 3 0 15 7];
  c=[1 10 29 25 29 40 21 0 13 4 0 7 21 9 7 22 6 2 29 11];
endfunction

function points()
  xbasc();
  [e,c]=donnees();
  plot2d(e,c,-4,rect=[-1,-1,24,41]);
  xtitle('Le nuage de points');
  legends(['donnees'],[-4],2);
endfunction

function [r,b,a]=lineaire()
  [e,c]=donnees();
  n=length(e);
  em=sum(e)/n;
  cm=sum(c)/n;
  b=((e-em)*(c-cm)')/((e-em)*(e-em)');
  a=cm-b*em;
  r=((e-em)*(c-cm)')/sqrt(((e-em)*(e-em)')*((c-cm)*(c-cm)'));
  xbasc();
  plot2d(e,c,-4,rect=[min(e)-1,min(c)-1,max(e)+1,max(c)+1]);
  plot2d(em,cm,-3);
  ra=[min(e)-1:0.1:max(e)+1]';
  plot2d(ra,a+b*ra,5);
  xtitle('La droite de regression');
  legends(['donnees','centre de gravite','droite de regression'],[-4,-3,5],2);
endfunction

function [r,b,a]=lineairegaussien()
  [e,c]=donnees();
  n=length(e);
  em=sum(e)/n;
  cm=sum(c)/n;
  b=((e-em)*(c-cm)')/((e-em)*(e-em)');
  a=cm-b*em;
  r=((e-em)*(c-cm)')/sqrt(((e-em)*(e-em)')*((c-cm)*(c-cm)'));
  s=sqrt((c-a-b*e)*(c-a-b*e)'/(n-2));
  t=cdft('T',n-2,0.975,0.025);
  xbasc();
  plot2d(e,c,-4,rect=[min(e)-1,min(c)-1,max(e)+1,max(c)+1]);
  ra=[min(e)-1:0.1:max(e)+1]';
  plot2d(ra,[a+b*ra a+b*ra+t*s*sqrt(1/n+(ra-em).^2/((e-em)*(e-em)')) ...
	     a+b*ra-t*s*sqrt(1/n+(ra-em).^2/((e-em)*(e-em)')) a+b*ra+t*s*sqrt(1+1/n+(ra-em).^2/((e-em)*(e-em)')) a+b*ra-t*s*sqrt(1+1/n+(ra-em).^2/((e-em)*(e-em)'))],[5 2 2 3 3]);
  xtitle('Droite de regression avec intervalles de confiance et de"+...
	 " prediction');
  legends(['donnees','droite de regression','intervalle de confiance', ...
	   'intervalle de prediction'],[-4,5,2,3],2);
endfunction


//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù


// La meme chose sauf que je tire des donnees gaussiennes
function [r,b,a,s]=lineairegaussientheo(n,al,be,si)
  e=rand(1,n);
  eps=grand(1,n,'nor',0,si);
  c=al+be*e+eps;
   em=sum(e)/n;
  cm=sum(c)/n;
  b=((e-em)*(c-cm)')/((e-em)*(e-em)');
  a=cm-b*em;
  r=((e-em)*(c-cm)')/sqrt(((e-em)*(e-em)')*((c-cm)*(c-cm)'));
  s=sqrt((c-a-b*e)*(c-a-b*e)'/(n-2));
  t=cdft('T',n-2,0.975,0.025);
  xbasc();
  plot2d(e,c,-4,rect=[0,min(c)-.1,1,max(c)+.1]);
  ra=[0:0.05:1]';
  plot2d(ra,[a+b*ra al+be*ra a+b*ra+t*s*sqrt(1/n+(ra-em).^2/((e-em)*(e-em)')) ...
	     a+b*ra-t*s*sqrt(1/n+(ra-em).^2/((e-em)*(e-em)')) a+b*ra+t* ...
	     s*sqrt(1+1/n+(ra-em).^2/((e-em)*(e-em)')) a+b*ra-t*s*sqrt(1+1/n+(ra-em).^2/((e-em)*(e-em)'))],[5 1 2 2 3 3]);
  xtitle('Droite de regression avec intervalles de confiance et de"+...
	 " prediction');
  legends(['donnees','droite de regression','droite des moyennes','intervalle de confiance', ...
	   'intervalle de prediction'],[-4,5,1,2,3],2);
endfunction


// NOUVEAU : Le test de Fisher pour la validité de l'explication linéaire (ou beta=0).

function [P]=testfisher(al,be,si,x,N,th)
  n=length(x);
  xm=sum(x)/n;
  vx=(x-xm)*(x-xm)';
  E=grand(N,n,'nor',0,si);
  Y=al+be*ones(N,1)*x+E;
  Ym=sum(Y,2)/n;
  B=Y*((x-xm)')/vx;
  Ystar=Ym*ones(1,n)+B*(x-xm);
  S2=diag((Y-Ystar)*(Y-Ystar)')/(n-2);
  tc=((B.^2)./S2)*vx;
  P=sum(tc>=th)/N;
endfunction


// La fonction de puissance du test de dépendance linéaire pour des variables gaussiennes en fonction de beta
function fonctionpuissance(bmax,N)
  x=rand(1,30);
  b=[-bmax:0.05:bmax];
  j=length(b);
  P=[];
  for i=1:j,
    P=[P testfisher(0,b(i),0.2,x,N,4.2)];
  end
  xbasc();
  plot2d(b',[P' 0.05*ones(b')],[5 1]);
  xtitle('La puissance du test d''explication lineaire en fonction de la valeur de b');
  legends(['puissance' 'seuil'],[5 1],2); 
endfunction


// Test de dépendance linéaire sur mon exemple
function [tc]=fisherexemple()
  [x,Y]=donnees();
  n=length(x);
  xm=sum(x)/n;
  vx=(x-xm)*(x-xm)'; 
  Ym=sum(Y)/n;
  B=Y*((x-xm)')/vx;
  Ystar=Ym*ones(1,n)+B*(x-xm);
  S2=(Y-Ystar)*(Y-Ystar)'/(n-2);
  tc=((B^2)/S2)*vx;
endfunction
// On obtient 32.2 à comparer aux valeurs de rejet pour une Fisher(1,18), qui sont 4.41 au niveau 0.05 et 8.29 au niveau 0.01. On rejette donc H0, ce qui veut dire qu'il n'est pas déraisonnable de donner une explication linéaire des Y par les x...