# Performance Evaluation and Dimensioning ofGIX/M/c/N Systems through Kernel Estimation

Authors: G. M. Gontijo; G. S. Atuncar; F. R. B. Cruz; L. Kerbache
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

We extend the analysis of queueing systems for real-life situations, where the arrival pattern of customers is unknown. In real systems, we must understand how the choice of a method of estimation influences the configuration of the system. Using kernel smoothing, we evaluate algorithms to estimate performance measures of aGIX/M/c/N system, including the invariant probability distribution of the number of customers in the system, the blocking probability, the average queue size, and the average client queue time. We successfully apply the method to the arrivals to a call center to plan and improve the performance of these important queueing systems.


---

## Body

## 1. Introduction

There is a large practical interest in investigating the behavior of general-arrival queueing systems, namely, those ofGIX/M/c/N type, because when managing real queueing systems, the behavior of the arrival process generally is not Markovian. In Kendall [1] notation, in these queueing systems, the interarrival times are independent and follow a general distribution (GI), the service times follow an exponential distribution (M), we have c identical servers working in parallel, and a maximum capacity of N users that are simultaneously allowed in the system, including those in service. Finally, X is a random variable representing the size of group (bulk) arrivals. Such queueing systems could be used in situations where we have relative control over how the servers work, but we do not know beforehand how the customers arrive at the system.Naturally, the mathematical model depends on the type of queueing system considered, and there are several methods to obtain such models. The most widely used methods are those that attempt to explain the density functions of interarrival and service times by means of parametric statistical models. Nevertheless, real data sometimes do not fit well into parametric models; instead, they produce intractable models. Exact results for performance evaluation of Markovian and some simple general queueing systems are known [2], but such systems sometimes are not found in real life.An alternative approach is using nonparametric methods to study queueing systems. Nonparametric methods that use kernel smoothing have received much attention lately [3]. Kernel estimators provide a simple way of finding structure in data sets without imposing a specific parametric model [4], which gives us flexibility to handle virtually any data set. There is extensive literature discussing queueing systems [5] and kernel smoothing [4] as separate concepts, but virtually no study has brought the two concepts together.The contribution of this paper is twofold. First, we implement algorithms to calculate performance measures of queueing systems, where the density of the interarrival time is estimated by the kernel method. Second, we evaluate the performance of these algorithms as a function of the kernel estimator, the smoothing window, the intensity rate, and the system size. We also present a case study in a call center that illustrates the usefulness of the method.First, we present a literature review about queueing systems and kernel estimators as well as the fundamental concepts required to understand the proposed queueing system model. We also discuss the use of kernel estimators, the chosen models, and the issue of selecting the smoothing parameter. Then, we describe how we estimated system performance and present the comparative results of simulations with the different methods discussed. We then apply our methods to a call center case and end with our main conclusions and some ideas for future work in the area.


## 2. Literature Review and Fundamental Concepts

### 2.1. Previous Works

There are many situations in real life where queues occur and queueing models may be helpful. Recently, queueing models have been used successfully in manufacturing processes [6, 7], transportation [8], airports, ports, and product distribution systems [9], computer and telecommunication systems [10], call center modeling [11], and the analysis of health systems [12]. Queues may cause the quality of the services or the prices of the goods to rise or fall, depending on their efficiency [8], which may be estimated by means of the mathematical tools developed in queueing theory.As we mentioned earlier, there is not much literature that incorporates both general arrival queueing systems and kernel smoothing. Takács [13] analyzed a closed solution for various systems that have nonspecific distributions, including some multiserver queues such as GI/M/c and M/G/c. Hokstad [14] established some closed form results to the GI/M/c/N system. Chaudhry and Templeton [15] analyzed various types of queue with bulk arrivals. Zhao [16] proposed a closed form solution for the GIX/M/c system, and Bareche and Aïssani [17] proposed a method to evaluate the proximity of GI/M/1 and M/M/1 systems when the density of the interarrival time is estimated by kernel estimators.Concerning kernel smoothing, Wand and Jones [4] introduced general fundamental concepts. Regarding the issues related to the asymmetry of the random variable under analysis (e.g., nonnegativity), Zhang et al. [18] proposed a boundary corrected kernel estimator based on pseudodata generation, transformation, and reflection around the Y-axis. Chen [19] proposed the use of a gamma kernel to avoid boundary problems present in certain situations. Scaillet [20] studied the application of other asymmetric kernels. Bouezmarni and Scaillet [21] were concerned about the consistency of these asymmetrical estimators. For recent developments in the area of kernel smoothing and a thorough literature review, see Atuncar et al. [22] and de Lima and Atuncar [3].


### 2.2. TheGIX/M/c/N Model with Partial Blocks

