# Interest of Boundary Kernel Density Techniques in Evaluating an Approximation Error of Queueing Systems Characteristics

Authors: Aïcha Bareche; Djamil Aïssani
Journal: International Journal of Mathematics and Mathematical Sciences
 (2014)

---

## 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 ofG
/
M
/
1 and M
/
M
/
1 queueing systems, when the density function of the general arrivals law G in the G
/
M
/
1 system is unknown and defined on a bounded support. To compute this error, we use two kinds of norms: the L
1 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.


---

## Body

## 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 ofG
/
M
/
1 and M
/
M
/
1 systems, when the density function of the general arrivals law G in the G
/
M
/
1 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 L
1 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 theG
/
M
/
1 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 G, 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 theG
/
M
/
1 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 M
/
M
/
1 queue.This paper is organized as follows: in Section2, 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.L
1 Norm

To determine the proximity of the stationary distributionsα
i and π
i of two queueing systems (G
/
M
/
1 and M
/
M
/
1 systems, resp.), Pedrono and Hellary [7] have defined the following metric (L
1 norm):
(1)
ɛ
=
max
⁡
i
≥
0
|
π
i
-
α
i
|
.
In our case, π
i is the stationary distribution of the M
/
M
/
1 system given by π
i
=
(
1
-
ρ
)
ρ
i, i
=
0,1
,
…, where ρ
=
λ
/
γ is the traffic intensity of the M
/
M
/
1 system, λ is the mean rate of the interarrival times, and γ is the mean service time. α
i is the stationary distribution of the G
/
M
/
1 system given by α
i
=
(
1
-
x
)
x
i, ∀
i
≥
0, where x is the unique solution (found numerically by the fixed point method) of the system:
(2)
x
=
∫
‍
e
-
γ
t
(
1
-
x
)
g
(
t
)
d
t
,
and g represents the density function of the general distribution G.


### 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(
N
,
B
(
N
)
), where B
(
N
) is the σ-algebra generated by all singletons {
j
}, j
∈
N.LetM
=
{
μ
j
} be the space of finite measures on B
(
N
) and N
=
{
f
(
j
)
} the space of bounded measurable functions on N. We associate with each transition kernel P the linear mapping (
μ
P
)
k
=
∑
j
≥
0
μ
j
P
j
k.
(3)
(
P
f
)
(
k
)
=
∑
i
≥
0
‍
f
(
i
)
P
k
i
.
Introduce on M the class of norms (weight norms) of the form
(4)
∥
μ
∥
υ
=
∑
j
≥
0
‍
υ
(
j
)
|
μ
j
|
,
where υ is an arbitrary measurable function (not necessarily finite) bounded below away from a positive constant. This norm induces in the space N the norm ∥
f
∥
υ
=
sup
⁡
k
≥
0
⁡
(
|
f
(
k
)
|
/
υ
(
k
)
).Let us considerB, the space of linear operators on the space {
μ
∈
M
:
∥
μ
∥
υ
<
∞
}, with norm ∥
Q
∥
υ
=
sup
⁡
k
≥
0
(
1
/
υ
(
k
)
)
∑
j
≥
0
υ
(
j
)
|
Q
k
j
|.Letα and μ be two invariant measures and suppose that these measures have finite υ-norm. Then
(5)
|
α
f
-
μ
f
|
≤
∥
α
-
μ
∥
υ
∥
f
∥
υ
inf
⁡
k
≥
0
υ
(
k
)
.


#### 2.2.2. Strong Stability Criterion

Definition 1 (see [1, 2]).
A Markov chainX with a transition kernel P and an invariant measure π is said to be υ-strongly stable with respect to the norm ∥
·
∥
υ defined in (4), if ∥
P
∥
υ
<
∞ and each stochastic kernel Q on the space (
N
,
B
(
N
)
) in some neighborhood {
Q
:
∥
Q
-
P
∥
υ
<
ϵ
} has a unique invariant measure μ
=
μ
(
Q
) and ∥
π
-
μ
∥
υ
→
0 as ∥
Q
-
P
∥
υ
→
0 in this neighborhood.


#### 2.2.3. Strong Stability ofM
/
M
/
1 System after Perturbation of the Arrival Flow

(a) Description of the Models. Let us consider a G
/
M
/
1
(
F
I
F
O
,
∞
) system where interarrival times are independently distributed with general distribution G
(
t
) and service times are distributed with E
γ
(
t
) (exponential with parameter γ).LetX
n
* be the number of customers left behind in the system by the nth departure. It is easy to prove that X
n
* forms a Markov chain [2] with a transition operator P
*
=
(
P
i
j
*
)
i
,
j
≥
0, where
(6)
P
i
j
*
=
{
d
i
+
1
-
j
*
=
∫
0
+
∞
‍
1
(
i
+
1
-
j
)
!
d
i
+
1
-
j
*
∫
0
+
∞
‍
×
e
-
γ
t
(
γ
t
)
i
+
1
-
j
d
G
(
t
)
,
if
1
≤
j
≤
i
+
1
,
1
-
∑
k
=
0
i
‍
d
k
*
,
if
j
=
0
,
0
,
otherwise
.
Consider also an M
/
M
/
1
(
F
I
F
O
,
∞
) system, which has Poisson arrivals with parameter λ and the same distribution of the service times as the preceding system. It is known that X
n (the number of customers left behind in the system by the nth departure) forms a Markov chain with a transition operator P
=
(
P
i
j
)
i
,
j
≥
0, where
(7)
P
i
j
=
{
d
i
+
1
-
j
=
λ
γ
i
+
1
-
j
(
λ
+
γ
)
i
+
2
-
j
,
if
1
≤
j
≤
i
+
1
,
1
-
∑
k
=
0
i
‍
d
k
=
(
γ
γ
+
λ
)
i
,
if
j
=
0
,
0
,
otherwise
.
Suppose that the arrival flow of the G
/
M
/
1 system is close to the Poisson one. This proximity is then characterized by the metric
(8)
w
=
w
(
G
,
E
λ
)
=
∫
0
∞
‍
|
G
-
E
λ
|
(
d
t
)
=
∫
0
∞
‍
|
g
-
e
λ
|
(
t
)
(
d
t
)
,
where |
a
| is the variation of the measure a and g and e
λ are the respective density functions of the distributions G and E
λ. The stationary distributions of the states of X
n
* and X
n are defined as follows:
(9)
α
k
=
lim
⁡
n
→
∞
Pr
(
X
n
*
=
k
)
,
k
=
0,1
,
2
,
…
,
π
k
=
lim
⁡
n
→
∞
Pr
(
X
n
=
k
)
,
k
=
0,1
,
2
,
…
,(b) Estimation of the Strong StabilityTheorem 2 (see [8]).
Suppose that the traffic intensity(
λ
/
γ
) of the M
/
M
/
1 system is smaller than 1. Therefore, for all β such that 1
<
β
<
γ
/
λ, the imbedded Markov chain X
n is υ-strongly stable, after a small perturbation of the interarrival times, for υ
(
k
)
=
β
k.
The margin between the transition operators is given by∥
P
*
-
P
∥
v
≤
(
1
+
β
)
w.
In addition, ifw
<
(
1
-
ρ
)
(
γ
-
λ
β
)
/
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
), we have the error on the stationary distributions as follows:
(10)
∥
α
-
π
∥
υ
≤
(
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
)
(
γ
-
λ
)
w
)
×
(
(
β
-
1
)
(
γ
-
λ
β
)
3
(
β
-
1
)
γ
+
λ
β
m
×
-
(
2
γ
-
λ
(
1
+
β
)
)
(
1
+
β
)
(
γ
-
λ
β
)
w
(
β
-
1
)
(
γ
-
λ
β
)
3
(
β
-
1
)
γ
+
λ
β
)
-
1
=
E
r
r
,
where w 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 x of the system (2) and the variation distance w expressed in (8), which both involve the density function of the unknown general distribution G. Statistical methods to do so will be discussed in the following.


