function kap8_3b()
% statistiscs with nQ(t)
tEnd = 10000;
seed = 4;

NQ = runSingleserver4B(seed, tEnd);
tIn = NQ.time;
nQIn = NQ.signals.values;

load_system('singleserver4B');

% get data at fixed times
t = (1:tEnd)';
nQ = interp1(tIn, nQIn, t, 'previous');
Y = mean(nQ)
S = std(nQ)

% compute total rho
rho = getRho(nQ)

% use repetitions
nRep = 100;
tEnd = 100;
t = (1:tEnd)';
Y = zeros(nRep, 1);
S = zeros(nRep, 1);
for I=1:nRep
  NQ = runSingleserver4B(2*I, tEnd);
  fprintf('  Runde %3d\n', I)
  tIn = NQ.time;
  nQIn = NQ.signals.values;
  nQ = interp1(tIn, nQIn, t, 'previous');
  Y(I) = mean(nQ);
  S(I) = std(nQ);
end

% analysis of sample means Y
Ytot = mean(Y)
Stot = std(Y);
S2tot = Stot^2
alpha = 5/100;
z = tinv(1-alpha/2, nRep-1);
dS = z*Stot/sqrt(nRep);
fprintf('Konfidenzintervall: [%6.4f, %6.4f]\n', Ytot-dS, Ytot+dS)

set(gcf(), 'Position',[1 1 700 400], 'Color',[1 1 1]);
subplot(1,2,1)
plot(1:nRep, Y, [0,nRep], [Ytot,Ytot])
title('Mittelwerte der Läufe')
xlim([0,100])

subplot(1,2,2)
plot(1:nRep,S, [0,nRep], [Stot,Stot])
title('Standardabweichungen der Läufe')
xlim([0,100])

annotation('textbox', [0.471 0.267 0.03 0.052],...
    'String','$\overline{Y}$',...
    'FitBoxToText','off', 'Interpreter', 'latex');
annotation('textbox', [0.914 0.267 0.03 0.052], ...
    'String','$S^2$',...
    'FitBoxToText','off', 'Interpreter', 'latex');

F = getframe(gcf(), [65, 20, 600, 370]);
imwrite(F.cdata, 'bild241.png');

%----------------------------------------------------------------
function rho = getRho(Xi)
% estimates correlation coefficient (Xi,X_i+1) of Xi
n = length(Xi);
Xbar = sum(Xi)/n;
S2 = (Xi - Xbar)'*(Xi - Xbar)/(n-1);
C = (Xi(1:end-1) - Xbar)'*(Xi(2:end) - Xbar)/(n-1);
rho = C/S2;
