Acessibilidade / Reportar erro

A DISTRIBUTION FOR THE SERVICE MODEL

ABSTRACT

In this paper, we propose a distribution that describes a specific system. The system has a heavy traffic, a fast service and the service rate depends on state of the system. This distribution we call the Maximum-Conway-Maxwell-Poisson-exponential distribution, denoted by MAXCOMPE distribution. The MAXCOMPE distribution is obtained by compound distributions in which we use the zero truncated Conway-Maxwell-Poisson distribution and the exponential distribution. This distribution has adjustment mechanism in order to re-establish the equilibrium of the system when the traffic flow increases and that is described by variations of the pressure parameter. Because of this, the MAXCOMPE distribution contains sub-models, such as, the Maximum-geometric-exponential distribution, the Maximum-Poisson-exponential distribution and the Maximum-Bernoulli-exponential distribution. The properties of the proposed distribution are discussed, including formal proof of its density function and explicit algebraic formulas for their reliability function and moments. The parameter estimation is based on the usual maximum likelihood method. Simulated and real data are shown to illustrate the applicability of the model.

Keywords:
MAXCOMPE distribution; service rate; server

1 INTRODUCTION

In this paper, the studied system is very specific. The system has a single-server, fast service and heavy traffic. Furthermore, the service rate depends on state of the system, that is, when the number of customers increases, the service rate increases. The arrivals of the customers in the system are attached to the service, in other words, when a service finishes, another customer arrives and he enters into service directly. Thus, the service time is the same as inter-arrival time. For this reason, this system has an adjustment mechanism in order to re-establish the equilibrium of the system and avoid congestion. The increase in service rate and the opening the new service channels are the adjustment mechanisms. Some examples of systems with this behavior are nonstop toll electronic, the traffic flow and the access on website.

We propose the Maximum-Conway-Maxwell-Poisson-exponential distribution, denoted byMAXCOMPE distribution to describe this system. Moreover, we are particularly interested on the maximum inter-arrival time or maximum service time. The methodology used to obtain this distribution is the compound distributions. We compose the zero truncated Conway-Maxwell-Poisson distribution and the exponential distribution.

The MAXCOMPE distribution contains sub-models that describe the variations of the system, such as, the Maximum-geometric-exponential distribution, denoted by MAXGE distribution, the Maximum-Poisson-exponential distribution, denoted by MAXPE distribution and Maximum-Bernoulli-exponential distribution, denoted by MAXBE distribution.

The compound distributions has been extensively studied. We mention the following articles: (Barreto-Souza et al., 20113 BARRETO-SOUZA W, DE MORAIS AL & CORDEIRO GM. 2011. The weibull-geometric distribution. Journal of Statistical Computation and Simulation, 81(5): 645-657.) introduced the Weibull-geometric distribution which generalizes the exponential-geometric distribution proposed by (Adamidis et al., 20051 ADAMIDIS K, DIMITRIKAPOULO T & LOUKAS S. 2005. On an extension of the exponential-geometric distribution. Statistics & probability letters, 73(3): 259-269.) and (Cancho et al., 20118 CANCHO VG, LOUZADA F & BARRIGA GD. 2011. The geometric-birnbaum-saunders regression model with cure rate. Journal of Statistical Planning and Inference.) proposed the Poisson-exponential distribution with two-parameters, (Cordeiro et al., 20129 CORDEIRO GM, RODRIGUES J & DE CASTRO M. 2012. The exponential com-poisson distribution. Statistical Papers, pages 1-12.) introduced a new three parameters distribution called exponential-Conway-Maxwell Poisson distribution, denoted by ECOMP distribution.

We are interested in systems which the service rate depends on state of the system. We would like to mention some articles with that characteristic. (Kimura, 199115 KIMURA T. 1991. Approximating the mean waiting time in the gi/g/s queue. Journal of the Operational Research Society, pages 959-970.) provided two distributions with an approximation for the mean waiting time in a GI /G /s queue. This approximation are weighted combinations of the exact mean waiting time for the GI /M /s and M /D /S queues. (Sharma & Sharma, 199424 SHARMA K & SHARMA G. 1994. A delay dependent queue without pre-emption with general linearly increasing priority function. Journal of the Operational Research Society, pages 948-953.) studied the delay dependent priority queueing discipline with non-preemptive priority classes and it obtained derived for the expected waiting time of a customer from the priority class. (Morabito & de Lima, 200021 MORABITO R & DE LIMA FC. 2000. Um modelo para analisar o problema de filas em caixas de supermercados: um estudo de caso. Pesquisa Operacional, 20(1): 59-71.) applied queueing theory to analyze the problem of congestion in supermarket check-outs. (Whitt, 200327 WHITT W. 2003. How multiserver queues scale with growing congestion-dependent demand. Operations Research, 51(4): 531-542.) studied how performance scales in the standard M /M /n queue in the presence of growing congestion-dependent customer demand and scaled the queue by letting the potential arrival rate be proportional to the number of servers. (Jongbloed & Koole, 200113 JONGBLOED G & KOOLE G. 2001. Managing uncertainty in call centres using poisson mixtures. Applied Stochastic Models in Business and Industry, 17(4): 307-318.) studied a call center as a queueing model with Poisson arrivals having an unknown varying arrival rate. (Bekker et al., 20045 BEKKER R, BORST S, BOXMA O & KELLA O. 2004. Queues with workload-dependent arrival and service rates. Queueing Systems, 46(3): 537-556.) considered the types of queues M /G /1 and G /G /1 using the workload-dependent arrival and a service speed. Some studies can be found of a single server queue where the service time is dependent linearly and randomly on the waiting time such as in (Boxma & Vlasiou, 20077 BOXMA O & VLASIOU M. 2007. On queues with service and interarrival times depending on waiting times. Queueing Systems, 56(3): 121-132.) and (Whitt, 199026 WHITT W. 1990. Queues with service times and interarrival times depending linearly and randomly upon waiting times. Queueing Systems, 6(1): 335-351.). (Iyer & Manjunath, 200612 IYER SK & MANJUNATH D. 2006. Queues with dependency between interarrival and service times using mixtures of bivariates. Stochastic models, 22(1): 3-20.) analyzed queueing models where the joint density of the inter-arrival time and the service time were described by a mixture of joint densities. (Droguett & Mosleh, 200710 DROGUETT EL & MOSLEH A. 2007. Time to failure assessment of products at service conditions from accelerated lifetime tests with stress-dependent spread in life. Pesquisa Operacional, 27(2): 209-233.) presented a Bayesian framework for the evaluate of product time to failure at service from ALT with stress dependent spread where the time to failure follows a Weibull distribution. (Marin et al., 200718 MARIN CV, DRURY CG, BATTA R & LIN L. 2007. Server adaptation in an airport security system queue. OR Insight, 20(4): 22-31.) studied the airport security queueing system for server behavior dependent to the queue length and implication on security. (Barth et al., 20104 BARTH W, MANITZ M & RAIK S. 2010. Analysis of two-level support systems with time-dependent overflow a banking application. Production and Operations Management, 19(6): 757-768.) evaluated the performance of call centers service with two levels of support and time-dependent mechanism. (Bekker et al., 20116 BEKKER RGK, NIELSEN B & BANG T. 2011. Queues with waiting time dependent service. Queueing Systems, 68(1): 61-78.) studied asteady-state waiting time distribution where the service depends on the waiting time of the first customer in line.