## 2.1.L
1 Norm

To determine the proximity of the stationary distributionsα
i and π
i of two queueing systems (G
/
M
/
1 and M
/
M
/
1 systems, resp.), Pedrono and Hellary [7] have defined the following metric (L
1 norm):
(1)
ɛ
=
max
⁡
i
≥
0
|
π
i
-
α
i
|
.
In our case, π
i is the stationary distribution of the M
/
M
/
1 system given by π
i
=
(
1
-
ρ
)
ρ
i, i
=
0,1
,
…, where ρ
=
λ
/
γ is the traffic intensity of the M
/
M
/
1 system, λ is the mean rate of the interarrival times, and γ is the mean service time. α
i is the stationary distribution of the G
/
M
/
1 system given by α
i
=
(
1
-
x
)
x
i, ∀
i
≥
0, where x is the unique solution (found numerically by the fixed point method) of the system:
(2)
x
=
∫
‍
e
-
γ
t
(
1
-
x
)
g
(
t
)
d
t
,
and g represents the density function of the general distribution G.


## 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(
N
,
B
(
N
)
), where B
(
N
) is the σ-algebra generated by all singletons {
j
}, j
∈
N.LetM
=
{
μ
j
} be the space of finite measures on B
(
N
) and N
=
{
f
(
j
)
} the space of bounded measurable functions on N. We associate with each transition kernel P the linear mapping (
μ
P
)
k
=
∑
j
≥
0
μ
j
P
j
k.
(3)
(
P
f
)
(
k
)
=
∑
i
≥
0
‍
f
(
i
)
P
k
i
.
Introduce on M the class of norms (weight norms) of the form
(4)
∥
μ
∥
υ
=
∑
j
≥
0
‍
υ
(
j
)
|
μ
j
|
,
where υ is an arbitrary measurable function (not necessarily finite) bounded below away from a positive constant. This norm induces in the space N the norm ∥
f
∥
υ
=
sup
⁡
k
≥
0
⁡
(
|
f
(
k
)
|
/
υ
(
k
)
).Let us considerB, the space of linear operators on the space {
μ
∈
M
:
∥
μ
∥
υ
<
∞
}, with norm ∥
Q
∥
υ
=
sup
⁡
k
≥
0
(
1
/
υ
(
k
)
)
∑
j
≥
0
υ
(
j
)
|
Q
k
j
|.Letα and μ be two invariant measures and suppose that these measures have finite υ-norm. Then
(5)
|
α
f
-
μ
f
|
≤
∥
α
-
μ
∥
υ
∥
f
∥
υ
inf
⁡
k
≥
0
υ
(
k
)
.


### 2.2.2. Strong Stability Criterion

Definition 1 (see [1, 2]).
A Markov chainX with a transition kernel P and an invariant measure π is said to be υ-strongly stable with respect to the norm ∥
·
∥
υ defined in (4), if ∥
P
∥
υ
<
∞ and each stochastic kernel Q on the space (
N
,
B
(
N
)
) in some neighborhood {
Q
:
∥
Q
-
P
∥
υ
<
ϵ
} has a unique invariant measure μ
=
μ
(
Q
) and ∥
π
-
μ
∥
υ
→
0 as ∥
Q
-
P
∥
υ
→
0 in this neighborhood.


### 2.2.3. Strong Stability ofM
/
M
/
1 System after Perturbation of the Arrival Flow

