Abstract
In many everyday situations in which a queue is formed, queueing models may play a key role. By using such models, which are idealizations of reality, accurate performance measures can be determined, such as traffic intensity (), which is defined as the ratio between the arrival rate and the service rate. An intermediate step in the process includes the statistical estimation of the parameters of the proper model. In this study, we are interested in investigating the finite-sample behavior of some well-known methods for the estimation of for single-server finite Markovian queues or, in Kendall notation, queues, namely, the maximum likelihood estimator, Bayesian methods, and bootstrap corrections. We performed extensive simulations to verify the quality of the estimators for samples up to 200. The computational results show that accurate estimates in terms of the lowest mean squared errors can be obtained for a broad range of values in the parametric space by using the Jeffreys’ prior. A numerical example is analyzed in detail, the limitations of the results are discussed, and notable topics to be further developed in this research area are presented.
1. Introduction
Queueing models are idealizations of many real systems. However, they enable the accurate determination of performance measures as long as a previous step has been fulfilled; that step is the statistical estimation of its parameters [1, 2]. It is impossible to discuss inference in queues without mentioning the pioneering work of Clarke in the 1950s, who describes maximum likelihood estimators for the arrival rates and service times in simple queues, and the work of Schruben and Kulkarni in the 1980s, who consider the problem of bias in queue estimation. It is also important to mention the Bayesian papers on the subject, including that of Muddapur in the 1970s, who published one of the first results that extended Clark’s methodology, the series of papers from Armero and Bayarri [3, 4], Armero and Conesa [5–8], Choudhury and Borthakur [9], and, more recently, Chowdhury and Mukherjee [10, 11], Cruz et al. [12], and Quinino and Cruz [13]. These are only a few examples of the papers in this important research area.
The purpose of this paper is to evaluate the behavior of a traffic intensity estimator (), which is defined as the ratio between the arrival rate () and the service rate (), specifically for single-server finite Markovian queues, which model various real systems [14]. The type of queue researched in this paper is typically seen in many service-oriented settings, where there is a finite queue in front of a server. Think of a gas station, where cars can queue up for a limited amount of space, and the traffic intensity should be not too high, as in this case cars might not join the queue. An alternative example is observed within a supermarket context where customers line up for the check out. Having a good understanding of the traffic intensity is crucial in these situations. Good estimations are needed to properly design the system (e.g., to install the gas pump buffer sizes needed) or to properly manage the system (e.g., to set waiting spaces for the store clerks at check out).
In summary, previous results obtained for infinite Markovian queues are extended here for finite Markovian queues. To reach this goal this paper combines some techniques and approaches from the work of Almeida and Cruz [15] (i.e., Bayesian inference and Monte Carlo simulation for evaluation of estimators under finite samples) with other classical tools (e.g., Gibbs sampling and bootstrapping).
The remainder of the paper is organized as follows. Section 2 details the queue equations and estimators for . The computational results are presented and discussed in Section 3, followed by Section 4, which concludes the text with some final remarks and topics for future research in the area.
2. Material and Methods
When you have Poisson arrivals, exponential service times, a single server, and limited waiting space, you have an queue; in Kendall notation, represents the number of customers simultaneously allowed in the queueing system. The probability of a number of users of the system, for , is given by [14]:where is the traffic intensity. Estimating traffic intensity is important as it is a key design parameter in production network design, routing of products, and so on.
2.1. Maximum Likelihood Estimator
Maximum likelihood estimation for the truncated geometric model is known since Thomasson and Kapadia [16]. Consider the stationary probability distribution given by (1). Next, consider a random sample of size , , where is the number of customers an outside observer finds in the system. In this case, a maximum of customers are allowed in the queueing system at once. Therefore, the likelihood function iswhere . Note that the likelihood is a function of traffic intensity and sample , although only its size, , and its sum, , which is a sufficient statistic for , are necessary.
Needless to say that, for the implementation of maximum likelihood estimator (MLE), any bounded optimization algorithm could be used. However for the sake of simplicity, the implementation used in this study was encoded in R [17] and can be seen in Listing 1. For convenience, the logarithm of the likelihood function was considered because it allows products to be turned into sums. Maximizing the log-likelihood is done numerically through an R internal function. However, tests (not shown) were conducted with the original likelihood function; they indicated that the results did not change significantly in terms of accuracy or computational effort.
| ||||||||||||||||
A nice feature about the MLEs is their invariance to transformations [18] in such a way that if is the MLE of and is a function of , then is the MLE of . Thus, the expected number of customers in the system and the average queue length have, respectively, the following MLEs [14]:in which is the MLE for .
2.2. Bayesian Inference
One of the alternatives for making inferences is the Bayesian method. One of its main differences from the classical method is that Bayesian inference allows the incorporation of some a priori information into the model of the unknown parameters. Unlike the classical method, the Bayesian method considers these parameters random variables, i.e., associates them with probability distributions. Therefore, the knowledge that the manager has regarding a given unknown parameter can be considered. Two different prior distributions are described as follows. However, other alternatives are possible, such as those proposed in the study by Armero and Bayarri [3]. Lingappaiah [19] also considered a Bayesian approach to this distribution.
2.2.1. Beta Prior
Therefore, the inference process for estimating starts with (2) and assumes an a priori beta distribution, that is, , which has been successfully used in inference in other Markovian queues [12, 13] and results in the following a posteriori distribution:
Because the a posteriori distribution of (5) is not a known distribution, it is necessary to use an approximation method to generate samples from the distribution. The implementation of the Bayesian estimator for beta prior is shown in Listing 2. Note that this is a bounded one-dimensional problem, and numerical integrations could be used as well to find the a posteriori distribution, which should be followed by another numerical integration to compute the estimator. An open research question is whether or not one method could be consistently superior to the other depending on the sample sizes in Gibbs sampling and the precision of the two numerical integrations.
| ||||||||||||||||||||
A sample is extracted from the a posteriori distribution by using the function arms (for size 5,000), which is available in R’s HI package [20]. The a posteriori distribution was represented by the logarithm of the probability density function (without the normalization constant). Because a quadratic loss function is considered, the Bayes point estimator is simply the average of the sample. Examples of a priori beta distributions are shown in Figure 1.
2.2.2. Jeffreys Prior
Someone might argue that it would be more natural to use a noninformative Jeffreys’ prior distribution, which is defined in terms of the Fisher information given by
Thus, the following prior distribution for can be obtained:
Thus, from (1), a can be found as follows. The logarithm of , from (1), is given by
First and second derivatives are given, respectively, by
It follows that
The expectation is given by [14]
Then
Thus, combining the likelihood, (2), and Jeffreys prior, given by (13), it is possible to find the following posterior probability distribution for :in which .
Note that, unlike the case shown earlier, (5), the posterior distribution given by (14) is fixed. Indeed, it is not possible to vary Jeffreys’ prior, which assumes a specific form. The implementation of the Bayesian estimator for Jeffreys’ prior is shown in Listing 3. Because the logarithm of the posterior distribution is used, variable Irho, as defined in Listing 3, must be such that 0 < Irho < ∞.
| ||||||||||||||||||||||||||||
2.3. Bootstrap Correction
Among the bias correction methods that are commonly used for estimators, the bootstrap method is widely used [21]. In its nonparametric version, the method consists of making several (usually ) resamplings (with replacement). The parameter of interest is reestimated by some (usually biased) method, . Then, the average is calculated, , and the bias is estimated usingwhere is the estimate obtained from the original sample. Then, the following corrected version of the estimator is obtained:
The procedure is illustrated in Figure 2. The R code [17] for the bootstrap corrected estimator is shown in Listing 4. The parameter is the maximum number of users simultaneously (in service and waiting) in the queue. Note that 1,000 bootstrap replicates are used and that the correction occurs as part of the MLE (see Listing 1).
| ||||||||||||||||||
Besides bias correction, the bootstrap method has been used by many researchers in the past with good results in confidence interval building and hypothesis testing [22]. As an example, an empirical bootstrap confidence interval is used in this work as a simple way of interval estimation for . If the distribution of was known then the critical values and could be found, where is its th percentile, and then which gives an confidence interval of
The bootstrap makes it possible to estimate the distribution of by the distribution of , where is the estimate obtained from an empirical bootstrap sample as explained earlier.
2.4. Simulation M/M/1/K of Queues
The number of users present in an queue follows the distribution given by (1). To efficiently generate random variables from a discrete distribution several methods are widely used in the literature, including function sample, from R [17]. The method used here is the discrete analog of the inverse transformation method, in which it is necessary to generate numbers , i.e., from a uniform distribution between 0 and 1, and to know the probabilities of interest, . Therefore, to simulate a discrete random variable with the probability function it is necessary to compute
because and follows the required probability distribution.
For queues and from (1), the following must hold:
Setting , it follows that
Therefore,where is the ceiling function; that is, its value is the least integer that is not inferior to .
Listing 5 presents the R implementation [17] used in this study with a sample call.
| ||||||||||||||||||||||
3. Computational Results
To analyze the performance of the estimators, 10,000 Monte Carlo replications were made using the R code [17] shown in Listing 6. Note that fESt(K, samp,…) can be any of the implemented estimation functions for (MLE, Bayesian, and bootstrap corrected MLE) and that the state of the random seed GlobalEnv$.Random.seed was stored immediately before the function fESt(K, samp,…), which can be stochastic, and was reloaded immediately after its call to ensure that the samples generated for the estimates were the same for all estimators.
| ||||||||||||||||||||||
3.1. Simulation Results
Samples were generated from (24) and Listing 5, for sizes ; and traffic intensities . In these scenarios, averages of 10,000 Monte Carlo replications were computed (via code from Listing 6), and point estimates of were computed by(i)the MLE via numerical maximization of the likelihood function, (2), and code from Listing 1;(ii)the Bayesian method via an a priori distribution and an average of 5,000 samples from the a posteriori distribution, (5), obtained by the function arms [23] from R’s HI package [20] via code from Listing 2;(iii)the Bayesian method via Jeffreys prior distribution and an average of 1,000 samples from the a posteriori distribution, (14), obtained by the function arms [23] from R’s HI package [20] via code from Listing 3;(iv)the bootstrap corrected MLE, (16), and code from Listing 4. Additionally, mean squared errors (MSEs) were calculated.
The results are shown in Tables 1, 2, and 3 and summarized in Figures 3, 4, and 5. In Figures 3, 4, and 5, the average estimation error and the mean squared error (MSE), defined as , are shown as functions of both the traffic intensity, (averaged over all sample sizes), and the sample size, (averaged over all traffic intensities).
(a) Average error for
(b) Average MSE for
(c) Average error for
(d) Average MSE for
(a) Average error for
(b) Average MSE for
(c) Average error for
(d) Average MSE for
(a) Average error for
(b) Average MSE for
(c) Average error for
(d) Average MSE for
For queues with capacity (Figure 3), an approximately constant average estimation error is observed for the MLE and the bootstrap corrected MLE, except when (Figure 3(a)). The Bayesian estimators (beta prior and Jeffreys’ prior) did not show equivalent performance; their estimates tended to overestimate the true value (positive error) when and to underestimate otherwise. Regarding the sample size, , all of the estimators showed a monotonic decrease in error (Figure 3(c)). From the MSE side, both Bayesian estimators presented themselves generally as the best alternative because they presented the lowest values most of the time. Although the bootstrap corrected MLE presented the lowest bias, the method achieves this performance at the cost of high variability, as reflected by its highest MSEs.
When the queue capacity was increased slightly to (Figure 4), a very similar behavior was noted. However, the difficulty of estimation for traffic intensities seemed to decrease, and the highest MSEs occurred when . The Bayesian estimators maintained the performance presented earlier, for ; that is, the estimates have positive bias for and negative bias otherwise. The errors of all of the estimators converged to zero as the sample size grew. The bootstrap corrected MLE presented an average estimation error near zero for samples . However, from the point of view of the MSE, the smaller values were again obtained using the Bayesian methods (beta and Jeffreys’ prior).
Finally, for queues with (Figure 5), the observed behavior could be considered, for practical purposes, as being equal to that of an infinite Markovian queue in terms of the average estimation error and the MSE. This behavior was observed for infinite Markovian queues ( queues) [15]. This finding is merely evidence of the correctness of our implementations and the quality of the computational results presented.
Additionally, computational experiments were performed for the estimators for , (3), and for , (4), and for their bootstrap corrected versions, with 1,000 resamplings, and . Table 4 and Figure 6 show the results obtained for , for sample sizes and averages of 10,000 Monte Carlo replications. Similarly for , the results are presented in Table 5 and Figure 7. In summary, at an extra cost of the bootstrap method and without inflation of the MSEs, the researcher may achieve with samples of size estimates for and with the same average error of the MLE for samples of size , a reduction that is relevant in practical terms because it may lead to reduction in time and cost to obtain the estimates. Note that the bootstrap method always provides smaller errors and MSEs than the MLE method for all estimates of and , even when . Also note the jump up and down in the errors and MSEs when transitions from to .
(a) Average error for
(b) Average MSE for
(c) Average error for
(d) Average MSE for
(a) Average error for
(b) Average MSE for
(c) Average error for
(d) Average MSE for
Finally, to illustrate the ease of use of the bootstrap in the interval estimation of the traffic intensity , computational experiments were performed. The length and coverage of empirical bootstrap intervals computed from (18) and from a normal distribution approximation (i.e., , where is the th percentile of the standard normal distribution and the standard deviation was estimated also by bootstrapping) were evaluated for , for sample sizes , averages of 10,000 Monte Carlo replications, and . The satisfactory performance of the bootstrap was demonstrated, as presented in Table 6, with the coverages approaching the nominal confidence of 95% (that is, ) as the sample sizes increase.
3.2. Numerical Example
To better illustrate an application of the method, a numerical application based on the data considered in Table 7, collected in a large supermarket network in a region of interest [12], is provided. The goal is to evaluate traffic intensity, which, for managerial reasons, should not exceed 87%; if it does, users may leave. The data comprise 200 random observations of the number of customers in the system at random times sufficiently spaced and previously defined by the person responsible for collecting the data to avoid correlation. The observed values (V) and frequency (F) are presented in Table 7. For instance, of the 200 observations, at 8 times, no customers were found in the system; at 21 times, only 1 customer was found, and so on.
The estimates are shown in Table 8; they were calculated using the MLE (Listing 1) and the Bayesian estimators (Listings 2 and 3) based on an a priori and 5,000 samples and Jeffreys’ distribution and 1,000 samples from the a posteriori distribution; the bootstrap corrected MLE (Listing 4) was based on 1,000 resamplings. The complete script is shown in Listing 7.
| ||||||||||||||||||||||||||||||||||||||
According to the results presented in the previous section, the Bayesian estimates should be the most reliable. It is notable that the system utilization seemed to be below the target (). Note that the analysis is based only on counts of the number of users in the system. It is not necessary to estimate the arrival and service rates separately to determine .
4. Conclusions and Final Observations
The problem of traffic intensity estimation in finite Markov queues ( queues) is presented as quite challenging. In fact, no estimator is absolutely superior to another in all parametric space. Although the estimates of the MLE and the bootstrap corrected MLE exhibit less bias, the Bayesian estimates (beta and Jeffreys’ prior) present the lowest MSE in general. Perhaps due to the skewness of the a posterior distribution, the Bayesian estimators do not present low bias. In general, for sample size and queues with , the average estimation error is less than 0.005; this value was only exceeded by the Bayesian estimator for queues with total capacity .
In regard to the behavior of the average estimation error and the average MSE as functions of the traffic intensity , the major errors are observed when the sample size is small () and the traffic intensities are , unlike in the case of the queues, which exhibit higher biases when . Perhaps due to the truncation of the number of users to the maximum queue length, , systems with high traffic intensities require more computational effort and are the most difficult to estimate.
Finally, it is important to note that, for queues with capacity , the system’s behavior is similar to that of an infinite Markovian queue (an queue), as expected. That is, the average estimation error is greater, and the MSE is highest when .
Future work in this area includes testing other Bayesian point estimators (e.g., the median, because of the asymmetry of the a posterior distribution), developing interval estimators, hypothesis testing methods, or even Kernel-based methods [24].
Data Availability
The data used to support the findings of this study are included within the article.
Disclosure
The Brazilian government funding agencies mentioned had no role in the study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
Special thanks are due to Gabriel and Carolina for helping with the algebra. This work was supported by the Brazilian agencies CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico of the Ministry for Science and Technology) [Grant nos. 304671/2014-2, 305841/2016-5] and FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais) [Grant nos. CEX-PPM-00564-17, APQ-02119-15, and BIP-00106-16].