mydata <- c(90, 53, 94, 303, 138, 116, 179, 82, 78, 64) #data in seconds
mydata <- mydata/60 #converting data from seconds to minutes
plotdist(mydata, histo = TRUE, demp = TRUE) #plotting using fitdistrplus
mean(mydata) #scale parameter- mean number of minutes until next customer
## [1] 1.995
1/(mean(mydata)) #MLE for rate parameter, ie. customers per minute
## [1] 0.5012531
fit_g <- fitdist(mydata, "gamma")
gofstat(fit_g) #we fail to reject the hypothesis that the data come from a gamma distribution
## Goodness-of-fit statistics
## 1-mle-gamma
## Kolmogorov-Smirnov statistic 0.21293574
## Cramer-von Mises statistic 0.06888454
## Anderson-Darling statistic 0.43708412
##
## Goodness-of-fit criteria
## 1-mle-gamma
## Akaike's Information Criterion 30.73512
## Bayesian Information Criterion 31.34029
summary(fit_g) #note that the package reports the rate, not the scale
## Fitting of the distribution ' gamma ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape 3.911542 1.6800310
## rate 1.960625 0.8985585
## Loglikelihood: -13.36756 AIC: 30.73512 BIC: 31.34029
## Correlation matrix:
## shape rate
## shape 1.000000 0.937168
## rate 0.937168 1.000000
Now that we know the MLE estimator for the rate and the number of customers who need to be seen (7), let’s use stochastic processes to generate a 50% CI.
set.seed(1234)
rdata <- rgamma(10000, shape=7, rate = 1/fit_g$estimate[[2]])
hist(rdata)
quantile(rdata, prob = c(0.25, 0.75))
## 25% 75%
## 10.04325 16.68872
quantile(rdata, prob = c(0.025, 0.975))
## 2.5% 97.5%
## 5.516821 25.492242