mu=1/5; % i.e. 5 minutes to shop
lambda=1/1; % arrival rate: 1 per min.
dt=0.01;
nsteps=1000000;
Q=zeros(1,nsteps);
Q(1)=0;
for i=2:nsteps
    Q(i)=Q(i-1);
    if rand()< lambda*dt   % arrival rate
        Q(i)=Q(i)+1;
    end;
    if (Q(i)>0)  
        if rand()< mu*dt*Q(i) % 'depature' rate
            Q(i)=Q(i)-1;  
        end;
    end;
end;

% Compute distributions


p=[];
for i=0:max(Q)
    p(end+1)=sum(Q==i)/nsteps;
end;


clf;
plot(0:max(Q), p, 'o');
title('Probabilities of having exactly k people in queue');
xlabel('k'); ylabel('\pi_k','rotation',0);


hold on;
rho=lambda/mu;
k=[0:max(Q)];
plot(k, rho.^k./factorial(k).*exp(-rho), '--');
legend('simulation', 'exact');


% Output pi_0: simulation vs. exact:
fprintf('pi_0 (simulation) = %g, pi_0 (exact)=%g \n', p(1), exp(-rho));