This paper is organized as follows. In Section 2, the MAXCOMPE distribution and some of its properties are presented. The expressions for the density function and r -th moment of theMAXCOMPE distribution are described. In Section 3, we discuss maximum likelihood estimation. In Section 4, some numerical results with simulation and real data are showed. Conclusions are presented in Section 5.

2 THE MAXCOMPE DISTRIBUTION

In this section, we present the MAXCOMPE distribution. Consider the arrivals of the customers occur at random, inter-arrival times are exponentially distributed with mean λ, the service time is exponentially distributed and its mean is dependent on the system state, given by μm = mΦμ, where Φ is pressure parameter which indicates the degree that the service rate is affected by the state of the system, m is the number of customers in the system and μ is the service rate. The service discipline is considered to be first-in-first-out (FIFO).

When the pressure parameter assumes the values zero, Φ =0, the system is not accelerated and an adjustment is not required, for Φ =1 the service rate is proportional to the state of the queuing system and the service channels are open, for Φ → ∞ the system is more accelerated and the service is very fast.

In Figure 1 we show the scheme of the system with the three variations of the pressure parameter.

Figure 1
The scheme of the System.

Let M be a random variable representing the number of customers in the system with a zero truncated COMP distribution (Cordeiro et al., 20129 CORDEIRO GM, RODRIGUES J & DE CASTRO M. 2012. The exponential com-poisson distribution. Statistical Papers, pages 1-12.). Its probability mass function is expressed as

where is a normalizing constant, ρ >0 is the traffic intensity and Ф ∈ [ ∞, ∞] is the pressure parameter.

Furthermore, Ф ∈ [- ∞,1) the work is not accelerated and the service center. For Ф ∈ (1, ∞] the work in the center is accelerated and some adjustments are needed. For Ф =0, we have the geometric distribution. For Ф =1, we have the Poisson distribution and for Φ → ∞ we have the Bernoulli distribution.

Let Yi , i =1,2,..., be independent random variables denoting the inter-arrival times. The random variable Yi is exponentially distributed with parameter λ >0 and its density function is given by

Consider the random variable Y = max(Y1 , ..., YM ) that represents the maximum inter-arrival time or the maximum service time.

In the following proposition, we provide the distribution that describes the studied system.

Proposition 1 Consider M, Y1, Y2, ... independent random variables defined above. The density function of the MAXCOMPE distribution is given by

where θ = (ρ, λ, Φ)T is the vector of parameters, λ >0 is the arrival rate, ρ <1 is the traffic intensity and Ф ∈ ( - ∞,) is the pressure parameter.

Proof.

• The conditional density function of the Y given M=m is provided by

where is the density function of the Y1 and is the probability function of the Y1 . Therefore, the density function of the MAXCOMPE distribution is given by

This completes the proof.

The distribution function of the MAXCOMPE distribution is given by

where .

The reliability function of the MAXCOMPE distribution is

The r-th moment of Y is given by

The Equation (6) can be calculated numerically by truncation (Gradshtein et al., 200711 GRADSHTEIN IS, RYZHIK IM, JEFFREY A & ZWILLINGER D. 2007. Table of integrals, series, and products. Academic press.). The truncation of this sum is described in Appendix APPENDIX In this appendix, we show the programs in software R used in this paper and the truncation of the sum in E (Y ). 1. The programs in software R used in this paper. • Programs in software R to simulate the M /M /1 model, Test Log-rank and Chi-squared tests. require(VGAM) rejlog=0 rejchi=0 for (iin 1:1000) n =10000 stop nat =1 number of the channels λ =0.9 the arrival rate mi =1 the service rate v=vector() customers.arrival=0 customes.service=0 customers.queue=0 system=numeric(n) wait=numeric(n) time = numeric(n) time-wait=numeric(n) clock1=numeric(nat) y=vector() time-service=numeric(n) alarm-clock=rep(-1,nat) for(i;in 2:n) { time.atual = (i -1)*0.1 time[i] = (i -1)*0.1 arrival = rpois (1,0.1) customers-arrival= customers-arrival + arrival customers-queue = customers-queue + arrival system[i] = system[i-1] + arrival if (sum (alarm -clock ==-1)>0 and customers -queue >0) { vagas = sum (alarm -clock ==-1) if (place !=sum (alarm -clock ==-1)) } stop("non conforming supositions! ") } minimum = min(place,customers-queue) for (j in 1:minimum) { clock =rexp (1,0.9) clock 1[which (alarm -clock ==-1)[1]]=clock alarm -clock [which (alarm -clock ==-1)[1]]=time .atual +clock customers -queue =customers -queue -1} } for (j in 1:nat ) { if (time .atual >=alarm -clock [j ] e alarm -clock [j ]!=-1) { system [i ]=system [i ]-1 alarm -clock [j ]=-1 customers -service =customers -service +1 y [customers -service ]=clock 1[j ] clock 1[j ]=0 } } Log-rank and Chi-squared tests time =c (y ,s 2) group =c (rep (0,length (y )),rep (1,length (s 2))) status =rep (1,length (y )+length (s 2)) mice =data .frame (time ,status ,group ) s =survdiff (Surv (time ,status ==1) group ,data =mice ) schi =chisq .test (y ,s 2) if (schisq >3.84) { rejlog =rejlog +1 } if (schip .value <0.05) { rejchi <-rejchi +1 } } • Program reliability function and estimated reliability function by Kaplan-Meiermethod. dado <-cbind (y ,cens =1) s 2=numeric () mu =0.9 traffic intensity sigma =0.9 service rate s 2=exp (mu *(1-exp (-sigma *y )))*exp (-sigma *y )/(exp (mu )-1) MAXPE distribution s 2=exp (-y *sigma )*(1-mu )*sigma /(1-mu *(1-exp (-y *sigma ))) MAXGE distribution require (survival ) tempos =dado [,1] cens =dado [,2] ekm =survfit (Surv (y ,cens ) 1) plot (ekm ,lty =1,lwd =2,conf .int =F ,xlab ="time ", ylab ="reliabilityfunction ") lines (y ,s 2,col ='blue ',lty =1,lwd =2) text (7,1,"trafficintensity =0.8") legend (3,0.8,c ('M /M /1','W (y )=0.876'),bty ='1',lty =c (1,4,3,2), lwd =c (2,2,2,2)) • Program optim for the MLEs. verosi <-function (param ) { mu <-param [1] sigma <-param [2] if (any(mu <1e - 200))return(.Machinedouble.xmax 5) if (any(sigma <1e - 200))return(.Machinedouble.xmax 5) f = sigma * (exp(-sigma * y)) * (1 - mu)/(1 - mu* (1 - exp(-sigma * y)))2 MAXGE distribution lv <-log (f ) sum (-lv ) } estima <-optim (c (0.888,0.9),verosi ,method ="BFGS ", hessian=T); estima estimahessian solve (estimahessian ) Criterions AIC and BIC lm 1<-lm (smax 1) AIC (lm 1) stopifnot (all .equal (AIC (lm 1), AIC (logLik (lm 1)))) BIC (lm 1) lm 2<-update (lm 1,...-Examination ) AIC (lm 1,lm 2) BIC (lm 1,lm 2) 2. The truncation of the sum in for all the models: The is We have an upper bound as The series τ has an upper bound given by the series that is convergent (Minka et al., 2003) because Therefore, the series τ is convergent. Furthermore, it is possible to truncate the series τ at some υ-th term such that where where R v is absolute truncated error. An upper bound can be found based in the series ρm /(m !)Φ, m =1,2,...,, decreases at a faster rate then a geometric series. Thus, there exists 0 < εv < 1 for all m > υ, so that ρm /(m + 1)Φ < εv Therefore, R v < εv , the error is small. .

