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

ident='DECIS.05';
if MSDOS .. 
  then cd("U:/promo/prenom.nom/decis") ;
  else cd("~/docs/Ensait/decis/scilab") ;
end//if
baserep=pwd();
myfig=baserep+'/../figures/';

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

global N dd vth1 vth2 
global N p1 p2 p3 p4 bor lesx gap alt gen
global obs hyp //test_chisquare

//*********************************************************
// les menus
//*********************************************************

function menu
  while %t do
  n=x_choose(['somme de plusieurs dés';'la population de tous les échantillons'
              'exercice_3_3_9';'générateur ŕ quatre choix';'exemple chisquare'],[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() ;
  case 5 then test_chisquare
  else return
  end//select
  end//while
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

//*********************************************************
// test_chisquare
//*********************************************************

function tmp=test_chisquare
global tmp tmq leso lesx obs hyp1 hyp2
  n=60; T=5;
  memo=grand('getsd'); grand('setsd',123456)
  tau=3 // la moyenne de la distrib exponentielle
  tmp=grand(1,ceil(n*1.5*1.2),'exp',tau);
  grand('setsd',memo);
  tmq=cumsum(tmp);
  [ind,leso,info]=dsearch(tmq,0:T:n*T)
  if info==0 then error 'not enough arrivals'; end//if
  lesx=0:4; borhist=[lesx-0.5, 4.5]; borobs=[lesx-0.5, 20.5];
  [ind,obs]=dsearch(leso,borobs);
  printf('\nE(a)= %f,  sd(a)= %f \n', mean(tmp), sqrt(variance(tmp)))
  moy=meanf(lesx, obs);
  printf(  'E(x)= %f, var(x)= %f \n', moy, variancef(lesx, obs) )
  deff('p=poiss1(n)',sprintf('p=%f^n/factorial(n)*%f', 1.2, exp(-1.2) ));
  hyp1=feval(lesx, poiss1)*n; hyp1($)=n-sum(hyp1(1:$-1)); 
  deff('p=poiss2(n)',sprintf('p=%f^n/factorial(n)*%f', moy, exp(-moy) ));
  hyp2=feval(lesx, poiss2)*n; hyp2($)=n-sum(hyp2(1:$-1)); 
  pears1=sum((obs-hyp1)^2 ./hyp1); pears2=sum((obs-hyp2)^2 ./hyp2);
  printf('\n obs = '); mprintf('%4d  '  , obs' )
  printf('\n hyp1= '); mprintf('%4.1f  ', hyp1'); printf('; Pearson= %7.4f, nu=%d, p=%5.3f'  , pears1, 4, cdfchi('PQ', pears1,4))
  printf('\n hyp2= '); mprintf('%4.1f  ', hyp2'); printf('; Pearson= %7.4f, nu=%d, p=%5.3f\n', pears2, 3, cdfchi('PQ', pears2,3))
  scf(1); clf(); 
  plot2d3(borhist, [hyp1,0], style=5); plot2d2(borhist, [hyp1,0], style=5)
  plot2d3(borhist, [hyp2,0], style=2); plot2d2(borhist, [hyp2,0], style=2)
  plot2d3(borhist, [ obs,0], style=1); plot2d2(borhist, [ obs,0], style=1)
  curax=gca();
  curax.sub_ticks=[0,0];
  curax.x_ticks=tlist(curax.x_ticks(1),lesx,string(lesx));
  curax.children(1).children.thickness=2;curax.children(2).children.thickness=2;
  curax.data_bounds(1)=-1;

endfunction

//*********************************************************
// somme de plusieurs dés
//*********************************************************

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','seed'],['500','0.15','0.4','0.25','59861689']);
  if sig==[] then return; end//if
  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
  seed=evstr(sig(5)); rand('seed', seed);
  // caveat: les seed de rand et grand sont indépendants...
  printf("N=%d, p1=%4.2f, p2=%4.2f, p3=%4.2f, p4=%4.2f, seed=%d\n",N,p1,p2,p3,p4,seed)
  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 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 exercice_3_3_9
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



menu