function []=contenu()
// table des matieres
// Ch. Suquet, 23/11/2010
//
tit="Bibliotheque de fonctions pour les lois binomiales  et divers :"
doccoefbinA = " - coefbinA(n), coefficients binomiaux C_n^k, k=0:n";
docloibin = " - loibin(n,p), vecteur des P(S=k), k=0:n, S de loi Bin(n,p)";
docfdrbin = " - fdrbin(n,p,x), vecteur des P(S<=x(i)), x vecteur ligne";
docinvfdrbin = " - invfdrbin(n,p,y) inverse generalise de la fdr, version scalaire"
docvinvfdrbin = " - vinvfdrbin(n,p,y) inverse generalise de la fdr, version vectorielle"
docbatonsbin = " - batonsbin(n,p),  diagramme en batons de loibin(n,p)";
doclimlocbin = " - limlocbin(n,p), diagramme en batons et gaussienne";
docuniform = " - uniform(x,c), diagramme en batons de la loi uniforme sur les valeurs de x, couleur c";
doccnctbin = " - cnctbin(n,p,seuil), intervalle de concentration de Bin(n,p)";
docfluctbin = " - fluctbin(n,p,seuil), intervalle de fluctuation de Bin(n,p)";
//
printf('%s',tit)
printf('\n%s',doccoefbinA)
printf('\n%s',docloibin)
printf('\n%s',docfdrbin)
printf('\n%s',docinvfdrbin)
printf('\n%s',docvinvfdrbin)
printf('\n%s',docbatonsbin)
printf('\n%s',doclimlocbin)
printf('\n%s',docuniform)
printf('\n%s',doccnctbin)
printf('\n%s',docfluctbin)
endfunction


function [y]=coefbinA(n)
// "A" comme additive
// fournit le vecteur ligne des C_n^k, k=0:n
// methode additive par le triangle de Pascal
A=[1]; // initialisation n=0
for i=1:n, A=[A 0]+[0 A];end
y=A
//
// Limitations: on peut aller jusqu'a n=1029
// pour le voir, taper
// y=coefbinA(1029);max(y)
// puis
// y=coefbinA(1030);max(y)
endfunction

function [y]=loibin(n,p)
// n entier, p probabilite
// retourne le vecteur ligne des masses ponctuelles de la loi Bin(n,p)
// p(k+1)=C_n^k p^k q^(n-k)
//
C=coefbinA(n);
q=1-p;
y=C.*p.^(0:n).*q.^(n:-1:0);
endfunction

function [y]=fdrbin(n,p,x)
// fonction de répartition de Bin(n,p), x est un vecteur de reels
// retourne le vecteur ligne des P(S<= x(i)), où S suit Bin(n,p)
B=loibin(n,p);
CB=cumsum(B);
R1=zeros(x(x<0)); // si x<0, y=0
R2=CB( floor( x((0<=x)&(x<=n)) ) + 1); 
// si 0<=x(i)<=n, on prend k=partie entiere de x(i) et y(i)=P(S<=k)=CB(k+1)
R3=ones( x(x>n) ); // si x(i)>n, y(i)=1
y=[R1 R2 R3]
endfunction

function [x]=invfdrbin(n,p,y)
// Inverse generalise de la fdr de Bin(n,p), 0<y<1
// x=Inf{t réel tel que F(t)>= y}
B=loibin(n,p);
CB=cumsum(B);
x=min( find(CB >= y) ) - 1
//
// message d'erreur si CB(n+1)<=y<1
if (CB(n+1)<=y)&(y<1) then
printf('erreur : le 3eme argument de invfdrbin vaut\n...
y= %.20f\n...
qui est trop proche de 1 car numeriquement\n...
fdrbin(%d,%f,%d)=%.20f',y,n,p,n,CB(n+1))
end
// essayer: y=1-10^(-16);invfdrbin(6,0.3,y)
// 
endfunction

function [x]=vinvfdrbin(n,p,y)
// Inverse généralise de la fdr de Bin(n,p), 0<y<1
// version vectorielle 
B=loibin(n,p);
CB=cumsum(B);
for i=1:length(y),
x(i)=min( find(CB >= y(i)) ) - 1;
end
if (min(x)<0) then
printf('erreur numerique : il y au moins une composante trop proche ...
de 1 dans le troisième argument de vinvfdrbin')
end
endfunction


//  Fonctions graphiques
//

function []=batonsbin(n,p,c)
// affiche le diagramme en batons de loibin(n,p)
// avec la couleur c  exemple c=color("red")  
x=0:n;y=loibin(n,p);
//clf();
plot2d3(x',y',c)
endfunction

function []=batonsbinpart(n,p,c,a,b)
// affiche le diagramme en batons de loibin(n,p)
// avec la couleur c  exemple c=color("red")
// affichage restreint a l'intervalle [a,b]  
y=loibin(n,p);
x=a:b;
z=y(a:b);
//clf();
plot2d3(x',z',c)
endfunction

function []=limlocbin(n,p)
// affiche le diagramme en batons de loibin(n,p) en bleu
// et la densite gaussienne approximante en rouge
x=0:n;y=loibin(n,p);
xx=0:0.1:n;moy=n.*p;var=n.*p.*(1-p);
yy=(2.*%pi.*var)^(-0.5).*exp(-(xx-moy).^2./(2.*var));
clf();
plot2d3(x',y',2);plot2d(xx',yy',5)
endfunction

function []=uniform(x,c)
// trace le diagramme en batons de la loi uniforme
// discrete sur les valeurs du vecteur x
// exemple uniform(0:50,color("green"))
y = ones(x)./length(x);
plot2d3(x',y',c)
endfunction

// Intervalle de concentration et intervalle de fluctuation

function [I,s] = cnctbin(n,p,seuil)
  
  // retourne le plus court intervalle dont la probabilite pour la loi
  // Bin(n,p) est superieure ou egale a s
  B=loibin(n,p) ; 
  [BT,L] = gsort(B,'g'); // tri par masse decroissante avec memoire dans L de l'indice d'origine
  j = min(find(cumsum(BT)>=seuil));
  a = min(L(1:j))-1 ; b = max(L(1:j))-1 ; // on enlève 1 à tous les indices car P est indexé à partir de 1 au lieu de 0
  I= [a,b];
  s= sum(B((a+1):(b+1))) // masse de I
   endfunction
  // NB cet algorithme  marche avec toutes les lois discretes telles que la suite des masses soit croissante jusqu'au max puis decroissante, ce qui assure que les indices recuperes dans J(1:m)forment bien un intervalle (pas de trou).   
  
function [I,s] = fluctbin(n,p,seuil)
    
    //retourne le plus court intervalle symetrique centre sur np
    // dont la probabilite pour la loi Bin(n,p) est >= s.
    // programme a completer pour traiter le cas ou l'intervalle
    // de fluctuation deborde de [0,n], cas np petit.
    B=loibin(n,p) ;
    m = floor(n*p) ;
    if m == n*p then
      a = m; b  = m ; s = B(m+1) ;  // decalage b_k=B(k+1) 
    else 
      a = m; b = m+1 ; s = B(a+1)+B(b+1) ; 
    end
    
    while s < seuil,
      s = s + B(a) + B(b+2) ; // b_{a-1} = B(a) et b_{b+1} = B(b+2)
      a = a-1 ; b = b+1 ;
    end 
    I = [a,b]; 
endfunction
 