The MAXCOMPE distribution has sub-models that are presented via the corollaries below.

Corollary 1 When the pressure parameter Ф assumes the value zero, the MAXCOMPE distribution becomes the Maximum-geometric-exponential distribution, denoted by MAXGE distribution. The density function of the MAXGE distribution is

where λ >0 and ρ <1.

In this case, the service rate is independent on state of the system.

The reliability function is given by

The r-th moment of Y is given by

The sum in Equation (9) is truncated the same way that in the Equation (6) (see Appendix APPENDIX In this appendix, we show the programs in software R used in this paper and the truncation of the sum in E (Y ). 1. The programs in software R used in this paper. • Programs in software R to simulate the M /M /1 model, Test Log-rank and Chi-squared tests. require(VGAM) rejlog=0 rejchi=0 for (iin 1:1000) n =10000 stop nat =1 number of the channels λ =0.9 the arrival rate mi =1 the service rate v=vector() customers.arrival=0 customes.service=0 customers.queue=0 system=numeric(n) wait=numeric(n) time = numeric(n) time-wait=numeric(n) clock1=numeric(nat) y=vector() time-service=numeric(n) alarm-clock=rep(-1,nat) for(i;in 2:n) { time.atual = (i -1)*0.1 time[i] = (i -1)*0.1 arrival = rpois (1,0.1) customers-arrival= customers-arrival + arrival customers-queue = customers-queue + arrival system[i] = system[i-1] + arrival if (sum (alarm -clock ==-1)>0 and customers -queue >0) { vagas = sum (alarm -clock ==-1) if (place !=sum (alarm -clock ==-1)) } stop("non conforming supositions! ") } minimum = min(place,customers-queue) for (j in 1:minimum) { clock =rexp (1,0.9) clock 1[which (alarm -clock ==-1)[1]]=clock alarm -clock [which (alarm -clock ==-1)[1]]=time .atual +clock customers -queue =customers -queue -1} } for (j in 1:nat ) { if (time .atual >=alarm -clock [j ] e alarm -clock [j ]!=-1) { system [i ]=system [i ]-1 alarm -clock [j ]=-1 customers -service =customers -service +1 y [customers -service ]=clock 1[j ] clock 1[j ]=0 } } Log-rank and Chi-squared tests time =c (y ,s 2) group =c (rep (0,length (y )),rep (1,length (s 2))) status =rep (1,length (y )+length (s 2)) mice =data .frame (time ,status ,group ) s =survdiff (Surv (time ,status ==1) group ,data =mice ) schi =chisq .test (y ,s 2) if (schisq >3.84) { rejlog =rejlog +1 } if (schip .value <0.05) { rejchi <-rejchi +1 } } • Program reliability function and estimated reliability function by Kaplan-Meiermethod. dado <-cbind (y ,cens =1) s 2=numeric () mu =0.9 traffic intensity sigma =0.9 service rate s 2=exp (mu *(1-exp (-sigma *y )))*exp (-sigma *y )/(exp (mu )-1) MAXPE distribution s 2=exp (-y *sigma )*(1-mu )*sigma /(1-mu *(1-exp (-y *sigma ))) MAXGE distribution require (survival ) tempos =dado [,1] cens =dado [,2] ekm =survfit (Surv (y ,cens ) 1) plot (ekm ,lty =1,lwd =2,conf .int =F ,xlab ="time ", ylab ="reliabilityfunction ") lines (y ,s 2,col ='blue ',lty =1,lwd =2) text (7,1,"trafficintensity =0.8") legend (3,0.8,c ('M /M /1','W (y )=0.876'),bty ='1',lty =c (1,4,3,2), lwd =c (2,2,2,2)) • Program optim for the MLEs. verosi <-function (param ) { mu <-param [1] sigma <-param [2] if (any(mu <1e - 200))return(.Machinedouble.xmax 5) if (any(sigma <1e - 200))return(.Machinedouble.xmax 5) f = sigma * (exp(-sigma * y)) * (1 - mu)/(1 - mu* (1 - exp(-sigma * y)))2 MAXGE distribution lv <-log (f ) sum (-lv ) } estima <-optim (c (0.888,0.9),verosi ,method ="BFGS ", hessian=T); estima estimahessian solve (estimahessian ) Criterions AIC and BIC lm 1<-lm (smax 1) AIC (lm 1) stopifnot (all .equal (AIC (lm 1), AIC (logLik (lm 1)))) BIC (lm 1) lm 2<-update (lm 1,...-Examination ) AIC (lm 1,lm 2) BIC (lm 1,lm 2) 2. The truncation of the sum in for all the models: The is We have an upper bound as The series τ has an upper bound given by the series that is convergent (Minka et al., 2003) because Therefore, the series τ is convergent. Furthermore, it is possible to truncate the series τ at some υ-th term such that where where R v is absolute truncated error. An upper bound can be found based in the series ρm /(m !)Φ, m =1,2,...,, decreases at a faster rate then a geometric series. Thus, there exists 0 < εv < 1 for all m > υ, so that ρm /(m + 1)Φ < εv Therefore, R v < εv , the error is small. ).

Corollary 2 When the pressure parameter Ф assumes the value one, Ф =1, the MAXCOMPE distribution becomes the Maximum-Poisson-exponential distribution, denoted by MAXPE distribution. The density function of the MAXPE distribution is given by

where λ >0 and ρ >0, in this case there is not restriction for ρ.

The service rate is directly proportional to the state of the system. Furthermore, new service channels are opened.

The reliability function is given by

The r-th moment of Y is given by