(a) Description of the Models. Let us consider a G
/
M
/
1
(
F
I
F
O
,
∞
) system where interarrival times are independently distributed with general distribution G
(
t
) and service times are distributed with E
γ
(
t
) (exponential with parameter γ).LetX
n
* be the number of customers left behind in the system by the nth departure. It is easy to prove that X
n
* forms a Markov chain [2] with a transition operator P
*
=
(
P
i
j
*
)
i
,
j
≥
0, where
(6)
P
i
j
*
=
{
d
i
+
1
-
j
*
=
∫
0
+
∞
‍
1
(
i
+
1
-
j
)
!
d
i
+
1
-
j
*
∫
0
+
∞
‍
×
e
-
γ
t
(
γ
t
)
i
+
1
-
j
d
G
(
t
)
,
if
1
≤
j
≤
i
+
1
,
1
-
∑
k
=
0
i
‍
d
k
*
,
if
j
=
0
,
0
,
otherwise
.
Consider also an M
/
M
/
1
(
F
I
F
O
,
∞
) system, which has Poisson arrivals with parameter λ and the same distribution of the service times as the preceding system. It is known that X
n (the number of customers left behind in the system by the nth departure) forms a Markov chain with a transition operator P
=
(
P
i
j
)
i
,
j
≥
0, where
(7)
P
i
j
=
{
d
i
+
1
-
j
=
λ
γ
i
+
1
-
j
(
λ
+
γ
)
i
+
2
-
j
,
if
1
≤
j
≤
i
+
1
,
1
-
∑
k
=
0
i
‍
d
k
=
(
γ
γ
+
λ
)
i
,
if
j
=
0
,
0
,
otherwise
.
Suppose that the arrival flow of the G
/
M
/
1 system is close to the Poisson one. This proximity is then characterized by the metric
(8)
w
=
w
(
G
,
E
λ
)
=
∫
0
∞
‍
|
G
-
E
λ
|
(
d
t
)
=
∫
0
∞
‍
|
g
-
e
λ
|
(
t
)
(
d
t
)
,
where |
a
| is the variation of the measure a and g and e
λ are the respective density functions of the distributions G and E
λ. The stationary distributions of the states of X
n
* and X
n are defined as follows:
(9)
α
k
=
lim
⁡
n
→
∞
Pr
(
X
n
*
=
k
)
,
k
=
0,1
,
2
,
…
,
π
k
=
lim
⁡
n
→
∞
Pr
(
X
n
=
k
)
,
k
=
0,1
,
2
,
…
,(b) Estimation of the Strong StabilityTheorem 2 (see [8]).
Suppose that the traffic intensity(
λ
/
γ
) of the M
/
M
/
1 system is smaller than 1. Therefore, for all β such that 1
<
β
<
γ
/
λ, the imbedded Markov chain X
n is υ-strongly stable, after a small perturbation of the interarrival times, for υ
(
k
)
=
β
k.
The margin between the transition operators is given by∥
P
*
-
P
∥
v
≤
(
1
+
β
)
w.
In addition, ifw
<
(
1
-
ρ
)
(
γ
-
λ
β
)
/
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
), we have the error on the stationary distributions as follows:
(10)
∥
α
-
π
∥
υ
≤
(
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
)
(
γ
-
λ
)
w
)
×
(
(
β
-
1
)
(
γ
-
λ
β
)
3
(
β
-
1
)
γ
+
λ
β
m
×
-
(
2
γ
-
λ
(
1
+
β
)
)
(
1
+
β
)
(
γ
-
λ
β
)
w
(
β
-
1
)
(
γ
-
λ
β
)
3
(
β
-
1
)
γ
+
λ
β
)
-
1
=
E
r
r
,
where w 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 x of the system (2) and the variation distance w expressed in (8), which both involve the density function of the unknown general distribution G. Statistical methods to do so will be discussed in the following.


## 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(
N
,
B
(
N
)
), where B
(
N
) is the σ-algebra generated by all singletons {
j
}, j
∈
N.LetM
=
{
μ
j
} be the space of finite measures on B
(
N
) and N
=
{
f
(
j
)
} the space of bounded measurable functions on N. We associate with each transition kernel P the linear mapping (
μ
P
)
k
=
∑
j
≥
0
μ
j
P
j
k.
(3)
(
P
f
)
(
k
)
=
∑
i
≥
0
‍
f
(
i
)
P
k
i
.
Introduce on M the class of norms (weight norms) of the form
(4)
∥
μ
∥
υ
=
∑
j
≥
0
‍
υ
(
j
)
|
μ
j
|
,
where υ is an arbitrary measurable function (not necessarily finite) bounded below away from a positive constant. This norm induces in the space N the norm ∥
f
∥
υ
=
sup
⁡
k
≥
0
⁡
(
|
f
(
k
)
|
/
υ
(
k
)
).Let us considerB, the space of linear operators on the space {
μ
∈
M
:
∥
μ
∥
υ
<
∞
}, with norm ∥
Q
∥
υ
=
sup
⁡
k
≥
0
(
1
/
υ
(
k
)
)
∑
j
≥
0
υ
(
j
)
|
Q
k
j
|.Letα and μ be two invariant measures and suppose that these measures have finite υ-norm. Then
(5)
|
α
f
-
μ
f
|
≤
∥
α
-
μ
∥
υ
∥
f
∥
υ
inf
⁡
k
≥
0
υ
(
k
)
.


## 2.2.2. Strong Stability Criterion

Definition 1 (see [1, 2]).
A Markov chainX with a transition kernel P and an invariant measure π is said to be υ-strongly stable with respect to the norm ∥
·
∥
υ defined in (4), if ∥
P
∥
υ
<
∞ and each stochastic kernel Q on the space (
N
,
B
(
N
)
) in some neighborhood {
Q
:
∥
Q
-
P
∥
υ
<
ϵ
} has a unique invariant measure μ
=
μ
(
Q
) and ∥
π
-
μ
∥
υ
→
0 as ∥
Q
-
P
∥
υ
→
0 in this neighborhood.


## 2.2.3. Strong Stability ofM
/
M
/
1 System after Perturbation of the Arrival Flow

(a) Description of the Models. Let us consider a G
/
M
/
1
(
F
I
F
O
,
∞
) system where interarrival times are independently distributed with general distribution G
(
t
) and service times are distributed with E
γ
(
t
) (exponential with parameter γ).LetX
n
* be the number of customers left behind in the system by the nth departure. It is easy to prove that X
n
* forms a Markov chain [2] with a transition operator P
*
=
(
P
i
j
*
)
i
,
j
≥
0, where
(6)
P
i
j
*
=
{
d
i
+
1
-
j
*
=
∫
0
+
∞
‍
1
(
i
+
1
-
j
)
!
d
i
+
1
-
j
*
∫
0
+
∞
‍
×
e
-
γ
t
(
γ
t
)
i
+
1
-
j
d
G
(
t
)
,
if
1
≤
j
≤
i
+
1
,
1
-
∑
k
=
0
i
‍
d
k
*
,
if
j
=
0
,
0
,
otherwise
.
Consider also an M
/
M
/
1
(
F
I
F
O
,
∞
) system, which has Poisson arrivals with parameter λ and the same distribution of the service times as the preceding system. It is known that X
n (the number of customers left behind in the system by the nth departure) forms a Markov chain with a transition operator P
=
(
P
i
j
)
i
,
j
≥
0, where
(7)
P
i
j
=
{
d
i
+
1
-
j
=
λ
γ
i
+
1
-
j
(
λ
+
γ
)
i
+
2
-
j
,
if
1
≤
j
≤
i
+
1
,
1
-
∑
k
=
0
i
‍
d
k
=
(
γ
γ
+
λ
)
i
,
if
j
=
0
,
0
,
otherwise
.
Suppose that the arrival flow of the G
/
M
/
1 system is close to the Poisson one. This proximity is then characterized by the metric
(8)
w
=
w
(
G
,
E
λ
)
=
∫
0
∞
‍
|
G
-
E
λ
|
(
d
t
)
=
∫
0
∞
‍
|
g
-
e
λ
|
(
t
)
(
d
t
)
,
where |
a
| is the variation of the measure a and g and e
λ are the respective density functions of the distributions G and E
λ. The stationary distributions of the states of X
n
* and X
n are defined as follows:
(9)
α
k
=
lim
⁡
n
→
∞
Pr
(
X
n
*
=
k
)
,
k
=
0,1
,
2
,
…
,
π
k
=
lim
⁡
n
→
∞
Pr
(
X
n
=
k
)
,
k
=
0,1
,
2
,
…
,(b) Estimation of the Strong StabilityTheorem 2 (see [8]).
Suppose that the traffic intensity(
λ
/
γ
) of the M
/
M
/
1 system is smaller than 1. Therefore, for all β such that 1
<
β
<
γ
/
λ, the imbedded Markov chain X
n is υ-strongly stable, after a small perturbation of the interarrival times, for υ
(
k
)
=
β
k.
The margin between the transition operators is given by∥
P
*
-
P
∥
v
≤
(
1
+
β
)
w.
In addition, ifw
<
(
1
-
ρ
)
(
γ
-
λ
β
)
/
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
), we have the error on the stationary distributions as follows:
(10)
∥
α
-
π
∥
υ
≤
(
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
)
(
γ
-
λ
)
w
)
×
(
(
β
-
1
)
(
γ
-
λ
β
)
3
(
β
-
1
)
γ
+
λ
β
m
×
-
(
2
γ
-
λ
(
1
+
β
)
)
(
1
+
β
)
(
γ
-
λ
β
)
w
(
β
-
1
)
(
γ
-
λ
β
)
3
(
β
-
1
)
γ
+
λ
β
)
-
1
=
E
r
r
,
where w 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 x of the system (2) and the variation distance w expressed in (8), which both involve the density function of the unknown general distribution G. Statistical methods to do so will be discussed in the following.


