# Read in NP times
times = sort(scan(file="NP_time.txt" , what = 0.0, n = -1, skip = 0))
# Fit a gaussian distribution to the given data
ftn = function(par) {
  # Return the least-squares difference between the measured and parameterised gamma cdfs.
  return(sum((ecdf(times)(times) - pgamma(times,shape=par[1],scale=par[2]))^2))
}
par = optim(c(1,1), ftn)
# plot the result as a histogram / pdf
pdf('pdf.pdf')
hist(times, freq=FALSE, breaks=seq(0,max(times),max(times)/10), 
     xlab="Patient time in system", ylab="pdf", main="")
lines(times, dgamma(times, shape=par$par[1], scale=par$par[2]))
dev.off()
# plot the result as a cdf
pdf("cdf.pdf")
plot(times, ecdf(times)(times), xlab = "Patient time in system", ylab="cdf", main="")
lines(times, pgamma(times, shape=par$par[1], scale=par$par[2]))
dev.off()
# print fitted parameters
cat(sprintf('Fitted scale: %f, shape: %f, (rate: %f)\n', par$par[1], par$par[2], 1/par$par[2]))