// Pierre L. Douillet
// www.douillet.info
// module stats
// version 07
// 2007-12-17

if MSDOS then 
  cd "c:/MapleV6/Users/Douillet/stats" ;
  mydir=pwd(); myfig=ans; 
else
  cd "~/docs/Ensait/stats" ;
  mydir=pwd()+"/datas/"
  myfig=pwd()+"/figures/"
end//if
ident='STATS';

global xx yy funy bord rep datas

//-------------------------------------------------------------------------------

function ceram
  nam = "nist-ceramic.txt"; 	shortnam="ceram";
  fich = mydir + nam ;
  hndl = mopen(fich, "r") ; rep = mgetl(hndl) ; mclose(hndl) ;
  datas = msscanf(-1, rep, "%s%s%s%s%s%s") ;
  lesx=eval(datas(2:$,$));
  curfig=scf(0); clf(); histplot (12, lesx)
  xtitle(sprintf("%d valeurs issues de ceramis-nist", size(lesx,'*')))
  curax=gca();
  curh=curax.children(1).children;
  curh.background=4;
  curh.foreground=2;
  curh.polyline_style=6;
  curh.fill_mode='on';
  curh.thickness=2;
  mea= mean(lesx);
  var=variance(lesx);
  hh=curax.data_bounds(2,2);
  plot2d([1,1]*mea,[0,hh*1.05], style=1)
  plot2d([1,1]*mea-sqrt(var),[0,hh*1.05], style=1)
  plot2d([1,1]*mea+sqrt(var),[0,hh*1.05], style=1)
  curax.sub_ticks=[0,0];
  curax.margins=[0.08,0.03,0.12,0.1] ;
  curax.font_size=4 ;
  curax.title.font_size=4 ;
  try
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(0, myfig+shortnam+"_histo")
  catch; end;
endfunction

//-------------------------------------------------------------------------------

function exo_regress
global xx yy funy bord
  xx=[4.11,5.73,5.47,5.16,2.44,6.98,2.94,4.34,2.47,6.16];
  yy=[12.6,12.8,11.2,11.9,9.58,16.9,7.59,12.5,8.14,16.8];
  n=size(xx,'*'); fre= eye(n,n);
  mx=mean(xx); vx=covar(xx,xx,fre);
  my=mean(yy); vy=covar(yy,yy,fre);
  cv=covar(xx,yy,fre);
  deff('y=funy(x)', sprintf("y=%16.12f*(x-%16.12f)+%16.12f",cv/vx,mx,my))
  dy= feval(xx,funy)-yy;
  vardy=dy*dy'/n;
  frv=vy/vardy;
  printf("\n(y-%6.3f)=%6.4f*(x-%6.3f) ; frv=%f\n",my,cv/vx,mx, frv)
  deff('y=bord(x)',   sprintf("y=%16.12f*sqrt(1+(x-%16.12f)^2/%16.12f)", sqrt(vardy*n/(n-2)),mx,vx ))
  deff('y=borsup(x)', sprintf("y=funy(x)+%16.12f*sqrt(1+(x-%16.12f)^2/%16.12f)", sqrt(vardy*n/(n-2)),mx,vx ))
  deff('y=borinf(x)', sprintf("y=funy(x)-%16.12f*sqrt(1+(x-%16.12f)^2/%16.12f)", sqrt(vardy*n/(n-2)),mx,vx ))

  curfig=scf(0); clf();
  plot(xx,yy,'o')
  lesxx=min(xx):(max(xx)-min(xx))/20:max(xx);
  lesyy=feval(lesxx,funy); lesdd=feval(lesxx,bord)
  plot2d(lesxx, lesyy, 1)
  plot2d(lesxx, lesyy+lesdd,5); plot2d(lesxx, lesyy-lesdd,5) ;

  curax=gca(); 
  plot(curax.data_bounds(:,1), my*[1;1])
  plot(mx*[1;1], curax.data_bounds(:,2))
  curax.sub_ticks=[0,0];
  xtitle(sprintf("\n(y -%6.3f) = %6.4f*(x -%6.3f) ; frv=%f\n",my,cv/vx,mx, frv))
  curax.title.font_size=4;
  
  try
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(0, myfig+"regression")
  catch; end;
endfunction

//-------------------------------------------------------------------------------

global N dd vth1 vth2 

function initdes(N_)
global N dd
  grand('setsd',300)
  N=N_
  function y=dd(S)
    y=sum(grand(1,S,'uin',1,6))
  endfunction
endfunction

function unseulde
  global vth1
  S=1; vth1=ones(1:6);
  lancers=feval(S*ones(1,N), dd); // ligne pour boucle for
  mu=mean(lancers) ; va=variance(lancers) ;
  expr=zeros(1:6*S); for j=lancers do expr(j)=expr(j)+1; end; // effectifs réels
  theo=vth1/6^S*N; // effectifs prévisionnels
  commesi=[]; for j=S:S*6 do commesi=[commesi, j*ones(1,vth1(j))]; end
  printf('%d lancers de %d dés ; moy=%9.6f, var=%9.6f \n',N,S,mu,va)

  curfig=scf(S); clf();
  bornes=S-0.5:1:S*6+0.5;
  histplot(bornes, lancers)
  histplot(bornes, commesi)
  xtitle(sprintf("%d lancers d''un seul dé", N))
  mise_en_forme(mu,va);
  try
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(S, myfig+"somme_un_de")
  catch; end;