The sum in Equation (12) is truncated the same way that in the Equation (6) (see Appendix APPENDIX In this appendix, we show the programs in software R used in this paper and the truncation of the sum in E (Y ). 1. The programs in software R used in this paper. • Programs in software R to simulate the M /M /1 model, Test Log-rank and Chi-squared tests. require(VGAM) rejlog=0 rejchi=0 for (iin 1:1000) n =10000 stop nat =1 number of the channels λ =0.9 the arrival rate mi =1 the service rate v=vector() customers.arrival=0 customes.service=0 customers.queue=0 system=numeric(n) wait=numeric(n) time = numeric(n) time-wait=numeric(n) clock1=numeric(nat) y=vector() time-service=numeric(n) alarm-clock=rep(-1,nat) for(i;in 2:n) { time.atual = (i -1)*0.1 time[i] = (i -1)*0.1 arrival = rpois (1,0.1) customers-arrival= customers-arrival + arrival customers-queue = customers-queue + arrival system[i] = system[i-1] + arrival if (sum (alarm -clock ==-1)>0 and customers -queue >0) { vagas = sum (alarm -clock ==-1) if (place !=sum (alarm -clock ==-1)) } stop("non conforming supositions! ") } minimum = min(place,customers-queue) for (j in 1:minimum) { clock =rexp (1,0.9) clock 1[which (alarm -clock ==-1)[1]]=clock alarm -clock [which (alarm -clock ==-1)[1]]=time .atual +clock customers -queue =customers -queue -1} } for (j in 1:nat ) { if (time .atual >=alarm -clock [j ] e alarm -clock [j ]!=-1) { system [i ]=system [i ]-1 alarm -clock [j ]=-1 customers -service =customers -service +1 y [customers -service ]=clock 1[j ] clock 1[j ]=0 } } Log-rank and Chi-squared tests time =c (y ,s 2) group =c (rep (0,length (y )),rep (1,length (s 2))) status =rep (1,length (y )+length (s 2)) mice =data .frame (time ,status ,group ) s =survdiff (Surv (time ,status ==1) group ,data =mice ) schi =chisq .test (y ,s 2) if (schisq >3.84) { rejlog =rejlog +1 } if (schip .value <0.05) { rejchi <-rejchi +1 } } • Program reliability function and estimated reliability function by Kaplan-Meiermethod. dado <-cbind (y ,cens =1) s 2=numeric () mu =0.9 traffic intensity sigma =0.9 service rate s 2=exp (mu *(1-exp (-sigma *y )))*exp (-sigma *y )/(exp (mu )-1) MAXPE distribution s 2=exp (-y *sigma )*(1-mu )*sigma /(1-mu *(1-exp (-y *sigma ))) MAXGE distribution require (survival ) tempos =dado [,1] cens =dado [,2] ekm =survfit (Surv (y ,cens ) 1) plot (ekm ,lty =1,lwd =2,conf .int =F ,xlab ="time ", ylab ="reliabilityfunction ") lines (y ,s 2,col ='blue ',lty =1,lwd =2) text (7,1,"trafficintensity =0.8") legend (3,0.8,c ('M /M /1','W (y )=0.876'),bty ='1',lty =c (1,4,3,2), lwd =c (2,2,2,2)) • Program optim for the MLEs. verosi <-function (param ) { mu <-param [1] sigma <-param [2] if (any(mu <1e - 200))return(.Machinedouble.xmax 5) if (any(sigma <1e - 200))return(.Machinedouble.xmax 5) f = sigma * (exp(-sigma * y)) * (1 - mu)/(1 - mu* (1 - exp(-sigma * y)))2 MAXGE distribution lv <-log (f ) sum (-lv ) } estima <-optim (c (0.888,0.9),verosi ,method ="BFGS ", hessian=T); estima estimahessian solve (estimahessian ) Criterions AIC and BIC lm 1<-lm (smax 1) AIC (lm 1) stopifnot (all .equal (AIC (lm 1), AIC (logLik (lm 1)))) BIC (lm 1) lm 2<-update (lm 1,...-Examination ) AIC (lm 1,lm 2) BIC (lm 1,lm 2) 2. The truncation of the sum in for all the models: The is We have an upper bound as The series τ has an upper bound given by the series that is convergent (Minka et al., 2003) because Therefore, the series τ is convergent. Furthermore, it is possible to truncate the series τ at some υ-th term such that where where R v is absolute truncated error. An upper bound can be found based in the series ρm /(m !)Φ, m =1,2,...,, decreases at a faster rate then a geometric series. Thus, there exists 0 < εv < 1 for all m > υ, so that ρm /(m + 1)Φ < εv Therefore, R v < εv , the error is small. )

Corollary 3 When Φ → ∞, the MAXCOMPE distribution becomes the Maximum-Bernoulli-exponential distribution, denoted by MAXBE distribution. The density function of the MAXBE is given by

where λ >0.

The reliability function is given by

The r-th moment of Y is given by

3 MAXIMUM LIKELIHOOD ESTIMATION

Let y= (y1, ..., yn) be a random sample of Y which it has the MAXCOMPE distribution with unknown parameter vector Θ = (ρ, λ, Ф)T. The log likelihood function for Θ is

The score function U(Θ) = (Uρ, Uλ, UФ)T has components

The maximum likelihood estimates (MLEs) of Θ are obtained by solving the Equations U (Θ)=0. We use the numerical maximization by optim package in software R (Rigby & Stasinopoulos, 200522 RIGBY R & STASINOPOULOS D. 2005. Generalized additive models for location, scale and shape. Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(3): 507-554.).

Under some regularity conditions, has asymptotic multivariate normal distribution with mean Θ and variance I-1(Θ)),

where I (Θ) is Fisher information matrix. Furthermore,

is the observed Fisher information matrix. This result can be used to construct confidence intervals for the parameters of the distribution.

The log likelihood function for the sub-models are shown below.

The log likelihood function of the MAXGE distribution is given by

The components of the score function are given by

The log likelihood function of the MAXPE distribution is

The components of the score function are given by

The log likelihood function of the MAXBE distribution is

The score function is

4 NUMERICAL RESULTS

