// Pierre L. Douillet
// www.douillet.info
// version 04
// module decis

ident='DECIS.04';
if MSDOS .. 
  then cd("c:/MapleV6/Users/Douillet/decis") ;
  else cd("~/docs/Ensait/decis/scilab") ;
end//if
baserep=pwd();
myfig=baserep+'/../figures/';

printf('\n'); printf('\n---------------------------------%s\n\n',baserep)

global N dd vth1 vth2 
global N p1 p2 p3 p4 bor lesx gap alt gen

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

  vth1=ones(1:6);
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;
  xpr2=feval(S*ones(1,N), dd);
  mu2=mean(xpr2) ; va2=variance(xpr2) ;

  theo=vth2/6^S*N;
  expr=zeros(1:6*S); for j=xpr2 do expr(j)=expr(j)+1; end;
  theo=theo(S:6*S); expr=expr(S:6*S);
  chi2=sum( (theo-expr)^2 ./ theo);
  printf('%d lancers de %d dés ; moy=%9.6f, var=%9.6f, chi2=%9.6f \n',N,S,mu2,va2, chi2)


  scf(S); clf();
  bar(vth2/6^S, width=0.4)
  bornes=S-0.5:1:S*6+0.5;
  histplot(bornes, xpr2)
  curax=gca();
  curax.children.children.background=12;
  curax.sub_ticks=[0,0];
  xtitle('somme de deux dés')
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;
  xpr3=feval(S*ones(1,N), dd);
  mu3=mean(xpr3) ; va3=variance(xpr3) ;

  theo=vth3/6^S*N;
  expr=zeros(1:6*S); for j=xpr3 do expr(j)=expr(j)+1; end;
  theo=theo(S:6*S); expr=expr(S:6*S);
  chi3=sum( (theo-expr)^2 ./ theo);
  printf('%d lancers de %d dés ; moy=%9.6f, var=%9.6f, chi2=%9.6f \n',N,S,mu3,va3, chi3)

  scf(S); clf();
  bar(vth3/6^S, width=0.4)
  bornes=S-0.5:1:S*6+0.5;
  histplot(bornes, xpr3)
  curax=gca();
  curax.children.children.background=12;
  curax.sub_ticks=[0,0];
  xtitle('somme de trois dés')
endfunction

function initpick
global N p1 p2 p3 p4 bor
  sig=x_mdialog('générateur à quatre choix', ['N','p1','p2','p3'],['500','0.5','0.3','0.1']);
  N=evstr(sig(1)) ; p1=max(0,evstr(sig(2))); p2=max(0,evstr(sig(3))); p3=max(0,evstr(sig(4))); p4=1-p1-p2-p3 ;
  if p4<0 then error ("bad values"); end
  printf("N=%d, p1=%4.2f, p2=%4.2f, p3=%4.2f, p4=%4.2f\n",N,p1,p2,p3,p4)
  bor=0.5:1:4.5
endfunction

function picka
global N p1 p2 p3 p4 bor lesx
  function y=gen
   z=rand(); if z<p1+p2 then if z<p1 then y=1; else y=2; end
   else if z<p1+p2+p3 then y=3; else y=4; end
   end
  endfunction
  lesx=feval(1:N,gen)
  scf(5); clf(); histplot(bor,lesx)
  xtitle(sprintf('%d générations par aiguillage',N))
endfunction

function pickb
global N p1 p2 p3 p4 bor lesx combien
  gap=[p1 p2 p3 p4]/max(p1,p2,p3,p4)
  combien=0
  function y=gen
  global combien
   while %t do z=rand(); y=ceil(z*4); combien=combien+1;
   if y-4*z<gap(y) then break ; end//if
   end//while
  endfunction
  lesx=feval(1:N,gen)
  scf(6); clf(); histplot(bor,lesx)
  xtitle(sprintf('%d générations par tir, en %d essais',N,combien))
endfunction