endfunction//unseulde

function mise_en_forme(mu,va)
  curax=gca();
  curax.sub_ticks=[0,0];
  curax.font_size=4 ;
  curax.title.font_size=4 ;
  curpoly=curax.children.children(2);
  curpoly.fill_mode="on"; curpoly.background=12; // comme il est dessous, il ne masque pas l'autre 
  hh=curax.data_bounds(2,2);
  plot2d([1,1]*mu,[0,hh*1.05], style=5)
  plot2d([1,1]*mu-sqrt(va),[0,hh*0.95], style=5)
  plot2d([1,1]*mu+sqrt(va),[0,hh*0.95], style=5)
endfunction

function deuxdes
  global vth2
  S=2; vth2=zeros(1:6*S);
  for i=1:6 do for j=1:size(vth1,'c') do
    vth2(i+j)=vth2(i+j)+vth1(j);
  end; end;
  lancers=feval(S*ones(1,N), dd); // ligne pour boucle for
  mu=mean(lancers) ; va=variance(lancers) ;
  expr=zeros(1:6*S); for j=lancers do expr(j)=expr(j)+1; end;
  theo=vth2/6^S*N;
  commesi=[]; for j=S:S*6 do commesi=[commesi, j*ones(1,vth2(j))]; end
  printf('%d lancers de %d dés ; moy=%9.6f, var=%9.6f \n',N,S,mu,va)

  curfig=scf(S); clf();
  bornes=S-0.5:1:S*6+0.5;
  histplot(bornes, lancers)
  histplot(bornes, commesi)
  xtitle(sprintf("%d lancers de deux dés", N))
  mise_en_forme(mu,va);
  try
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(S, myfig+"somme_deux_des")
  catch; end;
endfunction

function troisdes
  S=3; vth3=zeros(1:6*S);
  for i=1:6 do for j=1:size(vth2,'c') do
    vth3(i+j)=vth3(i+j)+vth2(j);
  end; end;

  lancers=feval(S*ones(1,N), dd); // ligne pour boucle for
  mu=mean(lancers) ; va=variance(lancers) ;
  expr=zeros(1:6*S); for j=lancers do expr(j)=expr(j)+1; end;
  theo=vth3/6^S*N;
  commesi=[]; for j=S:S*6 do commesi=[commesi, j*ones(1,vth3(j))]; end
  printf('%d lancers de %d dés ; moy=%9.6f, var=%9.6f \n',N,S,mu,va)

  curfig=scf(S); clf();
  bornes=S-0.5:1:S*6+0.5;
  histplot(bornes, lancers)
  histplot(bornes, commesi)
  xtitle(sprintf("%d lancers de trois dés", N))
  mise_en_forme(mu,va);
  try
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(S, myfig+"somme_trois_des")
  catch; end;
endfunction

function pearson=sumgamma
  n=1000; s=2.4; t=3.1;
  xx=grand(1,n,'gam',s,1); yy=grand(1,n,'gam',t,1); zz=xx+yy;
  curfig=scf(1); clf(); histplot(10,xx)
  deff('d=fx(x)',sprintf('d=x^%18.16f*exp(-x)*%18.16f', s-1, 1/gamma(s)))
  lesx=0:0.2:10; plot2d(lesx, feval(lesx, fx),5)
  xtitle('les X')

  curfig=scf(2); clf(); histplot(10,yy)
  deff('d=fy(x)',sprintf('d=x^%18.16f*exp(-x)*%18.16f', t-1, 1/gamma(t)))
  lesx=0:0.2:12; plot2d(lesx, feval(lesx, fy),5)
  xtitle('les Y')

  curfig=scf(3); clf();
  numbar=10; borp=[0:1/numbar:1-1/numbar,0.99999]; bor= cdfgam('X',(s+t)*ones(borp),ones(borp),borp,1-borp);
  histplot(bor,zz)
  deff('d=fz(x)',sprintf('d=x^%18.16f*exp(-x)*%18.16f', s+t-1, 1/gamma(s+t)))
  lesx=0:0.2:16; plot2d(lesx, feval(lesx, fz),5)
  xtitle('les Z')
  cumul=zeros(1:numbar+1);
  for j=2:numbar+1 do cumul(j)= size(find(zz<bor(j)),'*'); end;
  exper=cumul(2:$)-cumul(1:$-1); theor=ones(exper)*n/numbar;
  // printf("chi2=%f', sum((exper-theor)^2 ./ theor))
  pearson=sum((exper-theor)^2 ./ theor)
endfunction