In this section, we present simulated data using a program made in the software R (see Appendix APPENDIX In this appendix, we show the programs in software R used in this paper and the truncation of the sum in E (Y ). 1. The programs in software R used in this paper. • Programs in software R to simulate the M /M /1 model, Test Log-rank and Chi-squared tests. require(VGAM) rejlog=0 rejchi=0 for (iin 1:1000) n =10000 stop nat =1 number of the channels λ =0.9 the arrival rate mi =1 the service rate v=vector() customers.arrival=0 customes.service=0 customers.queue=0 system=numeric(n) wait=numeric(n) time = numeric(n) time-wait=numeric(n) clock1=numeric(nat) y=vector() time-service=numeric(n) alarm-clock=rep(-1,nat) for(i;in 2:n) { time.atual = (i -1)*0.1 time[i] = (i -1)*0.1 arrival = rpois (1,0.1) customers-arrival= customers-arrival + arrival customers-queue = customers-queue + arrival system[i] = system[i-1] + arrival if (sum (alarm -clock ==-1)>0 and customers -queue >0) { vagas = sum (alarm -clock ==-1) if (place !=sum (alarm -clock ==-1)) } stop("non conforming supositions! ") } minimum = min(place,customers-queue) for (j in 1:minimum) { clock =rexp (1,0.9) clock 1[which (alarm -clock ==-1)[1]]=clock alarm -clock [which (alarm -clock ==-1)[1]]=time .atual +clock customers -queue =customers -queue -1} } for (j in 1:nat ) { if (time .atual >=alarm -clock [j ] e alarm -clock [j ]!=-1) { system [i ]=system [i ]-1 alarm -clock [j ]=-1 customers -service =customers -service +1 y [customers -service ]=clock 1[j ] clock 1[j ]=0 } } Log-rank and Chi-squared tests time =c (y ,s 2) group =c (rep (0,length (y )),rep (1,length (s 2))) status =rep (1,length (y )+length (s 2)) mice =data .frame (time ,status ,group ) s =survdiff (Surv (time ,status ==1) group ,data =mice ) schi =chisq .test (y ,s 2) if (schisq >3.84) { rejlog =rejlog +1 } if (schip .value <0.05) { rejchi <-rejchi +1 } } • Program reliability function and estimated reliability function by Kaplan-Meiermethod. dado <-cbind (y ,cens =1) s 2=numeric () mu =0.9 traffic intensity sigma =0.9 service rate s 2=exp (mu *(1-exp (-sigma *y )))*exp (-sigma *y )/(exp (mu )-1) MAXPE distribution s 2=exp (-y *sigma )*(1-mu )*sigma /(1-mu *(1-exp (-y *sigma ))) MAXGE distribution require (survival ) tempos =dado [,1] cens =dado [,2] ekm =survfit (Surv (y ,cens ) 1) plot (ekm ,lty =1,lwd =2,conf .int =F ,xlab ="time ", ylab ="reliabilityfunction ") lines (y ,s 2,col ='blue ',lty =1,lwd =2) text (7,1,"trafficintensity =0.8") legend (3,0.8,c ('M /M /1','W (y )=0.876'),bty ='1',lty =c (1,4,3,2), lwd =c (2,2,2,2)) • Program optim for the MLEs. verosi <-function (param ) { mu <-param [1] sigma <-param [2] if (any(mu <1e - 200))return(.Machinedouble.xmax 5) if (any(sigma <1e - 200))return(.Machinedouble.xmax 5) f = sigma * (exp(-sigma * y)) * (1 - mu)/(1 - mu* (1 - exp(-sigma * y)))2 MAXGE distribution lv <-log (f ) sum (-lv ) } estima <-optim (c (0.888,0.9),verosi ,method ="BFGS ", hessian=T); estima estimahessian solve (estimahessian ) Criterions AIC and BIC lm 1<-lm (smax 1) AIC (lm 1) stopifnot (all .equal (AIC (lm 1), AIC (logLik (lm 1)))) BIC (lm 1) lm 2<-update (lm 1,...-Examination ) AIC (lm 1,lm 2) BIC (lm 1,lm 2) 2. The truncation of the sum in for all the models: The is We have an upper bound as The series τ has an upper bound given by the series that is convergent (Minka et al., 2003) because Therefore, the series τ is convergent. Furthermore, it is possible to truncate the series τ at some υ-th term such that where where R v is absolute truncated error. An upper bound can be found based in the series ρm /(m !)Φ, m =1,2,...,, decreases at a faster rate then a geometric series. Thus, there exists 0 < εv < 1 for all m > υ, so that ρm /(m + 1)Φ < εv Therefore, R v < εv , the error is small. ) and two real data.

4.1 Simulation

The studied system was simulated by M /M /1 model. This model can be applied to provide approximate of other models (Whitt, 198925 WHITT W. 1989. Planning queueing simulations. Management Science, 35(11): 1341-1366.). We simulated two M /M /1 models with two pairs of values (ρ, λ, Ф): (0.8,0,8,0) and (0.9,0.9,0). These parameters were chosen because the system has heavy traffic, thus, with ρ =0.8 and 0.9 the server has 80% and 90% of observed time occupied. We establish 10,000 arrivals in the system as the ending point of the simulation. The sample size in each simulation is 9.000 service times. We choose the simulations nearest of the studied system. In other words, the system with the queue that has the shortest length.

Therefore, the simulated models have single server, heavy traffic and the service rate independent on state of the system. Due to the fact, the MAXGE distribution was chosen to model the maximum service time.

In the Tables 1 and 2 we show the summary of the simulations chosen to represent the system. We observe that the standard-deviation increases when the traffic intensity increases.

Table 1
Summary of the simulated data with ρ =0.8.

Table 2
Summary of the simulated data with ρ =0.9.

The parameters of the model are estimated by maximum likelihood method. Numerical maximization of the log likelihood function is accomplished by using the RS method (Rigby & Stasinopoulos, 200522 RIGBY R & STASINOPOULOS D. 2005. Generalized additive models for location, scale and shape. Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(3): 507-554.). In Table 3 we show the MLEs of the MAXGE distribution, its to the respective confidence interval (CI) and standard error (sde).

Table 3
The MLEs, the confidence intervals (CI) 95% and sde for the parameters of the MAXGE distribution.

We simulate 1,000 samples sizes of the service times to evaluate the performance of the MLEs of the model for each used value ρ and λ. We consider different the sample sizes for observe the behavior of the estimates of the parameter, such as, n =50, n =250, n =500, n =700 and 1,000. In Tables 4 and 5 we show sample means, variances and estimated bias of the MLEs. We note that the bias and the variances of the MLEs decreases as n increases.

Table 4
Maximum likelihood simulation study of the MAXGE distribution for ρ =0.8 and λ =0.8.

Table 5
Maximum likelihood simulation study of the MAXGE distribution for ρ =0.9 and λ =0.9.

We can study the goodness of fit of the proposed model by the comparison between the observed and predicted values. The simplest way to make this comparison is graphically, which consists to compare the reliability function of the MAXGE distribution with the estimated reliability function by Kaplan-Meier method (Kaplan & Meier, 195814 KAPLAN EL & MEIER P. 1958. Nonparametric estimation from incomplete observations. Journal of the American statistical association, 53(282): 457-481.).

In the Figure 2 we compare the reliability function of the MAXGE distribution with the estimated reliability function. We observe that the reliability function of the MAXGE distribution has a close concordance with the estimated reliability function.

Figure 2
Comparison between the reliability function of the MAXGE distribution and the estimated reliability function with traffic intensity ρ =0.8 and ρ =0.9.

We use Log-rank and Chi-squared tests to evaluate if a closeness of the curves is significant or not (Kleinbaum & Klein, 201216 KLEINBAUM DG & KLEIN M. 2012. Kaplan-meier survival curves and the log-rank test. In Survival analysis, pages 55-96. Springer. and Massey Jr, 195119 MASSEY JR FJ. 1951. The kolmogorov-smirnov test for goodness of fit. Journal of the American statistical Association, 46(253): 68-78.). The null hypothesis of the tests cited is that there is no difference between the reliability functions.

We make 1,000 replications with different ending points to obtain the number of times that the null hypothesis is not rejected. The concordance percentages of the tests presented in Table 6 show that concordance increases as the traffic intensity, and ending point increases.

Table 6:
Concordance percentage of the Log-rank and Chi-squared tests.

