// Estimation de la longueur de la file pour la chantier routier

// Paramètres
N=100;     // nombre d'échantillons   (npmbre de simulations)
n=2000;    // taille des échantillons (durée de chaque simulation)
//
a=4/60;    // a= 4 sec. = 4/60 mn.
lambda1= 4;
lambda2= 3;
// d1 + d2 + 2L/v = 5mn ; L=.6km ; v= 0.5 km/mn ; d2 = 3/4 d1 ;
// d'où d1 = 52/35 mn ; d2 = 39/35 mn.
L=.6;
v=.5;
d1=52/35;
d2=39/35;
K=int(d1/a);
mu1= lambda1 * d1;
mu2= lambda1 * (d2+2*L/v);
// Simulation de la suite X(i) (1 <= i <= n)
Y=grand(N,n,'poi',mu1);
Z=grand(N,n,'poi',mu2);
X=zeros(N,n);
for i=1:n-1,
  X(1:N,i+1)= max(zeros(N,1), X(1:N,i)+Y(1:N,i)-K) + Z(1:N,i);
end
// Moyenne et écart-type empirique de X(n)
if N>1,
  Moy=sum(X(1:N,n)) / N ;
  Var=(sum((X(1:N,n)).^2)-Moy^2) / (N-1);
end;
printf('Moyenne  empirique  des  X(n) : %.2f\n', Moy);
printf('Ecart-type empirique des X(n) : %.2f\n', sqrt(Var));
// La longueur moyenne d'un véhicule étant de 8 mètres, on s'intéresse
// à la probabilité que X(n) > 31 (longueur 31 * 8 = 248 m).
// On la calcule sur *toutes* les trajectoires pour plus de précision
// (une seule --très longue-- suffirait d'après le théorème ergodique,
//  comme on peut le vérifier en faisant N=1 et n=200000).
Freq1=sum(bool2s(X > 31 * ones(X))) / N /n;
Freq2=sum(bool2s(X > 25 * ones(X))) / N /n;
printf('Fréquence de X > 31 (250 m) : %.3f\n', Freq1);
printf('Fréquence de X > 25 (200 m) : %.3f\n', Freq2);
// Tracé de l'histogramme des fréquences
xbasc();xselect();
histplot([0:50],X,rect=[0,0,50,0.10],nax=[10,6,5,3]);
xsegs([31;31],[0;0.05],style=5);
// La forme de l'histogramme n'incite pas à faire l'approximation normale,
// ni poissonienne car pour X Poisson E(X)=Var(X)...

// On notera que la précision sur les fréquences dépend du produit N*n.
// Pour N=1 et n=200000, on aurait la même précision, mais le temps de calcul
// serait beaucoup plus long qu'avec avec N=100 et n=2000.