function drawgamma(a,b)
  deff('d=fn(x)',sprintf('d=x^%18.16f*exp(x*%18.16f)*%18.16f', a-1, -1/b, 1/gamma(a)/b^a))
  curfig=scf(5); clf();
  lesx=[0:0.01:1]*a*b*4; lesy=feval(lesx, fn); lesY=(lesx .* lesy)/(a*b)
  plot2d(lesx,lesy,1)
  plot2d(lesx,lesY,1)
  xtitle(sprintf('loi gamma en nombre et en masse, forme=%4.2f, échelle=%4.2f',a,b))
  
  curax=gca();
  hh=curax.data_bounds(2,2);
  plot2d([1,1]*a*b,[0,hh*1.05], style=1)
  curax.sub_ticks=[0,0];
  curax.title.font_size=4 ; curax.title.font_style=8 ;
  try
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(5, myfig+"drawgamma")
  catch ; end;

endfunction



function n=menusommedes
  n=x_choose(['initdes';'unseuldé';'deuxdes';'troisdes';'fin'],[ident+' somme: ';'choisir un item'])
  select n
  case 1 then txt=['N']; sig=x_mdialog('enter N',txt,['500']); initdes(evstr(sig)); menusommedes()
  case 2 then unseulde() ; menusommedes()
  case 3 then deuxdes() ; menusommedes()
  case 4 then troisdes() ; menusommedes()
  case 5 then n=-1
  end//select
endfunction



function menu
  n=x_choose([ 'histogrammes : ceram.txt'; 'exo_regress';'somme de plusieurs dés'
    'tracé loi gamma'; 'somme de deux lois Gamma réduites'; 'fermer les figures'],[ident+': ';'choisir un item'])
  select n
  case 1 then ceram(); menu()
  case 2 then exo_regress(); menu()
  case 3 then nn=menusommedes(); if nn==-1 then break else menu() ; end
  case 4 then sig=x_mdialog('loi gamma', ['a','b'],['2.5','10']);
    a=max(0,evstr(sig(1))); b=max(0,evstr(sig(2))); 
    drawgamma(a,b) ; menu()
  case 5 then sumgamma() ; menu()
  case 6 then for j=1:0 do close(j); end; // ne fonctionne pas
  end//select
endfunction

function football
global rep datas
  col1=3; col2=4;
  nam = "football.txt"; 	shortnam="foot";
  fich = mydir + nam ;
  hndl = mopen(fich, "r") ; rep = mgetl(hndl) ; mclose(hndl) ;
  rep=strsubst(rep,","," ")
  datas = msscanf(-1, rep, "%s%s%s%s%s") ;
  xx=eval(datas(2:$,col1));
  yy=eval(datas(2:$,col2));

  n=size(xx,'*'); fre= eye(n,n);
  mx=mean(xx); vx=covar(xx,xx,fre);
  my=mean(yy); vy=covar(yy,yy,fre);
  cv=covar(xx,yy,fre);
  deff('y=funy(x)', sprintf("y=%16.12f*(x-%16.12f)+%16.12f",cv/vx,mx,my))
  dy= feval(xx,funy)-yy;
  vardy=dy'*dy/n; //caveat : vecteurs colonne !!!
  frv=vy/vardy;
  printf("\n(y-%6.3f)=%6.4f*(x-%6.3f) ; frv=%f\n",my,cv/vx,mx, frv)
  deff('y=bord(x)',   sprintf("y=%16.12f*sqrt(1+(x-%16.12f)^2/%16.12f)", sqrt(vardy*n/(n-2)),mx,vx ))
  deff('y=borsup(x)', sprintf("y=funy(x)+%16.12f*sqrt(1+(x-%16.12f)^2/%16.12f)", sqrt(vardy*n/(n-2)),mx,vx ))
  deff('y=borinf(x)', sprintf("y=funy(x)-%16.12f*sqrt(1+(x-%16.12f)^2/%16.12f)", sqrt(vardy*n/(n-2)),mx,vx ))

  curfig=scf(0); clf();
  plot(xx,yy,'o')
  lesxx=min(xx):(max(xx)-min(xx))/20:max(xx);
  lesyy=feval(lesxx,funy); lesdd=feval(lesxx,bord)
  plot2d(lesxx, lesyy, 1)
  plot2d(lesxx, lesyy+lesdd,5); plot2d(lesxx, lesyy-lesdd,5) ;

  curax=gca(); 
  plot(curax.data_bounds(:,1), my*[1;1])
  plot(mx*[1;1], curax.data_bounds(:,2))
  curax.sub_ticks=[0,0];
  xtitle(sprintf("football : [%s,%s]\n(y -%6.3f) = %6.4f*(x -%6.3f) ; frv=%f\n",datas(1,col1), datas(1,col2),my,cv/vx,mx, frv))
  curax.title.font_size=4;
  
  try
    set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2))
    xs2eps(0, myfig+shortnam+"-h-w")
  catch; end;


endfunction

menu
return


xx=[10;20;30;40]
yy=[15,25,35,45]
fre=[20,4,1,0;10,36,9,0;0,5,10,0;0,0,0,5]
N=sum(fre)
fy=sum(fre,'r'), fx=sum(fre,'c')
mx=sum(fx .* xx)/N, my=sum(fy .* yy)/N
vx=sum(fx .* xx^2)/N-mx^2, vy=sum(fy .* yy^2)/N-my^2
cv=covar(xx,yy,fre)