We obtain the E (Y ) by Equation (9) where r =1 and the sum was truncated (see Appendix APPENDIX In this appendix, we show the programs in software R used in this paper and the truncation of the sum in E (Y ). 1. The programs in software R used in this paper. • Programs in software R to simulate the M /M /1 model, Test Log-rank and Chi-squared tests. require(VGAM) rejlog=0 rejchi=0 for (iin 1:1000) n =10000 stop nat =1 number of the channels λ =0.9 the arrival rate mi =1 the service rate v=vector() customers.arrival=0 customes.service=0 customers.queue=0 system=numeric(n) wait=numeric(n) time = numeric(n) time-wait=numeric(n) clock1=numeric(nat) y=vector() time-service=numeric(n) alarm-clock=rep(-1,nat) for(i;in 2:n) { time.atual = (i -1)*0.1 time[i] = (i -1)*0.1 arrival = rpois (1,0.1) customers-arrival= customers-arrival + arrival customers-queue = customers-queue + arrival system[i] = system[i-1] + arrival if (sum (alarm -clock ==-1)>0 and customers -queue >0) { vagas = sum (alarm -clock ==-1) if (place !=sum (alarm -clock ==-1)) } stop("non conforming supositions! ") } minimum = min(place,customers-queue) for (j in 1:minimum) { clock =rexp (1,0.9) clock 1[which (alarm -clock ==-1)[1]]=clock alarm -clock [which (alarm -clock ==-1)[1]]=time .atual +clock customers -queue =customers -queue -1} } for (j in 1:nat ) { if (time .atual >=alarm -clock [j ] e alarm -clock [j ]!=-1) { system [i ]=system [i ]-1 alarm -clock [j ]=-1 customers -service =customers -service +1 y [customers -service ]=clock 1[j ] clock 1[j ]=0 } } Log-rank and Chi-squared tests time =c (y ,s 2) group =c (rep (0,length (y )),rep (1,length (s 2))) status =rep (1,length (y )+length (s 2)) mice =data .frame (time ,status ,group ) s =survdiff (Surv (time ,status ==1) group ,data =mice ) schi =chisq .test (y ,s 2) if (schisq >3.84) { rejlog =rejlog +1 } if (schip .value <0.05) { rejchi <-rejchi +1 } } • Program reliability function and estimated reliability function by Kaplan-Meiermethod. dado <-cbind (y ,cens =1) s 2=numeric () mu =0.9 traffic intensity sigma =0.9 service rate s 2=exp (mu *(1-exp (-sigma *y )))*exp (-sigma *y )/(exp (mu )-1) MAXPE distribution s 2=exp (-y *sigma )*(1-mu )*sigma /(1-mu *(1-exp (-y *sigma ))) MAXGE distribution require (survival ) tempos =dado [,1] cens =dado [,2] ekm =survfit (Surv (y ,cens ) 1) plot (ekm ,lty =1,lwd =2,conf .int =F ,xlab ="time ", ylab ="reliabilityfunction ") lines (y ,s 2,col ='blue ',lty =1,lwd =2) text (7,1,"trafficintensity =0.8") legend (3,0.8,c ('M /M /1','W (y )=0.876'),bty ='1',lty =c (1,4,3,2), lwd =c (2,2,2,2)) • Program optim for the MLEs. verosi <-function (param ) { mu <-param [1] sigma <-param [2] if (any(mu <1e - 200))return(.Machinedouble.xmax 5) if (any(sigma <1e - 200))return(.Machinedouble.xmax 5) f = sigma * (exp(-sigma * y)) * (1 - mu)/(1 - mu* (1 - exp(-sigma * y)))2 MAXGE distribution lv <-log (f ) sum (-lv ) } estima <-optim (c (0.888,0.9),verosi ,method ="BFGS ", hessian=T); estima estimahessian solve (estimahessian ) Criterions AIC and BIC lm 1<-lm (smax 1) AIC (lm 1) stopifnot (all .equal (AIC (lm 1), AIC (logLik (lm 1)))) BIC (lm 1) lm 2<-update (lm 1,...-Examination ) AIC (lm 1,lm 2) BIC (lm 1,lm 2) 2. The truncation of the sum in for all the models: The is We have an upper bound as The series τ has an upper bound given by the series that is convergent (Minka et al., 2003) because Therefore, the series τ is convergent. Furthermore, it is possible to truncate the series τ at some υ-th term such that where where R v is absolute truncated error. An upper bound can be found based in the series ρm /(m !)Φ, m =1,2,...,, decreases at a faster rate then a geometric series. Thus, there exists 0 < εv < 1 for all m > υ, so that ρm /(m + 1)Φ < εv Therefore, R v < εv , the error is small. ) and it is possible to calculate by MAPLE or in the software R, therefore, we obtained that E (Y )=3.87 unit time and E (Y )=2.51 unit time when the intensity traffic is ρ =0.8 and ρ =0.9, respectively.

In Figure 3 we compare the MAXGE distribution with the exponential-Conway-Maxwell Poisson distribution, denoted by ECOMP distribution (Cordeiro et al., 20129 CORDEIRO GM, RODRIGUES J & DE CASTRO M. 2012. The exponential com-poisson distribution. Statistical Papers, pages 1-12.), and the M /M /1 model. The ECOMP distribution models is for the minimum time.

Figure 3
Comparison between the reliability functions of the MAXGE distribution, ECOMP distribution and M/M/1 model with traffic intensity ρ =0.8 and ρ=0.9.

In order to select the best distribution to explain the data, we consider the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) which are defined by -2(ℓ()+q ) and (ℓ()+q log(n )), respectively, where q is the number of estimated parameters under the distribution and n is the sample size. The best distribution has lower values of AIC and BIC (Akaike, 19742 AKAIKE H. 1974. A new look at the statistical model identification. Automatic Control, IEEE Transactions on, 19(6): 716-723. and Schwarz et al., 197823 SCHWARZ G ET AL. 1978. Estimating the dimension of a model. The annals of statistics, 6(2): 461-464.).

In Table 7 we show the values of the criterions AIC and BIC for the MAXGE distribution, ECOMP distribution and M /M /1. We observe that the MAXGE distribution has the lower values of AIC and BIC.

Table 7
Values of the criterions AIC and BIC.

4.2 Real Data

Next, we analyze two real data: first of all, the express checkout of a Brazilian supermarket and secondly the accesses to the website.

4.2.1 The express checkout

A Brazilian supermarket was chosen to collect the data. It has ten normal checkouts and three express checkouts. In a particular day of collection of the data, one of the express checkouts had the same behavior of the studied system. The sample consists of the 85 service times in minutes.

In Table 8 we show the summary of the express checkout data. We observe that the maximum service time is 5.93 minutes and the minimum service time is 1.6 minutes.

Table 8
The summary of the express checkout data.

We use the Kolmogorov-Smirnov test to check whether a sample comes from a population with exponential distribution (H0 ). Therefore, with the significance level α =0.05 and p -value =0.3497 the null hypothesis is not rejected for the exponential distribution (Lilliefors, 196717 LILLIEFORS HW. 1967. On the kolmogorov-smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association, 62(318): 399-402.).

In the express checkout, the service rate is independent on state of the system, therefore, the pressure parameter is zero, Ф =0. For this reason, we choose the MAXGE distribution to model the maximum service time.

In Table 9 we show the MLEs of the parameters ρ and λ, the confidence intervals (CI) 95%and sde.

Table 9
The MLEs, confidence intervals (CI) 95% and sde.

Using the Equation (9) with r =1 and we truncated the sum in the sample size. We obtain E (Y )=2.72 minutes for the maximum service time.

In Figure 4 we compare the reliability function of the MAXGE distribution with the estimated reliability function. We observe that the reliability function of the MAXGE distribution fits with the estimated reliability function.

Figure 4
Comparison between the reliability function of the MAXGE distribution and the estimated reliability function.

In Figure 5 we compare the MAXGE distribution, ECOMP distribution and M /M /1 model.

Figure 5
Comparison between the reliability functions of the MAXGE distribution, ECOMP distribution and M/M/1 model.

In Table 10 we show that the MAXGE distribution has the lower values of the criterions AICand BIC.