## 3. Different Nonparametric Considered Estimates

LetX
1
,
…
,
X
n be a sample coming from a random variable X of density function f and distribution F. The Parzen-Rosenblatt classical kernel estimate (see [4, 5]) of the density f
(
x
) for each point x
∈
R is given by
(11)
f
n
(
x
)
=
1
n
h
n
∑
j
=
1
n
‍
K
(
x
-
X
j
h
n
)
,
where K is a symmetric density function called kernel and h
n is the smoothing parameter (or bandwidth).In practice, the critical step in the kernel density estimation is the choice of the bandwidthh
n 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 lineR [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. f
(
x
) is then estimated, for x
≥
0, as follows:
(12)
f
~
n
(
x
)
=
1
n
h
n
∑
j
=
1
n
‍
[
K
(
x
-
X
j
h
n
)
+
K
(
x
+
X
j
h
n
)
]
.
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:
(13)
f
^
b
(
x
)
=
1
n
∑
i
=
1
n
‍
K
(
x
,
b
)
(
X
i
)
,
where b is the bandwidth and the asymmetric kernel K can be taken as a Gamma density K
G with parameters (
x
/
b
+
1
,
b
) given by
(14)
K
G
(
x
b
+
1
,
b
)
(
t
)
=
t
x
/
b
e
-
t
/
b
b
x
/
b
+
1
Γ
(
x
/
b
+
1
)
.


## 4. Numerical Examples: Simulation

Recall that our object is the numerical evaluation of the approximation error between the stationary distributions ofG
/
M
/
1 and M
/
M
/
1 systems, when the density function of the general arrival law G in the G
/
M
/
1 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 L
1 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.Algorithm 1:Approximation error with L
1 norm.
(1)Begin
(2) Introduce the sample sizen and the traffic intensity ρ of the M
/
M
/
1 system;
(3) Determine the solutionx of the system (2);
(4)For  i  going from 1 to  n  do
π
i
⟵
(
1
-
ρ
)
ρ
i;
α
i
⟵
(
1
-
x
)
x
i;
q
i
⟵
|
π
i
-
α
i
|;
End for;
(5) Computation of the proximity error:
ε
⟵
max
⁡
i
=
1
,
n
¯
|
π
i
-
α
i
|;
(6)End.Algorithm 2:Approximation error with the weight norm.
(0)Begin
(1) Generate a sample of sizen of general probability distribution G with
theoretical densityg
(
x
);
(2) Estimate the theoretical densityg
(
x
) by a nonparametric density noted
in generalg
n
*
(
x
);
(3) Introduce the mean service rateγ of the G
/
M
/
1 system;
(4) Determine the mean arrival rateλ of the G
/
M
/
1 system: λ
⟵
1
∫
‍
x
g
n
*
(
x
)
d
x;
(5) Verify the stability: ifλ
/
γ
≥
1 then “the system is not stable” go to (9);
else{
* the system is v-strongly stable for 1
<
β
<
γ
/
λ
*
} go to (6);
(6) Determine the proximity ofg
n
*
(
x
) and e
λ
(
x
)
=
λ
e
-
λ
x:
w
⟵
∫
‍
|
g
n
*
(
x
)
-
e
λ
(
x
)
|
d
x;
(7) Determine the approximation domain (β
min
⁡
<
β
<
β
max
⁡
) of the system:
β
min
⁡
⟵
min
⁡
(
β
,
1
<
β
<
γ
λ
and
w
<
(
β
(
γ
+
λ
)
-
γ
-
λ
β
2
)
(
γ
-
λ
β
)
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
)
(
β
(
λ
+
γ
)
-
γ
)
);
β
max
⁡
⟵
max
⁡
(
β
,
1
<
β
<
γ
λ
and
w
<
(
β
(
γ
+
λ
)
-
γ
-
λ
β
2
)
(
γ
-
λ
β
)
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
)
(
β
(
λ
+
γ
)
-
γ
)
);
if(
β
min
⁡
≥
β
max
⁡
) then “the proximity is not sufficient”, go to (9);
else{
* the approximation is valid *
} go to (8);
(8) Determine the minimal errorErr on the stationary distribution:
Err
⟵
min
⁡
(
(
1
+
β
)
(
2
γ
-
λ
(
1
+
β
)
)
(
γ
-
λ
)
w
(
(
β
-
1
)
(
γ
-
λ
β
)
3
/
(
(
β
-
1
)
γ
+
λ
β
)
)
-
(
2
γ
-
λ
(
1
+
β
)
)
(
1
+
β
)
(
γ
-
λ
β
)
w
,
β
min
⁡
≤
β
≤
β
max
⁡
);
(9)End.


### 4.1. Algorithms

See Algorithms1 and 2.


### 4.2. Numerical Examples

#### 4.2.1. Approximation Error withL
1 Norm

