function bi=biID_gamma(a,lambda,alpha,beta)
//    Retourne le quantile d'ordre beta
//    pour la loi gamma de parametres a et lambda
//    
if beta<alpha,
   bi = cdfgam("X",a,lambda,beta,1-beta);
else
   bi = 0;
end;

function bs=bsID_gamma(a,lambda,alpha,beta)
//    Retourne le quantile d'ordre 1-alpha+beta
//    pour la loi gamma de parametres a et lambda
//    
bs = cdfgam("X",a,lambda,1-alpha+beta,alpha-beta);

function io=IDoptim_gamma(a,lambda,alpha)
//    Retourne l'intervalle de dispersion optimal
//    de niveau 1-alpha pour la loi gamma
//    de parametres a et lambda.
//
                     // longueur d'un intervalle de dispersion
deff("l=long(b)",..
"l=bsID_gamma(a,lambda,alpha,b)-biID_gamma(a,lambda,alpha,b)");
                     // densite de la loi gamma
deff("y=dens(x)","y=((lambda^a)/gamma(a))*(x^(a-1)).*exp(-x*lambda)");
                     // derivee de la longueur
deff("d=lder(b)",..
"d=(1)./dens(bsID_gamma(a,lambda,alpha,b))..
  -(1)./dens(biID_gamma(a,lambda,alpha,b))");
                     // preparation de la fonction cout
deff("[f,g,ind]=cout(b,ind)","f=long(b),g=lder(b)");
b0 = alpha/2;        // initialisation
eps = 0.000000001;   
                     // optimisation
[l,b] = optim(cout,"b",eps,alpha-eps,b0);
                     // borne inferieure
bi = biID_gamma(a,lambda,alpha,b)
                     // borne superieure
bs = bsID_gamma(a,lambda,alpha,b)
                     // intervalle
io = [bi,bs];