Table 10
Value of the criterions AIC and BIC.

4.2.2 The accesses to the website

We used collected data of the accesses on web site "Tendências Profissionais" where the main aim was to collect data to define the profile of the communication professionals in Brazil (the data are of the Midia group). For this purpose, a survey with 26 questions was allocated on this website. The survey was distributed through social networks (Facebook, Twitter, Orkut, MySpace, Youtube) and 1,000 emails were sent to the main communication agencies in Brazil. On the website "Tendências Profissionais", the arrivals and the outputs of the customers were registered. Data collection started on 20 October, 2010 and it was available for 20 days. We analyze in a particular day which we collected 385 accesses on website.

For access the Tendências Profissionais:

(http://tendenciasprofissionais.wordpress.com/grupodepesquisa-mid/).

In Table 11 we show the summary of the accesses on website. We observe that the maximum time inter-access is 15.2 seconds and the minimum inter-access is 0.0001 seconds.

Table 11
The summary of the accesses on webite in seconds.

The Kolmogorov-Smirnov test was used to check whether a sample comes from a population which is exponential distributed (H0 ), with the significance level α =0.05 and p -value =0.2497 the null hypothesis is not rejected for the exponential distribution.

On the website each access can be considered a service channel, therefore, the service rate dependent of the state of the system, in this case the pressure parameter assumes values one, Φ. For this reason, a possible model for the maximum inter-access is the MAXPE distribution.

In Table 12 we show the MLEs of the parameters ρ and λ, the confidence intervals 95%(CI) and sde for the parameters.

Table 12
The MLEs, the confidence intervals 95% (CI) and the sde for the parameters of the MAXPE distribution.

We calculate the E (Y ) by Equation (12) with r=1. We truncated the sum in the sample size and we obtain E (Y )=8.42 seconds.

In Figure 6 we compare the reliability function of the MAXPE distribution with the estimatedreliability function. We observe that the reliability function of the MAXPE distribution is a closely concordance with estimated reliability function.

Figure 6
Comparison between the reliability function of the MAXPE distribution and the estimated reliability function.

In Figure 7 we compare the MAXPE distribution with the M /M /∞ model.

Figure 7
Comparison between the reliability functions of the MAXPE distributionwith the M/M/∞ model.

In Table 13 we observe that the MAXPE distribution has the lower values criterions AIC and BIC, therefore, the MAXPE distribution is the best model.

Table 13
Value of the criterions AIC and BIC.

5 CONCLUSIONS

In this paper, we proposed the MAXCOMPE distribution for the maximum service time in which the service rate is dependent on state of the system and the arrival of customer is attached to the service. The MAXCOMPE distribution was obtained by compound distributions in which we used the zero truncated COMP distribution and the exponential distribution for the maximum service time. This distribution has an adjustment mechanism, described by the pressure parameter Ф, in order to re-establish the equilibrium of the system when the number of the customer increases. We studied three values for the pressure parameter of the MAXCOMPE distribution, Ф =0,1 and Φ → ∞, consequently, other distributions are obtained. For Ф =0, the service rate independent on state of the system and the server is able to absorb all the works, in this case, MAXCOMPE distribution becomes the MAXGE distribution. For Ф =1, the service rate increases proportionally of the state of the system, the new service channels are open and the MAXCOMPE distribution becomes the MAXPE distribution. Finally, for Φ → ∞ the service is very fast, we observe only one customer in the system and the MAXCOMPE distribution becomes the MAXBE distribution. The properties of the proposed distribution are discussed, including a formal proof of its density function and moments. To illustrate the applicability of the model the simulations and real date are used. We simulated two systems with different traffic intensities, such as, ρ=0.8 and 0.9. These values are chosen, because the studied system has heavy traffic, for this reason, the parameter of the traffic intensity has high values. A restriction is the used the small traffic intensity because the MAXCOMPE distribution is zero truncated, therefore, the system is never idle. Furthermore, we observed that the sample mean of the estimates of the parameters is near of the true value. In addition, the fit of the proposed model was better when compared with the models: ECOMP, M /M /1 and M /M /∞. We used the criterions AIC e BIC to select the best distribution to explain the data. We used two real data, data of the express checkout for illustrated the MAXGE distribution where there is single server. And data of the access on the website for illustrated the MAXPE distribution. To conclude, we believe that the MAXCOMPE distribution has a practical approach and it can be applied to various practical situations.

REFERENCES

  • 1
    ADAMIDIS K, DIMITRIKAPOULO T & LOUKAS S. 2005. On an extension of the exponential-geometric distribution. Statistics & probability letters, 73(3): 259-269.
  • 2
    AKAIKE H. 1974. A new look at the statistical model identification. Automatic Control, IEEE Transactions on, 19(6): 716-723.
  • 3
    BARRETO-SOUZA W, DE MORAIS AL & CORDEIRO GM. 2011. The weibull-geometric distribution. Journal of Statistical Computation and Simulation, 81(5): 645-657.
  • 4
    BARTH W, MANITZ M & RAIK S. 2010. Analysis of two-level support systems with time-dependent overflow a banking application. Production and Operations Management, 19(6): 757-768.
  • 5
    BEKKER R, BORST S, BOXMA O & KELLA O. 2004. Queues with workload-dependent arrival and service rates. Queueing Systems, 46(3): 537-556.
  • 6
    BEKKER RGK, NIELSEN B & BANG T. 2011. Queues with waiting time dependent service. Queueing Systems, 68(1): 61-78.
  • 7
    BOXMA O & VLASIOU M. 2007. On queues with service and interarrival times depending on waiting times. Queueing Systems, 56(3): 121-132.
  • 8
    CANCHO VG, LOUZADA F & BARRIGA GD. 2011. The geometric-birnbaum-saunders regression model with cure rate. Journal of Statistical Planning and Inference.
  • 9
    CORDEIRO GM, RODRIGUES J & DE CASTRO M. 2012. The exponential com-poisson distribution. Statistical Papers, pages 1-12.
  • 10
    DROGUETT EL & MOSLEH A. 2007. Time to failure assessment of products at service conditions from accelerated lifetime tests with stress-dependent spread in life. Pesquisa Operacional, 27(2): 209-233.
  • 11
    GRADSHTEIN IS, RYZHIK IM, JEFFREY A & ZWILLINGER D. 2007. Table of integrals, series, and products. Academic press.
  • 12
    IYER SK & MANJUNATH D. 2006. Queues with dependency between interarrival and service times using mixtures of bivariates. Stochastic models, 22(1): 3-20.
  • 13
    JONGBLOED G & KOOLE G. 2001. Managing uncertainty in call centres using poisson mixtures. Applied Stochastic Models in Business and Industry, 17(4): 307-318.
  • 14
    KAPLAN EL & MEIER P. 1958. Nonparametric estimation from incomplete observations. Journal of the American statistical association, 53(282): 457-481.
  • 15
    KIMURA T. 1991. Approximating the mean waiting time in the gi/g/s queue. Journal of the Operational Research Society, pages 959-970.
  • 16
    KLEINBAUM DG & KLEIN M. 2012. Kaplan-meier survival curves and the log-rank test. In Survival analysis, pages 55-96. Springer.
  • 17
    LILLIEFORS HW. 1967. On the kolmogorov-smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association, 62(318): 399-402.
  • 18
    MARIN CV, DRURY CG, BATTA R & LIN L. 2007. Server adaptation in an airport security system queue. OR Insight, 20(4): 22-31.
  • 19
    MASSEY JR FJ. 1951. The kolmogorov-smirnov test for goodness of fit. Journal of the American statistical Association, 46(253): 68-78.
  • 20
    MINKA TP, SHMUELI G, KADANE JB, BORLES S & BOATWRIGHT P. 2003. Computing with the com-poisson distribution. Pittsburgh, PA: Department of Statistics, Carnegie Mellon University.
  • 21
    MORABITO R & DE LIMA FC. 2000. Um modelo para analisar o problema de filas em caixas de supermercados: um estudo de caso. Pesquisa Operacional, 20(1): 59-71.
  • 22
    RIGBY R & STASINOPOULOS D. 2005. Generalized additive models for location, scale and shape. Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(3): 507-554.
  • 23
    SCHWARZ G ET AL. 1978. Estimating the dimension of a model. The annals of statistics, 6(2): 461-464.
  • 24
    SHARMA K & SHARMA G. 1994. A delay dependent queue without pre-emption with general linearly increasing priority function. Journal of the Operational Research Society, pages 948-953.
  • 25
    WHITT W. 1989. Planning queueing simulations. Management Science, 35(11): 1341-1366.
  • 26
    WHITT W. 1990. Queues with service times and interarrival times depending linearly and randomly upon waiting times. Queueing Systems, 6(1): 335-351.
  • 27
    WHITT W. 2003. How multiserver queues scale with growing congestion-dependent demand. Operations Research, 51(4): 531-542.

APPENDIX

In this appendix, we show the programs in software R used in this paper and the truncation of the sum in E (Y ).

1. The programs in software R used in this paper.

• Programs in software R to simulate the M /M /1 model, Test Log-rank and Chi-squared tests.

require(VGAM)

rejlog=0

rejchi=0

for (iin 1:1000) n =10000 stop

nat =1 number of the channels

λ =0.9 the arrival rate

mi =1 the service rate

v=vector()

customers.arrival=0

customes.service=0

customers.queue=0

system=numeric(n)

wait=numeric(n)

time = numeric(n)

time-wait=numeric(n)

clock1=numeric(nat)

y=vector()

time-service=numeric(n)

alarm-clock=rep(-1,nat)

for(i;in 2:n)

{

time.atual = (i -1)*0.1

time[i] = (i -1)*0.1

arrival = rpois (1,0.1)

customers-arrival= customers-arrival + arrival

customers-queue = customers-queue + arrival

system[i] = system[i-1] + arrival

if (sum (alarm -clock ==-1)>0 and customers -queue >0)

{

vagas = sum (alarm -clock ==-1)

if (place !=sum (alarm -clock ==-1))

} stop("non conforming supositions! ")

} minimum = min(place,customers-queue)

for (j in 1:minimum) { clock =rexp (1,0.9)

clock 1[which (alarm -clock ==-1)[1]]=clock

alarm -clock [which (alarm -clock ==-1)[1]]=time .atual +clock

customers -queue =customers -queue -1}

}