Example 1.
We consider the following four cases.
First Case. We consider a G
/
M
/
1 system such that the density function of the general law G is given by
(15)
g
(
x
)
=
{
1
2
e
-
x
+
e
-
2
x
,
if
x
≥
0
;
0
,
otherwise
.
By generating a sample coming from the general law G with the same density g
(
x
) defined above, we use the kernel density method to estimate the theoretical density g
(
x
) by using the different estimators given in the three following cases.
Second Case. We use the kernel estimate g
n
(
x
) defined in (11).
Third Case. We use the kernel estimate g
~
n
(
x
) defined in (12).
Fourth Case. We use the kernel estimate g
^
b
(
x
) defined in (13) with the Gamma kernel given in (14).
For the last three cases, we take the sample sizen
=
200 and the number of simulations R
=
100. In all the cases, we introduce the service mean time: γ
=
10.
We first determine the interarrival mean rateλ
=
1
/
∫
‍
x
g
(
x
)
d
x. We get some performance measures of the system: (i)
interarrival mean rate:λ
=
4
/
3;
(ii)
traffic intensity:ρ
=
2
/
15.
For calculation of the proximity error of theG
/
M
/
1 and M
/
M
/
1 systems, defined in (1), we use Algorithm 1. We obtain the results of Table 1.Table 1
Proximity errorε with different estimates.
g
(
x
)
g
n
(
x
)
g
~
n
(
x
)
g
^
b
(
x
)
x
0.1475
0.0668
0.1223
0.1364
ε
0.0142
0.0665
0.0220
0.0170Discussion. We note in Table 1 that the approximation errors found by the kernel method considering the correction of boundary effects (ɛ
=
0.0220 in case of mirror image estimate and ɛ
=
0.0170 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 (ɛ
=
0.0665). 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 theG
/
M
/
1 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 M
/
M
/
1 system is given by (see [12])
(16)
W
(
t
)
=
1
-
ρ
e
-
γ
(
1
-
ρ
)
t
,
t
≥
0
.
The same formula is valid to determine the PDF of the waiting time in the G
/
M
/
1 system, with replacing ρ by x (x is the unique solution of system (2)).
Using results of Table1 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 M
/
M
/
1 queue. For the G
/
M
/
1 system, we use the theoretical density g
(
t
) and its different estimates (
g
n
(
t
)
,
g
~
n
(
t
)
,
g
^
b
(
t
)
). An illustration is given in Figure 1.
Note that according to Figure1, curves of waiting times distributions for the G
/
M
/
1 system obtained using the theoretical density g
(
t
) and its estimates g
~
n
(
t
) and g
^
b
(
t
) are close to that of the waiting times distribution for the M
/
M
/
1 system, contrary to the curve of the waiting times distribution for the G
/
M
/
1 system obtained using the classical kernel estimate g
n
(
t
) which is far away from that of the waiting times distribution for the M
/
M
/
1 system. Once more, this is due to the boundary effects caused by using a symmetric fixed kernel estimate.Figure 1
Comparison of waiting times distributions ofM
/
M
/
1 and G
/
M
/
1 systems.


#### 4.2.2. Approximation Error with the Weight Norm

To realize this work, we use Algorithm2. The Epanechnikov kernel [10] is used throughout for estimators involving symmetric kernels. The bandwidth h
n is chosen to minimize the criterion of the “least squares cross-validation” [10]. The smoothing parameters b and k are chosen according to a bandwidth selection method which leads to an asymptotically optimal window in the sense of minimizing L
1 distance [13].Example 2.
We generate samples of sizen
=
50 issued from different laws. We take the number of simulations R
=
100. For each case of law, we replace the nonparametric density g
n
*
(
t
) (defined in Algorithm 2) by the density function g
n
(
t
) found by applying the Parzen-Rosenblatt classical kernel estimate defined in (11) to estimate the theoretical density g
(
t
) of each sample. Using Algorithm 2, we obtain the results of Table 2.Table 2
w and Err measures with different samples.
Exp(1)
Weibull(2, 0.5, 0)
Gamma(1, 3)
γ
=
10
γ
=
10
γ
=
2
Inter-arrival mean timeλ
0.9190
1.8244
0.2750
Traffic intensityρ of the system
0.0919
<
1
0.1824
<
1
0.1375
<
1
Stability domain
1
<
β
<
10.8811
1
<
β
<
5.4814
1
<
β
<
7.2740
Variationw
0.2444
0.3502
0.1615Discussion. First, notice that according to Table 2, it seems that the proposed method approximates certain values with some differences, for instance, the Exp
(
1
) with the Exp
(
0.919
). This may be explained by the error done when replacing the theoretical density g
(
x
) by the nonparametric classical kernel estimate g
n
(
x
) in the formula used to compute the interarrival mean rate λ; that is, λ
=
1
/
∫
‍
x
g
(
x
)
d
x. 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 Table2, the application of the Parzen-Rosenblatt classical kernel estimate (11) for the approximation of the G
/
M
/
1 system by the M
/
M
/
1 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 w (e.g., w
=
0.2444 for Exp
(
1
) and w
=
0.3502 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 M
/
M
/
1 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 w and the approximation error Err on the stationary distributions for the different estimates. Figure 2 describes the evolution of the error Err in function of β.Table 3
w and Err measures with different estimates.
g
(
x
)
g
n
(
x
)
g
~
n
(
x
)
g
^
b
(
x
)
Variation distancew
0.0711
0.2104
0.0895
0.0792
Error Err
0.21
0.35
0.26ErrorErr in function of β (a) with the theoretical density g
(
x
) and (b) with Gamma kernel estimate g
^
b
(
x
).
(a)
(b)Discussion. We note in Table 3 that the approximation error on the stationary distributions of the G
/
M
/
1 and M
/
M
/
1 systems was given when applying the kernel density method by considering the correction of boundary effects in the case of using mirror image estimate (Err
=
0.35) or in the case of using the asymmetric Gamma kernel estimate (Err
=
0.26). 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 (Err) on the stationary distributions of the quoted systems could not be given.Notice following Figure2(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 (1.07235
<
β
<
1.3) (resp., 1.08146
<
β
<
1.31). 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 (1.8
<
β
<
2.2) (resp., 1.9
<
β
<
2.2) (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 Err which corresponds in our case to β
=
1.3, Err
=
0.21 (resp., β
=
1.31, Err
=
0.26). 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 Corollary5 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 f.Corollary 5.
Suppose that the assumptions of Theorem2 hold, then, for any function f such that ∥
f
∥
υ
<
∞, one has
(17)
|
π
f
-
π
~
f
|
≤
E
r
r
∥
f
∥
υ
,
where E
r
r is defined in (10).Proof.
The proof follows from Theorem2 together with the fact that for any measure α it holds that |
α
f
|
≤
∥
α
∥
υ
∥
f
∥
υ (see formula (5)).Example for the Choice of the Individual Performance Measure
fExample 4 (case: mean number of the waiting customers in the system).
LetQ and Q
~ be the mean numbers of the waiting customers of the nominal (ideal) system and the perturbed one, respectively. To predict the deviation |
Q
-
Q
~
|, it suffices to take the individual performance measure f as follows:
(18)
f
(
k
)
=
{
0
,
if
k
=
0
;
k
-
1
,
if
k
≥
1
.
We have ∥
f
∥
υ
=
sup
⁡
k
≥
0
⁡
(
|
f
(
k
)
|
/
υ
(
k
)); then, for f
(
k
) defined in (18), it is easy to show that ∥
f
∥
υ
=
(
1
/
ln
⁡
(
β
)
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
); thus formula (17) can be written as follows:
(19)
|
π
f
-
π
~
f
|
=
|
Q
-
Q
~
|
≤
Err
1
ln
⁡
(
β
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
)
=
Q
Err
.By dividing expression (19) by λ and by using Little’s formulas, we get
(20)
|
Q
λ
-
Q
~
λ
|
=
|
W
-
W
~
|
≤
Err
1
λ
ln
⁡
(
β
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
)
=
W
Err
.
This last expression provides a bound for the deviation between the mean waiting times of the M
/
M
/
1 and G
/
M
/
1 systems.Applying this formula and results of Table3, we give in Table 4 some values for the error on the mean waiting times of the considered systems, using the theoretical density g
(
x
) defined in (15) and its asymmetric Gamma kernel estimate g
^
b
(
x
) for the G
/
M
/
1 system. We use values of β minimizing the error on the stationary distributions Err defined in (10).Table 4
Error on mean waiting times ofM
/
M
/
1 and G
/
M
/
1 systems.
g
(
x
)
g
^
b
(
x
)
β
1.3
1.31
Err
0.21
0.26
W
Err
0.1699
0.2028


## 4.1. Algorithms

See Algorithms1 and 2.


## 4.2. Numerical Examples

### 4.2.1. Approximation Error withL
1 Norm

Example 1.
We consider the following four cases.
First Case. We consider a G
/
M
/
1 system such that the density function of the general law G is given by
(15)
g
(
x
)
=
{
1
2
e
-
x
+
e
-
2
x
,
if
x
≥
0
;
0
,
otherwise
.
By generating a sample coming from the general law G with the same density g
(
x
) defined above, we use the kernel density method to estimate the theoretical density g
(
x
) by using the different estimators given in the three following cases.
Second Case. We use the kernel estimate g
n
(
x
) defined in (11).
Third Case. We use the kernel estimate g
~
n
(
x
) defined in (12).
Fourth Case. We use the kernel estimate g
^
b
(
x
) defined in (13) with the Gamma kernel given in (14).
For the last three cases, we take the sample sizen
=
200 and the number of simulations R
=
100. In all the cases, we introduce the service mean time: γ
=
10.
We first determine the interarrival mean rateλ
=
1
/
∫
‍
x
g
(
x
)
d
x. We get some performance measures of the system: (i)
interarrival mean rate:λ
=
4
/
3;
(ii)
traffic intensity:ρ
=
2
/
15.
For calculation of the proximity error of theG
/
M
/
1 and M
/
M
/
1 systems, defined in (1), we use Algorithm 1. We obtain the results of Table 1.Table 1
Proximity errorε with different estimates.
g
(
x
)
g
n
(
x
)
g
~
n
(
x
)
g
^
b
(
x
)
x
0.1475
0.0668
0.1223
0.1364
ε
0.0142
0.0665
0.0220
0.0170Discussion. We note in Table 1 that the approximation errors found by the kernel method considering the correction of boundary effects (ɛ
=
0.0220 in case of mirror image estimate and ɛ
=
0.0170 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 (ɛ
=
0.0665). 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 theG
/
M
/
1 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 M
/
M
/
1 system is given by (see [12])
(16)
W
(
t
)
=
1
-
ρ
e
-
γ
(
1
-
ρ
)
t
,
t
≥
0
.
The same formula is valid to determine the PDF of the waiting time in the G
/
M
/
1 system, with replacing ρ by x (x is the unique solution of system (2)).
Using results of Table1 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 M
/
M
/
1 queue. For the G
/
M
/
1 system, we use the theoretical density g
(
t
) and its different estimates (
g
n
(
t
)
,
g
~
n
(
t
)
,
g
^
b
(
t
)
). An illustration is given in Figure 1.
Note that according to Figure1, curves of waiting times distributions for the G
/
M
/
1 system obtained using the theoretical density g
(
t
) and its estimates g
~
n
(
t
) and g
^
b
(
t
) are close to that of the waiting times distribution for the M
/
M
/
1 system, contrary to the curve of the waiting times distribution for the G
/
M
/
1 system obtained using the classical kernel estimate g
n
(
t
) which is far away from that of the waiting times distribution for the M
/
M
/
1 system. Once more, this is due to the boundary effects caused by using a symmetric fixed kernel estimate.Figure 1
Comparison of waiting times distributions ofM
/
M
/
1 and G
/
M
/
1 systems.


### 4.2.2. Approximation Error with the Weight Norm

To realize this work, we use Algorithm2. The Epanechnikov kernel [10] is used throughout for estimators involving symmetric kernels. The bandwidth h
n is chosen to minimize the criterion of the “least squares cross-validation” [10]. The smoothing parameters b and k are chosen according to a bandwidth selection method which leads to an asymptotically optimal window in the sense of minimizing L
1 distance [13].Example 2.
We generate samples of sizen
=
50 issued from different laws. We take the number of simulations R
=
100. For each case of law, we replace the nonparametric density g
n
*
(
t
) (defined in Algorithm 2) by the density function g
n
(
t
) found by applying the Parzen-Rosenblatt classical kernel estimate defined in (11) to estimate the theoretical density g
(
t
) of each sample. Using Algorithm 2, we obtain the results of Table 2.Table 2
w and Err measures with different samples.
Exp(1)
Weibull(2, 0.5, 0)
Gamma(1, 3)
γ
=
10
γ
=
10
γ
=
2
Inter-arrival mean timeλ
0.9190
1.8244
0.2750
Traffic intensityρ of the system
0.0919
<
1
0.1824
<
1
0.1375
<
1
Stability domain
1
<
β
<
10.8811
1
<
β
<
5.4814
1
<
β
<
7.2740
Variationw
0.2444
0.3502
0.1615Discussion. First, notice that according to Table 2, it seems that the proposed method approximates certain values with some differences, for instance, the Exp
(
1
) with the Exp
(
0.919
). This may be explained by the error done when replacing the theoretical density g
(
x
) by the nonparametric classical kernel estimate g
n
(
x
) in the formula used to compute the interarrival mean rate λ; that is, λ
=
1
/
∫
‍
x
g
(
x
)
d
x. 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 Table2, the application of the Parzen-Rosenblatt classical kernel estimate (11) for the approximation of the G
/
M
/
1 system by the M
/
M
/
1 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 w (e.g., w
=
0.2444 for Exp
(
1
) and w
=
0.3502 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 M
/
M
/
1 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 w and the approximation error Err on the stationary distributions for the different estimates. Figure 2 describes the evolution of the error Err in function of β.Table 3
w and Err measures with different estimates.
g
(
x
)
g
n
(
x
)
g
~
n
(
x
)
g
^
b
(
x
)
Variation distancew
0.0711
0.2104
0.0895
0.0792
Error Err
0.21
0.35
0.26ErrorErr in function of β (a) with the theoretical density g
(
x
) and (b) with Gamma kernel estimate g
^
b
(
x
).
(a)
(b)Discussion. We note in Table 3 that the approximation error on the stationary distributions of the G
/
M
/
1 and M
/
M
/
1 systems was given when applying the kernel density method by considering the correction of boundary effects in the case of using mirror image estimate (Err
=
0.35) or in the case of using the asymmetric Gamma kernel estimate (Err
=
0.26). 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 (Err) on the stationary distributions of the quoted systems could not be given.Notice following Figure2(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 (1.07235
<
β
<
1.3) (resp., 1.08146
<
β
<
1.31). 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 (1.8
<
β
<
2.2) (resp., 1.9
<
β
<
2.2) (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 Err which corresponds in our case to β
=
1.3, Err
=
0.21 (resp., β
=
1.31, Err
=
0.26). 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 Corollary5 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 f.Corollary 5.
Suppose that the assumptions of Theorem2 hold, then, for any function f such that ∥
f
∥
υ
<
∞, one has
(17)
|
π
f
-
π
~
f
|
≤
E
r
r
∥
f
∥
υ
,
where E
r
r is defined in (10).Proof.
The proof follows from Theorem2 together with the fact that for any measure α it holds that |
α
f
|
≤
∥
α
∥
υ
∥
f
∥
υ (see formula (5)).Example for the Choice of the Individual Performance Measure
fExample 4 (case: mean number of the waiting customers in the system).
LetQ and Q
~ be the mean numbers of the waiting customers of the nominal (ideal) system and the perturbed one, respectively. To predict the deviation |
Q
-
Q
~
|, it suffices to take the individual performance measure f as follows:
(18)
f
(
k
)
=
{
0
,
if
k
=
0
;
k
-
1
,
if
k
≥
1
.
We have ∥
f
∥
υ
=
sup
⁡
k
≥
0
⁡
(
|
f
(
k
)
|
/
υ
(
k
)); then, for f
(
k
) defined in (18), it is easy to show that ∥
f
∥
υ
=
(
1
/
ln
⁡
(
β
)
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
); thus formula (17) can be written as follows:
(19)
|
π
f
-
π
~
f
|
=
|
Q
-
Q
~
|
≤
Err
1
ln
⁡
(
β
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
)
=
Q
Err
.By dividing expression (19) by λ and by using Little’s formulas, we get
(20)
|
Q
λ
-
Q
~
λ
|
=
|
W
-
W
~
|
≤
Err
1
λ
ln
⁡
(
β
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
)
=
W
Err
.
This last expression provides a bound for the deviation between the mean waiting times of the M
/
M
/
1 and G
/
M
/
1 systems.Applying this formula and results of Table3, we give in Table 4 some values for the error on the mean waiting times of the considered systems, using the theoretical density g
(
x
) defined in (15) and its asymmetric Gamma kernel estimate g
^
b
(
x
) for the G
/
M
/
1 system. We use values of β minimizing the error on the stationary distributions Err defined in (10).Table 4
Error on mean waiting times ofM
/
M
/
1 and G
/
M
/
1 systems.
g
(
x
)
g
^
b
(
x
)
β
1.3
1.31
Err
0.21
0.26
W
Err
0.1699
0.2028


## 4.2.1. Approximation Error withL
1 Norm

Example 1.
We consider the following four cases.
First Case. We consider a G
/
M
/
1 system such that the density function of the general law G is given by
(15)
g
(
x
)
=
{
1
2
e
-
x
+
e
-
2
x
,
if
x
≥
0
;
0
,
otherwise
.
By generating a sample coming from the general law G with the same density g
(
x
) defined above, we use the kernel density method to estimate the theoretical density g
(
x
) by using the different estimators given in the three following cases.
Second Case. We use the kernel estimate g
n
(
x
) defined in (11).
Third Case. We use the kernel estimate g
~
n
(
x
) defined in (12).
Fourth Case. We use the kernel estimate g
^
b
(
x
) defined in (13) with the Gamma kernel given in (14).
For the last three cases, we take the sample sizen
=
200 and the number of simulations R
=
100. In all the cases, we introduce the service mean time: γ
=
10.
We first determine the interarrival mean rateλ
=
1
/
∫
‍
x
g
(
x
)
d
x. We get some performance measures of the system: (i)
interarrival mean rate:λ
=
4
/
3;
(ii)
traffic intensity:ρ
=
2
/
15.
For calculation of the proximity error of theG
/
M
/
1 and M
/
M
/
1 systems, defined in (1), we use Algorithm 1. We obtain the results of Table 1.Table 1
Proximity errorε with different estimates.
g
(
x
)
g
n
(
x
)
g
~
n
(
x
)
g
^
b
(
x
)
x
0.1475
0.0668
0.1223
0.1364
ε
0.0142
0.0665
0.0220
0.0170Discussion. We note in Table 1 that the approximation errors found by the kernel method considering the correction of boundary effects (ɛ
=
0.0220 in case of mirror image estimate and ɛ
=
0.0170 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 (ɛ
=
0.0665). 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 theG
/
M
/
1 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 M
/
M
/
1 system is given by (see [12])
(16)
W
(
t
)
=
1
-
ρ
e
-
γ
(
1
-
ρ
)
t
,
t
≥
0
.
The same formula is valid to determine the PDF of the waiting time in the G
/
M
/
1 system, with replacing ρ by x (x is the unique solution of system (2)).
Using results of Table1 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 M
/
M
/
1 queue. For the G
/
M
/
1 system, we use the theoretical density g
(
t
) and its different estimates (
g
n
(
t
)
,
g
~
n
(
t
)
,
g
^
b
(
t
)
). An illustration is given in Figure 1.
Note that according to Figure1, curves of waiting times distributions for the G
/
M
/
1 system obtained using the theoretical density g
(
t
) and its estimates g
~
n
(
t
) and g
^
b
(
t
) are close to that of the waiting times distribution for the M
/
M
/
1 system, contrary to the curve of the waiting times distribution for the G
/
M
/
1 system obtained using the classical kernel estimate g
n
(
t
) which is far away from that of the waiting times distribution for the M
/
M
/
1 system. Once more, this is due to the boundary effects caused by using a symmetric fixed kernel estimate.Figure 1
Comparison of waiting times distributions ofM
/
M
/
1 and G
/
M
/
1 systems.


## 4.2.2. Approximation Error with the Weight Norm

To realize this work, we use Algorithm2. The Epanechnikov kernel [10] is used throughout for estimators involving symmetric kernels. The bandwidth h
n is chosen to minimize the criterion of the “least squares cross-validation” [10]. The smoothing parameters b and k are chosen according to a bandwidth selection method which leads to an asymptotically optimal window in the sense of minimizing L
1 distance [13].Example 2.
We generate samples of sizen
=
50 issued from different laws. We take the number of simulations R
=
100. For each case of law, we replace the nonparametric density g
n
*
(
t
) (defined in Algorithm 2) by the density function g
n
(
t
) found by applying the Parzen-Rosenblatt classical kernel estimate defined in (11) to estimate the theoretical density g
(
t
) of each sample. Using Algorithm 2, we obtain the results of Table 2.Table 2
w and Err measures with different samples.
Exp(1)
Weibull(2, 0.5, 0)
Gamma(1, 3)
γ
=
10
γ
=
10
γ
=
2
Inter-arrival mean timeλ
0.9190
1.8244
0.2750
Traffic intensityρ of the system
0.0919
<
1
0.1824
<
1
0.1375
<
1
Stability domain
1
<
β
<
10.8811
1
<
β
<
5.4814
1
<
β
<
7.2740
Variationw
0.2444
0.3502
0.1615Discussion. First, notice that according to Table 2, it seems that the proposed method approximates certain values with some differences, for instance, the Exp
(
1
) with the Exp
(
0.919
). This may be explained by the error done when replacing the theoretical density g
(
x
) by the nonparametric classical kernel estimate g
n
(
x
) in the formula used to compute the interarrival mean rate λ; that is, λ
=
1
/
∫
‍
x
g
(
x
)
d
x. 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 Table2, the application of the Parzen-Rosenblatt classical kernel estimate (11) for the approximation of the G
/
M
/
1 system by the M
/
M
/
1 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 w (e.g., w
=
0.2444 for Exp
(
1
) and w
=
0.3502 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 M
/
M
/
1 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 w and the approximation error Err on the stationary distributions for the different estimates. Figure 2 describes the evolution of the error Err in function of β.Table 3
w and Err measures with different estimates.
g
(
x
)
g
n
(
x
)
g
~
n
(
x
)
g
^
b
(
x
)
Variation distancew
0.0711
0.2104
0.0895
0.0792
Error Err
0.21
0.35
0.26ErrorErr in function of β (a) with the theoretical density g
(
x
) and (b) with Gamma kernel estimate g
^
b
(
x
).
(a)
(b)Discussion. We note in Table 3 that the approximation error on the stationary distributions of the G
/
M
/
1 and M
/
M
/
1 systems was given when applying the kernel density method by considering the correction of boundary effects in the case of using mirror image estimate (Err
=
0.35) or in the case of using the asymmetric Gamma kernel estimate (Err
=
0.26). 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 (Err) on the stationary distributions of the quoted systems could not be given.Notice following Figure2(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 (1.07235
<
β
<
1.3) (resp., 1.08146
<
β
<
1.31). 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 (1.8
<
β
<
2.2) (resp., 1.9
<
β
<
2.2) (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 Err which corresponds in our case to β
=
1.3, Err
=
0.21 (resp., β
=
1.31, Err
=
0.26). 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 Corollary5 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 f.Corollary 5.
Suppose that the assumptions of Theorem2 hold, then, for any function f such that ∥
f
∥
υ
<
∞, one has
(17)
|
π
f
-
π
~
f
|
≤
E
r
r
∥
f
∥
υ
,
where E
r
r is defined in (10).Proof.
The proof follows from Theorem2 together with the fact that for any measure α it holds that |
α
f
|
≤
∥
α
∥
υ
∥
f
∥
υ (see formula (5)).Example for the Choice of the Individual Performance Measure
fExample 4 (case: mean number of the waiting customers in the system).
LetQ and Q
~ be the mean numbers of the waiting customers of the nominal (ideal) system and the perturbed one, respectively. To predict the deviation |
Q
-
Q
~
|, it suffices to take the individual performance measure f as follows:
(18)
f
(
k
)
=
{
0
,
if
k
=
0
;
k
-
1
,
if
k
≥
1
.
We have ∥
f
∥
υ
=
sup
⁡
k
≥
0
⁡
(
|
f
(
k
)
|
/
υ
(
k
)); then, for f
(
k
) defined in (18), it is easy to show that ∥
f
∥
υ
=
(
1
/
ln
⁡
(
β
)
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
); thus formula (17) can be written as follows:
(19)
|
π
f
-
π
~
f
|
=
|
Q
-
Q
~
|
≤
Err
1
ln
⁡
(
β
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
)
=
Q
Err
.By dividing expression (19) by λ and by using Little’s formulas, we get
(20)
|
Q
λ
-
Q
~
λ
|
=
|
W
-
W
~
|
≤
Err
1
λ
ln
⁡
(
β
)
β
-
(
1
/
ln
⁡
(
β
)
+
1
)
=
W
Err
.
This last expression provides a bound for the deviation between the mean waiting times of the M
/
M
/
1 and G
/
M
/
1 systems.Applying this formula and results of Table3, we give in Table 4 some values for the error on the mean waiting times of the considered systems, using the theoretical density g
(
x
) defined in (15) and its asymmetric Gamma kernel estimate g
^
b
(
x
) for the G
/
M
/
1 system. We use values of β minimizing the error on the stationary distributions Err defined in (10).Table 4
Error on mean waiting times ofM
/
M
/
1 and G
/
M
/
1 systems.
g
(
x
)
g
^
b
(
x
)
β
1.3
1.31
Err
0.21
0.26
W
Err
0.1699
0.2028


## 5. Conclusion

By comparing the results of Tables1 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 Err (according to weight norm) is quite large compared to the error ɛ (according to L
1 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 Figure1 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 (G
/
M
/
1 and M
/
M
/
1).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 theG
/
G
/
1 system, risk and inventory models, and networks.



---

*Source: 871357-2014-08-14.xml*