function pickc
global N p1 p2 p3 p4 bor lesx gap alt combien gen
  chances=[p1 p2 p3 p4]*4 ; qui=[1:4] ; gap=zeros(qui); alt= zeros(qui);
  while size(qui,'*')>0 do
    [M,k]=max(chances); [m,j]=min(chances);
    ici=qui(j); gap(ici)=m; alt(ici)=qui(k); chances(k)=M-(1-m);
    chances=[chances(1:j-1),chances(j+1:$)]; qui=[qui(1:j-1),qui(j+1:$)]
  end
  function y=gen
   z=rand(); y=ceil(z*4);
   if y-4*z>gap(y) then y=alt(y) ; end//if
  endfunction
  lesx=feval(1:N,gen)
  scf(7); clf(); histplot(bor,lesx)
  xtitle(sprintf('%d générations par Walker',N))
endfunction

function touslesech(S)
  //S est le nombre de dés
  tous=[1:6]'; for s=2:S do
   lesval=tous; lesuns=ones(lesval(:,1)); tous=[];
   for d=1:6 do tous=[tous;[d*lesuns,lesval]]; end;
  end;
  N=size(tous,'r'); lesm= sum(tous,'c')/S;
  espm=sum(lesm)/N; varm= (lesm'*lesm)/N-espm^2;
  printf("N=%4d,   E(m)=%f,   V(m)=%10.7f\n",N,espm,varm)
try
  lesS=sum(tous .* tous,'c')/(S-1)-lesm^2*S/(S-1);
  espS=sum(lesS)/N; varS= (lesS'*lesS)/N-espS^2;
  printf("N=%4d, E(S^2)=%f, V(S^2)=%10.7f\n",N,espS,varS)
catch; end
endfunction


function test_student
lesx=[80,115,101,97,103]
nn=size(lesx,'*');
mx=mean(lesx); sx=sqrt(variance(lesx));
conf=0.95; p=(1+conf)/2; tt=cdft('T',nn-1,p,1-p);
dx=sx*tt/sqrt(nn);
disp(lesx)
printf("\n    mx=%5.2f, sx=%5.2f, tt=%5.2f, dx=%5.2f, min=%5.2f, max=%5.2f\n\n",mx,sx,tt,dx,mx-dx,mx+dx)

conf2=1-0.001; p=(1+conf2)/2; 
nn2=nn
for j=1:8 do
  tt2=cdft('T',nn2-1,p,1-p);
  nn3=(sx*tt2/dx)^2;
  printf("    nn2=%3d, tt2=%f, nn3=%5.2f\n', nn2, tt2, nn3)
  if abs(nn2-nn3)<1 then break ;end
  nn2=round((nn3+nn2)/2)
end
  nn2=ceil(max(nn2,nn3)); printf("n=%d\n",nn2);
endfunction


function menupick
  while %t do
  n=x_choose(['initpick';'picka';'pickb';'pickc'],[ident+' pick: ';'choisir un item'])
  select n
  case 1 then initpick(); 
  case 2 then picka();
  case 3 then pickb();
  case 4 then pickc();
  else return
  end//select
  end//while
endfunction

function menusomme
  while %t
  n=x_choose(['initdes';'deuxdes';'troisdes'],[ident+' somme: ';'choisir un item'])
  select n
  case 1 then txt=['N']; sig=x_mdialog('enter N',txt,['500']); initdes(evstr(sig));
  case 2 then deuxdes() ;
  case 3 then troisdes() ;
  else return
  end//select
  end//while
endfunction

function menu
  while %t do
  n=x_choose(['somme de plusieurs dés';'la population de tous les échantillons'
              'test de Student';'générateur à quatre choix'],[ident+': ';'choisir un item'])
  select n
  case 1 then menusomme(); 
  case 2 then for S=1:5 do touslesech(S); end;
  case 3 then test_student(); 
  case 4 then menupick() ; 
  else return
  end//select
  end//while
endfunction

menu