for (j in 1:nat )

{ if (time .atual >=alarm -clock [j ] e alarm -clock [j ]!=-1) { system [i ]=system [i ]-1

alarm -clock [j ]=-1

customers -service =customers -service +1

y [customers -service ]=clock 1[j ]

clock 1[j ]=0 }

}

Log-rank and Chi-squared tests time =c (y ,s 2)

group =c (rep (0,length (y )),rep (1,length (s 2)))

status =rep (1,length (y )+length (s 2))

mice =data .frame (time ,status ,group )

s =survdiff (Surv (time ,status ==1) group ,data =mice )

schi =chisq .test (y ,s 2)

if (schisq >3.84)

{

rejlog =rejlog +1

}

if (schip .value <0.05)

{

rejchi <-rejchi +1

}

}

• Program reliability function and estimated reliability function by Kaplan-Meiermethod.

dado <-cbind (y ,cens =1)

s 2=numeric ()

mu =0.9 traffic intensity

sigma =0.9 service rate

s 2=exp (mu *(1-exp (-sigma *y )))*exp (-sigma *y )/(exp (mu )-1) MAXPE distribution

s 2=exp (-y *sigma )*(1-mu )*sigma /(1-mu *(1-exp (-y *sigma ))) MAXGE distribution

require (survival )

tempos =dado [,1]

cens =dado [,2]

ekm =survfit (Surv (y ,cens ) 1)

plot (ekm ,lty =1,lwd =2,conf .int =F ,xlab ="time ",

ylab ="reliabilityfunction ")

lines (y ,s 2,col ='blue ',lty =1,lwd =2)

text (7,1,"trafficintensity =0.8")

legend (3,0.8,c ('M /M /1','W (y )=0.876'),bty ='1',lty =c (1,4,3,2),

lwd =c (2,2,2,2))

• Program optim for the MLEs.

verosi <-function (param ) {

mu <-param [1]

sigma <-param [2]

if (any(mu <1e - 200))return(.Machinedouble.xmax 5)

if (any(sigma <1e - 200))return(.Machinedouble.xmax 5)

f = sigma * (exp(-sigma * y)) * (1 - mu)/(1 - mu* (1 - exp(-sigma * y)))2 MAXGE distribution

lv <-log (f )

sum (-lv )

}

estima <-optim (c (0.888,0.9),verosi ,method ="BFGS ", hessian=T);

estima

estimahessian

solve (estimahessian )

Criterions AIC and BIC

lm 1<-lm (smax 1)

AIC (lm 1)

stopifnot (all .equal (AIC (lm 1),

AIC (logLik (lm 1))))

BIC (lm 1)

lm 2<-update (lm 1,...-Examination )

AIC (lm 1,lm 2)

BIC (lm 1,lm 2)

2. The truncation of the sum in for all the models:

The is

We have an upper bound as

The series τ has an upper bound given by the series that is convergent (Minka et al., 200320 MINKA TP, SHMUELI G, KADANE JB, BORLES S & BOATWRIGHT P. 2003. Computing with the com-poisson distribution. Pittsburgh, PA: Department of Statistics, Carnegie Mellon University.) because

Therefore, the series τ is convergent.

Furthermore, it is possible to truncate the series τ at some υ-th term such that

where

where R v is absolute truncated error.

An upper bound can be found based in the series ρm /(m !)Φ, m =1,2,...,, decreases at a faster rate then a geometric series. Thus, there exists 0 < εv < 1 for all m > υ, so that ρm /(m + 1)Φ < εv Therefore, R v < εv , the error is small.

Publication Dates

  • Publication in this collection
    Sep-Dec 2015

History

  • Received
    30 June 2014
  • Accepted
    03 Sept 2015
Sociedade Brasileira de Pesquisa Operacional Rua Mayrink Veiga, 32 - sala 601 - Centro, 20090-050 Rio de Janeiro RJ - Brasil, Tel.: +55 21 2263-0499, Fax: +55 21 2263-0501 - Rio de Janeiro - RJ - Brazil
E-mail: sobrapo@sobrapo.org.br