Vijaya Laxmi and Gupta [23] described a generalization of the system GI/M/c/N when customers arrive in groups of size X with P(X=i)=gi (i≥1) and mean E(X)=g̅. The GIX/M/c/N is a finite capacity system such that a customer that arrives to the saturated system is refused with a probability that we will call PBL, or partial blocks. This refers to the case in which an arrival group of a size greater than the remaining spots in the system is partially denied according to the number of the remaining vacancies until the system is complete. Let Bn be the number of clients who were served between the arrival of the nth customer and its successor. Therefore, the number of clients the nth customer finds in the system at the arrival, Yn, would depend on Xn and Bn such that Yn+1=[min(Yn+Xn,N)-Bn]+. Because Yn+1 depends only on Yn, Bn, and Xn but not on Yn-1, Yn-2, and so on, the stochastic process {Yn} is a first-order Markov Chain.Vijaya Laxmi and Gupta [23] report that when the traffic intensity rate ρ=(λg̅/cμ), where λ is the arrival rate and μ service rate, is smaller than 1, this Markov Chain has an invariant probability distribution of πk=limn→∞P(Yn=k), k=0,1,2,… associated with the number of clients an arbitrary customer finds in the system at arrival. The πk’s are often called prearrival probabilities.Prearrival probabilities can be determined by the following system of linear equations:(2.1)πk=∑j=0Npjkπj,
for k=0,1,2,…,N-1, and(2.2)∑j=0Nπj=1,
which gives(2.3)[(p0,0-1)p1,0⋯pN,0p0,1(p1,1-1)⋮⋮⋱⋮p0,N-1⋯⋯(pN,N-1-1)1⋯⋯1][π0π1⋮⋮πN]=[00⋮01],
where pjk are named transition probabilities such that(2.4)pjk=P{Yn+1=k∣Yn=j}=∑i=1∞P{Xn=i,Yn+1=k∣Yn=j}=∑i=1∞giP{Yn+1=k∣Yn=j,Xn=i}=∑i=1∞giP(Bn=min⁡(j+i,N)-k)=∑i=max⁡(1,j-k)N-jgiP(Bn=j+i-k)+∑i=N-j+1∞giP(Bn=N-k)=∑i=max⁡(1,j-k)N-jgi∫0∞πj+i,k(z)dτ(z)+∑i=N-j+1∞gi∫0∞πN,k(z)dτ(z),
where, j≤N, k≤N, πj+i,k(z) is the serving probability of (j-k+i) clients under the assumption of the interarrival time τn=z and  τ(z) is the interarrival time distribution.We will analyze howπj+i,k(z) behaves for all possibilities of (j+1) and (k) in regard to the number of servers c.(a)(j+i)≥c and k≥c
Whenj+i≥c and k≥c, there will be more clients than the servers can handle in the entire interval. Because the service process is Markovian, we can treat the server group as a single unit that serves customers at a rate cμz and a Poisson distributed transition probability
(2.5)πj+i,k(z)=e-cμz(cμz)j+i-k(j+i-k)!.(b)(j+i)<c
Whenj+i<c, all clients within the system are being served and only k customers will remain in the system to time z. Knowing that the probability of a service time greater than z is e-μz, we can describe this transition probability as a Binomial distribution
(2.6)πj+i,k(z)=(j+ik)e-kμz(1-e-μz)j+i-k.(c)(j+i)≥c and k<c
Whenj+i≥c and k<c, there will be (j+i-c) customers waiting and c customers being served at the beginning of the interval, but there will be (c-k) free servers at the end. Let y be the interval ending immediately before (j+i-c+1) clients are served. If each service time is exponentially distributed with a rate cμ, then y is gamma distributed with a shape (j+i-c+1) and rate a cμ. The other c customers will be served in a time (z-y), and only k will remain. The transition probability of this subinterval will follow a Binomial distribution with success probability e-μ(z-y). The values πj+i,k(z) can be obtained by the convolution of these two variables
(2.7)πj+i,k(z)=∫0ze-cμy(cμ)j+i-c+1yj+i-c(j+i-c)!(ck)e-kμ(z-y)(1-e-μ(z-y))c-kdy=(ck)e-kμz∫0z(cμy)j+i-c(j+i-c)!ekμye-cμy(1-eμye-μz)c-kcμdy=(ck)e-kμz∫0z(cμy)j+i-c(j+i-c)![e-μy(1-eμye-μz)]c-kcμdy=(ck)e-kμz[∫0z(cμy)j+i-c(j+i-c)!(e-μy-e-μz)c-kcμdy].
Therefore, from (2.4), we can get the transitions probabilities pjk given by (2.4) when we take each transition probability πj+i,k(z) as πN,k(z) when i>N-j(2.8)pjk={∑i=max⁡{1,k-j}N-jβj+i-kgi+βN-k∑k=N-j+1∞gi,k≥c,∑i=max⁡{1,k-j}N-jVj+i,kgi+VN,k∑k=N-j+1∞gi,0<k<c,1-∑r=1Npjr,k=0,
where Vj+i,k and βj+i-k are as follows:
(2.9)Vj+i,k={0,j<k<c,∫0∞(j+ik)e-kμz(1-e-μz)j+i-kdτ(z),k≤j≤c,∫0∞∫0z(ck)e-kμz(cμy)j-c(j-c)!cμ(e-μy-e-μz)c-kdydτ(z),k<c<j,βj+i-k=∫0∞e-cμz(cμz)j+i-k(j+i-k)!dτ(z),j+i-k≥0.
There is a relationship between the vector of prearrival probabilities π and the vector of arbitrary time probabilities P, related to the number of people that an outside observer finds in the system. Viajaya Laxmi and Gupta [23] established a method that proved the relationship
(2.10)Pk={ρcmin⁡{k,c}g̅∑i=0k-1πi∑j=k-i∞gj,0<k≤N,1-∑i=1NPi,k=0.
Performance measures (below) are used, including the average queue length Lq, blocking probability of an arbitrary customer PBL, and average waiting time in the queue Wq, to analyze the efficiency of queueing systems
(2.11)Lq=∑i=0N(i-c)Pi,PBL=∑i=0Nπi∑j=N-i∞g̅-1∑k=j+1∞gj,Wq=Lq[g̅λ(1-PBL)].


### 2.3. Kernel Estimators

Suppose that we have a sample of the interarrival times,X1,…,Xn, with an unknown density τ(t). The kernel estimator is an analytical tool that provides an effective way of revealing the structure behind such a sample.


#### 2.3.1. Gamma Kernel Estimator

Recently, Chen [19] suggested an asymmetric kernel with naturally varying shape as a way to avoid allocating weight for negative values. The gamma kernel estimators are always nonnegative, free of boundary bias, and achieve the optimal rate of convergence for the mean square error (MSE) in the nonnegative kernel estimator class. Bouezmarni and Scaillet [21] showed that this estimator is consistent and able to avoid boundary bias. Let KG(p,q) be the gamma density function with parameters p (shape) and q (rate). The gamma kernel considered is(2.12)KG(tb+1,b)(Xj)=Xjt/be-Xj/bb(t/b)+1Γ[(t/b)+1],
where b is a smoothing parameter that satisfies the condition b→0, nb→∞ as n→∞. The gamma kernel estimator is(2.13)τ̂(t;b)=n-1∑j=1nKG(tb+1,b)(Xj).
The smoothing parameter b is critical for the overall performance of the kernel estimator considered. A small b leads to a relatively bumpy density, while a large one results in a smooth density. There are several methods to determine the best fit, from a minimization of the mean integrated squared error (MISE) of τ̂(t;b) to the asymptotic behavior of the MISE (AMISE).(a) Least Squares Cross-Validation (LSCV) Method
The least squares cross-validation (LSCV) method starts from the MISE expansion(2.14)MISE{τ̂(x;b)}=E∫τ̂(x;b)2dx-2E∫τ̂(x;b)τ(x)dx+E∫τ(x)2dx.
The minimization of the first term is equivalent to the minimization of
(2.15)MISE{τ̂(x;h)}-E∫τ(x)2dx=E[∫τ̂(x;h)2dx-2∫τ̂(x;h)τ(x)dx].
The right-hand side is unknown, because it depends on τ. However, an unbiased estimator for this quantity is
(2.16)LSCV(h)=∫τ̂(x;h)2dx-2n-1∑i=1nτ̂-i(Xi;h),
where τ̂-i=(Xi;h) is the density estimate based on the sample with Xi deleted; this is often called the “leave-one-out” density estimator. A disadvantage of this method is that it suffers from high variation.(b) Asymptotic Behavior of the MISE (AMISE) Method
An alternative parameter selector is to consider the asymptotic behavior of the MISE of the gamma kernel estimator. Chen [19] uses some aspects of the gamma distribution and a Taylor expansion to determine the MISE as follows:
(2.17)MISE(τ̂)=b2∫0∞{xτ′(x)+12xτ′′(x)}2dx+(2nbπ)-1∫0∞x-1/2τ(x)dx+O(n-1b-1/2+b2).
The asymptotic MISE disregards the last term; therefore, the optimal b  that minimizes the leading terms above is
(2.18)bAMISE=[(2π)-1∫0∞x-1/2τ(x)dx4n∫0∞(xτ′(x)+2-1xτ′′(x))2dx]2/5,
where the functions τ, τ′, and τ′′ are unknown. These quantities are obtained from the fitted gamma density with parameters adjusted from the sample. This solution still requires further study, but our paper shows promising results.


#### 2.3.2. Zhang et al. [18] Estimator

Zhang et al. [18] submitted a model that works particularly well when τ(0)>0 and combines pseudodata creation, its transformation, and its reflection around the Y-axis in the following three steps.Step 1.
Transform the original dataX1,…,Xn to g(X1),…,g(Xn), while keeping the original data, where g is a nonnegative, continuous, and monotonically increasing function from [0,∞) to [0,∞).
Based on extensive simulations, the transformation that best suits a broad variety of densities is(2.19)g(x)=x+dx2+Adx3,
where A>1/3 and d=f′(0)/f(0).Step 2.
Reflect the pseudodata,g(X1),…,g(Xn), around the origin.Step 3.
Based on the enlarged data sample,-g(X1),…,-g(Xn),X1,…,Xn, define the new estimator as
(2.20)τ̂n(x,h)=1nh∑j=1n{K(x-Xjh)+K(x+g(Xj)h)},x≥0,
where h is a smoothing parameter and K is a symmetric probability function with support [-1,1] like the Epanechnikov kernel
(2.21)K(t)=34(1-t2)I[-1,1].
Notice that the transformation g defined above is not available in practice, because d is unknown. A good estimator can be obtained when d is written as (d/dx)logf(x)∣x=0,
(2.22)dn=log⁡fn(h)-log⁡fn(0)h,
where fn(h)=fn*(h)+1/n2, fn(0)=max(fn*(0),1/n2), fn*(h)=(1/nh)∑j=1nK((x-Xj)/h), fn*(0)=(1/nh0)∑j=1nK0(Xj/h0), and K0 is a so-called endpoint kernel, satisfying
(2.23)∫-10K0(t)dt=1,∫-10tK0(t)dt=0,∫-10t2K0(t)dt≠0,h0={[∫-11K(t)dt]2∫-10K0(t)2dt[∫-10t2K0(t)dt]2∫-11K(t)2dt}h.
Zhang et al. [18] proved that for t≥h, the effect of reflected pseudodata is insignificant, and the estimator can be reduced to the Parzen-Rosenblatt estimator
(2.24)τ̂n(t,h)=1nh∑j=1nK(t-Xjh).
They also established that
(2.25)∫0∞τ̂n(t)dt=1+1n∑i=1n∫-Xi/h-g(Xi)/hK(z)dz.
Thus, τ̂n(t) only integrates to 1 when dn=0, so gn(Xi)=Xi, or when Xi=0 for all Xi’s, because gn(0)=0. However, when n→∞, both limits of the second term will eventually converge to 0 and τ̂n(t) will integrate to 1 asymptotically.
Zhang and Karunamuni [24] used the endpoint kernel
(2.26)K0(t)=12(1+t)(12+t)I[-1,0]
and showed that this kernel minimizes the MSE when estimating τ(0). Therefore, h0=2h is approximately the optimal smoothing parameter for estimating τ(0) except when τ(0)=0.
Chiu [25] described a parameter selecting method that considered the optimal h that minimizes the asymptotic MISE when K is a symmetric probability function with up to the fourth moment being finite
(2.27)hAMISE=[∫K(x)2dxn(∫x2K(x)dx)2∫τ′′(x)2dx]1/5,
where the function ∫τ′′(x)2dx is unknown.
Chiu’s [25] “plug-in” method consists of estimating this quantity through the characteristic function of the sample
(2.28)φ̂(λ)=n-1∑j=1neiλXj,
and calculating the optimal h with the previous formula. The characteristic function of τ is
(2.29)φ(λ)=∫eiλxτ(x)dx.
By the inversion formula, we have
(2.30)τ(x)=(2π)-1∫eiλxφ(λ)dλ,
thus,
(2.31)τ′′(x)=(2π)-1∫λ2eiλxφ(λ)dλ,∫τ′′(x)2dx=∫[(2π)-1∫λ2eiλxφ(λ)dλ]2dx=∫λ4[(2π)-1∫eiλxφ(λ)dλ]2dx=∫λ4[τ(x)]2dx.
Using Parseval’s identity, we can show that
(2.32)∫λ4[τ(x)]2dx=(2π)-1∫λ4|φ(λ)|2dλ.
Chiu [25] introduced a cutoff value Λ for λ such that |φ̂(λ)|2<c/n. Extensive computational experiments from Bessegato et al. [26] show that c=3 is the value that minimizes the estimator variance. The final “plug in” estimator is then
(2.33)ĥAMISE=[∫K(x)2dxn(∫x2K(x)dx)2π-1∫0Λλ4[|φ̂(λ)|2-n-1]dλ]1/5.


## 2.1. Previous Works

There are many situations in real life where queues occur and queueing models may be helpful. Recently, queueing models have been used successfully in manufacturing processes [6, 7], transportation [8], airports, ports, and product distribution systems [9], computer and telecommunication systems [10], call center modeling [11], and the analysis of health systems [12]. Queues may cause the quality of the services or the prices of the goods to rise or fall, depending on their efficiency [8], which may be estimated by means of the mathematical tools developed in queueing theory.As we mentioned earlier, there is not much literature that incorporates both general arrival queueing systems and kernel smoothing. Takács [13] analyzed a closed solution for various systems that have nonspecific distributions, including some multiserver queues such as GI/M/c and M/G/c. Hokstad [14] established some closed form results to the GI/M/c/N system. Chaudhry and Templeton [15] analyzed various types of queue with bulk arrivals. Zhao [16] proposed a closed form solution for the GIX/M/c system, and Bareche and Aïssani [17] proposed a method to evaluate the proximity of GI/M/1 and M/M/1 systems when the density of the interarrival time is estimated by kernel estimators.Concerning kernel smoothing, Wand and Jones [4] introduced general fundamental concepts. Regarding the issues related to the asymmetry of the random variable under analysis (e.g., nonnegativity), Zhang et al. [18] proposed a boundary corrected kernel estimator based on pseudodata generation, transformation, and reflection around the Y-axis. Chen [19] proposed the use of a gamma kernel to avoid boundary problems present in certain situations. Scaillet [20] studied the application of other asymmetric kernels. Bouezmarni and Scaillet [21] were concerned about the consistency of these asymmetrical estimators. For recent developments in the area of kernel smoothing and a thorough literature review, see Atuncar et al. [22] and de Lima and Atuncar [3].


## 2.2. TheGIX/M/c/N Model with Partial Blocks

Vijaya Laxmi and Gupta [23] described a generalization of the system GI/M/c/N when customers arrive in groups of size X with P(X=i)=gi (i≥1) and mean E(X)=g̅. The GIX/M/c/N is a finite capacity system such that a customer that arrives to the saturated system is refused with a probability that we will call PBL, or partial blocks. This refers to the case in which an arrival group of a size greater than the remaining spots in the system is partially denied according to the number of the remaining vacancies until the system is complete. Let Bn be the number of clients who were served between the arrival of the nth customer and its successor. Therefore, the number of clients the nth customer finds in the system at the arrival, Yn, would depend on Xn and Bn such that Yn+1=[min(Yn+Xn,N)-Bn]+. Because Yn+1 depends only on Yn, Bn, and Xn but not on Yn-1, Yn-2, and so on, the stochastic process {Yn} is a first-order Markov Chain.Vijaya Laxmi and Gupta [23] report that when the traffic intensity rate ρ=(λg̅/cμ), where λ is the arrival rate and μ service rate, is smaller than 1, this Markov Chain has an invariant probability distribution of πk=limn→∞P(Yn=k), k=0,1,2,… associated with the number of clients an arbitrary customer finds in the system at arrival. The πk’s are often called prearrival probabilities.Prearrival probabilities can be determined by the following system of linear equations:(2.1)πk=∑j=0Npjkπj,
for k=0,1,2,…,N-1, and(2.2)∑j=0Nπj=1,
which gives(2.3)[(p0,0-1)p1,0⋯pN,0p0,1(p1,1-1)⋮⋮⋱⋮p0,N-1⋯⋯(pN,N-1-1)1⋯⋯1][π0π1⋮⋮πN]=[00⋮01],
where pjk are named transition probabilities such that(2.4)pjk=P{Yn+1=k∣Yn=j}=∑i=1∞P{Xn=i,Yn+1=k∣Yn=j}=∑i=1∞giP{Yn+1=k∣Yn=j,Xn=i}=∑i=1∞giP(Bn=min⁡(j+i,N)-k)=∑i=max⁡(1,j-k)N-jgiP(Bn=j+i-k)+∑i=N-j+1∞giP(Bn=N-k)=∑i=max⁡(1,j-k)N-jgi∫0∞πj+i,k(z)dτ(z)+∑i=N-j+1∞gi∫0∞πN,k(z)dτ(z),
where, j≤N, k≤N, πj+i,k(z) is the serving probability of (j-k+i) clients under the assumption of the interarrival time τn=z and  τ(z) is the interarrival time distribution.We will analyze howπj+i,k(z) behaves for all possibilities of (j+1) and (k) in regard to the number of servers c.(a)(j+i)≥c and k≥c
Whenj+i≥c and k≥c, there will be more clients than the servers can handle in the entire interval. Because the service process is Markovian, we can treat the server group as a single unit that serves customers at a rate cμz and a Poisson distributed transition probability
(2.5)πj+i,k(z)=e-cμz(cμz)j+i-k(j+i-k)!.(b)(j+i)<c
Whenj+i<c, all clients within the system are being served and only k customers will remain in the system to time z. Knowing that the probability of a service time greater than z is e-μz, we can describe this transition probability as a Binomial distribution
(2.6)πj+i,k(z)=(j+ik)e-kμz(1-e-μz)j+i-k.(c)(j+i)≥c and k<c
Whenj+i≥c and k<c, there will be (j+i-c) customers waiting and c customers being served at the beginning of the interval, but there will be (c-k) free servers at the end. Let y be the interval ending immediately before (j+i-c+1) clients are served. If each service time is exponentially distributed with a rate cμ, then y is gamma distributed with a shape (j+i-c+1) and rate a cμ. The other c customers will be served in a time (z-y), and only k will remain. The transition probability of this subinterval will follow a Binomial distribution with success probability e-μ(z-y). The values πj+i,k(z) can be obtained by the convolution of these two variables
(2.7)πj+i,k(z)=∫0ze-cμy(cμ)j+i-c+1yj+i-c(j+i-c)!(ck)e-kμ(z-y)(1-e-μ(z-y))c-kdy=(ck)e-kμz∫0z(cμy)j+i-c(j+i-c)!ekμye-cμy(1-eμye-μz)c-kcμdy=(ck)e-kμz∫0z(cμy)j+i-c(j+i-c)![e-μy(1-eμye-μz)]c-kcμdy=(ck)e-kμz[∫0z(cμy)j+i-c(j+i-c)!(e-μy-e-μz)c-kcμdy].
Therefore, from (2.4), we can get the transitions probabilities pjk given by (2.4) when we take each transition probability πj+i,k(z) as πN,k(z) when i>N-j(2.8)pjk={∑i=max⁡{1,k-j}N-jβj+i-kgi+βN-k∑k=N-j+1∞gi,k≥c,∑i=max⁡{1,k-j}N-jVj+i,kgi+VN,k∑k=N-j+1∞gi,0<k<c,1-∑r=1Npjr,k=0,
where Vj+i,k and βj+i-k are as follows:
(2.9)Vj+i,k={0,j<k<c,∫0∞(j+ik)e-kμz(1-e-μz)j+i-kdτ(z),k≤j≤c,∫0∞∫0z(ck)e-kμz(cμy)j-c(j-c)!cμ(e-μy-e-μz)c-kdydτ(z),k<c<j,βj+i-k=∫0∞e-cμz(cμz)j+i-k(j+i-k)!dτ(z),j+i-k≥0.
There is a relationship between the vector of prearrival probabilities π and the vector of arbitrary time probabilities P, related to the number of people that an outside observer finds in the system. Viajaya Laxmi and Gupta [23] established a method that proved the relationship
(2.10)Pk={ρcmin⁡{k,c}g̅∑i=0k-1πi∑j=k-i∞gj,0<k≤N,1-∑i=1NPi,k=0.
Performance measures (below) are used, including the average queue length Lq, blocking probability of an arbitrary customer PBL, and average waiting time in the queue Wq, to analyze the efficiency of queueing systems
(2.11)Lq=∑i=0N(i-c)Pi,PBL=∑i=0Nπi∑j=N-i∞g̅-1∑k=j+1∞gj,Wq=Lq[g̅λ(1-PBL)].


## 2.3. Kernel Estimators

Suppose that we have a sample of the interarrival times,X1,…,Xn, with an unknown density τ(t). The kernel estimator is an analytical tool that provides an effective way of revealing the structure behind such a sample.


### 2.3.1. Gamma Kernel Estimator

Recently, Chen [19] suggested an asymmetric kernel with naturally varying shape as a way to avoid allocating weight for negative values. The gamma kernel estimators are always nonnegative, free of boundary bias, and achieve the optimal rate of convergence for the mean square error (MSE) in the nonnegative kernel estimator class. Bouezmarni and Scaillet [21] showed that this estimator is consistent and able to avoid boundary bias. Let KG(p,q) be the gamma density function with parameters p (shape) and q (rate). The gamma kernel considered is(2.12)KG(tb+1,b)(Xj)=Xjt/be-Xj/bb(t/b)+1Γ[(t/b)+1],
where b is a smoothing parameter that satisfies the condition b→0, nb→∞ as n→∞. The gamma kernel estimator is(2.13)τ̂(t;b)=n-1∑j=1nKG(tb+1,b)(Xj).
The smoothing parameter b is critical for the overall performance of the kernel estimator considered. A small b leads to a relatively bumpy density, while a large one results in a smooth density. There are several methods to determine the best fit, from a minimization of the mean integrated squared error (MISE) of τ̂(t;b) to the asymptotic behavior of the MISE (AMISE).(a) Least Squares Cross-Validation (LSCV) Method
The least squares cross-validation (LSCV) method starts from the MISE expansion(2.14)MISE{τ̂(x;b)}=E∫τ̂(x;b)2dx-2E∫τ̂(x;b)τ(x)dx+E∫τ(x)2dx.
The minimization of the first term is equivalent to the minimization of
(2.15)MISE{τ̂(x;h)}-E∫τ(x)2dx=E[∫τ̂(x;h)2dx-2∫τ̂(x;h)τ(x)dx].
The right-hand side is unknown, because it depends on τ. However, an unbiased estimator for this quantity is
(2.16)LSCV(h)=∫τ̂(x;h)2dx-2n-1∑i=1nτ̂-i(Xi;h),
where τ̂-i=(Xi;h) is the density estimate based on the sample with Xi deleted; this is often called the “leave-one-out” density estimator. A disadvantage of this method is that it suffers from high variation.(b) Asymptotic Behavior of the MISE (AMISE) Method
An alternative parameter selector is to consider the asymptotic behavior of the MISE of the gamma kernel estimator. Chen [19] uses some aspects of the gamma distribution and a Taylor expansion to determine the MISE as follows:
(2.17)MISE(τ̂)=b2∫0∞{xτ′(x)+12xτ′′(x)}2dx+(2nbπ)-1∫0∞x-1/2τ(x)dx+O(n-1b-1/2+b2).
The asymptotic MISE disregards the last term; therefore, the optimal b  that minimizes the leading terms above is
(2.18)bAMISE=[(2π)-1∫0∞x-1/2τ(x)dx4n∫0∞(xτ′(x)+2-1xτ′′(x))2dx]2/5,
where the functions τ, τ′, and τ′′ are unknown. These quantities are obtained from the fitted gamma density with parameters adjusted from the sample. This solution still requires further study, but our paper shows promising results.


### 2.3.2. Zhang et al. [18] Estimator

Zhang et al. [18] submitted a model that works particularly well when τ(0)>0 and combines pseudodata creation, its transformation, and its reflection around the Y-axis in the following three steps.Step 1.
Transform the original dataX1,…,Xn to g(X1),…,g(Xn), while keeping the original data, where g is a nonnegative, continuous, and monotonically increasing function from [0,∞) to [0,∞).
Based on extensive simulations, the transformation that best suits a broad variety of densities is(2.19)g(x)=x+dx2+Adx3,
where A>1/3 and d=f′(0)/f(0).Step 2.
Reflect the pseudodata,g(X1),…,g(Xn), around the origin.Step 3.
Based on the enlarged data sample,-g(X1),…,-g(Xn),X1,…,Xn, define the new estimator as
(2.20)τ̂n(x,h)=1nh∑j=1n{K(x-Xjh)+K(x+g(Xj)h)},x≥0,
where h is a smoothing parameter and K is a symmetric probability function with support [-1,1] like the Epanechnikov kernel
(2.21)K(t)=34(1-t2)I[-1,1].
Notice that the transformation g defined above is not available in practice, because d is unknown. A good estimator can be obtained when d is written as (d/dx)logf(x)∣x=0,
(2.22)dn=log⁡fn(h)-log⁡fn(0)h,
where fn(h)=fn*(h)+1/n2, fn(0)=max(fn*(0),1/n2), fn*(h)=(1/nh)∑j=1nK((x-Xj)/h), fn*(0)=(1/nh0)∑j=1nK0(Xj/h0), and K0 is a so-called endpoint kernel, satisfying
(2.23)∫-10K0(t)dt=1,∫-10tK0(t)dt=0,∫-10t2K0(t)dt≠0,h0={[∫-11K(t)dt]2∫-10K0(t)2dt[∫-10t2K0(t)dt]2∫-11K(t)2dt}h.
Zhang et al. [18] proved that for t≥h, the effect of reflected pseudodata is insignificant, and the estimator can be reduced to the Parzen-Rosenblatt estimator
(2.24)τ̂n(t,h)=1nh∑j=1nK(t-Xjh).
They also established that
(2.25)∫0∞τ̂n(t)dt=1+1n∑i=1n∫-Xi/h-g(Xi)/hK(z)dz.
Thus, τ̂n(t) only integrates to 1 when dn=0, so gn(Xi)=Xi, or when Xi=0 for all Xi’s, because gn(0)=0. However, when n→∞, both limits of the second term will eventually converge to 0 and τ̂n(t) will integrate to 1 asymptotically.
Zhang and Karunamuni [24] used the endpoint kernel
(2.26)K0(t)=12(1+t)(12+t)I[-1,0]
and showed that this kernel minimizes the MSE when estimating τ(0). Therefore, h0=2h is approximately the optimal smoothing parameter for estimating τ(0) except when τ(0)=0.
Chiu [25] described a parameter selecting method that considered the optimal h that minimizes the asymptotic MISE when K is a symmetric probability function with up to the fourth moment being finite
(2.27)hAMISE=[∫K(x)2dxn(∫x2K(x)dx)2∫τ′′(x)2dx]1/5,
where the function ∫τ′′(x)2dx is unknown.
Chiu’s [25] “plug-in” method consists of estimating this quantity through the characteristic function of the sample
(2.28)φ̂(λ)=n-1∑j=1neiλXj,
and calculating the optimal h with the previous formula. The characteristic function of τ is
(2.29)φ(λ)=∫eiλxτ(x)dx.
By the inversion formula, we have
(2.30)τ(x)=(2π)-1∫eiλxφ(λ)dλ,
thus,
(2.31)τ′′(x)=(2π)-1∫λ2eiλxφ(λ)dλ,∫τ′′(x)2dx=∫[(2π)-1∫λ2eiλxφ(λ)dλ]2dx=∫λ4[(2π)-1∫eiλxφ(λ)dλ]2dx=∫λ4[τ(x)]2dx.
Using Parseval’s identity, we can show that
(2.32)∫λ4[τ(x)]2dx=(2π)-1∫λ4|φ(λ)|2dλ.
Chiu [25] introduced a cutoff value Λ for λ such that |φ̂(λ)|2<c/n. Extensive computational experiments from Bessegato et al. [26] show that c=3 is the value that minimizes the estimator variance. The final “plug in” estimator is then
(2.33)ĥAMISE=[∫K(x)2dxn(∫x2K(x)dx)2π-1∫0Λλ4[|φ̂(λ)|2-n-1]dλ]1/5.


## 2.3.1. Gamma Kernel Estimator

Recently, Chen [19] suggested an asymmetric kernel with naturally varying shape as a way to avoid allocating weight for negative values. The gamma kernel estimators are always nonnegative, free of boundary bias, and achieve the optimal rate of convergence for the mean square error (MSE) in the nonnegative kernel estimator class. Bouezmarni and Scaillet [21] showed that this estimator is consistent and able to avoid boundary bias. Let KG(p,q) be the gamma density function with parameters p (shape) and q (rate). The gamma kernel considered is(2.12)KG(tb+1,b)(Xj)=Xjt/be-Xj/bb(t/b)+1Γ[(t/b)+1],
where b is a smoothing parameter that satisfies the condition b→0, nb→∞ as n→∞. The gamma kernel estimator is(2.13)τ̂(t;b)=n-1∑j=1nKG(tb+1,b)(Xj).
The smoothing parameter b is critical for the overall performance of the kernel estimator considered. A small b leads to a relatively bumpy density, while a large one results in a smooth density. There are several methods to determine the best fit, from a minimization of the mean integrated squared error (MISE) of τ̂(t;b) to the asymptotic behavior of the MISE (AMISE).(a) Least Squares Cross-Validation (LSCV) Method
The least squares cross-validation (LSCV) method starts from the MISE expansion(2.14)MISE{τ̂(x;b)}=E∫τ̂(x;b)2dx-2E∫τ̂(x;b)τ(x)dx+E∫τ(x)2dx.
The minimization of the first term is equivalent to the minimization of
(2.15)MISE{τ̂(x;h)}-E∫τ(x)2dx=E[∫τ̂(x;h)2dx-2∫τ̂(x;h)τ(x)dx].
The right-hand side is unknown, because it depends on τ. However, an unbiased estimator for this quantity is
(2.16)LSCV(h)=∫τ̂(x;h)2dx-2n-1∑i=1nτ̂-i(Xi;h),
where τ̂-i=(Xi;h) is the density estimate based on the sample with Xi deleted; this is often called the “leave-one-out” density estimator. A disadvantage of this method is that it suffers from high variation.(b) Asymptotic Behavior of the MISE (AMISE) Method
An alternative parameter selector is to consider the asymptotic behavior of the MISE of the gamma kernel estimator. Chen [19] uses some aspects of the gamma distribution and a Taylor expansion to determine the MISE as follows:
(2.17)MISE(τ̂)=b2∫0∞{xτ′(x)+12xτ′′(x)}2dx+(2nbπ)-1∫0∞x-1/2τ(x)dx+O(n-1b-1/2+b2).
The asymptotic MISE disregards the last term; therefore, the optimal b  that minimizes the leading terms above is
(2.18)bAMISE=[(2π)-1∫0∞x-1/2τ(x)dx4n∫0∞(xτ′(x)+2-1xτ′′(x))2dx]2/5,
where the functions τ, τ′, and τ′′ are unknown. These quantities are obtained from the fitted gamma density with parameters adjusted from the sample. This solution still requires further study, but our paper shows promising results.


## 2.3.2. Zhang et al. [18] Estimator

Zhang et al. [18] submitted a model that works particularly well when τ(0)>0 and combines pseudodata creation, its transformation, and its reflection around the Y-axis in the following three steps.Step 1.
Transform the original dataX1,…,Xn to g(X1),…,g(Xn), while keeping the original data, where g is a nonnegative, continuous, and monotonically increasing function from [0,∞) to [0,∞).
Based on extensive simulations, the transformation that best suits a broad variety of densities is(2.19)g(x)=x+dx2+Adx3,
where A>1/3 and d=f′(0)/f(0).Step 2.
Reflect the pseudodata,g(X1),…,g(Xn), around the origin.Step 3.
Based on the enlarged data sample,-g(X1),…,-g(Xn),X1,…,Xn, define the new estimator as
(2.20)τ̂n(x,h)=1nh∑j=1n{K(x-Xjh)+K(x+g(Xj)h)},x≥0,
where h is a smoothing parameter and K is a symmetric probability function with support [-1,1] like the Epanechnikov kernel
(2.21)K(t)=34(1-t2)I[-1,1].
Notice that the transformation g defined above is not available in practice, because d is unknown. A good estimator can be obtained when d is written as (d/dx)logf(x)∣x=0,
(2.22)dn=log⁡fn(h)-log⁡fn(0)h,
where fn(h)=fn*(h)+1/n2, fn(0)=max(fn*(0),1/n2), fn*(h)=(1/nh)∑j=1nK((x-Xj)/h), fn*(0)=(1/nh0)∑j=1nK0(Xj/h0), and K0 is a so-called endpoint kernel, satisfying
(2.23)∫-10K0(t)dt=1,∫-10tK0(t)dt=0,∫-10t2K0(t)dt≠0,h0={[∫-11K(t)dt]2∫-10K0(t)2dt[∫-10t2K0(t)dt]2∫-11K(t)2dt}h.
Zhang et al. [18] proved that for t≥h, the effect of reflected pseudodata is insignificant, and the estimator can be reduced to the Parzen-Rosenblatt estimator
(2.24)τ̂n(t,h)=1nh∑j=1nK(t-Xjh).
They also established that
(2.25)∫0∞τ̂n(t)dt=1+1n∑i=1n∫-Xi/h-g(Xi)/hK(z)dz.
Thus, τ̂n(t) only integrates to 1 when dn=0, so gn(Xi)=Xi, or when Xi=0 for all Xi’s, because gn(0)=0. However, when n→∞, both limits of the second term will eventually converge to 0 and τ̂n(t) will integrate to 1 asymptotically.
Zhang and Karunamuni [24] used the endpoint kernel
(2.26)K0(t)=12(1+t)(12+t)I[-1,0]
and showed that this kernel minimizes the MSE when estimating τ(0). Therefore, h0=2h is approximately the optimal smoothing parameter for estimating τ(0) except when τ(0)=0.
Chiu [25] described a parameter selecting method that considered the optimal h that minimizes the asymptotic MISE when K is a symmetric probability function with up to the fourth moment being finite
(2.27)hAMISE=[∫K(x)2dxn(∫x2K(x)dx)2∫τ′′(x)2dx]1/5,
where the function ∫τ′′(x)2dx is unknown.
Chiu’s [25] “plug-in” method consists of estimating this quantity through the characteristic function of the sample
(2.28)φ̂(λ)=n-1∑j=1neiλXj,
and calculating the optimal h with the previous formula. The characteristic function of τ is
(2.29)φ(λ)=∫eiλxτ(x)dx.
By the inversion formula, we have
(2.30)τ(x)=(2π)-1∫eiλxφ(λ)dλ,
thus,
(2.31)τ′′(x)=(2π)-1∫λ2eiλxφ(λ)dλ,∫τ′′(x)2dx=∫[(2π)-1∫λ2eiλxφ(λ)dλ]2dx=∫λ4[(2π)-1∫eiλxφ(λ)dλ]2dx=∫λ4[τ(x)]2dx.
Using Parseval’s identity, we can show that
(2.32)∫λ4[τ(x)]2dx=(2π)-1∫λ4|φ(λ)|2dλ.
Chiu [25] introduced a cutoff value Λ for λ such that |φ̂(λ)|2<c/n. Extensive computational experiments from Bessegato et al. [26] show that c=3 is the value that minimizes the estimator variance. The final “plug in” estimator is then
(2.33)ĥAMISE=[∫K(x)2dxn(∫x2K(x)dx)2π-1∫0Λλ4[|φ̂(λ)|2-n-1]dλ]1/5.


## 3. Experimental Results

This section presents some results of simulations forGIX/M/c/N systems with partial blocks, where the interarrival time is estimated through the following kernel methods:(i)
gamma kernel estimator with LSCV method,(ii)
gamma kennel estimator with optimalbAMISE,(iii)
Zhang et al. [18] estimator with Chiu’s [25] “plug in” method.To evaluate the performance of the estimators above, we will compare the mean square error (MSE) of each estimated prearrival probability in the following way.Step 1.
Generate a sample of sizen of general interarrival distribution τ.Step 2.
Calculate the mean service rateμ=g̅/ρcE(τ).Step 3.
Estimate the optimal smoothing parameterh or b.Step 4.
Use the kernel density method to estimate the theoretical density functionτ(x).Step 5.
Find each estimated transition probability.Step 6.
Solve the linear system:(3.1)[(p̂0,0-1)p̂1,0⋯p̂N,0p̂0,1(p̂1,1-1)⋮⋮⋱⋮p̂0,N-1⋯⋯(p̂N,N-1-1)1⋯⋯1][π̂0π̂1⋮⋮π̂N]=[00⋮01].
The algorithm above was coded in R 2.8.0 (or earlier versions, see [27]) and is available upon request from the authors for educational and research purposes.
The experiments are based on random samples of interarrival times of sizen=100. Slightly larger and smaller samples were also tested, but the results (not shown) are similar. The theoretical interarrival distributions considered in this experiment were(i)
Weibull distribution with shape = 2 and rate = 20,(ii)
gamma distribution with shape = 10 and rate = 2,(iii)
gamma mixture distribution of 0.45 × gamma (5; 2) + 0.55 × gamma (30; 1).The Weibull density has τ′(0)≠0, the pure gamma has τ′(0)=0, and the gamma mixture is bimodal. We choose Weibull and gamma distributions because of their well-known flexibility for modeling real databases. For simplicity, the group size X is constant and equal to 1. Two different numbers of servers were considered, c=5 and 10, with maximum capacities of N=20 and 25, respectively, (which result in buffers of fixed sizes equal to 15). Figures 1 and 2 show the results, where πk is the theoretical distribution of the number of clients and MSE(π̂k) is the mean square error of the estimates.πk and MSE(π̂k) for c=5.
(a)
Weibull (2; 20)(b)
Gamma (10; 2)(c)
0.45 Gamma (5; 2) + 0.55 Gamma (30; 1)πk and MSE(π̂k) for c=10.
(a)
Weibull (2; 20)(b)
Gamma (10; 2)(c)
0.45 Gamma (5; 2) + 0.55 Gamma (30; 1)


## 4. Discussion

In general, the errors decreased as the number of serversc increased, and they were dependent on the theoretical distribution considered. We obtained the largest errors for the last distribution (gamma mixture). The LSCV method is better than bAMISE on the first and last distributions. Zhang et al.’s [18] estimator has the worst performance on the first and a competitive performance on the second distribution (gamma distribution). All the estimators performed well with the second distribution. The last distribution (gamma mixture) has a particular behavior. Although the errors are large for most of the estimates, the error is small for the blocking probability (i.e., π̂20, for the system with c=5, Figure 1, and π̂25, for c=10, Figure 2). The blocking probability is an important performance measure, because it indicates the fraction of costumers lost by the system.Finally, simulations showed that whenτ′(0)≠0 the gamma kernel method had the best performance. This suggests the Zhang et al. [18] method does not work well when τ(0)>0. At the same time, its behavior with the bimodal density showed a very low MSE for probabilities near the maximum state. This implies a good estimation of the blocking probability and other performance measures.The method used to select the smoothing parameter for the gamma kernel estimator had no effect. A better selector would have the function of its AMISE optimal parameter estimated like the “plug-in” method for symmetric kernels. The combined use of smoothing parameter selection and Bayesian techniques is promising (see de Lima and Atuncar [3]).


## 5. Application to a Real Call Center Data Set

We analyzed a real database of 7,761 phone calls made on April 7, 2006, from 8:00 am until 1:00 pm to a call center (data available from the authors upon request). We worked on a set of interarrival and service times to find the required minimum size of the server facility taking into account that the arrival rates are not homogeneous. In fact, larger interarrival times were detected toward the beginning of the day rather than throughout the day. This will lead to the adjustment of a different system for each hour of service that had an approximately homogeneous arrival rate. Using this data, we will model aGIX/M/c/N system with a kernel interarrival density estimate for each hour of service.The data set is presented in seconds, which is the precision given by the data acquisition system. This constraint will lead to ties in the arrival times (i.e., more than one call can arrive in the same second). Therefore, we treated same-time calls as part of a single arrival group and fit a discrete probability distribution. Table1 shows the observed frequency of group sizes X and the fitted distribution D by hour. These arrivals were observed in each second.Table 1
Observed frequency of group sizesX and fitted distributions D.
Time:8 am-9 am9 am-10 am10 am-11 am11 am-12 am12 am-1 pm1 arrival0.9150.8080.7940.8820.9162 arrivals0.0850.1580.1720.1080.0753 arrivals0.0000.0310.0300.0090.0094 arrivals0.0000.0030.0040.0010.000D~Poisson (0.082)Geom. (0.813)Geom. (0.804)Geom. (0.885)Geom. (0.914)We used the gamma kennel-estimator with optimalbAMISE as the kernel density estimation method. This method gave estimated densities bounded at τ=0 and illustrates the information we miss when data is rounded. We also considered the gamma kernel estimator with LSCV method, but it did not behave as well with discrete data. Table 2 shows the smoothing parameter calculated, and Figure 3 shows the estimated densities.Table 2
AMISE smoothing parameter for gamma kernel.
Time:8 am-9 am9 am-10 am10 am-11 am11 am-12 am12 am-1 pmbAMISE0.15090.05050.04760.06250.0932Estimated interarrival times by hour.
(a)
8 am-9 am(b)
9 am-10 am(c)
10 am-11 am(d)
11 am-12 am(e)
12 am-1 pmService time distribution fits an exponential distribution with parameterλ=0.003339 and standard error of 3.8×10-5. The behavior of the system then becomes restricted to the choice of the number of servers c and capacity of the system N. In Table 3, the minimum number of servers required is set to maintain the stability of the system at each hour considered.Table 3
Minimum number of servers.
Time:8 am-9 am9 am-10 am10 am-11 am11 am-12 am12 am-1 pmmin{c}19971106440To optimize system performance, it is necessary to establish some criterion function for the adopted design. If we only consider the effective arrival rate(λEF=λ[1-PBL]), we can describe this performance measure according to the number of servers and the maximum queue size. Figure 4 shows this relationship from 8 to 9 am. We saw similar results (not shown) for different periods. A good criterion would be the C=(C1/λEF+C2c+C3N), where C1, C2, and C3 are costs related to each parameter and defined by the system environment. As an example, if we have C1=500, C2=1 and C3=1, the state with minimal criterion would be c=19 and N=20. If cost C1 is raised to 600, the best configuration would be c=21 and N=21.Figure 4
Effective arrival rate as a function ofN and c.Figure5 shows the invariant distributions estimated for other periods with different queue and server facility sizes. The method seems helpful for adjusting the maximum size N, and it is clear that the period plays a key role.Estimated invariant distributions.
(a)
9 am-10 am(b)
10 am-11 am(c)
11 am-12 am(d)
12 am-1 pm


## 6. Conclusions

The main focus of this paper is on using kernels in the analysis of queuing systems. We studied the adequacy of the kernel estimator methods for calculating the invariant probability distribution and performance measures of queueing systems that have general interarrival distribution times with bulk arrivals. Discussions in simulations showed that the method is effective. Also, we successfully applied the method to the calls from a call center to plan and improve the performance of these important queueing systems.Future research may take other directions—for instance, determining how the estimate of the prearrival invariant distribution moves away from its real value. One approach is from the variance of the bias and variance of each probability estimated. For example, the variance of an estimate of transition probabilitypij could depend on the variance of term β̂0, so we would need to find(6.1)var(β̂0)=var(∫0∞e-cμzτ̂(z,h)dz).
This research can be applied to several areas of practical interest, including health and industry. For example, it could be critical in a medical emergency room that must optimize resource allocation.Another possible future direction for research in this field involves the development of models to deal with dependency for the interarrival times. Such a case is particularly true in computer networks. Thus, in computer communications networks, people are interested in finding proper bounds of arrivals (see Li and Zhao [28]), with the assumption that arrivals are dependent (see Li and Zhao [29]), which are recently published results based on the so-called deterministic queuing theory (for details, see Le Boudec and Thiran [30]).



---

*Source: 348262-2011-10-26.xml*