Abstract
We show the interest of nonparametric methods taking into account the boundary correction techniques for a numerical evaluation of an approximation error between the stationary distributions of and queueing systems, when the density function of the general arrivals law in the system is unknown and defined on a bounded support. To compute this error, we use two kinds of norms: the norm and the weight norm. Numerical examples based on simulation studies are presented for the two cases of considered norms. A comparative study of the results has been provided.
1. Introduction
When modeling practical problems, one may often replace a real system by another one which is close to it in some sense but simpler in structure and/or components. This approximation is necessary because real systems are generally very complicated, so their analysis cannot lead to analytical results or it leads to complicated results which are not useful in practice.
To overcome the difficulties encountered in obtaining exact and interpretable solutions for many queueing systems, analysts use approximation methods. Use of these methods allows approaching the characteristics of a complex model by those of a simpler one. It is interesting in this case to measure the resulting approximation error.
One of these approximation methods is the strong stability [1, 2]. This technique is applicable to all operations research models which can be represented by a Markov chain. According to this approach, we suppose that the perturbation is small with respect to a certain operators norm (weight norm). Such a strict condition allows us to obtain better estimations on the characteristics of the perturbed chain, for instance, the perturbed stationary distributions.
In this paper, we focus on the evaluation of the approximation error between the stationary distributions of and systems, when the density function of the general arrivals law in the system is unknown and defined on a bounded support and must first be assessed by an appropriate nonparametric method [3–6]. To determine this error, different norms are used, namely, the norm defined in [7] and the weight norm of the approximation strong stability method [1, 2].
Moreover, as the strong stability method assumes that the perturbation is small, then we suppose that the arrivals law of the system is close to the Poisson one with parameter . This conducts us to consider the problem of boundary bias correction when performing nonparametric estimation of the unknown density of the law , since the density function of the exponential law is defined on the positive real line [3, 6].
On the other hand, in practice, we are often more interested in the deviation between the average characteristics (e.g., the mean waiting time) of the ideal model and the perturbed one than in the difference between stationary probabilities. Indeed, in most of the cases the model is used for the calculation of distributions of waiting times. We may be able to compare the resulting distributions of waiting times (which is the most common goal of the model) given by solving the queue.
This paper is organized as follows: in Section 2, we present the two norms under consideration. Some techniques for the correction of boundary effects used in kernel density estimation are discussed in Section 3. The main results of this work are presented in Section 4 which is devoted to the numerical evaluation of the approximation error on the stationary distributions of the two considered systems, by combining between a norm (approximate method) and nonparametric methods. Illustrative numerical examples based on simulation results are exposed, allowing effectuating a comparative study.
2. Different Considered Norms
2.1. Norm
To determine the proximity of the stationary distributions and of two queueing systems ( and systems, resp.), Pedrono and Hellary [7] have defined the following metric ( norm): In our case, is the stationary distribution of the system given by , , where is the traffic intensity of the system, is the mean rate of the interarrival times, and is the mean service time. is the stationary distribution of the system given by , , where is the unique solution (found numerically by the fixed point method) of the system: and represents the density function of the general distribution .
2.2. Operators Norm of the Strong Stability Method (Weight Norm)
2.2.1. Preliminaries and Notations
In this section, we introduce some necessary notations appropriate for our case study and recall the basic definition of the strong stability method. For a general framework, see [1, 2].
Consider the measurable space , where is the -algebra generated by all singletons , .
Let be the space of finite measures on and the space of bounded measurable functions on . We associate with each transition kernel the linear mapping . Introduce on the class of norms (weight norms) of the form where is an arbitrary measurable function (not necessarily finite) bounded below away from a positive constant. This norm induces in the space the norm .
Let us consider , the space of linear operators on the space , with norm .
Let and be two invariant measures and suppose that these measures have finite -norm. Then
2.2.2. Strong Stability Criterion
Definition 1 (see [1, 2]). A Markov chain with a transition kernel and an invariant measure is said to be -strongly stable with respect to the norm defined in (4), if and each stochastic kernel on the space in some neighborhood has a unique invariant measure and as in this neighborhood.
2.2.3. Strong Stability of System after Perturbation of the Arrival Flow
(a) Description of the Models. Let us consider a system where interarrival times are independently distributed with general distribution and service times are distributed with (exponential with parameter ).
Let be the number of customers left behind in the system by the th departure. It is easy to prove that forms a Markov chain [2] with a transition operator , where Consider also an system, which has Poisson arrivals with parameter and the same distribution of the service times as the preceding system. It is known that (the number of customers left behind in the system by the th departure) forms a Markov chain with a transition operator , where Suppose that the arrival flow of the system is close to the Poisson one. This proximity is then characterized by the metric where is the variation of the measure and and are the respective density functions of the distributions and . The stationary distributions of the states of and are defined as follows:
(b) Estimation of the Strong Stability
Theorem 2 (see [8]). Suppose that the traffic intensity of the system is smaller than 1. Therefore, for all such that , the imbedded Markov chain is -strongly stable, after a small perturbation of the interarrival times, for .
The margin between the transition operators is given by .
In addition, if , we have the error on the stationary distributions as follows:
where is defined in (8), and are defined in (9), and .
Given the error defined in (1) and the bound in formula (10) of Theorem 2, it remains to be computed the solution of the system (2) and the variation distance expressed in (8), which both involve the density function of the unknown general distribution . Statistical methods to do so will be discussed in the following.
3. Different Nonparametric Considered Estimates
Let be a sample coming from a random variable of density function and distribution . The Parzen-Rosenblatt classical kernel estimate (see [4, 5]) of the density for each point is given by where is a symmetric density function called kernel and is the smoothing parameter (or bandwidth).
In practice, the critical step in the kernel density estimation is the choice of the bandwidth which controls the smoothness of the kernel estimate (11). This problem has been widely studied and many methods have been proposed (see the review paper [9] and the monograph of Silverman [10]).
Several results are known in the literature when the density function is defined on the real line [4, 10]. In the case of a density function defined on a bounded support, the boundary effects are present since there is a bias near the border [3, 6, 10]. This problem is due to usage of a fixed kernel which assigns a weight outside the support when smoothing is carried out near the boundary. To resolve this problem, many recent methods have been elaborated.
Schuster [6] suggests creating the mirror image of the data in the other side of the boundary and then applying the estimator (11) for the set of the initial data and their reflection. is then estimated, for , as follows: Another simple idea to avoid the problem of boundary effects is the use of a flexible kernel, which never assigns a weight out of the support of the density function and which corrects automatically and implicitly the boundary effects. We can cite the asymmetric kernels [3] given by the following form: where is the bandwidth and the asymmetric kernel can be taken as a Gamma density with parameters given by
4. Numerical Examples: Simulation
Recall that our object is the numerical evaluation of the approximation error between the stationary distributions of and systems, when the density function of the general arrival law in the system is unknown and must be estimated by nonparametric methods. In this section, we focus on developing algorithms that take into account these nonparametric methods and which are adapted to the case of the two norms mentioned in Section 2. Indeed, to precise the proximity error between the two systems according to the norm (resp., to the weight norm), we give Algorithm 1 (resp., Algorithm 2; see [11] for the general steps), detailed in Section 4.1 below. Numerical illustrative examples based on simulation under Matlab 7.1 environment are presented in Section 4.2.
| ||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||
4.1. Algorithms
4.2. Numerical Examples
4.2.1. Approximation Error with Norm
Example 1. We consider the following four cases.
First Case. We consider a system such that the density function of the general law is given by
By generating a sample coming from the general law with the same density defined above, we use the kernel density method to estimate the theoretical density by using the different estimators given in the three following cases.
Second Case. We use the kernel estimate defined in (11).
Third Case. We use the kernel estimate defined in (12).
Fourth Case. We use the kernel estimate defined in (13) with the Gamma kernel given in (14).
For the last three cases, we take the sample size and the number of simulations . In all the cases, we introduce the service mean time: .
We first determine the interarrival mean rate . We get some performance measures of the system: (i)interarrival mean rate: ;(ii)traffic intensity: .
For calculation of the proximity error of the and systems, defined in (1), we use Algorithm 1. We obtain the results of Table 1.
Discussion. We note in Table 1 that the approximation errors found by the kernel method considering the correction of boundary effects ( in case of mirror image estimate and in case of asymmetric Gamma kernel estimate) are better than the one found by the classical kernel method without taking into consideration the correction of boundary effects (). The kernel density method is so a good tool for the numerical evaluation of the proximity of the stationary distributions of two systems to approximate, in order to have an idea of the proximity of the characteristics of these systems and the possibility of substituting ones by others.
Remark 3. In most of the cases the model is used for the calculation of distributions of waiting times. It can be solved numerically, by solving (2). The probability density function (PDF) of the waiting time in the system is given by (see [12])
The same formula is valid to determine the PDF of the waiting time in the system, with replacing by ( is the unique solution of system (2)).
Using results of Table 1 together with (16), we may compare the resulting distributions of waiting times (which is the most common goal of the model) given by solving the queue. For the system, we use the theoretical density and its different estimates . An illustration is given in Figure 1.
Note that according to Figure 1, curves of waiting times distributions for the system obtained using the theoretical density and its estimates and are close to that of the waiting times distribution for the system, contrary to the curve of the waiting times distribution for the system obtained using the classical kernel estimate which is far away from that of the waiting times distribution for the system. Once more, this is due to the boundary effects caused by using a symmetric fixed kernel estimate.
4.2.2. Approximation Error with the Weight Norm
To realize this work, we use Algorithm 2. The Epanechnikov kernel [10] is used throughout for estimators involving symmetric kernels. The bandwidth is chosen to minimize the criterion of the “least squares cross-validation” [10]. The smoothing parameters and are chosen according to a bandwidth selection method which leads to an asymptotically optimal window in the sense of minimizing distance [13].
Example 2. We generate samples of size issued from different laws. We take the number of simulations . For each case of law, we replace the nonparametric density (defined in Algorithm 2) by the density function found by applying the Parzen-Rosenblatt classical kernel estimate defined in (11) to estimate the theoretical density of each sample. Using Algorithm 2, we obtain the results of Table 2.
Discussion. First, notice that according to Table 2, it seems that the proposed method approximates certain values with some differences, for instance, the with the . This may be explained by the error done when replacing the theoretical density by the nonparametric classical kernel estimate in the formula used to compute the interarrival mean rate ; that is, . Add to this the errors committed when using approximative numerical methods in the programming process (e.g., the trapezes method used for computing the integral in the above formula).
Note also that according to Table 2, the application of the Parzen-Rosenblatt classical kernel estimate (11) for the approximation of the system by the one when using the strong stability method does not allow us to determine the error on the stationary distributions between the two systems. This is due to the importance of the value of the variation distance (e.g., for and for Weibull(2,0.5,0)). Therefore, the kernel density method applied for the study of the strong stability of the queueing systems (e.g., the system) does not precise the proximity error of the laws of the two systems (real and ideal) but affirms and reinforces the order of the importance of the smallness of the perturbation done in the study of the strong stability of the systems.
Example 3. Following the previous example, the classical kernel estimate (Parzen-Rosenblatt estimate) has shown its insufficiency for determining the approximation error on the stationary distributions of the corresponding systems. It is why we consider again in this example the study of this problem by using the same hyperexponential law defined in (15) and by applying the classical kernel estimate and the other boundary correction techniques for an eventual comparison. We follow the same process used in the simulation study in Section 4.2.1. We use again Algorithm 2 and we list in Table 3 the variation distance and the approximation error on the stationary distributions for the different estimates. Figure 2 describes the evolution of the error in function of .
(a)
(b)
Discussion. We note in Table 3 that the approximation error on the stationary distributions of the and systems was given when applying the kernel density method by considering the correction of boundary effects in the case of using mirror image estimate () or in the case of using the asymmetric Gamma kernel estimate (). But when applying the kernel density method without taking into consideration the correction of boundary effects in the case of using the Parzen-Rosenblatt classical kernel estimate, the approximation error () on the stationary distributions of the quoted systems could not be given.
Notice following Figure 2(a) (resp., Figure 2(b)) that the error, being important at the start, decreases speedily for the values of in the neighborhood of the lower bound () (resp., ). This may be explained by the fact that they are at the boundary of the stability domain (critical region). We notice also that the error increases speedily in the neighborhood of the upper bound () (resp., ) (critical region). In contrast, everywhere else, the error increases reasonably with the values of (favorable region). Nevertheless, it would be necessary to consider the minimal error which corresponds in our case to , (resp., , ). Results obtained by considering the asymmetric Gamma kernel estimate meet up with those obtained when using the real theoretical density.
Remark 4. In practice, we are often more interested in the deviation between the average characteristics (e.g., the mean queue length) of the nominal chain and the perturbed one than in the difference between stationary probabilities. For this purpose, we give Corollary 5 that allows us to predict the characteristics deviation of both systems using the results of Theorem 2 and that for an appropriate choice of the individual performance measure .
Corollary 5. Suppose that the assumptions of Theorem 2 hold, then, for any function such that , one has where is defined in (10).
Proof. The proof follows from Theorem 2 together with the fact that for any measure it holds that (see formula (5)).
Example for the Choice of the Individual Performance Measure
Example 4 (case: mean number of the waiting customers in the system). Let and be the mean numbers of the waiting customers of the nominal (ideal) system and the perturbed one, respectively. To predict the deviation , it suffices to take the individual performance measure as follows: We have ); then, for defined in (18), it is easy to show that ; thus formula (17) can be written as follows:
By dividing expression (19) by and by using Little’s formulas, we get This last expression provides a bound for the deviation between the mean waiting times of the and systems.
Applying this formula and results of Table 3, we give in Table 4 some values for the error on the mean waiting times of the considered systems, using the theoretical density defined in (15) and its asymmetric Gamma kernel estimate for the system. We use values of minimizing the error on the stationary distributions defined in (10).
5. Conclusion
By comparing the results of Tables 1 and 3, we note that the conclusions converge for the two cases of considered norms in the sense that the classical kernel method has shortcomings and boundary correction techniques are more appropriate in our case study. However, we note that the error (according to weight norm) is quite large compared to the error (according to norm). This is firstly due to the difference between the two norms. In addition, during the application of the strong stability method, we add to the error committed by the nonparametric estimation the error resulting from the perturbation performed in this case.
Again, following the results of Figure 1 and Table 4, the same finding is observed concerning the waiting time characteristic.
To summarize, the comparative study between the results obtained by applying the different nonparametric methods to both considered norms and for some specific characteristics shows the impact and interest of those that take into account the correction of boundary effects to determine an approximation error between the considered systems ( and ).
Systems used in this paper are relatively simple. They serve more as an illustrative support for a good comprehension of the techniques used to solve the posed problem. It would be interesting to consider the results of this work for the approximation of more complex systems, such as the system, risk and inventory models, and networks.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.