#!/usr/bin/R --no-save < 
times = scan(file="emergency-1-2-2.txt" , what = 0.0, n = -1, skip = 1)

plotPatients = function(points) {
  # plot the given array using lines on a graph.
  L = length(points)
  plot(c(1, L), 
       c(0, max(points)), type="n", 
       xlab="patient number", ylab="Patient time in system (hours)")
  lines(1:L, points)
}

pdf('short.pdf')
plotPatients(times[1:100])
dev.off()
pdf('long.pdf')
plotPatients(times)
dev.off()

# Now calculate the mean and variance using a blocking procedure
burn_in = 1000
block_size = 10000
# for each set of simulated data...
for (fname in list.files(path='.', pattern='emergency-.*txt')) {
  times = scan(file=fname, what = 0.0, n = -1, skip = 1, quiet=TRUE)
  times = times[burn_in:length(times)] # throw away initial burn_in data
  # calculate number of blocks
  blocks.num = floor(length(times) / block_size)
  means = rep(0, blocks.num)
  # calculate mean for each block
  for (i in 1:(blocks.num)) {
    range = (block_size * (i-1)):(block_size*i)
    means[i] = mean(times[range])
  }
  # calculate mean of means and error margin
  mean.best_guess = mean(means)
  mean.std = sd(means)      
  cat(sprintf('%s: Mean: %6.4f; 95%% confidence interval [%6.4f, %6.4f]\n', 
              fname, mean.best_guess, mean.best_guess - mean.std*2,
              mean.best_guess + mean.std * 2))
  cat(sprintf('%s\t%.3f\t[%.3f, %.3f]\n', 
              fname, mean.best_guess, mean.best_guess - mean.std*2,
              mean.best_guess + mean.std * 2));
}