getf("IDgamma.sci");
//
// Representation d'un intervalle de dispersion optimal pour la loi 
// Gamma G(5,0.5) (loi de khi-deux a 10 ddl).
//
a= 5; l=0.5;
deff("y=dens(x)","y=((l^a)/gamma(a))*(x^(a-1)).*exp(-x*l)");
xset("window",0); xbasc();
plotframe([0,0,30,0.1],[0,10,0,10],[%f,%f],..
["loi gamma G(5,0.5)","",""]);
inter=IDoptim_gamma(a,l,0.05);
bi = inter(1); bs = inter(2);
xset("use color",0);
xset("pattern",13);
xpol =linspace(0,bi,20);
ypol = dens(xpol);
xpol = [xpol,[bi,0,0]];
ypol = [ypol,[0,0,dens(0)]];
xfpoly(xpol,ypol);
xpol =linspace(bs,30,20);
ypol = dens(xpol);
xpol = [xpol,[30,bs,bs]];
ypol = [ypol,[dens(30),0,dens(bs)]];
xfpoly(xpol,ypol);
xset("pattern",1);
x = [0:0.1:30];
y= dens(x);
plot2d(x,y,1,"000");
plot2d([bi,bi],[0,dens(bi)],1,"000");
plot2d([bs,bs],[0,dens(bs)],1,"000");
//
// Comparaison des intervalles de dispersion 
// de niveau 0.9 pour la loi gamma G(10,1). 
// 
beta = linspace(0.001,0.099,100);
a = 10*ones(beta); lambda = ones(beta);
alpha=0.1*ones(beta);
bi = cdfgam("X",a,lambda,beta,1-beta);
bs = cdfgam("X",a,lambda,1-alpha+beta,alpha-beta);
xset("window",1); xbasc();
plotframe([0,0,0.1,30],[0,10,0,10],[%f,%f],..
["Intervalles de dispersion loi gamma","",""]);
plot2d([beta',beta'],[bi',bs'],[1,1],"000");
//
// Representation d'intervalles de dispersion optimaux 
// pour la loi Gamma G(10,lambda) en fonction de lambda
//
a = 10; alpha = 0.1;
lambda=linspace(0.1,10,500);
inter=[];
for l=lambda;
  inter=[inter;IDoptim_gamma(a,l,alpha)];
end;
lambda=lambda'*[1,1];
xset("window",2); xbasc();
plotframe([0,0,10,50],[0,10,0,10],[%f,%f],..
["Intervalles de dispersion loi gamma","",""]);
plot2d2("gnn",lambda,inter,[1,1],"000");



