// Simulation d'une file M/GI/1.
// Aux instants de départ d'un client :
// Y_{n+1} = Y_n + 1_{Y_n>0} + K_{n+1}  où K_n est le nombre de clients
// arrivés pendant le service du n-ième client.

// Durée de la simulation (en nombre d'arrivées observées).
T=10000;

// La durée de service est une v.a. D. :
// on simule d'abord un échantillon de taille T de v.a. de même loi que D.

// 1) Dans le cas le plus simple, D est constante : 
//D=2;
//D=D * ones(T,1);
//ED=D(1);        // espérance

// 2) D uniforme sur [a,b] :
//a=1;b=3;
//D=a * ones(T,1)+ (b-a) * rand(T,1,'uniform');
//ED= (a+b)/2; // espérance

// 3) D exponentielle de paramètre a
a=1/2;
D=-log(rand(T,1,'uniform')) ./ a; 
ED= 1/a;    // espérance

// Paramètre du processus de Poisson :
// lambda = nombre moyen d'arrivées par unité de temps.
lambda=0.45;   // cas stable   (rho=0.9  si ED=2)
//lambda=0.5;  // cas limite   (rho=1    si ED=2)
//lambda=0.52; // cas instable (rho=1.04 si ED=2)

// Initialisation de Y 
Y= zeros(T+1,1);
// Calcul de rho et affichage (file stable ou instable)
rho= lambda * ED;
if rho < 1 then
   printf('rho = %f < 1 : file stable.\n', rho)
else
   printf('rho = %f >= 1 : file instable.\n', rho)
 end

 // Boucle de calcul des K(n) et Y(n)
for n=1:T,
  // K(n)= nombre de clients arrivant dans un intervalle de largeur D(n),
  // K(n) suit une loi de Poisson de paramètre lambda * D(n)
  // N.B. : grand ne permet pas de prendre un paramètre vectoriel (lambda * D),
  // ce qui impose de calculer les K(n) à l'intérieur de la boucle...
  K(n)=grand(1,1,'poi', lambda * D(n));
  if Y(n)==0 then
     Y(n+1)=K(n);
  else
     Y(n+1)=Y(n) -1 + K(n);
  end
end
printf('Nombre moyen empirique de clients présents = %f', sum(Y) ./ T)

// Affichage trajectoire (si rho>=1) ou histogramme (si rho<1).
xbasc();xselect();
if rho < 1 then
  // Si stable : histogramme (approxime le loi limite).
  Y=sort(Y); // tri de Y par ordre décroissant
  // on élimine les + grandes valeurs : celles dont la fréquence cumulée
  // est inférieure à< 0.2%.
  h =Y(max(1,int(T / 500)));
  histplot([-1:h], Y) // la tranche [-1,0] contient la fréquence de 'Y=0',
                      // voir help('histplot') : classes = ]n-1,n]
  // Tracé de la loi invariante pour les services exponentiels
  // (loi géométrique de paramètre rho) :
  pi=(1-rho) * (rho * ones([0:h])).^[0:h];
  plot2d([0:h]',pi', style=5);
else
  // Si instable : affichage de la trajectoire
  plot2d([0:T]', Y);
end

// Si rho<1, l'histogramme approche la probabilité invariante 
// (théorème ergodique des chaînes de Markov).
// Vérifier dans le cas rho<1,
// -- que la fréquence de 'Y=0' est proche de pi_0=(1-rho),
// -- que, dans le cas où les services  sont exponentiels la moyenne
//    empirique du nombre de clients présents est proche de rho / (1-rho).
// Dans le cas de services exponentiels on connait la loi limite
// (géométrique de paramètre rho), on peut la faire tracer cf. ci-dessus.