%
%  macm 202 -- 03 mar 03 -- djm
%
%  w9qstats.m:  single server queue statistics
%
%     - single-queue statistics

%  ordered (nonzero) wait time data

ord_t = wsort(begin:end)';
Nt = size(ord_t,2);
mu = 1/mean(ord_t);
ord_z  = 1-exp(-mu*ord_t);

disp([num2str(Nt) ' = # who queued'])
disp([num2str(begin-1) ' = # no queue'])

%  EDF points

Fplus = (1:Nt)/Nt;  Fminus = (0:Nt-1)/Nt;
edfx = [ord_t ; ord_t];  edfy = [Fminus ; Fplus];

figure(2);  clf
plot(ord_t,Fplus,'.')
hold on;
plot(edfx(:),edfy(:),'k')
plot([0 ord_t],[0 ord_z],'r')
title('\bf EDF of ordered (nonzero) waiting times')
xlabel('\bf ordered wait times, x_i (i=1...n)')
ylabel('\bf empirical (sample) distribution function, F_n(x_i) & F(x)')
axis([0 ord_t(end) 0 1]);

edfx = [ord_z ; ord_z];

figure(3);  clf
plot(ord_z,Fplus,'.')
hold on;
plot(edfx(:),edfy(:),'k')
plot([0 ; 1],[0 ; 1],'r')
title('\bf uniform EDF of (nonzero) waiting times')
xlabel('\bf transformed ordered wait times, z_i = F(x_i) = 1-e^{-x_i/\theta} (\theta=mean)')
ylabel('\bf EDF, F(z_i) & line of slope 1')

%  edf statistics

disp(' ')
Dplus = max(abs(Fplus-ord_z));
Dminus  = max(abs(ord_z-Fminus));
Dmax   = max([Dplus Dminus]);  disp(['D = ' num2str(Dmax)])
Dmax   = max([Dplus Dminus]);  disp(['stat = ' num2str(sqrt(Nt)*Dmax)])

Wsquared = sum((ord_z - (Fplus+Fminus)/2).^2) + 1/(12*Nt) ;
disp(['W^2 = ' num2str(Wsquared)])
disp(' ')

%  expected values of ordered statistics

ord_m = cumsum(1./(Nt:-1:1));
figure(4);  clf
plot(ord_m,ord_t,'.')
hold on;
plot([0 ; mu*ord_t(end)],[0 ; ord_t(end)],'r')
title('\bf expected value correlation line')
xlabel('\bf expected value of ordered times, m_i')
ylabel('\bf transformed & ordered times, t_i')

t_mean = mean(ord_t);  disp(['mean t_i = ' num2str(t_mean)])
m_mean = mean(ord_m);  disp(['mean m_i = ' num2str(m_mean)])

rsquared = (sum((ord_m-m_mean).*(ord_t-t_mean)))^2 / ...
    (sum((ord_m-m_mean).^2) * sum((ord_t-t_mean).^2));
disp(['r^2 = ' num2str(rsquared)])
disp(' ')