- Email: [email protected]
Bayesian hierarchical models in manufacturing bulk service queues
Journal of Statistical Planning and Inference 136 (2006) 335 – 354 www.elsevier.com/locate/jspi
Bayesian hierarchical models in manufacturing bulk service queues Carmen Armero∗ , David Conesa Department of Statistics and Operations Research, Universitat de València, C Dr. Moliner 50, 46100 Burjassot (Valencia), Spain Received 23 April 2003; received in revised form 21 July 2003; accepted 12 July 2004 Available online 15 September 2004
Abstract In this paper, Queueing Theory and Bayesian statistical tools are used to analyze the congestion of various manufacturing bulk service queues with the same characteristics that are working independently of one another and in equilibrium. Hierarchical models are discussed in order to develop the whole inferential process for the parameters governing the system. Markov Chain Monte Carlo methods and numerical inversion of transforms are addressed to compute the posterior predictive distributions of the usual measures of performance in practice. © 2004 Elsevier B.V. All rights reserved. MSC: 60K25; 62F15; 62M05 Keywords: Hierarchical Bayesian analysis; Bulk service queues; Congestion in manufacturing systems; Markov Chain Monte Carlo integration; Numerical inversion of transforms; Simulation
1. Introduction Nowadays, manufacturers are forced to continuously improve productivity so as to lower their manufacturing costs, while maintaining certain quality standards. In this context, there are usually two possible situations: whilst some are interested in improving the quality and efficiency of a particular manufacturing system that is currently in use, others are ∗ Corresponding author. Tel.: 34-96-3544-309; fax: 34-96-3544-735.
E-mail addresses: [email protected] (C. Armero), [email protected] (D. Conesa).. 0378-3758/$ - see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2004.07.007
336
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
interested in the possible behavior of a new plant or facility to be installed having the same characteristics as the one already running. In both cases, in order to make strategic and tactical decisions, it becomes necessary to have some previous knowledge of the performance (productivity) of the former system. In this paper, we use Queueing Theory and Bayesian statistical tools to analyze the congestion of a single-stage system with various independent bulk service queues (models in which the service facility can serve various customers simultaneously) with the same characteristics and running in equilibrium. This situation describes, for instance, a factory or a warehouse with various equal (or similar) machines (ovens, stamping presses, etc.) with or without human operators, working independently (that is, generating their own lines) and, moreover, operating correctly for a long time (the machines do not necessarily have to be in the same factory, they could equally be in different geographical localizations). There are a wide variety of bulk service models, depending on the specific arrival process of the customers (pieces, jobs, etc.), the service mechanism of the machine, the capacity of the whole system and the queue discipline (see Chaudhry and Templeton (1983) for a detailed description of these models). Nevertheless, in this study, we will only concentrate on the most standard one, usually denoted as M/M(a, b)/1, where the customers arrive singly according to a Poisson process with arrival rate, ; the service distribution is exponential with service rate, ; the single server of the system operates according to the so-called general bulk service rule, in which the batch service size is always at least the quorum, a, and can never exceed the maximum capacity of the server, that is b; the arrival and service processes are independent; the queue discipline for the groups is first-in, first-out (FIFO); and no restrictions on the capacity of the system are considered. Bulk service queues are also very popular in many other practical situations such as transportation systems (buses, guided tours, medical evacuation systems, etc.) and computing (information-transmission systems handling messages, etc.). Some references to the usefulness of these models are given by Griffiths (1995), who analyzes in depth the congestion of the Suez Canal by modeling the throughput of ships as a bulk service queue, and Abellán et al. (2003) and Armero et al. (2003), who deal with the waiting list for a renal transplant in the País Valencià (Spain). Our main objective is to statistically analyze the congestion of a system formed by a set of K independent work stations, each of them generating an M/M(a, b)/1 queueing system in the steady-state (such as, for example, a system of elevators in a skyscraper, a set of furnaces for the curing process of the brake pads in the automotive friction industry, or even a service company offering tourist guides for sightseeing tours around a city). The remainder of the paper is organized as follows. In Section 2, we use Bayesian hierarchical models to develop the whole inferential process for the parameters of the system. Section 3 concentrates on the study and practical computation of the joint congestion of the system, through the posterior predictive distribution of the most relevant measures of performance for each queue (number of customers or jobs in the queue, waiting times in the queue or in the system, busy and idle times, batch service size, etc.). Finally, in order to illustrate the results obtained, Section 4 discusses a fictitious numerical example that considers a plant with various machines.
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
337
2. The Bayesian hierarchical model for estimation In order to learn about the congestion of the system, we first make inference about the parameters governing the whole system, = (1 , . . . , K ) with i = (i , i ), where i and i are, respectively, the arrival and the service rate of the ith queue (i = 1, . . . , K). From a statistical point of view, we subjectively assume that the fact that all queues are consolidated and operating for a long time indicates that they are working at equilibrium and so, the ergodic condition applies for each one of the K queues. Note also that this assumption implies that the parameters can only move freely in the reduced parametric space e = {(1 , . . . , K ); i < bi , i = 1, . . . , K}. The problem of analyzing a queueing system statistically is not a new one, indeed, as Bhat et al. (1997) state, it constitutes a very important field of research in Queueing Theory. Most of the studies in this area are in the framework of the frequentist approach, but lately there has been an increasing interest in Bayesian methods (see Armero and Bayarri (1999) and Conesa (2000) for an exhaustive list of up to date Bayesian references). It is well known that in general the two basic sources of information in Bayesian inference are the prior beliefs and the data. Therefore, in order to estimate , our first step consists of assigning a prior distribution that suitably describes our prior knowledge. The simplest approaches to this problem are either to analyze each queue individually or to carry out a joint analysis without any distinction among the different queues. Nevertheless, we will not pursue these analyses here because, even though machines are identical, there is always variability in their performance (they are usually managed by human operators); and moreover, because it is reasonable to assume that the arrival and service rates of all the queues in the system are related. For this reason, we are going to construct a specific model that allows us to include all these characteristics in the analysis. We start by assuming that the K pairs {i =(i , i ); i =1, . . . , K} corresponding to the K queues are a random sample from a common prior distribution p( | ) that depends on an unknown hyperparametric vector, . This kind of modeling requires a second level in which we must express our knowledge about through another prior distribution, p(), usually known as hyperprior distribution. The resulting structure is usually known as a Bayesian hierarchical model (for a detailed description of these, see, for example, Berger (1985), Gelman et al. (1995) and more recently, Draper (2001)). Note that this model considers the queues exchangeable and they are only independent conditional on . Once we have specified all our prior knowledge of (, ) through the joint prior distribution, p(, ) = p(|)p(), the following stage is to get experimental information. There are various methods to obtain sample information on queues (see Armero and Bayarri (1996) for a revision on this topic). In this paper, we use a simple experiment that consists of observing, simultaneously or not, the arrival and the service process of each one of the queues. In particular, for the ith queue, we record the interarrival times of nai consecutive customers and the times spent in serving nsi groups (where nai and nsi are fixed in advance). It follows that in accordance with the general hypothesis of the M/M(a, b)/1 queue, the nai interarrival times Yi = (Yi1 , . . . , Yinai ) are a random sample from Ex(i ) and the nsi recorded service times Si = (Si1 , . . . , Sinsi ) constitute a random sample from Ex(i ).
338
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
The resulting likelihood function of for all the data z = (z1 , . . . , zK ) (where zi = (yi , si ) represent the data of each queue) is () =
K
i (i , i ) ∝
K i=1
i=1
ni ai e−i tai ni si e−i tsi ,
i < bi , i = 1, . . . , K, (2.1)
ai where tai = nj =1 yij is the total time until the arrival of the nai customer considered in si the ith queue, and tsi = nl=1 sil is the total time taken by the server of the ith queue to complete all the services under consideration. Note that the restrictions in the domain of the likelihood in (2.1) correspond to the ergodic condition of each queue. Next, we consider the problem of assigning a prior distribution to the unknown parameters of the model. As mentioned above, we specify the joint prior distribution of (, ) in two stages. In the first, we assign a prior distribution, p(|), for given the hyperparameters. Then, we select a hyperprior distribution p() for those hyperparameters. For simplicity, we select the conditional prior distribution for each pair i = (i , i ) proportional to the product of two independent Gamma distributions (a member of the common natural conjugate prior distribution that would result from an individual analysis of each queue): p(|) =
K
p(i , i |)
i=1
= C()
K
K
Ga(i |, ) Ga(i |, ),
i < b i , i = 1, . . . , K,
(2.2)
i=1
where the hyperparameters are = (, , , ) and the common proportionality constant corresponding to each queue is ( + ) b (b) −1 C() = F 1, + ; + 1; , ()() (b + )(+) b + where F (a, b; c; z) is the Gauss’s hypergeometric function (see Gradshteyn and Ryzhik (1965) for a description of this function). Note that i and i given are not independent due to the restriction i < b i , i = 1, . . . , K. A hyperprior distribution, p(), has to be chosen. In general, the specification of this type of distribution is a very challenging problem in Bayesian statistics (see, for instance, Berger (1985)). In spite of this, it is worth noting that in situations when the prior knowledge of is vague, we can always use non-informative priors (never forgetting that when using improper prior densities, we must always check if the corresponding posterior distribution is proper). Nevertheless, when dealing with real problems, the experts in the subject should normally have enough information (maybe from historical data) in order to construct a nonvague hyperprior distribution, or at least, to constrain the hyperparameters (appropriately reparameterized) into a finite region (see García-Donato (2001) and Gelman et al. (1995) for examples of this situation). An alternative way of estimating the hyperparameters is to use the so-called empirical Bayes (EB) methods, which basically consist of using the experimental information to
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
339
ˆ of the unknown hyperparameters, , and then to proceed in the produce an estimate, , ˆ as the prior distribution of the parametric usual Bayesian way, but working with p(|) vector . Lehoczky (1990), Thiruvaiyaru and Basawa (1992), Sohn (1996), Armero and Conesa (2000) and Sohn (2000) have applied EB methods to the statistical analysis of different queueing systems. There has been a lot written about the advantages and disadvantages of both procedures and so we will not pursue this question here. Nevertheless, as Berger (1985) states, EB methods work fine when the number of queues is high because the i ’s are a random sample from the general population p(|), and so, the bigger the sample, the better the information about . For small samples, however, these methods do not reflect the errors introduced in the estimation of hyperparameters, and only in those particular situations in which all the ˆ the error would possible values of are, on the whole, tightly grouped in their estimation , be small leading to a good performance of the method. In general, EB methods require less computational effort than the full hierarchical Bayesian model. In this case, as the ergodic condition needs to be considered in the general EB procedure, the marginal distribution of the data dependent on results in terms of Gauss hypergeometric functions, and so, new computational problems appear in the estimation of . Robbins (1983), Maritz and Lwin (1989) and Carlin and Louis (2000) are very good references on the subject. Once the elicitation of p() is completed, the joint prior distribution that contains all our prior knowledge of the parameters and hyperparameters is p(, ) = p(|)p() = p()
K
p(i |)
i=1 K
= C p(, , , )
K
Ga(i |, )Ga(i |, ), i < b i , i = 1, . . . , K.
i=1
(2.3) This prior information is then updated, via Bayes’ theorem, with the experimental information contained in the likelihood () in (2.1), to get the joint posterior distribution: p(, |z) ∝ ()p(|)p() = C K p(, , , )
K i=1
ni ai e−i tai ni si e−i tsi Ga(i |, )Ga(i |, ), (2.4)
where i < bi , for i = 1, . . . , K. The distribution in (2.4) contains all our knowledge of the system, but it does not have a known analytical expression. Therefore, in order to provide more intuitive information about and , numerical methods will be needed to evaluate, in practice, the main characteristics of this posterior distribution. Among others, one feasible (indeed the most popular) possibility is to use Markov chain Monte Carlo (MCMC) methods, that draw samples from any intractable posterior by running a cleverly constructed Markov chain over a long period, the stationary distribution of which, being the one we want to simulate from.
340
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
Among the different ways of building these chains, the most popular are Gibbs sampling and Metropolis–Hastings algorithm (see Gilks et al. (1996) for a detailed discussion of these). When, as in our case, the distribution of interest is multivariate, these methods produce a random sample from the joint posterior distribution p(, |z) in (2.4). This means we can use Monte Carlo integration methods to approximate the main features of the marginal posterior distribution of the arrival, service rates and, consequently, of the traffic intensity of each of the K queues. These marginal distributions will provide us with the first intuitive impression of the congestion and performance of the queues we were looking for. MCMC methods have also been successfully used by Rios Insua et al. (1998) in order to deal with the posterior distribution of the parameters in an M/Hk /1 queueing system.
3. Predicting the whole congestion of the system As we have just commented, MCMC methods allow us to extract almost all the information contained in the joint posterior distribution in (2.4) regarding the parameters and hyperparameters (, ), managing this system. Nevertheless, it is well known that, in Queueing Theory, when a system is working in steady state (like the one we are dealing with), its congestion is better described through the so-called measures of performance of the queue (number of customers in the system, waiting times in the system and in the queue, idle and busy periods, etc.). In general, the posterior predictive distribution of a congestion variable can be obtained through its steady-state distribution given the parameters of the model and their marginal posterior distribution. However, in our system, this general procedure gets very complicated in practice: the conditional distribution of these variables is usually expressed is terms of transforms without analytical inversion, and since we have no closed expression for the corresponding marginal posterior distribution of the parameters we can only work with simulated samples from it. In this section, we will discuss some of the practical questions involved when computing the posterior predictive distribution of the measures most often used to describe congestion in bulk service queues, especially those presenting computational problems.All these measures are, of course, interconnected, and their utility in a real problem will depend on its particular characteristics and the objectives of the manufacturers who are analyzing (or paying for) it. 3.1. The state of the system Let (Hi , N q i ) denote the state of the ith queue where Hi = 1 means that the server is busy processing a batch of s items, a s b, and Hi = 0 indicates that the server is idle because the number of items in the system is less than the quorum a. In both cases, N q i is the number of customers waiting in queue. Given the arrival and service rate, its probability distribution is (Medhi, 1991)
r a+1 − rib+1 a + i p((0, 0)|ri ) = 1 − ri (1 − ri )2
−1 ,
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
p((0, n)|ri ) = p((1, n)|ri ) =
1 − rin+1 1 − ri 1 − rib 1 − ri
341
p((0, 0)|ri ),
n = 1, . . . , a − 1
p((0, 0)|ri ) rin+1 ,
n = 0, 1, . . . ,
where, for simplicity, p((hi , ni )|ri ) denotes the probability p((Hi , N q i ) = (hi , ni )|ri ), and ri is the unique real root in (0, 1) of the equation: i = ri + ri2 + · · · + rib . i
(3.1)
Taking into account that the K queues work independently, given the parameters of the system, the joint conditional distribution that describes the state of the whole system, denoted as p((h1 , n1 ), . . . , (hK , nK )|r1 , . . . , rK ) from now on, is the product of its corresponding marginal distributions, in particular, K i=1 p((hi , ni )|ri ). Therefore, the joint posterior predictive distribution of the state of all the queues is p((h1 , n1 ), . . . , (hK , nK )|z) = E (r1 ,...,rK |z) [p((h1 , n1 ), . . . , (hK , nK )|r1 , . . . , rK )]
K = p((hi , ni )|ri )p(r1 , . . . , rK |z)d(r1 , . . . , rK ).
(3.2)
i=1
As we do not have a closed expression for p(r1 , . . . , rK |z) nor for the joint posterior predictive, we can approximate the probabilities of this distribution, by using the following procedure based on MCMC methods: (1) From the random sample of size m of the joint posterior distribution in (2.4) obtained us(j ) (j ) ing MCMC techniques, consider {(1 , . . . , K ); j = 1, . . . , m}, as a random sample from the marginal posterior distribution of the parameters, p(|z). (j ) (2) Generate a sample from p(r1 , . . . , rK |z), where each generated component ri is ob(j ) (j ) (j ) tained by solving Eq. (3.1) for each simulated pair i = (i , i ), ∀i. (3) Approximate each posterior predictive probability as p((h1 , n1 ), . . . , (hK , nK )|z) m 1 (j ) (j ) ≈ p((h1 , n1 ), . . . , (hK , nK )|r1 , . . . , rK ) m j =1
m K 1 (j ) ≈ p((hi , ni )|ri ). m
(3.3)
j =1 i=1
This procedure is similar to the one used in Wiper (1998) to obtain the posterior predictive distribution of the number of customers in Er /M/1 queueing systems. It is worth
342
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
mentioning that the computational effort of this procedure (and the next two presented in this section) depends mainly on two aspects: the MCMC simulation from the marginal posterior distribution of the parameters and the approximation of the posterior predictive distribution. The MCMC simulation is reasonably easy to handle with any programming language (C + +, Java, etc.) or by using specific software such as BUGS䉸 (Gilks et al., 1994). In this procedure, the latter aspect depends only on the difficulty of generating a random sample from the posterior distribution p(r1 , . . . , rK |z), which it turns into a root finding problem that can be implemented in most programming languages without too much difficulty. In Section 4, we discuss our computational experience with all the procedures presented in this section in more detail. From (3.3) it is easy to compute the posterior distribution of the situation of idleness or busyness of all the K machines. Note also that (3.3) does not provide an intuitive visualization of the congestion of each queue. To do so, it would be more suitable to analyze each one of them separately through their corresponding marginal distribution, p((hi , ni )|z), that could be computed directly, marginalizing (3.3), or using a similar procedure to the one presented above, also based on MCMC methods: m 1 (j ) p((hi , ni )|ri ), p((hi , ni )|z) ≈ m
(3.4)
j =1
(j )
where {ri ; j = 1, . . . , m} is a random sample from the marginal distribution p(ri |z). It would be very useful for manufacturing purposes (to learn about the performance of the system, improve the system behavior, or even for planning and building integrated decision projects, for example) to increase our information about the system. Two ways of doing this will be to study the service batch size (a large value of the expected size would indicate that machines always work at their maximum capacity), or/and the waiting time in the queue for an arriving job (of special importance when dealing with pieces that can be damaged by a long wait). Obviously these special features could make the managers think about the possibility of increasing the quorum of the machines, buying new more powerful machines, etc. For both congestion measures, their conditional distribution depends on the root of Eq. (3.1), and so, in order to approximate their corresponding posterior predictive distribution, we could just reproduce procedures similar to the ones presented above. 3.2. Idle periods As the interest of managers in charge of the system is usually to maximize its output (always with the final objective of maximizing benefits), it could be interesting to study how long the idle periods of the machines that make up the whole system are. Due to the fact that the K queues work independently, given all the parameters of the system, the joint probability density function (p.d.f.) of the duration D of an idle period in each queue will be the product of its corresponding marginal distributions: d(t1 , . . . , tK |1 , . . . , K , r1 , . . . , rK ) =
K i=1
di (ti |i , ri ),
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
343
where for the ith queue, ri is as in (3.1) and the conditional distribution of Di is expressed in terms of the Laplace transform (LT) of their corresponding p.d.f. (Medhi, 1991):
∞ di∗ (u|i , ri ) = e−ut i di (ti |i , ri ) dti 0
ai − [ri (u + i )]a 1 i (1 − ri ) = . (1 − ria ) (u + i )a i − ri (u + i )
(3.5)
This LT cannot be analytically inverted, and so numerical inversion becomes necessary. This problem is not new in Queueing Theory. Nowadays, computational analysis plays a fundamental role when analyzing practical congestion (Neuts (1981), Abate and Whitt (1992) and Tijms (1997) are classical references in the subject). We could approximate the posterior predictive distribution of all the idle periods: K ((r1 ,1 ),...,(rK ,K )|z) d(t1 , . . . , tK |z) = E di (ti |ri , i ) (3.6) i=1
using a procedure based, not only on MCMC methods, but also on algorithms for numerically inverting the LTs involved: (1) Use the results in steps 1 and 2 of the general procedure given in Section 3.1 to compute (j ) (j ) (j ) (j ) (3.2), and consider {(1 , r1 ), . . . , (K , rK ), j = 1, . . . , m} as a random sample from p((1 , r1 ), . . . , (K , rK )|z). (j ) (j ) (j ) (j ) (2) For each simulated vector (i , ri ), compute di (ti |i , ri ) by numerically inverting (j ) (j ) its LT, di∗ (u|i , ri ), given in (3.5). (3) Approximate the p.d.f. of the joint posterior predictive distribution as d(t1 , . . . , tK |z) ≈
m K 1 (j ) (j ) di (ti |i , ri ). m
(3.7)
j =1 i=1
As in the previous case, the computational effort of this procedure depends on the MCMC simulation from the marginal posterior distribution of the parameters and the approximation of the marginal conditional distribution di (ti |i , ri ), which now incorporates the additional computational threat of inverting LTs mentioned above. There are many algorithms for numerically inverting LTs and they are not too difficult to implement in any programming language. Abate and Whitt (1992) provide a detailed discussion of most of them and plenty of references on the subject. It is worth mentioning, that although these are very fast algorithms, the biggest computational cost for approximating the marginal posterior predictive distribution comes from the fact that we have to numerically invert an LT for each simulated vector obtained via MCMC. It is clear that the corresponding marginal posterior predictive distribution for the idle period of each machine, more easily computed as di (ti |z) ≈
m 1 (j ) (j ) di (ti |i , ri ), m j =1
(3.8)
344
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354 (j )
(j )
(j )
(j )
where di (ti |i , ri ) is obtained by numerically inverting di∗ (u|i , ri ) as in step 2 above, will provide a convenient visualization of the congestion of each machine. Working in a similar way, we could also analyze the number of items in the queue when an idle period begins. Given all the parameters of the system, the joint distribution of the number of items waiting to be processed by each machine when an idle period begins, p(N d 1 = n1 , . . . , Nd K = nK |r1 , . . . , rK ), with 0 ni < a, i = 1, . . . , K, is the product of their corresponding marginal distributions, p(N d i = ni |ri ), now expressed in terms of their probability generating function (p.g.f.) (Medhi, 1991):
G(Nd i |ri ) (u) =
∞
p(N d i = ni |ri )uni =
ni =0
(1 − ri )(1 − (ri u)a ) . (1 − ria )(1 − ri u)
(3.9)
In order to compute the corresponding joint posterior predictive distribution, we could use a similar procedure to the one presented in this section to deal with idle periods in which, instead of inverting LTs, we would need methods to numerically invert p.g.f.s.
3.3. Busy periods Another complementary way of describing the congestion of our system is to analyze the duration, B, of the busy periods. As usual, its joint conditional (on the parameters) distribution is the product of its corresponding marginal distributions, the p.d.f.s of which are given in terms of their LTs (Medhi, 1991): bi∗ (u|i , Ri ) =
i (1 − Ria ) , u + i − i Ria
(3.10)
where here Ri is the unique real root in (0, 1) of the equation: i Rib+1 − (u + i + i )Ri + i = 0.
(3.11)
In order to approximate the joint posterior predictive distribution of the duration of a busy period for each machine, b(t1 , . . . , tK |z), we could use a similar procedure (with similar computational behavior) to the one presented in Section 3.2 for the idle periods: (1) From the generated sample of the joint posterior distribution p(, ) in (2.4) consider (j ) (j ) {(1 , . . . , K ); j = 1, . . . , m} as a random sample from the marginal posterior distribution p(|z). (j ) (j ) (j ) (2) For each simulated pair i , compute bi (ti |i , Ri ) by numerical inversion of its (j ) (j ) corresponding LT bi∗ (u|i , Ri ), taking into account that Eq. (3.11) must be solved at every evaluation of the LT necessary for its numerical inversion.
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
345
(3) Approximate the p.d.f. of the joint posterior predictive distribution as m
b(t1 , . . . , tK |z) ≈
1 (j ) (j ) (j ) (j ) b(t1 , . . . , tK |1 , . . . , K , R1 , . . . , RK ) m j =1
=
m K 1 (j ) (j ) bi (ti |i , Ri ). m
(3.12)
j =1 i=1
As in the previous measures, in order to have a more convenient visualization of the congestion of each queue, we can consider their corresponding marginal posterior predictive distribution, bi (ti |z), which could be approximated as bi (ti |z) ≈
m 1 (j ) (j ) bi (ti |i , Ri ), m
(3.13)
j =1
(j )
(j )
where each bi (ti |i , Ri ) is obtained by numerically inverting its corresponding LT (j ) (j ) bi∗ (u|i , Ri ) as stated in stage 2. 4. A numerical example with machines and facilities In this final section, with the main objective of illustrating the previously presented results, we analyze the congestion in a manufacturing scenario formed by K = 6 independent machines that work processing items in packets according to an M/M(9, 20)/1 queue model. For each machine i in the system, i = 1, . . . , 6, Table 1 displays the results of a simulated experiment like the one proposed in Section 2. In particular, it contains the interarrival time of nai = 20 consecutive items, the time needed to process nsi = 20 packets of pieces and the maximum likelihood estimate (m.l.e.) ˆ i of the traffic intensity, obtained by maximizing the likelihood function in (2.1) subject to the ergodic condition i < 20i (in our case using MATLAB䉸 ). Note that, although we are considering all the machines to be similar and working at equilibrium, they produce quite different results. In particular, we observe a high level of congestion in the performance of machines 4 and 5 (in both cases, the corresponding Table 1 Sufficient statistics for the arrival and service rate for each machine corresponding to the total time (tai ) to achieve the arrival of 20 consecutive pieces and the processing time (tsi ) for 20 packets of pieces, respectively, and the m.l.e. ( ˆ i ) of the traffic intensity for each machine Machine i
1
2
3
4
5
6
tai tsi ˆ i
13.9 158.8 0.5712
9.4 186.7 0.9931
11.6 135.9 0.5858
8.4 195.5 0.9999
10.4 302.4 0.9999
9.2 169.9 0.9233
346
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354 6 Machine 1 Machine 2 Machine 3 Machine 4 Machine 5 Machine 6
5 4 3 2 1 0
(a)
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4 0.6 Unique machine
0.8
1.0
6 5 4 3 2 1 0
(b)
Fig. 1. Marginal posterior distribution of the traffic intensity of each machine when they are analyzed: (a) separately; (b) together without any distinction among them.
m.l.e.s, ˆ 4 and ˆ 5 , are practically one), while machines 1 and 3 seem to have a very good performance ( ˆ 1 = 0.5712 and ˆ 3 = 0.5858). It is important to point out that if we did not distinguish among the different machines, the m.l.e. of the common would be ˆ = 0.9135, an intermediate level between the lower values for machines 1 and 3 and the higher ones of machines 4 and 5. With the sole idea of providing a reference framework, we start the statistical analysis of the congestion of this system by considering the two most extreme possible situations: each machine is analyzed completely separated, sharing no common information with the others, or, all machines are identical (i = ) with customers arriving at each machine with identical rates (i = ). Firstly, we present the unpooled analysis. Fig. 1(a) shows the marginal posterior distribution of the traffic intensity for each machine. In particular, the posterior distribution corresponding to machine i has been obtained assigning a Jeffrey’s non-informative prior for the corresponding parametric vector i and combining, via Bayes’ theorem, the data collected only on that machine. In the case of analyzing all the machines together,
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
347
Fig. 1(b) displays the posterior distribution of the traffic intensity, computed from the posterior distribution of the common arrival and the service rate, both computed using a Jeffrey’s non-informative prior and the data collected in all the machines (a total of 120 interarrival times and 120 service times). In both graphics the assumption of equilibrium means that the right tail of all the densities does not approach zero when the traffic intensity reaches to the critical value 1 (this special feature is discussed in depth in Armero and Bayarri, 1994). Note also, that the different sample sizes considered when computing both posterior distributions are evident in the smaller variability of the posterior distribution of traffic intensity in the pooled analysis. After these two extreme analyses, we will analyze the congestion of our manufacturing system by using a full Bayesian hierarchical model that allows us to maintain the particular features of each machine (and/or its human operator) and also incorporates those common to all of them. To do so, we assume that the parameters {i = (i , i ); i = 1, . . . , 6} of each machine are a random sample from the common population distribution in (2.2): p(i , i |) ∝ Ga(i |, )Ga(i |, ),
i < 20i ,
(4.1)
where = (, , , ) represents the vector of unknown hyperparameters. We now express our previous knowledge about . At this stage, with the sole objective of better understanding how Bayesian hierarchical models work in queues, we consider two possible situations: one in which our knowledge about could be considered as scarce, being our selection for the hyperprior distribution, p(), the product of four (one for each hyperparameter) exponential distributions with mean 1/0.001, Ex(0.001), and a second, in which we express a greater knowledge of with the product of four independent Ex(0.1). From this point onwards, we will refer to this first situation as LIH (least informative hyperprior), and to the second as MIH (most informative hyperprior). In both cases, we have assumed independence among the hyperparameters, and also chosen proper distributions in order to obtain proper joint posterior distributions. We have selected the exponential density for the hyperprior for all parameters, because, in addition to the uniform and the Gamma distributions, this is very common when working with exponential data. In both situations, the joint posterior distribution of the parameters and hyperparameters, p(, |z), although proper, is analytically intractable. Therefore, we will approximate these distributions by simulation using Markov Chain techniques. One way to do this is via the Gibbs sampler (see, for instance Casella and George (1992), for a detailed description of this method), implemented by Gilks et al. (1994) in the BUGS 䉷 (Bayesian inference Using Gibbs Sampling) program. In particular, we work with the Windows 䉷 version, WinBUGS, developed by Spiegelhalter et al. (1999) (see http://www.mrc-bsu.cam.ac. uk/bugs/ for more information about this project). (0) First we choose a vector of starting values for the parameters (for instance, i =0.17 and (0) i =0.01, i=1, . . . , 6) and the hyperparameters (for instance, (0) =(0) =(0) =(0) =1000 for the LIH case, and (0) = (0) = (0) = (0) = 10 for the MIH case). Then, we decide the duration of the burn-in period (the instant when we think that convergence to the stationary distribution has been reached) to be the first 5000 simulations and to work with the next 10,000 ones. In those situations in which convergence could not be guaranteed, we select every 50th iteration of each chain to contribute to the analysis (usually called a thin value of
348
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
6 Machine 1 Machine 2 Machine 3 Machine 4 Machine 5 Machine 6
5 4 3 2 1 0 0.0
0.2
(a)
0.4 0.6 Product of Ex(0.1)
0.8
1.0
0.4 0.6 Product of Ex(0.001)
0.8
1.0
6 Machine 1 Machine 2 Machine 3 Machine 4 Machine 5 Machine 6
5 4 3 2 1 0 0.0 (b)
0.2
Fig. 2. Estimated marginal posterior distribution of the traffic intensity of each machine in the MIH case (at the top) and in the LIH (at the bottom).
50). The results obtained with WinBUGS are a random sample (each element is a vector with 16 components, 2 × 6 parameters + 4 hyperparameters) from the joint posterior distribution p(, |z). Convergence diagnosis of the resulting Markov chain is a major part of these kind of simulation techniques. To do so, we use CODA, a menu-driven set of S-Plus 䉷 functions which serve as an output processor for the BUGS 䉷 software (Best et al., 1995). Following the authors’recommendation, we have used a combination of diagnostic methods plus visual inspection of the trace plot and summary statistics to asses the convergence of the random samples generated by the WinBUGS software for all the parameters of this model with a reasonable degree of confidence. From the obtained random sample of the joint posterior distribution p(, |z), we get a (j ) random sample { i , j = 1, . . . , 10, 000} from the marginal posterior distribution of the traffic intensity for machine i, p( i |z), i = 1, . . . , 6. Fig. 2 shows a kernel estimation of each one of these marginal distributions in both situations of prior knowledge about
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
349
considered. Note that, in the LIH case, the resulting densities are quite similar, indicating that the hierarchical structure emphasises to the common characteristics of the six machines than the particular differences among them. Nevertheless, when our knowledge about is more accurate (the MIH case), the corresponding densities show greater variety because now the hierarchical model recognizes and distinguishes with greater intensity the individual performance of each machine. However, if we compare the results in Fig. 1, where all the machines were analyzed separately, with the ones in the MIH situation, we observe how the model takes into account the common characteristics of the machines by greatly smoothing the high level of congestion on machine 5 and displacing the values of the traffic intensity on machines 1 and 3 toward more congested regions. We continue with our analysis of the congestion by considering some of the measures of performance studied in Section 3. For simplicity, we only deal here with the marginal posterior distribution of the state of the system, the waiting time in the queue and the duration of the busy periods. All of them are analyzed in the unpooled case and also in the two applications of the Bayesian hierarchical model considered. In these two latter situations, the marginal posterior distributions have been approximated using a random sample of the joint posterior distribution obtained with WinBUGS and also, when needed, using methods for numerically inverting the transforms involved. In particular, we work with the Lattice-Poisson algorithm for inverting probability generating functions and the Euler algorithm to deal with Laplace transforms. Both methods are based on the Fourier series method and have been extensively discussed by Abate and Whitt (1992). Following these authors’ recommendation, we also always use the Gaver–Stehfest (1966,1970) algorithm (based on the Post–Widder (1930,1934) formula and also discussed in depth in Abate and Whitt, 1992) to check the results given by the Euler algorithm. As an exception, because the Laplace transform of the conditioned distribution of the busy period is given in terms of a real root of an equation, we only use the Gaver–Stefesth algorithm (which works with real variable) instead of the Euler algorithm (which works with complex variable) to invert this particular transform. It is worth mentioning that these algorithms consist of evaluating the expression of the transform several times and so, they do not present many difficulties with high values of a and b. In the unpooled analysis of the machines, the calculus of the relevant predictive distribution for each machine has been performed using similar procedures to those in the last section based on Monte Carlo integration. Nevertheless, the expression of the joint posterior distribution of the arrival and service rate usually has a nice form (for instance, as a result of a conjugate analysis) from which it is easy to simulate without the need of MCMC. In relation to the methods used to find roots of the equations involved in these procedures, good results with moderate values of the quorum b motivated us to choose a hybrid algorithm based on a combination of Bisection and Newton–Raphson algorithms (see Press et al. (1986) for a description of this routine). Obviously, this is our particular choice here but other alternatives could have been used. For instance, in order to compute steady-state probabilities of M/G(a, b)/1 queues, Chaudhry et al. (1987) use the Muller’s method, adapted from the version in Conte and de Boor (1972), over any other for its good behavior with big values of b. Fig. 3 displays, the posterior predictive distribution of the state of each machine, for all three situations of prior knowledge considered. The bars in each graphic are the probabilities
350
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
0.06 0.04
0.04
0.02
0.04
0.04
0.1724
0.04
0.04
0 2 4 6 8
0.00 0 2 4 6 8
1 3 5 7 9 11 14
0.0540
0.04
0.04
0 2 4 6 8
0.06
0.00 0 2 4 6 8
1 3 5 7 9 11 14
0.04
0.02
0.02
0.00
0.00 0 2 4 6 8
0.04
0.0785
0.04
0.00
0 2 4 6 8
1 3 5 7 9 11 14
1 3 5 7 9 11 14
0.06 0.0836
0.04
0.0549
0.02
0.00
0.00
0 2 4 6 8
1 3 5 7 9 11 14
0.02
0.02
0.0373
0.02
0.06
0.06
1 3 5 7 9 11 14
0.06 0.0345
0 2 4 6 8
1 3 5 7 9 11 14
0 2 4 6 8
1 3 5 7 9 11 14
0.06 0.0402
0.0505
0.02
0.00
0.00
1 3 5 7 9 11 14
0.06 0.0599
0.02
0.02
0 2 4 6 8
1 3 5 7 9 11 14
0.06
0.06
0.0621
0.02
0.00
0.00
1 3 5 7 9 11 14
0.06 0.1565
0.02
0.02
0 2 4 6 8
1 3 5 7 9 11 14
0.06
0.06
0.0524
0.00 0 2 4 6 8
1 3 5 7 9 11 14
1 3 5 7 9 11 14
0.02
0.00 0 2 4 6 8
0 2 4 6 8
0.06 0.0715
0.02
0.00
0.04
0.00 1 3 5 7 9 11 14
0.06 0.0695
0.03
0.0591
0.02 0 2 4 6 8
1 3 5 7 9 11 14
0.06
0.04
0.04
0.00 0 2 4 6 8
0.04
0.1495
0.02
0.00
0.04
0.06
0.06 0.1806
0.00 0 2 4 6 8
1 3 5 7 9 11 14
0 2 4 6 8
1 3 5 7 9 11 14
Fig. 3. Marginal posterior distribution of the state of each of the six machines working in the system. The bar charts in the column on the left correspond to the separate analysis; in the central column to the MIH Bayesian hierarchical model case; and the column on the right to the LIH Bayesian hierarchical model case. In each graphic, the number in the top left corner is the probability that the machine is idle.
of the number, Nq , of pieces waiting in the queue to be processed when the machine remains idle (probabilities on the left, computed from 0 to 8 pieces, one unit less than the quorum a = 9) and when the machine is busy (probabilities on the right, computed for simplicity until the value Nq = 14 pieces). The probability that the machine is idle is also presented inside the frame in the top left corner. Note that graphics in Fig. 3 show the same general behavior and relationships as those previously observed when dealing with the traffic intensity. In the unpooled analysis (left column), machines 1 and 3 seem to be the least congested (idle about 18.06% and 17.24% of the time, respectively), while machines 4 and 5 are the most congested (with approximate values of idleness of 5.40% and 4.02%, respectively). Obviously, if this was a real problem
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
351
these results would make the people in charge of this facility investigate the particular characteristics of machines 1 and 3 and think over future actions to increase their efficiency. The partial pooling situation, which incorporates the common characteristics among the machines through the hierarchical Bayes modeling, is also presented in Fig. 3: bar charts in the central column correspond to the MIH case, and the ones in the right column correspond to the less informative hyperprior (LIH) distribution. If we compare them with the results of the separate analysis, we observe that the hierarchical model tends to make the performance of the machines more similar. Note that the hierarchical model provides an analysis that jointly integrates the particular peculiarities of each machine (human operator, geographical localization, reliability, ageing, etc.) with all their common characteristics. The statistical analysis is computationally more complex but also more flexible and accurate for the current operation of the system. In the LIH case, this effect is more intense in machines 1 and 3. They are idle for about 5.91% and 6.21% of the available time, quantities that are only slightly bigger than machines 4 and 5, the two most congested ones (with approximate percentages of 3.73% and 3.99%, respectively). In the same way, if we compare the behavior of each machine in the three situations considered, the differences in machines 1 and 3 are spectacular. For instance, the predictive probability that machine 1 is idle is 0.1806 in the separate analysis, 0.1495 in the MIH hierarchical model and 0.0591 in the LIH. But note that we can hardly observe any differences among the rest of the machines. The probability of being idle for the two most congested ones differs by less than 1%. For instance, in the most congested machine 5, in the unpooled analysis the probability of being idle is 0.0402, while in the LIH analysis it is 0.037.3 and 0.0345 in the MIH. For all cases and machines, Fig. 4 displays the marginal posterior distribution (expressed through the cumulative distribution function) of the waiting time in the queue of a piece before being processed. In order to clarify their behavior, we have also included the median values of these distributions. In these graphics, we can observe in general the same type of results and relationships as the ones previously commented on relation to the state of the system. Nevertheless, through basic observation of the graphics corresponding to machines 1, 3 and 5, it seems that the waiting time in the queue, as a measure of performance of a queueing system, is less sensitive to the attempts of the Bayesian hierarchical model (in the LIH when the information about the hyperparameters is scarce) to make the machines uniform than it was in the case of the state of the system previously discussed. Another relevant thing to point out is that the posterior probability of a piece being processed immediately, without waiting in the queue is so small that in all the graphics it seems to be zero. This predictive probability is the posterior probability that the machine is idle and, at the same time, there are eight (the quorum of the machine minus one) pieces waiting in the queue for the arrival of a new piece that allows the machine to process a packet with nine pieces, that is p((0, 8)|z). We finish our predictive analysis of this system with the study of the duration of the busy periods of each machine. Fig. 5 displays the cumulative distribution function of the marginal posterior distribution of this variable in all the situations considered, including in each case the median busy time. We can observe small differences among the six machines, although now machine 5 looks a bit different from the others. It is also worth noting that the separate analysis and the two applications of the Bayesian hierarchical model considered
352
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4
MIH LIH Separately
0.2
Med. = 7.7 Med. = 17.0 Med. = 7.2
0.0
0.4 0.2
Med. = 15.9 Med. = 19.3 Med. = 16.9
20
30 40 Machine 1
50
20
30 40 Machine 2
50
Med. = 18.9 Med. = 20.1 Separately Med. = 21.7 MIH
LIH
0.8
10
60
1.0
1.0
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0 20
30 40 Machine 4
50
60
30 40 Machine 3
50
60
0.8
0.6
10
20
1.0 Med. = 42.0 MIH Med. = 26.2 LIH Separately Med. = 40.5
0.6
0.0
Med. = 6.8 MIH Med. = 16.1 LIH Separately Med. = 6.5
0.0 10
60
0.4 0.2
0.0 10
0.8
MIH LIH Separately
Med. = 13.0 MIH Med. = 18.4 LIH Separately Med. = 14.1
0.0 10
20
30 40 Machine 5
50
60
10
20
30 40 Machine 6
50
60
Fig. 4. Posterior cumulative distribution function of the waiting time in the queue for a piece and its median values for each of the six machines in the separate analysis and in the LIH and MIH applications of the Bayesian hierarchical model considered. 1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4
Med. = 6.9 MIH Med. = 10.6 LIH Separately Med. = 6.6
0.2 0.0
0.8
Med. = 9.6 MIH Med. = 11.2 LIH Separately Med. = 9.6
0.2
20
30 40 Machine 1
50
60
10
1.0
Med. = 10.2 MIH Med. = 11.2 LIH Separately Med. = 10.3
0.8
0.4
20
30 40 Machine 2
50
60
10
0.6
0.4
0.4
0.4
0.2
0.2
0.2
30 40 Machine 4
50
60
50
60
Med. = 8.7 MIH Med. = 10.8 LIH Separately Med. = 8.7
0.0
0.0 20
30 40 Machine 3
0.8
0.6
10
20
1.0 Med. = 15.0 MIH Med. = 12.6 LIH Separately Med. = 15.8
0.6
0.0
Med. = 6.0 MIH Med. = 10.3 LIH Separately Med. = 5.7
0.2 0.0
0.0 10
1.0
0.4
10
20
30 40 50 Machine 5
60
10
20
30 40 Machine 6
50
60
Fig. 5. Marginal posterior distribution of the duration of a busy period and the busy median time for each of the six machines in the separate analysis and in the LIH and MIH applications of the Bayesian hierarchical model considered.
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
353
provide very similar results in the different machines of the manufacturing system. These features could indicate that the duration of a busy period is a robust measure of performance in relation to the use of an unpooled or a partially pooled analysis as explored here using the Bayesian hierarchical model.
Acknowledgements The authors would like to thank two anonymous referees for their interest and careful reading of the paper. Their helpful comments have lead to this improved final version.
References Abate, J., Whitt, W., 1992. The Fourier-series method for inverting transforms of probability distributions. Queueing Systems 10, 5–87. Abellán, J.J., Armero, C., Conesa, D., Pérez-Panadés, J., Zurriaga, O., Martínez-Beneito, M.A., Vanaclocha, H., García-Blasco, M.J., 2003. Analysis of the Renal Transplant Waiting List in the País Valencià (Spain). Technical Report TR06-2003, Departament d’Estadística i Investigació Operativa, Universitat de València. Available online at the following url adress http://matheron.uv.es/investigar/technicals.html. Armero, C., Bayarri, M.J., 1994. Prior assessments for prediction in queues. The Statistician 43, 139–153. Armero, C., Bayarri, M.J., 1996. Bayesian questions and answers in queues. Bernardo, J.M., Berger, J.O., Dawid, A.P., Smith, A.F.M. (Eds.), Bayesian Statistics, vol. 5. Oxford University Press, Oxford, pp. 3–23. Armero, C., Bayarri, M.J., 1999. Dealing with uncertainties in queues and networks of queues: a Bayesian approach. In: Ghosh, S. (Ed.), Multivariate, Design and Sampling. Marcel Dekker, New York, pp. 579–608. Armero, C., Conesa, D., 2000. Prediction in Markovian bulk arrival queues. Queueing Systems 34, 327–350. Armero, C., Abellán, J.J., Conesa, D., Pérez-Panadés, J., Martínez-Beneito, M.A., Zurriaga, O., Vanaclocha, H., García-Blasco, M.J., 2003. Waiting for a kidney in the País Valencià (Spain). Technical Report TR052003, Departament d’Estadística i Inv. Operativa, Universitat de València. Available online at the url adress http://matheron.uv.es/investigar/technicals.html. Berger, J.O., 1985. Statistical Decision Theory and Bayesian Analysis, Springer, New York. Best, N.G., Cowles, M.K., Vines, S.K., 1995. CODA Convergence Diagnosis and Output Analysis Software for Gibbs Sampler Output: Version 0.3. Bhat, U.N., Miller, G.K., Rao, S.S., 1997. Statistical analysis of queueing systems. In: Dshalalow, J.H. (Ed.), Frontiers in Queueing. CRC Press, Boca Raton, FL (Chapter 13). Carlin, B.P., Louis, T.A., 2000. Bayes and Empirical Bayes Methods for Data Analysis, Chapman & Hall/CRC Press, London. Casella, G., George, E.I., 1992. Explaining the Gibbs sampler. Amer. Statist. 46 (3), 167–174. Chaudhry, M.L., Templeton, J.G.C., 1983. A First Course in Bulk Queues, Wiley, New York. Chaudhry, M.L., Madill, B.R., Brière, G., 1987. Computational analysis of steady-state probabilities of M/Ga,b /1 and related nonbulk queues. Queueing Systems 2, 93–114. Conesa, D., 2000. Inferencia y predicción en colas con ingresos o llegadas en grupos. Ph.D. Thesis, Universitat de València (in Spanish). Conte, S.D., de Boor, C., 1972. Elementary Numerical Analysis, McGraw-Hill, New York. Draper, D., 2001. Bayesian Hierarchical Modelling, Preliminary draft available in the url address http://www.bath.ac.uk/m˜asdd. García-Donato, G., 2001. Control de procesos. Alternativas Empírico-Bayes y Bayesianas a los métodos clásicos, Master Degree work, Universitat de València (in Spanish). Gaver Jr., D.P., 1966. Observing stochastic processes and approximate transform inversion. Oper. Res. 14, 444–459. Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 1995. Bayesian Data Analysis, Chapman & Hall, London.
354
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
Gilks, W.R., Thomas, A., Spiegelhalter, D.J., 1994. A language and program for complex Bayesian modelling. The Statistician 43, 169–178. Gilks, W.R., Richardson, S., Spiegelhalter, D.J., 1996. Markov Chain Monte Carlo in Practice, Chapman & Hall, London. Gradshteyn, I.S., Ryzhik, I.M., 1965. Table of Integrals, Series and Products, Academic Press, Inc., Boston. Griffiths, J.D., 1995. Queueing at the Suez Canal. J. Opl. Res. Soc. 46, 1299–1309. Lehoczky, J., 1990. Statistical methods. In: Heyman, D.P., Sobel, M.J. (Eds.), Handbooks in O.R. and M.S. Elsevier Science Publishers B.V. (North-Holland), Amsterdam. (Chapter 6). Maritz, J.S., Lwin, T., 1989. Empirical Bayes Methods, Chapman & Hall, London. Medhi, J., 1991. Stochastic Models in Queueing Theory, Academic Press, Inc., Boston. Neuts, M.F., 1981. Matrix-Geometric Solutions in Stochastic Models—An Algorithmic Approach, The John Hopkins University Press, Baltimore. Post, E.L., 1930. Generalized differentiation. Trans. Amer. Math. Soc. 32, 723–781. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1986. Numerical Recipes, Cambridge University Press, Cambridge. Rios Insua, D., Wiper, M., Ruggieri, F., 1998. Bayesian analysis of M/Er /1 and M/Hk /1 queues. Queueing Systems 30, 289–308. Robbins, H., 1983. Some thoughts on empirical Bayes estimation. Ann. Stat. 713–723. Spiegelhalter, D.J., Thomas, A., Best, N.G., 1999. Winbugs Version 1.2 User Manual. MRC Biostatistics Unit. Stehfest, H., 1970. Algorithm 368. Numerical inversion of Laplace transforms. Commun. ACM 13, 479–491 (erratum 13, 624). Thiruvaiyaru, D., Basawa, I.V., 1992. Empirical Bayes estimation for queueing systems and networks. Queueing Systems 11, 179–202. Tijms, H.C., 1997. Computational methods in queueing. In: Dshalalow, J.H. (Ed.), Frontiers in Queueing. CRC Press, Boca Raton, FL. (Chapter 12). Widder, D.V., 1934. Numerical inversion of the Laplace integral and the related moment problem. Trans. Amer. Math. Soc. 36, 107–200. Wiper, M., 1998. Bayesian analysis of Er /M/1 and Er /M/c queues. J. Statist. Planning Inference 69, 65–79. Young Sohn, S., 1996. Empirical Bayesian analysis for traffic intensity: M/M/1 queues with covariates. Queueing Systems 22, 383–401. Young Sohn, S., 2000. Robust design of server capability in M/M/1 queues with both partly random arrival and service queues. Comput. Oper. Res. 29, 433–440.
Bayesian hierarchical models in manufacturing bulk service queues Carmen Armero∗ , David Conesa Department of Statistics and Operations Research, Universitat de València, C Dr. Moliner 50, 46100 Burjassot (Valencia), Spain Received 23 April 2003; received in revised form 21 July 2003; accepted 12 July 2004 Available online 15 September 2004
Abstract In this paper, Queueing Theory and Bayesian statistical tools are used to analyze the congestion of various manufacturing bulk service queues with the same characteristics that are working independently of one another and in equilibrium. Hierarchical models are discussed in order to develop the whole inferential process for the parameters governing the system. Markov Chain Monte Carlo methods and numerical inversion of transforms are addressed to compute the posterior predictive distributions of the usual measures of performance in practice. © 2004 Elsevier B.V. All rights reserved. MSC: 60K25; 62F15; 62M05 Keywords: Hierarchical Bayesian analysis; Bulk service queues; Congestion in manufacturing systems; Markov Chain Monte Carlo integration; Numerical inversion of transforms; Simulation
1. Introduction Nowadays, manufacturers are forced to continuously improve productivity so as to lower their manufacturing costs, while maintaining certain quality standards. In this context, there are usually two possible situations: whilst some are interested in improving the quality and efficiency of a particular manufacturing system that is currently in use, others are ∗ Corresponding author. Tel.: 34-96-3544-309; fax: 34-96-3544-735.
E-mail addresses: [email protected] (C. Armero), [email protected] (D. Conesa).. 0378-3758/$ - see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2004.07.007
336
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
interested in the possible behavior of a new plant or facility to be installed having the same characteristics as the one already running. In both cases, in order to make strategic and tactical decisions, it becomes necessary to have some previous knowledge of the performance (productivity) of the former system. In this paper, we use Queueing Theory and Bayesian statistical tools to analyze the congestion of a single-stage system with various independent bulk service queues (models in which the service facility can serve various customers simultaneously) with the same characteristics and running in equilibrium. This situation describes, for instance, a factory or a warehouse with various equal (or similar) machines (ovens, stamping presses, etc.) with or without human operators, working independently (that is, generating their own lines) and, moreover, operating correctly for a long time (the machines do not necessarily have to be in the same factory, they could equally be in different geographical localizations). There are a wide variety of bulk service models, depending on the specific arrival process of the customers (pieces, jobs, etc.), the service mechanism of the machine, the capacity of the whole system and the queue discipline (see Chaudhry and Templeton (1983) for a detailed description of these models). Nevertheless, in this study, we will only concentrate on the most standard one, usually denoted as M/M(a, b)/1, where the customers arrive singly according to a Poisson process with arrival rate, ; the service distribution is exponential with service rate, ; the single server of the system operates according to the so-called general bulk service rule, in which the batch service size is always at least the quorum, a, and can never exceed the maximum capacity of the server, that is b; the arrival and service processes are independent; the queue discipline for the groups is first-in, first-out (FIFO); and no restrictions on the capacity of the system are considered. Bulk service queues are also very popular in many other practical situations such as transportation systems (buses, guided tours, medical evacuation systems, etc.) and computing (information-transmission systems handling messages, etc.). Some references to the usefulness of these models are given by Griffiths (1995), who analyzes in depth the congestion of the Suez Canal by modeling the throughput of ships as a bulk service queue, and Abellán et al. (2003) and Armero et al. (2003), who deal with the waiting list for a renal transplant in the País Valencià (Spain). Our main objective is to statistically analyze the congestion of a system formed by a set of K independent work stations, each of them generating an M/M(a, b)/1 queueing system in the steady-state (such as, for example, a system of elevators in a skyscraper, a set of furnaces for the curing process of the brake pads in the automotive friction industry, or even a service company offering tourist guides for sightseeing tours around a city). The remainder of the paper is organized as follows. In Section 2, we use Bayesian hierarchical models to develop the whole inferential process for the parameters of the system. Section 3 concentrates on the study and practical computation of the joint congestion of the system, through the posterior predictive distribution of the most relevant measures of performance for each queue (number of customers or jobs in the queue, waiting times in the queue or in the system, busy and idle times, batch service size, etc.). Finally, in order to illustrate the results obtained, Section 4 discusses a fictitious numerical example that considers a plant with various machines.
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
337
2. The Bayesian hierarchical model for estimation In order to learn about the congestion of the system, we first make inference about the parameters governing the whole system, = (1 , . . . , K ) with i = (i , i ), where i and i are, respectively, the arrival and the service rate of the ith queue (i = 1, . . . , K). From a statistical point of view, we subjectively assume that the fact that all queues are consolidated and operating for a long time indicates that they are working at equilibrium and so, the ergodic condition applies for each one of the K queues. Note also that this assumption implies that the parameters can only move freely in the reduced parametric space e = {(1 , . . . , K ); i < bi , i = 1, . . . , K}. The problem of analyzing a queueing system statistically is not a new one, indeed, as Bhat et al. (1997) state, it constitutes a very important field of research in Queueing Theory. Most of the studies in this area are in the framework of the frequentist approach, but lately there has been an increasing interest in Bayesian methods (see Armero and Bayarri (1999) and Conesa (2000) for an exhaustive list of up to date Bayesian references). It is well known that in general the two basic sources of information in Bayesian inference are the prior beliefs and the data. Therefore, in order to estimate , our first step consists of assigning a prior distribution that suitably describes our prior knowledge. The simplest approaches to this problem are either to analyze each queue individually or to carry out a joint analysis without any distinction among the different queues. Nevertheless, we will not pursue these analyses here because, even though machines are identical, there is always variability in their performance (they are usually managed by human operators); and moreover, because it is reasonable to assume that the arrival and service rates of all the queues in the system are related. For this reason, we are going to construct a specific model that allows us to include all these characteristics in the analysis. We start by assuming that the K pairs {i =(i , i ); i =1, . . . , K} corresponding to the K queues are a random sample from a common prior distribution p( | ) that depends on an unknown hyperparametric vector, . This kind of modeling requires a second level in which we must express our knowledge about through another prior distribution, p(), usually known as hyperprior distribution. The resulting structure is usually known as a Bayesian hierarchical model (for a detailed description of these, see, for example, Berger (1985), Gelman et al. (1995) and more recently, Draper (2001)). Note that this model considers the queues exchangeable and they are only independent conditional on . Once we have specified all our prior knowledge of (, ) through the joint prior distribution, p(, ) = p(|)p(), the following stage is to get experimental information. There are various methods to obtain sample information on queues (see Armero and Bayarri (1996) for a revision on this topic). In this paper, we use a simple experiment that consists of observing, simultaneously or not, the arrival and the service process of each one of the queues. In particular, for the ith queue, we record the interarrival times of nai consecutive customers and the times spent in serving nsi groups (where nai and nsi are fixed in advance). It follows that in accordance with the general hypothesis of the M/M(a, b)/1 queue, the nai interarrival times Yi = (Yi1 , . . . , Yinai ) are a random sample from Ex(i ) and the nsi recorded service times Si = (Si1 , . . . , Sinsi ) constitute a random sample from Ex(i ).
338
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
The resulting likelihood function of for all the data z = (z1 , . . . , zK ) (where zi = (yi , si ) represent the data of each queue) is () =
K
i (i , i ) ∝
K i=1
i=1
ni ai e−i tai ni si e−i tsi ,
i < bi , i = 1, . . . , K, (2.1)
ai where tai = nj =1 yij is the total time until the arrival of the nai customer considered in si the ith queue, and tsi = nl=1 sil is the total time taken by the server of the ith queue to complete all the services under consideration. Note that the restrictions in the domain of the likelihood in (2.1) correspond to the ergodic condition of each queue. Next, we consider the problem of assigning a prior distribution to the unknown parameters of the model. As mentioned above, we specify the joint prior distribution of (, ) in two stages. In the first, we assign a prior distribution, p(|), for given the hyperparameters. Then, we select a hyperprior distribution p() for those hyperparameters. For simplicity, we select the conditional prior distribution for each pair i = (i , i ) proportional to the product of two independent Gamma distributions (a member of the common natural conjugate prior distribution that would result from an individual analysis of each queue): p(|) =
K
p(i , i |)
i=1
= C()
K
K
Ga(i |, ) Ga(i |, ),
i < b i , i = 1, . . . , K,
(2.2)
i=1
where the hyperparameters are = (, , , ) and the common proportionality constant corresponding to each queue is ( + ) b (b) −1 C() = F 1, + ; + 1; , ()() (b + )(+) b + where F (a, b; c; z) is the Gauss’s hypergeometric function (see Gradshteyn and Ryzhik (1965) for a description of this function). Note that i and i given are not independent due to the restriction i < b i , i = 1, . . . , K. A hyperprior distribution, p(), has to be chosen. In general, the specification of this type of distribution is a very challenging problem in Bayesian statistics (see, for instance, Berger (1985)). In spite of this, it is worth noting that in situations when the prior knowledge of is vague, we can always use non-informative priors (never forgetting that when using improper prior densities, we must always check if the corresponding posterior distribution is proper). Nevertheless, when dealing with real problems, the experts in the subject should normally have enough information (maybe from historical data) in order to construct a nonvague hyperprior distribution, or at least, to constrain the hyperparameters (appropriately reparameterized) into a finite region (see García-Donato (2001) and Gelman et al. (1995) for examples of this situation). An alternative way of estimating the hyperparameters is to use the so-called empirical Bayes (EB) methods, which basically consist of using the experimental information to
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
339
ˆ of the unknown hyperparameters, , and then to proceed in the produce an estimate, , ˆ as the prior distribution of the parametric usual Bayesian way, but working with p(|) vector . Lehoczky (1990), Thiruvaiyaru and Basawa (1992), Sohn (1996), Armero and Conesa (2000) and Sohn (2000) have applied EB methods to the statistical analysis of different queueing systems. There has been a lot written about the advantages and disadvantages of both procedures and so we will not pursue this question here. Nevertheless, as Berger (1985) states, EB methods work fine when the number of queues is high because the i ’s are a random sample from the general population p(|), and so, the bigger the sample, the better the information about . For small samples, however, these methods do not reflect the errors introduced in the estimation of hyperparameters, and only in those particular situations in which all the ˆ the error would possible values of are, on the whole, tightly grouped in their estimation , be small leading to a good performance of the method. In general, EB methods require less computational effort than the full hierarchical Bayesian model. In this case, as the ergodic condition needs to be considered in the general EB procedure, the marginal distribution of the data dependent on results in terms of Gauss hypergeometric functions, and so, new computational problems appear in the estimation of . Robbins (1983), Maritz and Lwin (1989) and Carlin and Louis (2000) are very good references on the subject. Once the elicitation of p() is completed, the joint prior distribution that contains all our prior knowledge of the parameters and hyperparameters is p(, ) = p(|)p() = p()
K
p(i |)
i=1 K
= C p(, , , )
K
Ga(i |, )Ga(i |, ), i < b i , i = 1, . . . , K.
i=1
(2.3) This prior information is then updated, via Bayes’ theorem, with the experimental information contained in the likelihood () in (2.1), to get the joint posterior distribution: p(, |z) ∝ ()p(|)p() = C K p(, , , )
K i=1
ni ai e−i tai ni si e−i tsi Ga(i |, )Ga(i |, ), (2.4)
where i < bi , for i = 1, . . . , K. The distribution in (2.4) contains all our knowledge of the system, but it does not have a known analytical expression. Therefore, in order to provide more intuitive information about and , numerical methods will be needed to evaluate, in practice, the main characteristics of this posterior distribution. Among others, one feasible (indeed the most popular) possibility is to use Markov chain Monte Carlo (MCMC) methods, that draw samples from any intractable posterior by running a cleverly constructed Markov chain over a long period, the stationary distribution of which, being the one we want to simulate from.
340
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
Among the different ways of building these chains, the most popular are Gibbs sampling and Metropolis–Hastings algorithm (see Gilks et al. (1996) for a detailed discussion of these). When, as in our case, the distribution of interest is multivariate, these methods produce a random sample from the joint posterior distribution p(, |z) in (2.4). This means we can use Monte Carlo integration methods to approximate the main features of the marginal posterior distribution of the arrival, service rates and, consequently, of the traffic intensity of each of the K queues. These marginal distributions will provide us with the first intuitive impression of the congestion and performance of the queues we were looking for. MCMC methods have also been successfully used by Rios Insua et al. (1998) in order to deal with the posterior distribution of the parameters in an M/Hk /1 queueing system.
3. Predicting the whole congestion of the system As we have just commented, MCMC methods allow us to extract almost all the information contained in the joint posterior distribution in (2.4) regarding the parameters and hyperparameters (, ), managing this system. Nevertheless, it is well known that, in Queueing Theory, when a system is working in steady state (like the one we are dealing with), its congestion is better described through the so-called measures of performance of the queue (number of customers in the system, waiting times in the system and in the queue, idle and busy periods, etc.). In general, the posterior predictive distribution of a congestion variable can be obtained through its steady-state distribution given the parameters of the model and their marginal posterior distribution. However, in our system, this general procedure gets very complicated in practice: the conditional distribution of these variables is usually expressed is terms of transforms without analytical inversion, and since we have no closed expression for the corresponding marginal posterior distribution of the parameters we can only work with simulated samples from it. In this section, we will discuss some of the practical questions involved when computing the posterior predictive distribution of the measures most often used to describe congestion in bulk service queues, especially those presenting computational problems.All these measures are, of course, interconnected, and their utility in a real problem will depend on its particular characteristics and the objectives of the manufacturers who are analyzing (or paying for) it. 3.1. The state of the system Let (Hi , N q i ) denote the state of the ith queue where Hi = 1 means that the server is busy processing a batch of s items, a s b, and Hi = 0 indicates that the server is idle because the number of items in the system is less than the quorum a. In both cases, N q i is the number of customers waiting in queue. Given the arrival and service rate, its probability distribution is (Medhi, 1991)
r a+1 − rib+1 a + i p((0, 0)|ri ) = 1 − ri (1 − ri )2
−1 ,
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
p((0, n)|ri ) = p((1, n)|ri ) =
1 − rin+1 1 − ri 1 − rib 1 − ri
341
p((0, 0)|ri ),
n = 1, . . . , a − 1
p((0, 0)|ri ) rin+1 ,
n = 0, 1, . . . ,
where, for simplicity, p((hi , ni )|ri ) denotes the probability p((Hi , N q i ) = (hi , ni )|ri ), and ri is the unique real root in (0, 1) of the equation: i = ri + ri2 + · · · + rib . i
(3.1)
Taking into account that the K queues work independently, given the parameters of the system, the joint conditional distribution that describes the state of the whole system, denoted as p((h1 , n1 ), . . . , (hK , nK )|r1 , . . . , rK ) from now on, is the product of its corresponding marginal distributions, in particular, K i=1 p((hi , ni )|ri ). Therefore, the joint posterior predictive distribution of the state of all the queues is p((h1 , n1 ), . . . , (hK , nK )|z) = E (r1 ,...,rK |z) [p((h1 , n1 ), . . . , (hK , nK )|r1 , . . . , rK )]
K = p((hi , ni )|ri )p(r1 , . . . , rK |z)d(r1 , . . . , rK ).
(3.2)
i=1
As we do not have a closed expression for p(r1 , . . . , rK |z) nor for the joint posterior predictive, we can approximate the probabilities of this distribution, by using the following procedure based on MCMC methods: (1) From the random sample of size m of the joint posterior distribution in (2.4) obtained us(j ) (j ) ing MCMC techniques, consider {(1 , . . . , K ); j = 1, . . . , m}, as a random sample from the marginal posterior distribution of the parameters, p(|z). (j ) (2) Generate a sample from p(r1 , . . . , rK |z), where each generated component ri is ob(j ) (j ) (j ) tained by solving Eq. (3.1) for each simulated pair i = (i , i ), ∀i. (3) Approximate each posterior predictive probability as p((h1 , n1 ), . . . , (hK , nK )|z) m 1 (j ) (j ) ≈ p((h1 , n1 ), . . . , (hK , nK )|r1 , . . . , rK ) m j =1
m K 1 (j ) ≈ p((hi , ni )|ri ). m
(3.3)
j =1 i=1
This procedure is similar to the one used in Wiper (1998) to obtain the posterior predictive distribution of the number of customers in Er /M/1 queueing systems. It is worth
342
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
mentioning that the computational effort of this procedure (and the next two presented in this section) depends mainly on two aspects: the MCMC simulation from the marginal posterior distribution of the parameters and the approximation of the posterior predictive distribution. The MCMC simulation is reasonably easy to handle with any programming language (C + +, Java, etc.) or by using specific software such as BUGS䉸 (Gilks et al., 1994). In this procedure, the latter aspect depends only on the difficulty of generating a random sample from the posterior distribution p(r1 , . . . , rK |z), which it turns into a root finding problem that can be implemented in most programming languages without too much difficulty. In Section 4, we discuss our computational experience with all the procedures presented in this section in more detail. From (3.3) it is easy to compute the posterior distribution of the situation of idleness or busyness of all the K machines. Note also that (3.3) does not provide an intuitive visualization of the congestion of each queue. To do so, it would be more suitable to analyze each one of them separately through their corresponding marginal distribution, p((hi , ni )|z), that could be computed directly, marginalizing (3.3), or using a similar procedure to the one presented above, also based on MCMC methods: m 1 (j ) p((hi , ni )|ri ), p((hi , ni )|z) ≈ m
(3.4)
j =1
(j )
where {ri ; j = 1, . . . , m} is a random sample from the marginal distribution p(ri |z). It would be very useful for manufacturing purposes (to learn about the performance of the system, improve the system behavior, or even for planning and building integrated decision projects, for example) to increase our information about the system. Two ways of doing this will be to study the service batch size (a large value of the expected size would indicate that machines always work at their maximum capacity), or/and the waiting time in the queue for an arriving job (of special importance when dealing with pieces that can be damaged by a long wait). Obviously these special features could make the managers think about the possibility of increasing the quorum of the machines, buying new more powerful machines, etc. For both congestion measures, their conditional distribution depends on the root of Eq. (3.1), and so, in order to approximate their corresponding posterior predictive distribution, we could just reproduce procedures similar to the ones presented above. 3.2. Idle periods As the interest of managers in charge of the system is usually to maximize its output (always with the final objective of maximizing benefits), it could be interesting to study how long the idle periods of the machines that make up the whole system are. Due to the fact that the K queues work independently, given all the parameters of the system, the joint probability density function (p.d.f.) of the duration D of an idle period in each queue will be the product of its corresponding marginal distributions: d(t1 , . . . , tK |1 , . . . , K , r1 , . . . , rK ) =
K i=1
di (ti |i , ri ),
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
343
where for the ith queue, ri is as in (3.1) and the conditional distribution of Di is expressed in terms of the Laplace transform (LT) of their corresponding p.d.f. (Medhi, 1991):
∞ di∗ (u|i , ri ) = e−ut i di (ti |i , ri ) dti 0
ai − [ri (u + i )]a 1 i (1 − ri ) = . (1 − ria ) (u + i )a i − ri (u + i )
(3.5)
This LT cannot be analytically inverted, and so numerical inversion becomes necessary. This problem is not new in Queueing Theory. Nowadays, computational analysis plays a fundamental role when analyzing practical congestion (Neuts (1981), Abate and Whitt (1992) and Tijms (1997) are classical references in the subject). We could approximate the posterior predictive distribution of all the idle periods: K ((r1 ,1 ),...,(rK ,K )|z) d(t1 , . . . , tK |z) = E di (ti |ri , i ) (3.6) i=1
using a procedure based, not only on MCMC methods, but also on algorithms for numerically inverting the LTs involved: (1) Use the results in steps 1 and 2 of the general procedure given in Section 3.1 to compute (j ) (j ) (j ) (j ) (3.2), and consider {(1 , r1 ), . . . , (K , rK ), j = 1, . . . , m} as a random sample from p((1 , r1 ), . . . , (K , rK )|z). (j ) (j ) (j ) (j ) (2) For each simulated vector (i , ri ), compute di (ti |i , ri ) by numerically inverting (j ) (j ) its LT, di∗ (u|i , ri ), given in (3.5). (3) Approximate the p.d.f. of the joint posterior predictive distribution as d(t1 , . . . , tK |z) ≈
m K 1 (j ) (j ) di (ti |i , ri ). m
(3.7)
j =1 i=1
As in the previous case, the computational effort of this procedure depends on the MCMC simulation from the marginal posterior distribution of the parameters and the approximation of the marginal conditional distribution di (ti |i , ri ), which now incorporates the additional computational threat of inverting LTs mentioned above. There are many algorithms for numerically inverting LTs and they are not too difficult to implement in any programming language. Abate and Whitt (1992) provide a detailed discussion of most of them and plenty of references on the subject. It is worth mentioning, that although these are very fast algorithms, the biggest computational cost for approximating the marginal posterior predictive distribution comes from the fact that we have to numerically invert an LT for each simulated vector obtained via MCMC. It is clear that the corresponding marginal posterior predictive distribution for the idle period of each machine, more easily computed as di (ti |z) ≈
m 1 (j ) (j ) di (ti |i , ri ), m j =1
(3.8)
344
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354 (j )
(j )
(j )
(j )
where di (ti |i , ri ) is obtained by numerically inverting di∗ (u|i , ri ) as in step 2 above, will provide a convenient visualization of the congestion of each machine. Working in a similar way, we could also analyze the number of items in the queue when an idle period begins. Given all the parameters of the system, the joint distribution of the number of items waiting to be processed by each machine when an idle period begins, p(N d 1 = n1 , . . . , Nd K = nK |r1 , . . . , rK ), with 0 ni < a, i = 1, . . . , K, is the product of their corresponding marginal distributions, p(N d i = ni |ri ), now expressed in terms of their probability generating function (p.g.f.) (Medhi, 1991):
G(Nd i |ri ) (u) =
∞
p(N d i = ni |ri )uni =
ni =0
(1 − ri )(1 − (ri u)a ) . (1 − ria )(1 − ri u)
(3.9)
In order to compute the corresponding joint posterior predictive distribution, we could use a similar procedure to the one presented in this section to deal with idle periods in which, instead of inverting LTs, we would need methods to numerically invert p.g.f.s.
3.3. Busy periods Another complementary way of describing the congestion of our system is to analyze the duration, B, of the busy periods. As usual, its joint conditional (on the parameters) distribution is the product of its corresponding marginal distributions, the p.d.f.s of which are given in terms of their LTs (Medhi, 1991): bi∗ (u|i , Ri ) =
i (1 − Ria ) , u + i − i Ria
(3.10)
where here Ri is the unique real root in (0, 1) of the equation: i Rib+1 − (u + i + i )Ri + i = 0.
(3.11)
In order to approximate the joint posterior predictive distribution of the duration of a busy period for each machine, b(t1 , . . . , tK |z), we could use a similar procedure (with similar computational behavior) to the one presented in Section 3.2 for the idle periods: (1) From the generated sample of the joint posterior distribution p(, ) in (2.4) consider (j ) (j ) {(1 , . . . , K ); j = 1, . . . , m} as a random sample from the marginal posterior distribution p(|z). (j ) (j ) (j ) (2) For each simulated pair i , compute bi (ti |i , Ri ) by numerical inversion of its (j ) (j ) corresponding LT bi∗ (u|i , Ri ), taking into account that Eq. (3.11) must be solved at every evaluation of the LT necessary for its numerical inversion.
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
345
(3) Approximate the p.d.f. of the joint posterior predictive distribution as m
b(t1 , . . . , tK |z) ≈
1 (j ) (j ) (j ) (j ) b(t1 , . . . , tK |1 , . . . , K , R1 , . . . , RK ) m j =1
=
m K 1 (j ) (j ) bi (ti |i , Ri ). m
(3.12)
j =1 i=1
As in the previous measures, in order to have a more convenient visualization of the congestion of each queue, we can consider their corresponding marginal posterior predictive distribution, bi (ti |z), which could be approximated as bi (ti |z) ≈
m 1 (j ) (j ) bi (ti |i , Ri ), m
(3.13)
j =1
(j )
(j )
where each bi (ti |i , Ri ) is obtained by numerically inverting its corresponding LT (j ) (j ) bi∗ (u|i , Ri ) as stated in stage 2. 4. A numerical example with machines and facilities In this final section, with the main objective of illustrating the previously presented results, we analyze the congestion in a manufacturing scenario formed by K = 6 independent machines that work processing items in packets according to an M/M(9, 20)/1 queue model. For each machine i in the system, i = 1, . . . , 6, Table 1 displays the results of a simulated experiment like the one proposed in Section 2. In particular, it contains the interarrival time of nai = 20 consecutive items, the time needed to process nsi = 20 packets of pieces and the maximum likelihood estimate (m.l.e.) ˆ i of the traffic intensity, obtained by maximizing the likelihood function in (2.1) subject to the ergodic condition i < 20i (in our case using MATLAB䉸 ). Note that, although we are considering all the machines to be similar and working at equilibrium, they produce quite different results. In particular, we observe a high level of congestion in the performance of machines 4 and 5 (in both cases, the corresponding Table 1 Sufficient statistics for the arrival and service rate for each machine corresponding to the total time (tai ) to achieve the arrival of 20 consecutive pieces and the processing time (tsi ) for 20 packets of pieces, respectively, and the m.l.e. ( ˆ i ) of the traffic intensity for each machine Machine i
1
2
3
4
5
6
tai tsi ˆ i
13.9 158.8 0.5712
9.4 186.7 0.9931
11.6 135.9 0.5858
8.4 195.5 0.9999
10.4 302.4 0.9999
9.2 169.9 0.9233
346
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354 6 Machine 1 Machine 2 Machine 3 Machine 4 Machine 5 Machine 6
5 4 3 2 1 0
(a)
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4 0.6 Unique machine
0.8
1.0
6 5 4 3 2 1 0
(b)
Fig. 1. Marginal posterior distribution of the traffic intensity of each machine when they are analyzed: (a) separately; (b) together without any distinction among them.
m.l.e.s, ˆ 4 and ˆ 5 , are practically one), while machines 1 and 3 seem to have a very good performance ( ˆ 1 = 0.5712 and ˆ 3 = 0.5858). It is important to point out that if we did not distinguish among the different machines, the m.l.e. of the common would be ˆ = 0.9135, an intermediate level between the lower values for machines 1 and 3 and the higher ones of machines 4 and 5. With the sole idea of providing a reference framework, we start the statistical analysis of the congestion of this system by considering the two most extreme possible situations: each machine is analyzed completely separated, sharing no common information with the others, or, all machines are identical (i = ) with customers arriving at each machine with identical rates (i = ). Firstly, we present the unpooled analysis. Fig. 1(a) shows the marginal posterior distribution of the traffic intensity for each machine. In particular, the posterior distribution corresponding to machine i has been obtained assigning a Jeffrey’s non-informative prior for the corresponding parametric vector i and combining, via Bayes’ theorem, the data collected only on that machine. In the case of analyzing all the machines together,
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
347
Fig. 1(b) displays the posterior distribution of the traffic intensity, computed from the posterior distribution of the common arrival and the service rate, both computed using a Jeffrey’s non-informative prior and the data collected in all the machines (a total of 120 interarrival times and 120 service times). In both graphics the assumption of equilibrium means that the right tail of all the densities does not approach zero when the traffic intensity reaches to the critical value 1 (this special feature is discussed in depth in Armero and Bayarri, 1994). Note also, that the different sample sizes considered when computing both posterior distributions are evident in the smaller variability of the posterior distribution of traffic intensity in the pooled analysis. After these two extreme analyses, we will analyze the congestion of our manufacturing system by using a full Bayesian hierarchical model that allows us to maintain the particular features of each machine (and/or its human operator) and also incorporates those common to all of them. To do so, we assume that the parameters {i = (i , i ); i = 1, . . . , 6} of each machine are a random sample from the common population distribution in (2.2): p(i , i |) ∝ Ga(i |, )Ga(i |, ),
i < 20i ,
(4.1)
where = (, , , ) represents the vector of unknown hyperparameters. We now express our previous knowledge about . At this stage, with the sole objective of better understanding how Bayesian hierarchical models work in queues, we consider two possible situations: one in which our knowledge about could be considered as scarce, being our selection for the hyperprior distribution, p(), the product of four (one for each hyperparameter) exponential distributions with mean 1/0.001, Ex(0.001), and a second, in which we express a greater knowledge of with the product of four independent Ex(0.1). From this point onwards, we will refer to this first situation as LIH (least informative hyperprior), and to the second as MIH (most informative hyperprior). In both cases, we have assumed independence among the hyperparameters, and also chosen proper distributions in order to obtain proper joint posterior distributions. We have selected the exponential density for the hyperprior for all parameters, because, in addition to the uniform and the Gamma distributions, this is very common when working with exponential data. In both situations, the joint posterior distribution of the parameters and hyperparameters, p(, |z), although proper, is analytically intractable. Therefore, we will approximate these distributions by simulation using Markov Chain techniques. One way to do this is via the Gibbs sampler (see, for instance Casella and George (1992), for a detailed description of this method), implemented by Gilks et al. (1994) in the BUGS 䉷 (Bayesian inference Using Gibbs Sampling) program. In particular, we work with the Windows 䉷 version, WinBUGS, developed by Spiegelhalter et al. (1999) (see http://www.mrc-bsu.cam.ac. uk/bugs/ for more information about this project). (0) First we choose a vector of starting values for the parameters (for instance, i =0.17 and (0) i =0.01, i=1, . . . , 6) and the hyperparameters (for instance, (0) =(0) =(0) =(0) =1000 for the LIH case, and (0) = (0) = (0) = (0) = 10 for the MIH case). Then, we decide the duration of the burn-in period (the instant when we think that convergence to the stationary distribution has been reached) to be the first 5000 simulations and to work with the next 10,000 ones. In those situations in which convergence could not be guaranteed, we select every 50th iteration of each chain to contribute to the analysis (usually called a thin value of
348
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
6 Machine 1 Machine 2 Machine 3 Machine 4 Machine 5 Machine 6
5 4 3 2 1 0 0.0
0.2
(a)
0.4 0.6 Product of Ex(0.1)
0.8
1.0
0.4 0.6 Product of Ex(0.001)
0.8
1.0
6 Machine 1 Machine 2 Machine 3 Machine 4 Machine 5 Machine 6
5 4 3 2 1 0 0.0 (b)
0.2
Fig. 2. Estimated marginal posterior distribution of the traffic intensity of each machine in the MIH case (at the top) and in the LIH (at the bottom).
50). The results obtained with WinBUGS are a random sample (each element is a vector with 16 components, 2 × 6 parameters + 4 hyperparameters) from the joint posterior distribution p(, |z). Convergence diagnosis of the resulting Markov chain is a major part of these kind of simulation techniques. To do so, we use CODA, a menu-driven set of S-Plus 䉷 functions which serve as an output processor for the BUGS 䉷 software (Best et al., 1995). Following the authors’recommendation, we have used a combination of diagnostic methods plus visual inspection of the trace plot and summary statistics to asses the convergence of the random samples generated by the WinBUGS software for all the parameters of this model with a reasonable degree of confidence. From the obtained random sample of the joint posterior distribution p(, |z), we get a (j ) random sample { i , j = 1, . . . , 10, 000} from the marginal posterior distribution of the traffic intensity for machine i, p( i |z), i = 1, . . . , 6. Fig. 2 shows a kernel estimation of each one of these marginal distributions in both situations of prior knowledge about
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
349
considered. Note that, in the LIH case, the resulting densities are quite similar, indicating that the hierarchical structure emphasises to the common characteristics of the six machines than the particular differences among them. Nevertheless, when our knowledge about is more accurate (the MIH case), the corresponding densities show greater variety because now the hierarchical model recognizes and distinguishes with greater intensity the individual performance of each machine. However, if we compare the results in Fig. 1, where all the machines were analyzed separately, with the ones in the MIH situation, we observe how the model takes into account the common characteristics of the machines by greatly smoothing the high level of congestion on machine 5 and displacing the values of the traffic intensity on machines 1 and 3 toward more congested regions. We continue with our analysis of the congestion by considering some of the measures of performance studied in Section 3. For simplicity, we only deal here with the marginal posterior distribution of the state of the system, the waiting time in the queue and the duration of the busy periods. All of them are analyzed in the unpooled case and also in the two applications of the Bayesian hierarchical model considered. In these two latter situations, the marginal posterior distributions have been approximated using a random sample of the joint posterior distribution obtained with WinBUGS and also, when needed, using methods for numerically inverting the transforms involved. In particular, we work with the Lattice-Poisson algorithm for inverting probability generating functions and the Euler algorithm to deal with Laplace transforms. Both methods are based on the Fourier series method and have been extensively discussed by Abate and Whitt (1992). Following these authors’ recommendation, we also always use the Gaver–Stehfest (1966,1970) algorithm (based on the Post–Widder (1930,1934) formula and also discussed in depth in Abate and Whitt, 1992) to check the results given by the Euler algorithm. As an exception, because the Laplace transform of the conditioned distribution of the busy period is given in terms of a real root of an equation, we only use the Gaver–Stefesth algorithm (which works with real variable) instead of the Euler algorithm (which works with complex variable) to invert this particular transform. It is worth mentioning that these algorithms consist of evaluating the expression of the transform several times and so, they do not present many difficulties with high values of a and b. In the unpooled analysis of the machines, the calculus of the relevant predictive distribution for each machine has been performed using similar procedures to those in the last section based on Monte Carlo integration. Nevertheless, the expression of the joint posterior distribution of the arrival and service rate usually has a nice form (for instance, as a result of a conjugate analysis) from which it is easy to simulate without the need of MCMC. In relation to the methods used to find roots of the equations involved in these procedures, good results with moderate values of the quorum b motivated us to choose a hybrid algorithm based on a combination of Bisection and Newton–Raphson algorithms (see Press et al. (1986) for a description of this routine). Obviously, this is our particular choice here but other alternatives could have been used. For instance, in order to compute steady-state probabilities of M/G(a, b)/1 queues, Chaudhry et al. (1987) use the Muller’s method, adapted from the version in Conte and de Boor (1972), over any other for its good behavior with big values of b. Fig. 3 displays, the posterior predictive distribution of the state of each machine, for all three situations of prior knowledge considered. The bars in each graphic are the probabilities
350
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
0.06 0.04
0.04
0.02
0.04
0.04
0.1724
0.04
0.04
0 2 4 6 8
0.00 0 2 4 6 8
1 3 5 7 9 11 14
0.0540
0.04
0.04
0 2 4 6 8
0.06
0.00 0 2 4 6 8
1 3 5 7 9 11 14
0.04
0.02
0.02
0.00
0.00 0 2 4 6 8
0.04
0.0785
0.04
0.00
0 2 4 6 8
1 3 5 7 9 11 14
1 3 5 7 9 11 14
0.06 0.0836
0.04
0.0549
0.02
0.00
0.00
0 2 4 6 8
1 3 5 7 9 11 14
0.02
0.02
0.0373
0.02
0.06
0.06
1 3 5 7 9 11 14
0.06 0.0345
0 2 4 6 8
1 3 5 7 9 11 14
0 2 4 6 8
1 3 5 7 9 11 14
0.06 0.0402
0.0505
0.02
0.00
0.00
1 3 5 7 9 11 14
0.06 0.0599
0.02
0.02
0 2 4 6 8
1 3 5 7 9 11 14
0.06
0.06
0.0621
0.02
0.00
0.00
1 3 5 7 9 11 14
0.06 0.1565
0.02
0.02
0 2 4 6 8
1 3 5 7 9 11 14
0.06
0.06
0.0524
0.00 0 2 4 6 8
1 3 5 7 9 11 14
1 3 5 7 9 11 14
0.02
0.00 0 2 4 6 8
0 2 4 6 8
0.06 0.0715
0.02
0.00
0.04
0.00 1 3 5 7 9 11 14
0.06 0.0695
0.03
0.0591
0.02 0 2 4 6 8
1 3 5 7 9 11 14
0.06
0.04
0.04
0.00 0 2 4 6 8
0.04
0.1495
0.02
0.00
0.04
0.06
0.06 0.1806
0.00 0 2 4 6 8
1 3 5 7 9 11 14
0 2 4 6 8
1 3 5 7 9 11 14
Fig. 3. Marginal posterior distribution of the state of each of the six machines working in the system. The bar charts in the column on the left correspond to the separate analysis; in the central column to the MIH Bayesian hierarchical model case; and the column on the right to the LIH Bayesian hierarchical model case. In each graphic, the number in the top left corner is the probability that the machine is idle.
of the number, Nq , of pieces waiting in the queue to be processed when the machine remains idle (probabilities on the left, computed from 0 to 8 pieces, one unit less than the quorum a = 9) and when the machine is busy (probabilities on the right, computed for simplicity until the value Nq = 14 pieces). The probability that the machine is idle is also presented inside the frame in the top left corner. Note that graphics in Fig. 3 show the same general behavior and relationships as those previously observed when dealing with the traffic intensity. In the unpooled analysis (left column), machines 1 and 3 seem to be the least congested (idle about 18.06% and 17.24% of the time, respectively), while machines 4 and 5 are the most congested (with approximate values of idleness of 5.40% and 4.02%, respectively). Obviously, if this was a real problem
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
351
these results would make the people in charge of this facility investigate the particular characteristics of machines 1 and 3 and think over future actions to increase their efficiency. The partial pooling situation, which incorporates the common characteristics among the machines through the hierarchical Bayes modeling, is also presented in Fig. 3: bar charts in the central column correspond to the MIH case, and the ones in the right column correspond to the less informative hyperprior (LIH) distribution. If we compare them with the results of the separate analysis, we observe that the hierarchical model tends to make the performance of the machines more similar. Note that the hierarchical model provides an analysis that jointly integrates the particular peculiarities of each machine (human operator, geographical localization, reliability, ageing, etc.) with all their common characteristics. The statistical analysis is computationally more complex but also more flexible and accurate for the current operation of the system. In the LIH case, this effect is more intense in machines 1 and 3. They are idle for about 5.91% and 6.21% of the available time, quantities that are only slightly bigger than machines 4 and 5, the two most congested ones (with approximate percentages of 3.73% and 3.99%, respectively). In the same way, if we compare the behavior of each machine in the three situations considered, the differences in machines 1 and 3 are spectacular. For instance, the predictive probability that machine 1 is idle is 0.1806 in the separate analysis, 0.1495 in the MIH hierarchical model and 0.0591 in the LIH. But note that we can hardly observe any differences among the rest of the machines. The probability of being idle for the two most congested ones differs by less than 1%. For instance, in the most congested machine 5, in the unpooled analysis the probability of being idle is 0.0402, while in the LIH analysis it is 0.037.3 and 0.0345 in the MIH. For all cases and machines, Fig. 4 displays the marginal posterior distribution (expressed through the cumulative distribution function) of the waiting time in the queue of a piece before being processed. In order to clarify their behavior, we have also included the median values of these distributions. In these graphics, we can observe in general the same type of results and relationships as the ones previously commented on relation to the state of the system. Nevertheless, through basic observation of the graphics corresponding to machines 1, 3 and 5, it seems that the waiting time in the queue, as a measure of performance of a queueing system, is less sensitive to the attempts of the Bayesian hierarchical model (in the LIH when the information about the hyperparameters is scarce) to make the machines uniform than it was in the case of the state of the system previously discussed. Another relevant thing to point out is that the posterior probability of a piece being processed immediately, without waiting in the queue is so small that in all the graphics it seems to be zero. This predictive probability is the posterior probability that the machine is idle and, at the same time, there are eight (the quorum of the machine minus one) pieces waiting in the queue for the arrival of a new piece that allows the machine to process a packet with nine pieces, that is p((0, 8)|z). We finish our predictive analysis of this system with the study of the duration of the busy periods of each machine. Fig. 5 displays the cumulative distribution function of the marginal posterior distribution of this variable in all the situations considered, including in each case the median busy time. We can observe small differences among the six machines, although now machine 5 looks a bit different from the others. It is also worth noting that the separate analysis and the two applications of the Bayesian hierarchical model considered
352
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4
MIH LIH Separately
0.2
Med. = 7.7 Med. = 17.0 Med. = 7.2
0.0
0.4 0.2
Med. = 15.9 Med. = 19.3 Med. = 16.9
20
30 40 Machine 1
50
20
30 40 Machine 2
50
Med. = 18.9 Med. = 20.1 Separately Med. = 21.7 MIH
LIH
0.8
10
60
1.0
1.0
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0 20
30 40 Machine 4
50
60
30 40 Machine 3
50
60
0.8
0.6
10
20
1.0 Med. = 42.0 MIH Med. = 26.2 LIH Separately Med. = 40.5
0.6
0.0
Med. = 6.8 MIH Med. = 16.1 LIH Separately Med. = 6.5
0.0 10
60
0.4 0.2
0.0 10
0.8
MIH LIH Separately
Med. = 13.0 MIH Med. = 18.4 LIH Separately Med. = 14.1
0.0 10
20
30 40 Machine 5
50
60
10
20
30 40 Machine 6
50
60
Fig. 4. Posterior cumulative distribution function of the waiting time in the queue for a piece and its median values for each of the six machines in the separate analysis and in the LIH and MIH applications of the Bayesian hierarchical model considered. 1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4
Med. = 6.9 MIH Med. = 10.6 LIH Separately Med. = 6.6
0.2 0.0
0.8
Med. = 9.6 MIH Med. = 11.2 LIH Separately Med. = 9.6
0.2
20
30 40 Machine 1
50
60
10
1.0
Med. = 10.2 MIH Med. = 11.2 LIH Separately Med. = 10.3
0.8
0.4
20
30 40 Machine 2
50
60
10
0.6
0.4
0.4
0.4
0.2
0.2
0.2
30 40 Machine 4
50
60
50
60
Med. = 8.7 MIH Med. = 10.8 LIH Separately Med. = 8.7
0.0
0.0 20
30 40 Machine 3
0.8
0.6
10
20
1.0 Med. = 15.0 MIH Med. = 12.6 LIH Separately Med. = 15.8
0.6
0.0
Med. = 6.0 MIH Med. = 10.3 LIH Separately Med. = 5.7
0.2 0.0
0.0 10
1.0
0.4
10
20
30 40 50 Machine 5
60
10
20
30 40 Machine 6
50
60
Fig. 5. Marginal posterior distribution of the duration of a busy period and the busy median time for each of the six machines in the separate analysis and in the LIH and MIH applications of the Bayesian hierarchical model considered.
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
353
provide very similar results in the different machines of the manufacturing system. These features could indicate that the duration of a busy period is a robust measure of performance in relation to the use of an unpooled or a partially pooled analysis as explored here using the Bayesian hierarchical model.
Acknowledgements The authors would like to thank two anonymous referees for their interest and careful reading of the paper. Their helpful comments have lead to this improved final version.
References Abate, J., Whitt, W., 1992. The Fourier-series method for inverting transforms of probability distributions. Queueing Systems 10, 5–87. Abellán, J.J., Armero, C., Conesa, D., Pérez-Panadés, J., Zurriaga, O., Martínez-Beneito, M.A., Vanaclocha, H., García-Blasco, M.J., 2003. Analysis of the Renal Transplant Waiting List in the País Valencià (Spain). Technical Report TR06-2003, Departament d’Estadística i Investigació Operativa, Universitat de València. Available online at the following url adress http://matheron.uv.es/investigar/technicals.html. Armero, C., Bayarri, M.J., 1994. Prior assessments for prediction in queues. The Statistician 43, 139–153. Armero, C., Bayarri, M.J., 1996. Bayesian questions and answers in queues. Bernardo, J.M., Berger, J.O., Dawid, A.P., Smith, A.F.M. (Eds.), Bayesian Statistics, vol. 5. Oxford University Press, Oxford, pp. 3–23. Armero, C., Bayarri, M.J., 1999. Dealing with uncertainties in queues and networks of queues: a Bayesian approach. In: Ghosh, S. (Ed.), Multivariate, Design and Sampling. Marcel Dekker, New York, pp. 579–608. Armero, C., Conesa, D., 2000. Prediction in Markovian bulk arrival queues. Queueing Systems 34, 327–350. Armero, C., Abellán, J.J., Conesa, D., Pérez-Panadés, J., Martínez-Beneito, M.A., Zurriaga, O., Vanaclocha, H., García-Blasco, M.J., 2003. Waiting for a kidney in the País Valencià (Spain). Technical Report TR052003, Departament d’Estadística i Inv. Operativa, Universitat de València. Available online at the url adress http://matheron.uv.es/investigar/technicals.html. Berger, J.O., 1985. Statistical Decision Theory and Bayesian Analysis, Springer, New York. Best, N.G., Cowles, M.K., Vines, S.K., 1995. CODA Convergence Diagnosis and Output Analysis Software for Gibbs Sampler Output: Version 0.3. Bhat, U.N., Miller, G.K., Rao, S.S., 1997. Statistical analysis of queueing systems. In: Dshalalow, J.H. (Ed.), Frontiers in Queueing. CRC Press, Boca Raton, FL (Chapter 13). Carlin, B.P., Louis, T.A., 2000. Bayes and Empirical Bayes Methods for Data Analysis, Chapman & Hall/CRC Press, London. Casella, G., George, E.I., 1992. Explaining the Gibbs sampler. Amer. Statist. 46 (3), 167–174. Chaudhry, M.L., Templeton, J.G.C., 1983. A First Course in Bulk Queues, Wiley, New York. Chaudhry, M.L., Madill, B.R., Brière, G., 1987. Computational analysis of steady-state probabilities of M/Ga,b /1 and related nonbulk queues. Queueing Systems 2, 93–114. Conesa, D., 2000. Inferencia y predicción en colas con ingresos o llegadas en grupos. Ph.D. Thesis, Universitat de València (in Spanish). Conte, S.D., de Boor, C., 1972. Elementary Numerical Analysis, McGraw-Hill, New York. Draper, D., 2001. Bayesian Hierarchical Modelling, Preliminary draft available in the url address http://www.bath.ac.uk/m˜asdd. García-Donato, G., 2001. Control de procesos. Alternativas Empírico-Bayes y Bayesianas a los métodos clásicos, Master Degree work, Universitat de València (in Spanish). Gaver Jr., D.P., 1966. Observing stochastic processes and approximate transform inversion. Oper. Res. 14, 444–459. Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 1995. Bayesian Data Analysis, Chapman & Hall, London.
354
C. Armero, D. Conesa / Journal of Statistical Planning and Inference 136 (2006) 335 – 354
Gilks, W.R., Thomas, A., Spiegelhalter, D.J., 1994. A language and program for complex Bayesian modelling. The Statistician 43, 169–178. Gilks, W.R., Richardson, S., Spiegelhalter, D.J., 1996. Markov Chain Monte Carlo in Practice, Chapman & Hall, London. Gradshteyn, I.S., Ryzhik, I.M., 1965. Table of Integrals, Series and Products, Academic Press, Inc., Boston. Griffiths, J.D., 1995. Queueing at the Suez Canal. J. Opl. Res. Soc. 46, 1299–1309. Lehoczky, J., 1990. Statistical methods. In: Heyman, D.P., Sobel, M.J. (Eds.), Handbooks in O.R. and M.S. Elsevier Science Publishers B.V. (North-Holland), Amsterdam. (Chapter 6). Maritz, J.S., Lwin, T., 1989. Empirical Bayes Methods, Chapman & Hall, London. Medhi, J., 1991. Stochastic Models in Queueing Theory, Academic Press, Inc., Boston. Neuts, M.F., 1981. Matrix-Geometric Solutions in Stochastic Models—An Algorithmic Approach, The John Hopkins University Press, Baltimore. Post, E.L., 1930. Generalized differentiation. Trans. Amer. Math. Soc. 32, 723–781. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1986. Numerical Recipes, Cambridge University Press, Cambridge. Rios Insua, D., Wiper, M., Ruggieri, F., 1998. Bayesian analysis of M/Er /1 and M/Hk /1 queues. Queueing Systems 30, 289–308. Robbins, H., 1983. Some thoughts on empirical Bayes estimation. Ann. Stat. 713–723. Spiegelhalter, D.J., Thomas, A., Best, N.G., 1999. Winbugs Version 1.2 User Manual. MRC Biostatistics Unit. Stehfest, H., 1970. Algorithm 368. Numerical inversion of Laplace transforms. Commun. ACM 13, 479–491 (erratum 13, 624). Thiruvaiyaru, D., Basawa, I.V., 1992. Empirical Bayes estimation for queueing systems and networks. Queueing Systems 11, 179–202. Tijms, H.C., 1997. Computational methods in queueing. In: Dshalalow, J.H. (Ed.), Frontiers in Queueing. CRC Press, Boca Raton, FL. (Chapter 12). Widder, D.V., 1934. Numerical inversion of the Laplace integral and the related moment problem. Trans. Amer. Math. Soc. 36, 107–200. Wiper, M., 1998. Bayesian analysis of Er /M/1 and Er /M/c queues. J. Statist. Planning Inference 69, 65–79. Young Sohn, S., 1996. Empirical Bayesian analysis for traffic intensity: M/M/1 queues with covariates. Queueing Systems 22, 383–401. Young Sohn, S., 2000. Robust design of server capability in M/M/1 queues with both partly random arrival and service queues. Comput. Oper. Res. 29, 433–440.