Bayesian Nonparametric Inference for M/G/1 Queueing
Systems
Moritz von Rohrscheidt∗1 and Cornelia Wichelhaus

†2

arXiv:1703.07072v1 [] 21 Mar 2017

1 Ruprecht-Karls
2 Technische

Universität Heidelberg
Universität Darmstadt

March 22, 2017

Abstract
In this work, nonparametric statistical inference is provided for the continuous-time
M/G/1 queueing model from a Bayesian point of view. The inference is based on observations of the inter-arrival and service times. Beside other characteristics of the system,
particular interest is in the waiting time distribution which is not accessible in closed
form. Thus, we use an indirect statistical approach by exploiting the Pollaczek-Khinchine
transform formula for the Laplace transform of the waiting time distribution. Due to this,
an estimator is defined and its frequentist validation in terms of posterior consistency and
posterior normality is studied. It will turn out that we can hereby make inference for the
observables separately and compose the results subsequently by suitable techniques.

Keywords: Bayesian Statistics, Nonparametric Inference, Queues in Continuous Time

1

Introduction

Much attention has been drawn to the analysis and the statistics of queueing systems. Beside approaches to statistical inference for queueing systems in the classical frequentist sense,
during the last three decades there has taken place mentionable research in what is frequently
called Bayesian queues. By that we mean inference for queueing systems from the Bayesian
point of view. Bayesian statistical approaches are often feasible and useful in the area of
Operations Research since one is able to express prior knowledge about the system which is
often present, see also e.g. the foreword by Dennis Lindley in De Finetti (1974). It is worth
to mention here some of the recent works on Bayesian queues. Armero and Bayarri (1994)
deals with Markovian queues in the Bayesian sense, i.e. queues with a homogeneous Poisson
process as input stream and exponential service times with varying number of service stations.
In this parametric model, the parametrically conjugate Gamma prior is put on the parameter
space which contains the arrival and service parameter with the aim to evolve a Bayesian test
for the stability of the queueing system.
∗
†

rohrscheidt@uni-heidelberg.de
wichelhaus@mathematik.tu-darmstadt.de

1

This work forms the base for further works on Bayesian inference for systems with Markovian
characters under several generalizations cf. e.g. Armero and Conesa (1998, 2004, 2006) and
Ausín et al. (2007, 2008).
Generalizations in many ways were performed as for example a non-Markovian setting for
the inter-arrival times, cf. Wiper (1998) and Ausín et al. (2007, 2008). Assuming the arrival
stream of customers to be Poisson is often well-suited and ensures desirable properties of the
stochastic process in the background of the system, meaning that there is a description in
terms of an embedded Markov chain. In contrast, the service times often have to be modeled
to be governed by a more general distribution. Moreover, more flexible classes of service time
distributions often preserve the aforementioned property. Thus, it is reasonable to generalize
the distribution of the service times which, in the classical parametric treatment, is taken to
be exponential. Models with generalizations of the service time distribution can be found in
Insua et al. (1998), where Erlang and hyperexponential distributions are employed to model
the service times more flexibly. Ausín et al. (2004) use the well known probabilistic result
that the class of phase-type distributions forms a dense set in the space of all probability
distributions [cf. Asmussen (2008, p.84)]. They model the service time distribution semiparametrically, assign prior distributions to its parameters for eventually inferring the system
by MCMC procedures.
However, one is often interested in a even more general approach for modeling the service time
distribution which leads to a non-parametric Bayesian approach to estimate the unknown service time distribution. The first work in a discrete-time framework concerning this was given
by Conti (1999) for a Geo/G/1 system which can be seen to be the discrete-time analog of the
continuous-time M /G/1. Geo/G/1 is frequently used to model communication systems where
information encoded in packages of fixed size is transmitted by a serving unit. The server is
assumed to be able to transmit one package during one time slot. Hence, the randomness of
the service times can rather be seen to be given in form of random batch sizes, i.e. marks of
the point process governing the arrival stream. Employing a Dirichlet process prior for the distribution of the magnitude of these marks, Conti obtains estimators for various characteristics
of the queue by a well known functional relationship [see also e.g. Grübel and Pitts (1992)]
of the inter-arrival and the service time distribution with the waiting time distribution. Since
the target of Conti’s work was to obtain estimates for the waiting times of customers, he took
this rather indirect approach mainly because assigning a prior distribution to the waiting time
distribution and updating it by data consisting of the marked arrival stream is infeasible. He
obtained large sample properties for the estimator of the probability generating function of
the waiting time distribution. These were uniform posterior consistency and a Bernstein-von
Mises type result. The former roughly says that the posterior distribution, i.e. the prior updated by the data, will center around the true value while the latter forms a Bayesian analog
to the central limit theorem and gives an idea how this centering looks like. For the Bernsteinvon Mises theorem results in Freedman (1963) were employed.
The aim of the present work is to take a similar way of Bayesian statistical inference for the
continuous-time M /G/1 system which is the more appropriate model in many situations where
time-continuity is more reasonable. For example one could think of customers arriving at a
cashier in a supermarket, cars arriving at a traffic jam, goods arriving at a storing center and
many others. Therefor, a continuous-time analog for the functional relationship of observables
and objects of interest is used which is known as the functional Pollaczek-Khinchine formula
in honor of Felix Pollaczek and Aleksandr Khinchine who derived the steady state behaviour
2

of the M /G/1 system in the 1930s. Thereby, the philosophy is analog to that of Conti (1999),
meaning that the main target will be the nonparametric estimation of the waiting time distribution on basis of observations of the arrival stream and the service times. However, a richer
class of prior distributions is used which allows to express prior knowledge more flexibly. We
feel that this is more reasonable for a system evolving in continuous time. The combination of
this larger class of prior distributions together with the assumption of continuous time leads
to more intricate proofs of results related to the discrete-time analog. As usual, the Bayesian
approach is appreciable if prior knowledge and only few data is available. However, since
typically large data samples are accessible, large sample properties as posterior consistency
and posterior normality are investigated.
The paper is organized the following way. In section 2 the underlying queueing model is
introduced and the functional relationship between observables and characteristics of interest
is briefly described. Section 3 is devoted to the assignment of prior distributions to the
random distributions governing the observables. Therefor, we review some facts of Bayesian
statistics and call to mind a large family of non-parametric priors for eventually obtaining
suitable estimators. Subsequently, frequentist validations of the suggested estimators is given
in section 4 and section 5. While section 4 deals with the concentration of the posterior law
of the random quantities in form of posterior consistency results, section 5 is devoted to the
question how these concentrations take place yielding Bernstein-von Mises type results.

3

2

The Queueing Model

The present section is devoted to the underlying queueing model. The model is chosen, in
Kendall’s notation [Kendall (1953)], to be the M /G/1/∞-FIFO system. This notation means
that (M) indistinguishable customers arrive consecutively according to a homogeneous Poisson
process. Here, the M stands for "Markovian" to depict the memorylessness property of the
exponential distribution which governs the independent inter-arrival times (An )n∈Z . Note that
most of the appreciable properties of the M /G/1 system are due to this assumption. Moreover,
customers requiring service are served according to a general service time distribution (G)
concentrated on R+ . The service times are assumed to form a sequence of i.i.d. random
variables (Sn )n∈Z . The consecutive services are accomplished by one (1) reliable service station
in a first-in-first-out (FIFO) manner. Customers who cannot be served immediately are stored
in an infinitely large waiting-room (∞) and thus form a queue which, provided the queue is
stable and in steady state, is distributed according to a stationary distribution. For the sake
of clarity of what follows, we introduce the following notation. We hereby assume that all
random variables are defined on a common underlying probability space (Ω, A, P).
• Qt ≡ queue length at time t ∈ R
• Nt ≡ number of customers in the entire system at time t ∈ R,
i.e. Nt = Qt + 1 in case Qt > 0 for t ∈ R
• Wn ≡ waiting time of customer n ∈ Z
• Dn ≡ sojourn time in the system of customer n ∈ Z
• An ≡ time between arrival of n-th and (n − 1)-st customer n ∈ Z
• Sn ≡ service time of customer n ∈ Z
Letting λ denote the arrival rate of the Poisson arrival stream, µ the mean of the amount of
time the server needs to complete the service of a customer and defining the so-called traffic
intensity ρ ∶= E[S]/E[A] = λµ, a well known assumption for the existence of a stationary
distribution for N is ρ < 1. If the queue is stable, by stationarity it is common to denote above
quantities without an index tacitly assuming that the system has already run for an infinitely
large time-horizon. If additional one assumes E[S 2 ] < ∞ it can be shown that
E[N ] = ρ +

ρ2 + λVar[S]
.
2(1 − ρ)

This formula is known as the Pollaczek-Khinchin mean formula and relates the mean of the
distribution of N to that of A, S and to the variance of S, Var[S]. However, this formula only
relates parameters of aforementioned distributions to each other instead of the entire distributions. However, in Bayesian statistics one is typically interested in the entire distribution of a
parameter one is uncertain in. A functional relationship is given by the Pollaczek-Khinchine
k
transform formula. To state it, define the following objects. Let n(z) = ∑∞
k=0 z P(N = k) be
∞
the probability generating function (p.g.f.) of the distribution of N and g(z) = ∫0 e−zs dG(s)
the Laplace-Stieltjes transform (LST) of the service time distribution. Then the following
functional relationship can be obtained, see e.g. Haigh (2004) or Nelson (2013).
n(z) = g(λ(1 − z))

(1 − z)(1 − ρ)
,
g(λ(1 − z)) − z

z ∈ [0, 1].

4

Manipulations of this formula give
q(z) =

(1 − z)(1 − ρ)
,
g(λ(1 − z)) − z

and

z ∈ [0, 1],

w(z) =

z(1 − ρ)
,
z − λ(1 − g(z))

z ∈ R+ ,

where q(z) denotes the p.g.f. of the queue-length distribution and w(z) the LST of the waitingtime distribution. The latter two quantities are of special statistical interest since they provide
the essential information about the development of the queue.

5

3

Prior Assignments and Estimators

All the transforms n(z), q(z) and w(z) introduced in section 2 depend on the arrival rate λ
and on the LST of the service time distribution g(z). These values are typically assumed to
be unknown and thus need to be inferred. For this inference, we choose a Bayesian approach
which is further introduced in the present section. Since it does not make sense to a subjectivist statistician to talk about "fixed but unknown" parameters [see De Finetti (1974)],
another philosophical concept is used. That is one interprets the "fixed but unknown" parameter as a random quantity itself. This approach leads to the concept of exchangeability which
provides a "meaningful, observable character" of the data. This approach is briefly reviewed.
Let (Ω, A, P) be an abstract probability space. Call an infinite sequence of random variables
(Xi )∞
i=1 with Xi ∶ Ω → R, i ∈ N, exchangeable if for any m ∈ N and any permutation π of m
elements it holds that
(E)

L [X1 , . . . , Xm ] = L [Xπ(1) , . . . , Xπ(m) ] ,
where L denotes the joint law of the respective random object.

Now, assume that we observe n inter-arrival times An1 = (A1 , . . . , An ) between the first
(n + 1) consecutive customers as well as the service times S1n = (S1 , . . . , Sn ) of the first n
customers. The data An1 and S1n are assumed to be the first n projections of two independent
∞
sequences of exchangeable random variables A∞
1 and S1 which, in turn, are assumed to be
independent of each other, i.e. L(Si ; i ∈ I∣Aj ; j ∈ J) = L(Si ; i ∈ I) for all finite subsets I, J ⊂ N.
We turn to de Finetti’s theorem for Polish spaces [Hewitt and Savage (1955)]. Let P(S)
denote the space of all probability measures on some Polish space S and consider in particular
P(R) equipped with the topology of weak convergence of measures. This leads to a measurable
space (P(R), BP(R) ) which is itself Polish [e.g. Kechris (1995)]. By the de Finetti theorem
for exchangeables, it holds for all n ∈ N and for all measurable subsets Ai ⊂ R, i = 1, . . . , n that
X1∞ = (X1 , X2 . . . ) is exchangeable if and only if there is a unique mixing measure ν ∈ P(P(R))
such that
n

P(Xi ∈ Ai ; i = 1, . . . , n) = ∫
P(R)

∏ P (Xi ∈ Ai )ν(dP ).
i=1

The right-hand side reflects the equivalent property of the sequence X1∞ being exchangeable
the following way. The data X1∞ are conditionally i.i.d. given some probability measure
P ∈ P(R), in symbols write X1∞ ∣P ∼ ⊗N P . The probability measure P itself is random and
distributed according to the mixing measure ν which is called the prior distribution from a
Bayesian statistical point of view. Moreover, notice that unconditionally the data X1∞ are
in general not independent. Indeed, by a result from Kingman (1978), it is easily seen that
in general data are positively correlated, a fact that makes Bayesian statistics a theory of
statistical prediction and thus convenient for other theories as e.g. machine learning.
For actual applications, however, the prior ν and the integration in de Finetti’s theorem
becomes infeasible in general. This is mainly due to the fact that one has merely a nonconstructive proof of the existence of ν. Another difficulty is that the set P(R) is quite large.
One has mainly two ways to choose to circumvent this problem in applications. The first
is to shrink the support of ν to a reasonable subset of P(R) by an additional "symmetry
6

constraint" on the sequence of random variables, thus, to put additional information on the
data, to boil down above integration to a finite dimensional parameter space. The second is
to model ν non-parametrically by genius probabilistic tools. We use the first approach here
to get a prior distribution for the random arrival rate λ ∶ Ω → R+ and the second to get a prior
distribution for the random service time distribution G ∶ Ω → P(R+ ).

3.1

Arrival Rate

In case of the inter-arrivals an additional symmetry assumption on the law of the infinite
exchangeable sequence A∞
1 will give rise to a mixing measure which is supported by all exponential distributions. This additional symmetry condition is stated as follows. Let n ∈ N
and for i = 1, ..., n, Bi ∈ BR+ , the Borel sigma-field on R+ . Furthermore, let ci ∈ R be real
constants such that ∑ni=1 ci = 0 and Ci = ci + Bi = {ci + r ∶ r ∈ Bi } ⊂ R+ , i = 1, . . . , n. The law of
the exchangeable sequence A∞
1 fulfills the symmetry condition
P(Ai ∈ Bi ; i = 1, ..., n) = P(Ai ∈ Ci ; i = 1, ..., n),

(S)

for all n, Bi , Ci as above if and only if it is a mixture of exponential distributions [Diaconis
and Ylvisaker (1985)], in symbols
n

P(Ai ∈ Bi ; i = 1, ..., n) = ∫ ∏ P (Ai ∈ Bi )ν(dP )
E

i=1
n

= ∫ ∏∫
R+ i=1

Bi

λe−λxi dxi ν̃(dλ),

where E denotes the space of all exponential distributions on R+ and ν̃ denotes the pushforward measure of ν along the natural parametrization ˜∶ E → R+ ; P ↦ λ of the exponential
distributions. Moreover, if in above situation additionally it holds E[A2 ∣A1 ] = αA1 + β for real
constants α, β > 0, then the mixing measure ν̃ can be shown to be a Gamma distribution.
Since the M /G/1 queueing model implies a mixture of exponential distributions for the
joint law of the inter-arrivals, from a Bayesian point of view, we assume A∞
1 to meet constraints
(E), (S) and that E[A2 ∣A1 ] = αA1 + β holds as well (note that by exchangeability, this extends
to all random inter-arrival times). The latter directly leads to a conjugate prior for the arrival
rate. To be more precise, if we assume λ to be a random variable distributed according to a
Gamma distribution with hyper parameters a, b > 0, then the posterior distribution given the
data X1n is a Gamma distribution as well with updated hyper-parameters (a + n, b + ∑ni=1 Ai ),
i.e. the family of gamma distributions is closed with respect to exponential sampling. In
summary, we assume the following sampling scheme
λ∣(a, b) ∼ Γ(a, b)
A∞
1 ∣λ ∼ ⊗ E(λ),
N

which leads, through Bayes theorem, to the posterior distribution, the Bayes estimate for

7

squared error loss and the posterior predictive distribution density, respectively, given by
n

λ∣(a, b), An1 ∼ Γ (a + n, b + ∑ Ai )
i=1

a+n
=
b + ∑ni=1 Ai
(a + n) (b + ∑ni=1 Ai )a+n
f (an+1 ∣An1 , (a, b)) =
.
(b + ∑ni=1 Ai + an+1 )a+n+1
EΓ [λ∣An1 , (a, b)]

The latter leads to the predictive value for the next observation
E[An+1 ∣An1 , (a, b)] =

n
1
b
.
∑ Ai +
a + n − 1 i=1
a+n−1

Note that the latter equation again reflects the learning process which does not exist in the frequentistic approach in such an explicit form and that for n = 1 it is given by E[A2 ∣A1 , (a, b)] =
1/aA1 + b/a.

3.2

Service Time Distribution

Since the M /G/1 model does not imply a parametric mixture for the service time random
variables as for the inter-arrivals, assigning a suitable prior is a more difficult task. Not being
able to shrink the support of the mixing measure ν in the de Finetti theorem to a finitedimensional set, we have to choose a prior that supports most of P(R+ ). The common way is
to parametrize P(R+ ) by a reasonable dense subset. This is taken to be the set of all discrete
∞
distributions on R+ , Pd (R+ ) = {P ∈ P(R+ ) ∶ P (⋅) = ∑∞
i=1 wi δxi (⋅); wi ∈ [0, 1], ∑i=1 wi = 1, xi ∈
R+ ; ∀i ∈ N}. It is well-known that the most famous non-parametric prior in Bayesian statistics,
namely the Dirichlet process prior obtained in Ferguson (1973), samples discrete probability
measures with probability one. Moreover, it is known that it has full weak support of all of
P(R+ ) if its base measure has support all of R+ . However, here we will use a slightly more
general family of prior distributions that enables us to model more flexibly prior beliefs of the
true data generating distribution G0 . This larger family will be a subclass of so called neutral
to the right prior processes, namely the beta-Stacy processes. Although these priors sample
discrete probabilities with probability one, too, an analogue of the result concerning the weak
support is known as well. The class of neutral to the right priors is now briefly introduced
and its most important properties will be stated.
3.2.1

Neutral to the right Priors

Let F(R+ ) denote the space of all cumulative distribution functions (c.d.f.) on R+ . Then, a
random distribution function F ∈ ΩF (R+ ) is said to be neutral to the right (NTR) if for each
k > 1 and 0<t1 < t2 < ... < tk the normalized increments
F (t1 ),

F (t2 ) − F (t1 )
F (tk ) − F (tk−1 )
, ...,
1 − F (t1 )
1 − F (tk−1 )

are independent assuming F (t) < 1, ∀t ∈ R+ . That is, for all i ∈ N, F̄ (ti )/F̄ (ti−1 ) is independent
of the sigma-field generated by F up to time ti−1 , σ({F (t) ∶ t < ti−1 ), where F̄ (⋅) ∶= 1 − F (⋅) is
the survival function associated to F . This essentially asserts that the proportion of mass that
8

F assigns to (ti , ∞) with respect to (ti−1 , ∞) does not depend on how F behaves left of ti−1 .
This property coined the name neutral to the right. Doksum (1974) has shown that F (⋅) ∈ ΩF
is NTR if and only if L(F (⋅)) = L(1 − exp[−A(⋅)]) for some independent increment process
A(⋅) which is almost surely non-decreasing, right continuous and such that limt→−∞ A(t) = 0
and limt→∞ A(t) = ∞. Such objects are called increasing additive processes, see e.g. Sato
(1999). For more details on the construction of NTR priors see e.g. Phadia (2015). Since
independent increment processes are well understood, the definition of NTR priors leads to
a rich class of non-parametric priors which are analytically tractable. Another nice feature
of NTR priors is that this family is conjugate with respect to (right-censored) exchangeable
data. A fact that makes NTR priors appreciable in statistical survival analysis. Notice, that
a Dirichlet process prior updated by right censored data is not longer Dirichlet.
The next proposition makes a statement about the weak support of a NTR prior. Recall
that the topological support of a measure is the smallest closed set with full measure. In the
following theorem F(R+ ) is identified with the space of all probability measures P(R+ ) which
is equipped with the weak topology.
Proposition 3.1. Dey et al. (2003)
Let F be a random distribution function which is governed by a NTR prior Π ∈ P(F), i.e.
F ∼ Π, and let A(⋅) = − log(1 − F (⋅)) be the corresponding positive increasing additive process
with Lévy measure L. Then, Π has full support if L has full support. The assertion remains
true if F(R+ ) is equipped with the sup-norm.
Proof. It has to be shown that all weak neighborhoods of any probability measure have positive Π-mass. Since continuous distributions are dense in P(R+ ) with respect to the weak
topology, it suffices to show the assertion for all weak neighborhoods of continuous distributions. Now, choose some random distribution function F0 and some  > 0 and consider
U ∶= {F ∶ sup ∣F (t) − F0 (t)∣ < }. Then [e.g. Ghosh and Ramamoorthi (2003))] W ⊆ U for
0≤t<∞

some weak neighborhood W of F0 . Moreover, there is δ > 0, m ∈ N and 0 < t1 < t2 < ⋯ < tm
such that {F ∶ ∣F (ti , ti+1 ] − F0 (ti , ti+1 ]∣ < δ, i = 1, . . . , m} ⊆ W . Hence, it is enough to show
that sup-neighborhoods restricted to compact sets possess positive prior probability. By the
homeomorphism φ ∶ F (⋅) ↦ − log[1 − F (⋅)] the problem can be translated to the analogous one
for the corresponding Lévy process. That is to show that L(A) gives positive probability to
sets of the form C ∶= {A ∶ sup ∣A(t) − φ−1 F0 (t)∣ < γ} for fixed r ∈ Q. Now, take a partition
0≤t≤r

ρ = ∪ki=1 (ai , ai+1 ] of (0, r] such that sup φ−1 F0 (aj , aj+1 ] < γ and define
1≤j≤k

Bi ∶= (ai , ai+1 ] × (φ−1 F0 (ai , ai+1 ] − γ/k, φ−1 F0 (ai , ai+1 ] + γ/k),
B ∶= ∩ki=1 {A ∶ #{(t, A({t}) ∈ Bi } = 1} .
It follows that B ⊆ C and, since L was assumed to have full support,
k

φ−1 Π(C) ≥ φ−1 Π(B) = ∏ L(Bi )e−L(Bi ) > 0.
i=1

9

3.2.2

Beta Processes and Beta-Stacy Processes

For our purposes, we choose a NTR prior with corresponding process Y (⋅) being driven by a
certain class of Lévy measures. These were studied by Hjort (1990) and Walker and Muliere
(1997). Hjort studied beta-processes from a survival analysis viewpoint and therefor elicited
a non-parametric prior for the cumulative hazard function (c.h.f.) given for F ∈ F(R+ ) by
t

H(t) = ∫

0

dF (s)
, t > 0.
F̄ (s)

Walker and Muliere (1997) give the definition of the analog of the beta-process as a prior for
the c.d.f. directly as follows. F is said to be distributed according to a beta-Stacy process
+
with parameters (c(⋅), H(⋅)) ∈ RR
+ × F(R+ ) (for short F ∼ BS(c, H)) if for all t ≥ 0 the
corresponding process Y (⋅) fulfilling F (⋅) = 1 − exp[−Y (⋅)] has Lévy measure
t

dx
−xc(s)(1−H(s))
dLt (x) =
c(s)dHc (s),
∫ e
1 − e−x
0

for all x > 0, where Hc (t) = H(t) − ∑ H(tk ) is the continuous part of H with ti as the fixed
k∶tk <t

points of discontinuity of H. Since we are in a continuous-time framework, we will always
choose H to be continuous. However, discontinuities appear in the Lévy measure governing
the posterior law. Note that EBS [F (⋅)] = H(⋅) is the prior guess on the c.d.f. F and the
function c(⋅) acts like a tuning parameter affecting the magnitude of the increments and is
sometimes interpreted as the "flexible" belief in the prior guess. For the sake of clarity, note
that the Dirichlet process with finite measure α as parameter admits a similar representation
with Lévy-measure
t

dx
−xc(1−ᾱ(s))
dDt (x) =
cᾱ(ds),
∫ e
1 − e−x
0

α(⋅)

where c = α(R) and ᾱ(⋅) denotes the c.d.f. corresponding to the probability measure c . In
the case of Dirichlet processes it is well-known that the prior guess on the random probability
measure equals 1/cα(⋅) and c itself is often interpreted as the strength of belief in the prior
guess. Thus, a Dirichlet process is a beta-Stacy process whose parameters are determined by
α alone.
As already mentioned, the beta-Stacy process is a parametrically conjugate prior, meaning
that the posterior law of F given exchangeable possibly right-censored data is a beta-Stacy
process as well. To be more precise, let S1n be the first n projections of infinite exchangeable
data S1∞ . Let S1∞ be conditional iid given the c.d.f. F and let F ∼ BS(c, H).
n

n

i=1

i=1

Further, define M (t) = ∑ 1[t,∞) (Si ) and N (t) = ∑ 1[0,t] (Si ). Then it holds [Walker and
Muliere (1997, Theorem 4)] that F ∣S1n ∼ BS(c∗n , Hn∗ ), where c∗n (⋅) and Hn∗ (⋅) are given by
Hn∗ (t) = 1 − ∏ (1 −
s∈[0,t]

c∗n (t) =

c(s)dH(s) + dN (s)
),
c(s)H̄(s) + M (s)

c(t)H̄(t) + M (t) − N (t)
.
H¯∗ (t)
10

Thereby, Hn∗ (t) is defined by means of the product integral, see Gill and Johansen (1990).
Note, that the posterior process possesses fixed points of discontinuity at the observations and
that the posterior guess on F , i.e. the Bayes estimate with respect to squared error loss is
given by F̂n (⋅) = Hn∗ (⋅).
Having a look at the Pollazcek-Khinchine transform formulas in section 2, one may ask if
they are well defined for almost all λ and F drawn from their respective prior and posterior
distributions. While the posterior of the mean inter-arrival times is straightforward by definition, the mean and the second moment of the random c.d.f. describing the general service
times deserves more attention. The existence of functionals of c.d.f.s drawn according to a
beta-Stacy process was studied in Epifani et al. (2003). Relating A(⋅) and c(⋅) to the existence
of a certain functional, they obtained sufficient conditions for moments of order m to exist
[equation (10) in their article]. We will assume throughout that this condition holds at least
for the moment of second order. Moreover, note that they obtained an explicit formula for
the prior moments [equation (11)] as well as for the posterior moments [equation (13)].
In summary, we assume the following sample scheme for the service times in the M /G/1
system.
G∣(c, H) ∼ BS(c, H)
S1∞ ∣G ∼ ⊗ G.
N

This leads to the posterior distribution
G∣(c, H), S1n ∼ BS(c∗n , Hn∗ ).
Notice again that this updating continues to hold for right-censored observations. Hence,
the model can be enlarged in that one does not have to keep exact track of the customers
and is still able to make reasonable inference even if e.g. it is solely known that the service of
customers has exceeded a certain threshold.

3.3

Estimators for Queueing Characteristics

In this subsection, we study the Bayes estimators, i.e. the posterior means, of several characteristics of the M /G/1 system. Not all posterior laws of each single characteristic are obtainable
explicitly. Hence, for those characteristics whose posterior laws are not obtainable in closed
form, we define suitable estimators by replacing Bayes estimators for the corresponding values.
Assume for the general service time distribution G ∼ BS(c, H) and let
Ĝn = EBS [G∣S1n ] ,
µ̂n = EBS [∫
λ̂n =

∞

0
EΓ [λ∣An1 ]

tdG(t)∣S1n ] and

be the Bayes estimates with respect to squared error loss. Note that with Mn (s) = ∑nj=1 δXj [s, ∞),
as in Epifani et al. (2003), µ̂n can be given as
µ̂n = ∫

0

∞

t

exp [− ∫

0

t β(s) + M (s) − 1
α(ds)
n
] exp [− ∫
Nn (ds)] dt,
β(s) + Mn (s)
β(s) + Mn (s)
0

11

where if H is continuous [see Phadia (2015)]
• β(t) = c(t)[1 − H(t)]
t

• α(t) = ∫0 c(s)dH(s).
Observe that ρ̂n = EBS⊗Γ [λµ∣S1n , An1 ] = µ̂n λ̂n by above independence assumption. However, in
Bayesian statistics, one is interested in the entire posterior law rather than merely in a certain
functional. The posterior law of ρ is obtainable explicitly in the following form
PBS⊗Γ (ρ ≤ t∣An1 , S1n ) = PBS⊗Γ (µλ ≤ t∣An1 , S1n )
∞

=∫

0

=∫
=

0

∞

PBS⊗Γ (µ ≤ t/λ∣An1 , S1n , λ)PBS⊗Γ (λ∣An1 , S1n )dλ
PBS (µ ≤ t/λ∣S1n , λ)PΓ (λ∣An1 )dλ

∞
(b + ∑ni=1 Ai )a+n
n
a+n−1 −λ(b+∑n
i=1 Ai ) dλ.
e
∫ PBS (µ ≤ t/λ∣S1 , λ)λ
Γ(a + n)
0

Note that PBS (µ ≤ t/λ∣S1n , λ) can be stated more explicitly in form of a density by means of
Proposition 4 in Regazzini et al. (2003). This explicit form is omitted here due to technical
reasons and left to the interested reader. If one observes the queueing system and is, as is
usually the case, not aware of the condition ρ0 < 1, one can thus obtain the posterior probability that the system is stable PBS⊗Γ (ρ < 1∣An1 , S1n ).
We close this section by defining estimates for:
• mean service time: µ̂n
• traffic intensity: ρ̂n = λ̂n µ̂n
∞

• LST of service-time distribution: gn∗ (z) = ∫0 e−zt dĜn (t)
• LST of waiting-time distribution: wn∗ (z) =
• pgf of queue-size distribution: qn∗ (z) =

z(1−ρ̂n )
∗ (z))
z−λ̂n (1−gn

(1−ρ̂n )(1−z)
∗ (λ̂ (1−z))−z
gn
n

• system-size distribution: n∗n (z) = qn∗ (z)gn∗ (λ̂n (1 − z)).
Notice that for the random arrival- and service rate and the random distribution function
of the service times the natural Bayes estimator is used, i.e. the minimizer with respect to
squared error loss. They are analytically tractable and obtainable in closed form. For the
remaining queueing characteristics obvious plug-in estimators are used. The reason therefor
is that a closed form of the push-forward law under the mapping
(λ, G(⋅)) ↦ f (⋅),
where f ∈ {n, g, q}, of the prior ΠBS⊗Γ is not easy to obtain. Neither is known how to update such a pushed prior by exchangeable data (S, A)n1 which we have access to in a natural
Bayesian way.
However, the goodness of these estimates in a rigorous mathematical sense is established
in the next sections.
12

4

Consistency Results

Posterior consistency provides a tool for validation of Bayesian procedures. Roughly speaking, it is defined to be the property of the posterior law to center around the true parameter
when the number of observed data increases. As a consequence of posterior consistency two
different priors will asymptotically lead to the same prediction, a fact often called merging of
prior opinions. So, prior information is consecutively washed away by the new state of information provided by the data. Posterior consistency is the most desired property of Bayesian
procedures since it states that one can recover the true measure that generates the data. For
examples of inconsistent Bayes procedures see e.g. Diaconis and Freedman (1986) and Kim
and Lee (2001) and references therein.
We will have to deal with two different posterior consistency issues, a parametric one for
the posterior of the arrival rate and a non-parametric one for the service time distribution
and values depending on it. The first was considered by Doob (1949) in a rather general way
by the use of martingale theory. It clarifies that, under very weak constraints (the state- and
the parameter spaces are assumed to be Polish spaces and the likelihood is identifiable), there
is a subset of the parameter space with full prior mass such that the sequence of posterior
laws is consistent at any parameter of that subset that is taken as the true one. However,
one problem is that in high-dimensional parameter spaces the set with full prior mass might
become topologically small. Therefore, especially nonparametric Bayesian statistical problems
deserve a deeper study of posterior consistency since the set of possible likelihoods can not
longer be assumed to be parametrized finitely.
First of all, the definition of posterior consistency in a rather general framework is given as
it can be found e.g. in Schervish (1995). Notice that for any Polish state space S the space of
all probability measures, P(S), can be equipped with the topology induced by weak convergence which has neighborhood bases for P ∈ P(S) given by the collection of sets of the form
UP ; = {Q ∈ P(R) ∶ ∣∫ fi dP − ∫ fi dQ∣ < , fi ∈ Cb (R), for all i = 1, .., k}. This topology makes
P(S) itself a Polish space [see Kechris (1995)] and the topology can be metrized e.g. by the
Prohorov metric. The Borel σ-field, B(P), with respect to the weak topology serves as a
natural measure-theoretical structure to turn P(S) into a measurable space (P(S), B(P)).
B(P) turns out to be the smallest σ-field on P(S) that makes the mappings P ↦ P (A)
(B(P), B([0, 1]))-measurable for all A ∈ B(S). In the sequel Π, occasionally with an appropriate index, shall denote a prior distribution and Πn the corresponding posterior distribution
after having seen the first n projections of the exchangeable data. A distribution (or a value
that parametrizes a distribution) indexed by 0 will always denote the true data generating
measure.
Definition 1. Let Π ∈ P(P(R)) be a prior distribution, i.e. the distribution of some random
probability measure P ∈ P(R) and let X1∞ be a sequence of exchangeable data which is conditionally i.i.d. given P . Moreover, let P0 ∈ P(R) be the true data generating distribution and
(Πn )n≥0 = (Π(⋅∣X1n ))n≥0 be the sequence of posterior laws. Then call (Πn )n≥0 weakly consistent
at P0 if for all weak neighborhoods UP0 , of P0 it holds
n→∞

Πn (UP0 , ) Ð→ 1,
for P0∞ -almost all data sequences X1∞ .
13

Here, P0∞ denotes the true joint law governing the sequence X1∞ .
Needless to say that in the case of the arrival rate λ things become easier because, due
to the additional judgment concerning the symmetry of the distribution of the exchangeable
sequence of inter-arrival times, one can reduce the problem to that of parametric consistency
that is rather easy to handle by conjugacy.
Proposition 4.1. Let ΠΓ stand for the Gamma prior of the random arrival rate λ and
(ΠΓ;n (⋅))n≥1 = (ΠΓ (⋅∣An1 ))n≥1 be the sequence of posterior laws. Then, for all  > 0 and for Pλ∞0
almost all sequences A∞
1 it holds
n→∞

ΠΓ;n (∣λ − λ0 ∣ < ) Ð→ 1.
Proof. Taking into account that the gamma prior is conjugate for exponentially distributed
data, the posterior density can be shown to be that of a Γ(a + n, b + ∑ni=1 Ai ) distribution. By
well known properties of the Gamma distribution one has EΓ;n [λ] = (a + n)(b + ∑ni=1 Ai )−1 and
VarΓ;n [λ] = (a + n) (b + ∑ni=1 Ai )−2 . Applying the continuous mapping theorem the assertion
follows from the strong law of large numbers and the Markov inequality.
Let us now turn to the nonparametric part concerning the random c.d.f. G. The following
two results are due to Dey et al. (2003) and give the consistency of the Bayes estimator and
the random c.d.f. G.
Lemma 4.2 (Dey et al. (2003)). Let G0 ∈ F(R+ ) be the true c.d.f. of the service times.
Furthermore suppose that G0 is continuous and that G0 (t) < 1 for all t ∈ R+ . Let G ∼ BS(c, H)
and ΠBS;n denote the posterior distribution of G. Then, for G∞
0 -almost all sequences of data
S1∞ and all t ∈ R+ it holds
n→∞

EBS;n [G(t)] Ð→ G0 (t).
Next, the posterior consistency of the random service time distribution G induced by a
beta-Stacy prior will be stated. Since the service time distribution function G(⋅) is seen to be
a random function rather than a random variable, we investigate deviations in the sup-norm.
Theorem 4.3 (Dey et al. (2003)). Let G0 ∈ F(R+ ) be the true c.d.f. of the service times.
Again, suppose that G0 is continuous and that G0 (t) < 1 for all t ∈ R+ . Let G ∼ BS(c, H) and
ΠBS;n denote the posterior distribution of G. Then, for all  > 0 it holds
n→∞

ΠBS;n ( sup ∣G(t) − G0 (t)∣ < ) Ð→ 1
0≤t<∞

for

G∞
0 -almost

all sequences S1∞ .

Remark 1. We stress some peculiarity of a certain class of neutral to the right measures. For
a full formal treatment see Dey et al. (2003). As already mentioned in the previous section,
neighborhoods w.r.t. the sup-norm contain some weak neighborhoods. These, in turn, are
given as finite intersections of sets of the form {G ∶ ∣G(t) − G0 (t)∣ < γ}. Thus, it is enough to
show that the posterior mass of such sets converges to one. But since the Lévy measure that
corresponds to the beta-Stacy prior is of the form L(dt, ds) = a(t, s)dsK(dt) for suitable a and
K, the convergence of the expected value of the posterior law to the true c.d.f. already ensures
that the posterior variance vanishes with increasing data size.
14

As an immediate consequence one has the uniform consistency of the Bayes estimator.
Corollary 4.4. Let the conditions of Lemma 4.2 be fulfilled. Then the Bayes estimate
EBS;n [G(⋅)] of the service time distribution is uniformly consistent at the true continuous
service time distribution G0 (⋅), that is
n→∞

sup ∣EBS;n [G(t)] − G0 (t)∣ Ð→ 0,
0≤t<∞
∞
for G∞
0 -almost all sequences S1 .

Proof. By above theorem one has
RRR
RRR
R
RRR
sup ∣EBS;n [G(t)] − G0 (t)∣ = sup RRR∫ (G(t) − G0 (t)) ΠBS;n (dG)RRRRR
RRR
0≤t<∞
0≤t<∞ RR
RRF
R
≤
sup ∣G(t) − G0 (t)∣ ΠBS;n (dG)
∫
{G∶sup∣G(t)−G0 (t)∣≥}

0≤t<∞

+

sup ∣G(t) − G0 (t)∣ ΠBS;n (dG)

∫

{G∶sup∣G(t)−G0 (t)∣<}

0≤t<∞

n→∞

< ΠBS;n ( sup ∣G(t) − G0 (t)∣ ≥ ) +  Ð→ .
0≤t<∞

Since  > 0 can be chosen arbitrarily small, the proof is completed and the assertion follows.
The next lemma establishes the uniform consistency on R+ of the service time LST in
posterior law.
Proposition 4.5. Let g(z) = ∫ e−sz dG(s) denote the LST of the random service time distribution G possessing a beta-Stacy process prior and g0 (z) = ∫ e−sz dG0 (s) the LST of the
corresponding true data generating distribution. Then, if the constraints of Theorem 4.3 are
fulfilled, it holds
n→∞

ΠBS;n ( sup ∣g(z) − g0 (z)∣ < ) Ð→ 1,
0≤z<∞
∞
for all  > 0 and for G∞
0 -almost all data sequences S1 .

Proof. The assertion of the proposition follows from Theorem 4.3 and the continuous mapping
theorem applied to the mapping G(⋅) ↦ g(⋅) which is continuous w.r.t. the sup-norm. Indeed,
take a δ > 0 and let UG0 ,δ be a uniform δ-neighborhood of the true service time distribution
and let G ∈ UG0 ,δ with corresponding LST g. Then, by integration by parts of the Riemann-

15

Stieltjes integral it holds that
sup ∣g(z) − g0 (z)∣ = sup ∣∫

0≤z<∞

0≤z<∞

∞
0

e−sz [G − G0 ](ds)∣

= sup ∣[e−sz [G(s) − G0 (s)]∣s=∞ − [e−sz [G(s) − G0 (s)]∣s=0 − ∫

0

0≤z<∞

∞

= sup ∣∫
0≤z<∞

≤ sup ∫
0≤z<∞

0
∞

0

< δ sup ∫

[G(s) − G0 (s)]de−zs ∣ = sup ∣∫
0≤z<∞

∞
0

∞

[G(s) − G0 (s)]de−zs ∣

[G(s) − G0 (s)]ze−zs ds∣

ze−zs ∣G(s) − G0 (s)∣ ds

∞

ze−zs ds = δ,

0

0≤z<∞

which completes the proof and shows the claim.
Next, the consistency of the estimator gn∗ (z) of the LST of the service time distribution
is studied. On basis of the previous result, the following lemma establishes the uniform
consistency on R+ of the estimator gn∗ (z).
Lemma 4.6. Let gn∗ (⋅), g0 (⋅) as above. Then, under the constraints of Theorem 4.3, it holds
sup ∣gn∗ (z) − g0 (z)∣ Ð→ 0,
n→∞

0≤z<∞

∞
for G∞
0 -almost all sequences of data S1 .

Proof. Let  > 0 be arbitrary and let Ĝn (⋅) ∶= EBS;n [G(⋅)]. Then by properties of the RiemannStieltjes integral, one has
sup ∣gn∗ (z) − g0 (z)∣

0≤z<∞

∞

= sup ∣∫
0≤z<∞

0

= sup ∣ [e

e−zs dĜn (s) − ∫

∞

0

−zs

Ĝn (s)]∣

s=∞

0≤z<∞

− [e

−zs

e−zs dG0 (s)∣
s=0

−zs

− [e
= sup ∣∫
0≤z<∞

0

= sup ∣∫
0≤z<∞

∞

0

∞

Ĝn (s)d(e

−sz

∞

Ĝn (s)]∣

∞

)−∫

0

−∫

0

Ĝn (s)de−zs

G0 (s)]∣s=∞ + [e−zs G0 (s)]∣s=0 + ∫

0

−sz

G0 (s)d(e

∞

G0 (s)de−zs ∣

)∣

[Ĝn (s) − G0 (s)]ze−zs ds∣

< sup ∣Ĝn (s) − G0 (s)∣ .
0≤s<∞

Hence, the assertion follows from Corollary 4.4 and the proof is completed.
∞
Since the Bayesian estimate of the random mean of G, i.e. EBS;n [∫0 tdG(t)], is used
to define estimators of several queueing characteristics, its posterior consistency is examined
next. For the most prominent prior process, namely the Dirichlet process a rather general
result is known. This is reviewed briefly. So, let P ∼ Dα be a random probability measure

16

that is distributed according to a Dirichlet prior with finite measure α as parameter. Let
f ∶ S → R be measurable. Then [c.f. Feigin and Tweedie (1989)] ∫ ∣f ∣ dα < ∞ ⇒ ∫ ∣f ∣ dP < ∞
with Dα probability one and EDα [∫ f dP ] = ∫ f dEDα [P ] = ∫ f dα. This fact in combination
with an assumption that the state space is countable makes it easy to show that posterior
consistency of the random measure P induces posterior consistency of the random mean of P .
However, here neither the state space is assumed to be countable nor the random measure is
assumed to possess a Dirichlet prior. Thus the posterior consistency of the random mean of G
which is drawn according to a beta-Stacy process is studied in more depth and an affirmative
result is given.
Lemma 4.7. Suppose G(⋅) is drawn according to a beta-Stacy process with parameters (α, β),
where α is a measure on R+ and β(s) ≥ 1 such that
(1.) ∫R+ [β(s)]−1 α(ds) = ∞
(2.) ∫R+ exp [− ∫0 [β(s)]−1 α(ds)] dt < ∞.
t

Then the posterior expectation of the random mean of G converges G∞
0 -a.s. to the true one.
∞
it
holds
-almost
all
sequences
of
data
X
That is for G∞
1
0
EBS;n [∫

tdG(t)] ∶= EBS [∫

R+

R+

n→∞

tdG(t)∣X1n ] Ð→ ∫

tdG0 (t).
R+

Proof. First of all note that, due to Epifani et al. (2003), the first moment of the random
mean under the prior does exist and is given as
∞

EBS;n [∫

tdG(t)] = ∫ exp [− ∫

R+

0

t
0

⎡
⎤
⎢
α(ds)
β(Xi ) + Mn (Xi ) − 1 ⎥⎥
⎢
]×⎢ ∏
⎥ dt,
β(s) + Mn (s)
⎢{i∶Xi ≤t} β(Xi ) + Mn (Xi ) ⎥
⎣
⎦

where Mn (s) = ∑nj=1 δXj [s, ∞). Thus, letting Nn (s) = ∑ni=1 δXi [0, s), one obtains
EBS;n [∫
∞

tdG(t)]

R+
t

= ∫ exp [− ∫
0

0

t
β(s) + Mn (s)
α(ds)
] exp [− ∫ log (
) dNn (s)] dt.
β(s) + Mn (s)
β(s) + Mn (s) − 1
0

Since β(s) ≥ 1, by elementary properties of the logarithm, it follows
1
β(s) + Mn (s)
1
≤ log (
)≤
,
β(s) + Mn (s)
β(s) + Mn (s) − 1
β(s) + Mn (s) − 1
which in turn implies
∞

t

∫ exp [− ∫
0

0

tdG(t)]

≤ EBS;n [∫
∞

R+

≤ ∫ exp [− ∫
0

t
Nn (ds)
α(ds)
] × exp [− ∫
] dt
β(s) + Mn (s)
0 β(s) + Mn (s) − 1

t
0

α(ds) + Nn (ds)
] dt.
β(s) + Mn (s)
17

such that it remains to show that the bounding terms converge to the mean of the true c.d.f..
From a straight-forward application of the monotone convergence theorem and by continuity
of the exponential function, it follows
∞

lim EBS;n [∫

n→∞

R+

∞

t

= ∫ exp [− ∫
0

0

lim

0 n→∞

0

1/nNn (ds)
] dt
1/nβ(s) + 1/nMn (s)

G0 (ds)
] dt.
1 − G0 (s)

G0 (ds)
1−G0 (s)

Now, note that

t

tdG(t)] = ∫ exp [− ∫

= λ(s)ds = P(s < S ≤ s + ∆s∣S > s), where λ(s) = lim∆s→0
t

t G (ds)

P(s<S≤s+∆s)
∆s[1−G0 (s)]

0
is the cumulative
usually denotes the hazard function. Hence Λ(t) = ∫0 λ(s)ds = ∫0 1−G
0 (s)
hazard until t ≥ 0. On the other hand Λ(t) = − log(1 − G0 (t)). Thus,

t

exp [− ∫

0

t
G0 (ds)
] = exp [− ∫ λ(s)ds] = exp [−Λ(t)] = exp [log(1 − G0 (t))] = 1 − G0 (t),
1 − G0 (s)
0

which in turn completes the proof.
The next lemma establishes the contraction of the mass of the posterior law of the random
mean.
Lemma 4.8. If, in addition to the assumptions of the previous lemma, the parameter (α, β)
meets the condition
√

∫

exp [− ∫

R+

t

[β(s)]−1 α(ds)] dt < ∞

0

then the posterior variance of the random mean ∫ tG(dt) vanishes G∞
0 -a.s. as the size of the
∞
∞
data increases. That is, for G0 -almost all sequences of data X1 it holds
VBS;n [∫
R+

n→∞

tG(dt)] ∶= VBS [∫
R+

tG(dt)∣X1n ] Ð→ 0.

Proof. Due to Lemma 4.7 and the continuous mapping theorem the claim follows if the second
moment of the random mean under the posterior converges a.s. to the square of the mean of
the true distribution function G0 (⋅). From Epifani et al. (2003) the second moment under the
prior law exists and the second moment under the posterior law of the random mean is given

18

by
∞

EBS;n [(∫

2

tdG(t)) ]

0

∞ ∞

s
α(dx)
α(dx)
] exp [− ∫
]
β(x) + Mn (x) + 1
0 β(x) + Mn (x)

r

= 2 ∫ ∫ exp [− ∫

0

r

0

×
∞

r

= 2 ∫ exp [− ∫

0

0

⎛
β(Xj ) + Mn (Xj ) − 1 ⎞
β(Xi ) + Mn (Xi ) ⎞ ⎛
dsdr
∏
∏
⎝{i∶Xi ≤r} β(Xi ) + Mn (Xi ) + 1 ⎠ ⎝{j∶Xj ≤s} β(Xj ) + Mn (Xj ) ⎠

r
α(dx)
β(x) + Mn (x)
] exp [− ∫ log (
) dNn (x)]
β(x) + Mn (x) + 1
β(x) + Mn (x) + 1
0

∞

s

∫ exp [− ∫
r

0

s
α(dy)
β(y) + Mn (y)
] exp [− ∫ log (
) dNn (y)] dsdr
β(y) + Mn (y)
β(y) + Mn (y) − 1
0

Hence, by similar arguments as in the proof of the previous lemma and using Fubini’s
theorem, it follows that for G∞
0 −almost all sequences of data
∞

EBS;n [(∫

0

2

∞

n→∞

r

tdG(t)) ] Ð→ 2 ∫ exp [− ∫
0

0

∞

s dG (y)
dG0 (x)
0
] ∫ exp [− ∫
] dsdr
1 − G0 (x)
0 1 − G0 (y)
r

∞ ∞

r

= 2 ∫ ∫ exp [− ∫
0

0

r

s dG (y)
dG0 (x)
0
] exp [− ∫
] dsdr.
1 − G0 (x)
0 1 − G0 (y)

(#)

Furthermore, a straight-forward application of Fubini’s theorem yields
∞

∞ ∞

s

∫ ∫ g(r, s)drds = ∫ ∫ g(r, s)dsdr,
0

0

r

0

where g is an arbitrary integrable function g ∶ [0, ∞)2 → R+ . Now, let f (r, s) denote the
integrand of the last integral of equation (#). Since f is symmetric, i.e. f (r, s) = f (s, r), it
follows from above equality and further application of Fubini’s theorem
∞ ∞

∞

∞ ∞

s

2 ∫ ∫ f (r, s)dsdr = ∫ ∫ f (r, s)drds + ∫ ∫ f (r, s)dsdr
r
∞ r

0

0

0
∞ ∞

r
∞ ∞

0

= ∫ ∫ f (s, r)dsdr + ∫ ∫ f (r, s)dsdr = ∫ ∫ f (r, s)drds.
0

0

0

r

0

0

Therefore,
∞

EBS;n [(∫

0

2

n→∞

tdG(t)) ] Ð→ (∫

0

∞

2

tdG0 (t)) ,

which completes the proof of the assertion.
19

Theorem 4.9. Under the assumptions of Lemma 4.7 and Lemma 4.8 the mean of the random
c.d.f. G possesses the property of posterior consistency. That is for all  > 0 and G∞
0 -almost
all sequences of data it holds
ΠBS (∣∫

n→∞

tG(dt) − ∫
R+

R+

tG0 (dt)∣ > ∣X1n ) Ð→ 0.

Proof. This is a direct consequence of the two previous lemmas and the Markov inequality.
Remark 2. A brief discussion of the assumptions of above lemmas and theorem, respectively,
is given. Assumption (1.) of Lemma 4.7 is an artifact of survival analysis and ensures that the
prior (α, β), that is to be chosen, actually leads to a cumulative hazard function. Assumption
(2.) of the first lemma and the additional assumption of the second lemma, respectively,
ensures the existence of the first, resp. second, moment of the posterior distribution of the
random mean. These conditions are given in Epifani et al. (2003) Proposition 4.
Moreover, this work extends results which concern the existence of certain random functionals w.r.t. to a random c.d.f. drawn according to a NTR prior. Those results can be
understood as a generalization of a work by Feigin and Tweedie (1989) which gives conditions
under which certain functionals of a random measure drawn according to a Dirichlet prior
exist in terms of the base measure (the prior parameter). At this place it should be emphasized
that they investigated this problem by creating a new approach to the Dirichlet process which is
often not mentioned in the literature. To be more precise, they show that the Dirichlet process
can be extracted as the invariant distribution of a measure-valued Markov chain and exploit
this theory to show sufficiency of the conditions of their existence theorem. Since Dirichlet
priors are as well included in the family of NTR priors the extension by Epifani et al. seems
natural.
The last assumption, i.e. that β(s) ≥ 1 ensures posterior consistency of the random mean
of the c.d.f. governed by the beta-Stacy prior and is used in the proofs of both lemmas. Hence,
a prior for the random c.d.f. G is to be chosen, in terms of its parameter, such that the
sequence of posteriors of the G-mean is consistent as long as this property is desired. It is
interesting to see how advanced information on the posterior consistency of the random mean
influences the prior knowledge of the random c.d.f. itself. To put it another way, information
concerning the posterior consistency of the random mean is indeed prior information on G.
However, the other way around is not true in general, i.e. posterior consistency of the random
c.d.f. does not generally imply posterior consistency of the random mean. That is because this
functional is not continuous. However, posterior consistency of G implies that of its truncated
mean. That means integration is restricted to a compact set and would lead to a weaker form
of consistency in the sense of compact convergence.
As an immediate consequence of the previous results one obtains the consistency of the
traffic intensity of the M /G/1 queue.
Corollary 4.10. Let µ ∶= ∫ sG(ds), µ0 ∶= ∫ sG0 (ds) be the random and true mean of the
service-time distribution and ρ0 = λ0 µ0 the true traffic intensity. Further, let ΠBS⊗Γ =
ΠBS ⊗ ΠΓ be the prior on (P(R+ ), R+ ) that is formed by taking the product-measure of ΠBS
and ΠΓ , ΠBS⊗Γ;n the posterior, respectively, and define λ̂n ∶= Eλ;n [λ] ∶= EΓ [λ∣An1 ], µ̂n ∶=
∞
∞
EBS;n [∫0 tdG(t)] ∶= EBS [∫0 tdG(t)∣S1n ] and ρ̂n ∶= EBS⊗Γ;n [ρ] ∶= EBS⊗Γ [ρ∣An1 , S1n ].
Then, under the assumptions of Proposition 4.1 and Theorem 4.9, one has
20

(i) for all  > 0
n→∞

ΠBS⊗Γ (∣ρ − ρ0 ∣ ≥ ∣S1n , An1 ) Ð→ 0,
∞
∞
for Pλ∞0 ⊗ G∞
0 -almost all data sequences (A1 , S1 ),

(ii)
n→∞

∣ρ̂n − ρ0 ∣ Ð→ 0,
∞
∞
for Pλ∞0 ⊗ G∞
0 -almost all data sequences (A1 , S1 ).

Proof. (i) Recalling the posterior consistency of the arrival rate λ in Proposition 4.1 and the
service rate µ in Theorem 4.9 one has
ΠBS⊗Γ;n (∣ρ − ρ0 ∣ ≥ ) = ΠBS⊗Γ;n (∣µλ − µλ0 + µλ0 − λ0 µ0 ∣ ≥ )
≤ ΠBS⊗Γ;n (∣µλ − µλ0 ∣ + ∣µλ0 − λ0 µ0 ∣ ≥ )
≤ ΠBS⊗Γ;n (∣µλ − µλ0 ∣ ≥ /2) + ΠBS⊗Γ;n (∣µ − µ0 ∣ ≥


)
2λ0

≤ ΠBS⊗Γ;n ((µ0 + δ)∣λ − λ0 ∣ ≥ /2, ∣µ − µ0 ∣ < δ)
+ ΠBS;n (∣µ − µ0 ∣ ≥ δ) + ΠBS;n (∣µ − µ0 ∣ ≥
≤ ΠΓ;n (∣λ − λ0 ∣ ≥


2(µ0 + δ)


)
2λ0

) + ΠBS;n (∣µ − µ0 ∣ ≥ δ) + ΠBS;n (∣µ − µ0 ∣ ≥


),
2λ0

from which the assertion of (i) follows.
(ii) Straightforward one has a.s.
∣ρ̂n − ρ0 ∣ ≤ ∣µ̂n λ0 − µ̂n λ̂n ∣ + λ0 ∣µ̂n − µ0 ∣
≤ (µ0 + γ)∣λ0 − λ̂n ∣ + λ0 ∣µ̂n − µ0 ∣
by selecting γ > 0 such that a.s. ∣µ̂n − µ0 ∣ < γ, for all n sufficiently large. The assertion of (ii)
then follows from the proof of Proposition 4.1 and Lemma 4.7 completing the proof.
We are now in a position to state the main theorem of this section, i.e. the consistency of
the estimators for the waiting time LST, the queue length p.g.f. and the system size p.g.f.,
defined in section 3.3.
Theorem 4.11. Under the assumptions of Proposition 4.1, Proposition 4.5 and Theorem 4.9
one has
(i) for all  > 0
n→∞

ΠBS⊗Γ ( sup ∣f (z) − f0 (z)∣ ≥ ∣S1n , An1 ) Ð→ 0.
0≤z<∞

(ii)
sup ∣fn∗ (z) − f0 (z)∣ Ð→ 0,
n→∞

0≤z<∞

for

Pλ∞0

∞
∞
⊗ G∞
0 -almost all data sequences (A1 , S1 ), where f (⋅) ∈ {n(⋅), q(⋅), w(⋅)}.

21

Proof. We show the result for f = w, that is in the case of the p.g.f. of the waiting-time
distribution. The other cases are treated similarly and we omit the details. We intent to
apply the continuous mapping theorem. Therefor note that the mapping (g(⋅), λ, µ) ↦ w(⋅)
is continuous with respect to the suitable topologies induced by the sup-norm. Hence, the
assertion of the theorem is implied by the results of the present section.
We shall stress that above theorem is important, especially for applications, since it states
that one can reduce the difficult task of making Bayesian inference for queueing characteristics
to that of making inference for the observables separately in order to compose it subsequently
to that of the objects of interest. Qualitatively this will lead to reasonable asymptotic inference. This is appreciable because it is a non-feasible task to place a prior on the distributions of
aforementioned characteristics in a way that it can be updated with data given by observations
that we have access to.

22

5

Posterior Normality

Asymptotic normality results, so called Bernstein-von Mises theorems, for the posterior law
serve as an additional validation of Bayesian procedures. Especially in situations where the
exact posterior law is not available or hard to compute they are useful from a rather applied
viewpoint to get an approximation of the posterior law. While the previous section has shown
that the posterior concentrates around the true value when data increases, normality results
give an idea how, asymptotically speaking, it does concentrate and how fluctuations around
this centering appear. Often it is true that the posterior law, if centered and rescaled appropriately, resembles a centered Gaussian distribution with a certain covariance structure. For
applications, special interest lies in this limiting covariance structure.
The earliest result dates back to de Laplace (1774) who approximates the posterior of
a beta distribution by a normal integral. Sergei Bernstein and Richard von Mises gave a
rather modern version of this approach by applying a general result about infinite products
of functions to the sequence of posterior densities, see e.g. Johnson (1967) for a review and
references. These approaches of expanding posterior densities coined the name. Nowadays,
the parametric case is well understood and results including centering with the MLE or the
Bayes estimate as well can be found in Schervish (1995) or Ghosh and Ramamoorthi (2003).
However, the results needed for an asymptotic behavior of the posterior of infinite-dimensional
parameters, i.e. in the nonparametric case, go deeper and deserve separate investigation. Posterior normality results involving Dirichlet process priors or NTR processes can be found in
Conti (1999) who employs a result by Freedman (1963) on asymptotic normality in the finitedimensional case.
We use a general result by Kim and Lee (2004) on the asymptotic normality of NTR
processes to obtain the asymptotic behavior of the posterior law of the waiting time LST.
This result is stated next. For the sake of clarity, it is stated for the non-censored case.
However, it can be enlarged to the situation where data is right-censored. For a positive real
number τ let (D[0, τ ], ∥⋅∥τ ) denote the space of all cadlag functions on [0, τ ] equipped with
the sup-norm and L[X] the law of a random object X.
Theorem 5.1 (Kim and Lee (2004)). Suppose that G(⋅) = 1 − exp[A(⋅)] is a random c.d.f.
drawn according to a NTR prior Π with corresponding increasing additive process A(⋅) whose
Lévy measure L is given by
t

L([0, t], B) = ∫ ∫
0

B

gs (x)
dxλ(s)ds,
x

1

where ∫0 gt (x)dx = 1 for all t ∈ R+ . Let the true c.d.f. G0 be continuous and such that
G0 (t) < 1 for all t ∈ R+ . Moreover let the following conditions be fulfilled for a τ > 0,
•

sup

(1 − x)gt (x) < ∞,

t∈[0,τ ],x∈[0,1]

• there is a function q(t) ∶ R+ → R such that 0 < inf q(t) < sup q(t) < ∞ and, for some
t∈[0,τ ]

23

t∈[0,τ ]

α > 1/2 and  > 0,
sup

t∈[0,τ ],x∈[0,]

∣

gt (x) − q(t)
∣ < ∞,
xα

• λ(t) is bounded and positive on (0, ∞).
∞
Then it holds that for G∞
0 -almost all sequences of data S1
√
lim Π ( n[G(⋅) − EΠ;n [G(⋅)]]∣S1n ) = L[A(⋅)]
n→∞

weakly on (D[0, τ ], ∥⋅∥τ ), where A(⋅) denotes a centered Gaussian process given by A(t) =
(G0 (t) − 1)W (A0 (t)) with a Brownian motion W. The covariance structure then is h(u, v) ∶=
Cov[A(u), A(v)] = (1 − G0 (u))(1 − G0 (v))min(A0 (u), A0 (v)).
Proof. See Kim and Lee (2004) for the proof as well as for a discussion of the constraints of
the theorem.
Above theorem can well be stated as follows. The posterior distribution of the scaled
and centered random process looks more and more like a Gaussian process as the sample
size increases. The limiting process is centered, i.e. the expectation function is the constant
function taking only the value zero. Next result ensures that above assertions hold for betaStacy processes.
Corollary 5.2 (Kim and Lee (2004)). The assertion of Theorem 5.1 for the random c.d.f. G
holds, with G being governed by a beta-Stacy process prior with parameters (c(⋅), H(⋅)), where
H is a continuous distribution function with continuous density a with respect to the Lebesgue
measure such that 0 < inf t c(t)[1 − H(t)] < supt c(t)[1 − H(t)] < ∞, i.e. for G∞
0 -almost all
sequences of data S1∞ it holds
√
lim ΠBS;n ( n[G(⋅) − Ĝn (⋅)]) = L[A(⋅)]
n→∞

weakly on (D(1) [0, ∞), ∥⋅∥∞ ) with A(⋅) as in Theorem 5.1, where D(1) [0, ∞) denotes the space
of cadlag functions bounded by one on [0, ∞).
Proof. By Theorem 5 of Dey et al. (2003), a beta-Stacy process is a transformed beta process.
More precisely, a process Λ is a beta-Stacy process with parameters (c(⋅), H(⋅)) if and only
⋅ dH(s)
if it is a beta process with parameters (c(⋅)[1 − H(⋅)], ∫0 1−H(s) ). Further on, Kim and Lee
(2004) showed that a beta process fulfills the constraints of the theorem as long as
positive and continuous on R+ and c(t)[1 − H(t)] is as required.

a(s)
1−H(s)

is

Next, we turn to the asymptotic behavior of the LST g of the service time distribution
G and its random mean µ ∶= ∫ xdG(x). Since a plug-in estimator for g is employed, we will
use the functional delta method to provide the asymptotic normality result. Recall that in
the previous section a posterior consistency result was given for the random mean µ, i.e. it
was shown that the posterior law of µ centers a.s. around its true value µ0 ∶= ∫ xG0 (x). This
suggests the conjecture that a normality result might be present in this case as well. However,
even if Regazzini et al. (2003) provide results to approximate the density of the random mean
of random measures, we were not able to show the density in our case to be approximated
24

by a normal density if the sample size increases. The main reason for this is that it does not
seem to be straight-forward to obtain a suitable expansion of this approximated density.
Hence, in the following we make the rather technical assumption that there is prior knowledge of the kind that service times can not exceed a certain sufficiently large threshold M ∈ R+ .
From a practical point of view this is a rather gentle constraint. Let FM ∶= FM (R+ ) ∶= {F ∈
F(R+ ) ∶ F (t) = 1, ∀t ≥ M } be the space of all c.d.f.’s whose corresponding probability measure has support [0, M ]. Since it is well known that the prior guess on the c.d.f. G under a
beta-Stacy prior BS(c, H) is given as EBS;n [G] = H, in the following we take H ∈ FM such
that H is continuous on [0, M ]. Recall that G(⋅) = 1 − exp[−A(⋅)], where A is a non-negative
increasing additive process with Lévy measure
t

dx
−xc(s)[1−H(s)]
dNt (x) =
c(s)H(ds).
∫ e
1 − e−x
0

From the theory of increasing additive processes it is well known [see e.g. Sato (1999)] that
for all t > 0 A(t) = At is a random variable governed by a infinitely divisible distribution φt .
Let φ̂t (ξ) denote the characteristic function of φ, i.e.
φ̂t (ξ) = ∫ exp[iξs]φ(ds) = exp [− ∫

∞
0

(1 − eiξx )dNt (x)] .

Furthermore, since (At ) has independent increments one has φ̂t (ξ) = φ̂s (ξ)φ̂s,t (ξ) for all s < t,
where φ̂s,t (ξ) denotes the characteristic function of the difference At − As . Now, since H ∈ FM
it follows for t > M,

φ̂t (ξ) = exp [− ∫
= exp [− ∫

∞
0

∞
0

(1 − eiξx ) dNt (x)] = exp [− ∫

∞
0

t
1 − eiξx
−xc(s)[1−H(s)]
c(s)H(ds)dx]
∫ e
−x
1−e
0

iξx

M
1−e
e−xc(s)[1−H(s)] c(s)H(ds)dx]
∫
−x
1−e
0

⋅ exp [− ∫

0

∞

t
1 − eiξx
−xc(s)[1−H(s)]
c(s)H(ds)dx]
∫ e
−x
1−e
M

= φ̂M (ξ),
the process (At )t≥M is a.s. constant. That implies that the corresponding c.d.f. G is constant
from M onwards. In order to ensure that it is indeed a distribution function, we set G(t) ∶= 1,
(M )
for all t ≥ M . The support of the truncated prior ΠBS law will still be all of FM .
In order to achieve a posterior normality result for the random LST g, we show that the
∞
mapping Φ ∶ D(1) [0, M ] → C (1) [0, ∞); G ↦ ∫0 e−sz dG(s) is Hadamard differentiable, where
D(1) [0, M ] and C (1) [0, ∞) denote the space of cadlag and continuous functions, bounded by
one, respectively. Moreover, since we are solely interested in distribution functions on R+ ,
w.l.o.g. it is assumed that D(1) [0, M ] consists only of functions starting at zero. For a good
reference on Hadamard differentiability and applications including the functional delta method
see Kosorok (2008).

25

Lemma 5.3. The mapping
Φ ∶ (D(1) [0, M ], ∥.∥M ) → (C (1) [0, ∞), ∥.∥∞ )
M

G ↦ Φ[G](●) ∶= ∫

e−●s G(ds)

0

is Hadamard differentiable.
t↘0

Proof. Let t ↘ 0 and ht ∈ D(1) [0, M ], such that ht Ð→ h ∈ D(1) [0, M ] w.r.t. the sup-norm.
Then, by properties of the Riemann-Stieltjes integral, one has
∥

M
Φ[G + tht ](z) − Φ[G](z)
−∫
e−zs dh(s)∥
t
0
∞
M

[ht (s) − h(s)](−z)e−zs ds∣

M

( sup ∣ht (x) − h(x)∣) (−z)e−zs ds∣

= sup ∣∫
0≤z<∞

0

≤ sup ∣∫
0≤z<∞

0

0≤x<M
M

≤ ∥ht − h∥M sup ∫
0≤z<∞
t↘0

ze−zs ds

0

= ∥ht − h∥M Ð→ 0.
Define
Φ′G ∶ D(1) [0, M ] → C (1) [0, ∞)
M

h ↦ Φ′G [h](z) = ∫

e−zs dh(s).

0

Since the Riemann-Stieltjes integral is linear in the integrator, the mapping Φ′G is linear.
Moreover, if sup ∣F (s) − G(s)∣ < δ, it follows
0≤s<M

RRR
RRR M
M
M
RRR
RRR
−zs
−zs
RRR∫ e dF (s) − ∫ e dG(s)RRR ≤ ∥F − G∥M sup ∫ ze−zs ds < δ,
RRR
RRR
0≤z<∞
0
0
R
R0
′
thus the continuity of ΦG . Hence the mapping Φ is Hadamard differentiable with derivative
Φ′G .
Now, the lemma will be applied in combination with the functional delta method to obtain
the posterior normality of the service time LST centered suitably at its respective Bayes
(M )
estimator. Write ΠBS;n for the posterior law induced by the M -truncated beta-Stacy process.
Corollary 5.4. Let gn∗M (z) = ∫0 e−zs dEBS;n [G](s). Under the assumptions of Theorem 5.1
∞
and Corollary 5.2, it holds for G∞
0 -almost all sequences S1
(M ) √
lim ΠBS;n ( n [g(⋅) − gn∗M (⋅)]) = L[G(⋅)]
M

n→∞

on C (1) ([0, ∞), ∣∣.∣∣∞ ), where G(z) is a centered Gaussian process with covariance structure
γ(⋅, ⋅) given by
M

γ(u, v) = Cov [G(u), G(v)] = uv ∫

0

M

∫

e−(us+vt) h(u, v)dudv,

0

where h(⋅, ⋅) is defined in Theorem 5.1.
26

Proof. By the functional delta method applied to the mapping in the previous lemma one has
G(z) = Φ′G [(G0 (⋅) − 1)W (A0 (⋅))] (z) = ∫

M
0

e−zs d [(G0 (s) − 1)W (A0 (s)] .

Using a Riemann sum approximation for above integral, one concludes that the process G(⋅)
is a Gaussian process. Further, by well-known properties of the Riemann-Stieltjes integral,
one gets
M

G(z) = z ∫

0

(1 − G0 (s))W (A0 (s)) e−sz ds.

Using Fubini’s theorem it is immediately seen that E[G(z)] = 0 for any z ∈ R+ . Again using
Fubini’s theorem, the covariance structure of G(⋅) is obtained as
Cov [G(u), G(v)] = E [G(u)G(v)]
M

= uvE [∫

0
M

= uv ∫

0
M

∫

0

0

M

= uv ∫

0

M

∫

M

∫

(1 − G0 (s))W (A0 (s))e−us (1 − G0 (t))W (A0 (t))e−vt dsdt]

e−(us+vt) (1 − G0 (s))(1 − G0 (t))E [W (A0 (s))W (A0 (t))] dsdt
e−(us+vt) h(s, t)dsdt.

0

Next, we investigate the posterior normality of the mean of the random c.d.f. G. Since
M
no exact results seem obtainable, we use the plug-in estimator µ∗n ∶= ∫0 [1 − EBS;n [G(t)]]dt
M
which, in general, does not equal EBS;n [∫0 [1 − G(t)]dt]. This estimator in combination with
the M -truncated c.d.f.’s enables us to use the functional delta method for obtaining normality
results.
Lemma 5.5. Let M be an arbitrary positive real number. Then, the mapping
Ψ ∶ (FM , ∥⋅∥M ) → ([0, M ], ∣⋅∣)
M

G ↦ Ψ[G] ∶= ∫

[1 − G(s)]ds
0

is Hadamard-differentiable with derivative Ψ′ [h] = − ∫0 h(s)ds.
M

t↘0

Proof. Take t ↘ 0 and ht ∈ FM , such that ht Ð→ h ∈ FM . Then
∣

M
M
Ψ[G + tht ] − Ψ[G]
t↘0
+∫
h(s)ds∣ = ∣∫
h(s) − ht (s)ds∣ ≤ M ∥h − ht ∥M Ð→ 0.
t
0
0

Obviously, the derivative of Ψ is linear and continuous w.r.t. to the considered topologies.
Corollary 5.6. Under the assumptions of Theorem 5.1 and Corollary 5.2 it holds for G∞
0
almost all data S1∞ that
(M ) √
lim ΠBS;n ( n[µ − µ∗n ]) = L[H],

n→∞

M

M

where H is a centered Gaussian random variable with variance η ∶= Var[H] = ∫0 ∫0 h(s, t)dsdt.
27

Proof. By the previous lemma and the functional delta method the limiting variable is given by
M
Ψ′ [(G0 (⋅) − 1)W (A0 (⋅))] = ∫0 (1 − G0 (s))W (A0 (s))ds which is seen to be centered Gaussian
by a Riemann sum approximation in combination with Fubini’s theorem. Moreover, again by
Fubini’ theorem
Var[H] = E [H2 ] = E [∫
M

=∫

0

M
0

M

∫

0

(1 − G0 (s))W (A0 (s))(1 − G0 (t))W (A0 (t))dsdt]

M

∫

h(s, t)dsdt.
0

Next we consider the asymptotic normality of the arrival-rate when centering with its
Bayes estimate.
Proposition 5.7. Let λ̂n ∶= EΓ;n [λ]. Then, for Pλ∞0 -almost all sequences of data A∞
1 one has
√
lim ΠΓ ( n [λ − λ̂n ] ∣An1 ) = N ,
n→∞

where N is a centered Gaussian random variable with precision λ20 .
Proof. By Theorem 1.4.3. in Ghosh and Ramamoorthi (2003) the convergence of the posterior
√
distribution of n [λ − λ̂n ] to a centered normal distribution directly follows. Moreover, the
variance of this limiting Gaussian variable is given by the inverse Fisher information at the
true arrival rate. Checking the necessary conditions for interchanging integral and derivative
is left to the interested reader. The Fisher information is obtained as
I(λ0 ) = Eλ0 [(

2
∂
∂2
= λ−2
log (λe−λA )) ] = −Eλ0 [ 2 [log(λ) − λA]]
0 .
∂λ
∂λ
∣λ=λ
0

We are now in a position to formulate the posterior normality of the waiting-time LST.
However, the same techniques can be applied to show posterior normality of several other
queueing characteristics like e.g. for the queue length p.g.f. or the sojourn time LST. The
waiting-time distribution is of special interest since it gives a qualitative idea about the loss of
information that can occur in a M /G/1 system and thus helps to ensure a well-working system.
However, since the exact posterior law of the waiting time distribution is not obtainable in
closed form, asymptotic approximations are given. These results extend results of the previous
section where it was shown that the plug-in estimator is reasonable to make inference. Roughly,
we prove that the posterior law of the LST follows, asymptotically speaking, a Gaussian
(M )
quantity that is centered around the plug-in estimator. Let ΠBS⊗Γ = ΠBS ⊗ ΠΓ denote the
prior on the parameter space R+ × FM .
z(1−ρ)
z−λ(1−g(z))
z(1−λ̂n µ∗n )
∗M (z))
z−λ̂n (1−gn

Theorem 5.8. Let w(z) =

be the LST of the waiting time distribution as given in

section 2 and wn∗M (z) =

be its plug-in estimator. Then, under the assumptions

∞
∞
of Theorem 5.1 and Corollary 5.2, for G∞
0 ⊗ Pλ0 -almost all sequences of data (S, A)1 it holds
that
√
lim ΠBS⊗Γ ( n [w(⋅) − wn∗M (⋅)] ∣An1 , S1n ) = L[Z(⋅)]
n→∞

28

weakly on C([0, ∞), ∣∣.∣∣∞ ), where Z(⋅) is a centered Gaussian process with covariance structure
ζ(u, v) = Cov[Z(u), Z(v)] given by
ζ(u, v) =

µ20
w0 (u)w0 (v) 2
w0 (u)(1 − g0 (u)) w0 (v)(1 − g0 (v))
[λ
η
+
×
] + λ−2
0
0
2
2
(1 − ρ0 )
λ0 (1 − g0 (u)) − u λ0 (1 − g0 (v)) − v
λ0

w0 (u)w0 (v)λ20
γ(u, v)
[λ0 (1 − g0 (u)) − u] [λ0 (1 − g0 (v)) − v]
µ0 w0 (u)w0 (v) w0 (u)(1 − g0 (u)) w0 (v)(1 − g0 (v))
−
[
+
]
[λ0 (1 − ρ0 )]2
u
v
⎡
⎤
⎥
λ0 w0 (u)w0 (v) ⎢⎢
−us
−tv
⎥,
w
(u)
e
h(s,
t)d(s,
t)
+
w
(v)
e
h(s,
t)d(s,
t)
−
0
0
∫
∫
⎢
⎥
2
(1 − ρ0 )
⎢
⎥
⎣
⎦
[0,M ]2
[0,M ]2

+

where γ(⋅, ⋅) is given in Corollary 5.4.
Proof. Taking a similar route of proving as in the proof of Theorem 3 in Conti (1999), we
√
begin with a decomposition of n [w(z) − wn∗M (z)]. The decomposition yields
√
−w(z) √
n [w(z) − wn∗M (z)] =
n [λµ − λ̂n µ∗n ]
1−ρ
√
1
1
+ z(1 − λ̂n µ∗n ) n [
]
−
z − λ(1 − g(z)) z − λ̂n (1 − g(z))
√
1
1
)
+ z(1 − λ̂n µ∗n ) n (
−
z − λ̂n (1 − g(z)) z − λ̂n (1 − gn∗M (z))
=∶ Z1;n (z) + Z2;n (z) + Z3;n (z)
Now, the three terms of the sum are investigated separately and it will be shown that they
possess the same asymptotic distribution as objects whose asymptotic is easier to obtain.
These objects will be tagged by an additional ∗ superscript. First write
Z1;n (z) =

−w(z) √
n [µ∗n (λ − λ̂n ) + (µ − µ∗n )(λ − λ0 ) + λ0 (µ − µ∗n )]
1−ρ

and note that by the uniform posterior consistency results of the previous section and the
M
continuity of the mapping G ↦ ∫0 (1 − G(s))ds one has
lim ΠBS⊗Γ [Z1;n (z)∣An1 , S1n ] = lim ΠBS⊗Γ [

n→∞

=∶

n→∞

∗
lim ΠBS⊗Γ [Z1;n
(z)∣An1 , S1n ]
n→∞

−w0 (z) √
n(µ0 (λ − λ̂n ) + λ0 (µ − µ∗n ))∣An1 , S1n ]
1 − ρ0

=∶ L[Z1 (z)].

For Z2;n (z), note that the mapping λ ↦ [z − λ(1 − g(z))]−1 is analytic in a suitably chosen
neighborhood of λ0 . Its derivative is given by
λ↦

1 − g(z)
[z − λ(1 − g(z))]2

.

29

Thus, a Taylor expansion of that mapping yields
Z2;n (z) = (1 − λ̂n µ∗n )z

1 − g(z)

√
2

n (λ − λ̂n ) ,

[z − λ̄(1 − g(z))]

for a suitably chosen λ̄ ∈ [λ0 , λ̂n ]. Therefore, using consistency results and continuous mapping, one gets
lim ΠBS⊗Γ [Z2;n (z)∣S1n , An1 ] = lim ΠΓ [(1 − ρ0 )z

n→∞

n→∞

1 − g0 (z)
[z − λ0 (1 − g0 (z))]

√
2

n (λ − λ̂n ) ∣An1 ]

∗
∣An1 ] =∶ L[Z2 (z)].
=∶ lim ΠΓ [Z2;n
n→∞

−1

Another Taylor expansion for the mapping x ↦ [z − λ̂n (1 − x)]
as before yields
(M )

lim ΠBS⊗Γ [Z3;n (z)∣S1n , An1 ] = lim ΠBS [−
n→∞

n→∞

and analogous reasoning

√
w0 (z)λ0
n [g(z) − gn∗M (z)] ∣S1n ]
λ0 (1 − g0 (z)) − z

(M )

∗
∣S1n ] =∶ L[Z3 (z)].
=∶ lim ΠBS [Z3;n
n→∞

Now, the convergence of the posterior law of the waiting time LST follows from the previous
results of the present section. What remains is the calculation of the covariance structure.
This, in turn, is easily obtained by above decomposition and the assumed independence of the
arrivals and services or their prior laws, respectively.
3

3

3

Cov [∑ Zi (u), ∑ Zi (v)] = ∑ Cov [Zi (u), Zi (v)] + ∑ Cov [Zi (u), Zj (v)]
i=1

i=1

i=1

i≠j

3

= ∑ Cov [Zi (u), Zi (v)]
i=1

+ Cov [Z1 (u), Z2 (v)] + Cov [Z2 (u), Z1 (v)]
+ Cov [Z1 (u), Z3 (v)] + Cov [Z3 (u), Z1 (v)]
By the previous results of this section it follows
3

µ20
w0 (u)w0 (v) 2
[λ
η
+
]
0
(1 − ρ0 )2
λ20
w0 (u)(1 − g0 (u)) w0 (v)(1 − g0 (v))
+ λ−2
×
0
λ0 (1 − g0 (u)) − u λ0 (1 − g0 (v)) − v
w0 (u)w0 (v)λ20
+
γ(u, v).
[λ0 (1 − g0 (u)) − u] [λ0 (1 − g0 (v)) − v]

∑ Cov [Zi (u), Zi (v)] =
i=1

Furthermore, by the independence assumption of the prior laws of the inter-arrival rate and
the service time distribution, one has
Cov [Z1 (u), Z2 (v)] =

µ0
w0 (u)w0 (v)(g0 (v) − 1)
×
.
v − λ0 (1 − g0 (v))
λ20 (1 − ρ0 )
30

has

Furthermore, using the previous results of the present section and Fubini’s theorem, one

w0 (u)w0 (v)λ20
E[H × G(v)]
(1 − ρ0 )(λ0 (1 − g0 (v)) − v)
w0 (u)w0 (v)λ20
=
(1 − ρ0 )(λ0 (1 − g0 (v)) − v)

Cov [Z1 (u), Z3 (v)] =

M

× E [∫

0

M

(1 − G0 (s))W (A0 (s))ds × v ∫

0

(G0 (t) − 1)W (A0 (t))e−tv dt]

M
M
vw0 (u)w0 (v)λ20
=
e−tv h(s, t)dsdt
∫
∫
(1 − ρ0 )(λ0 (1 − g0 (v)) − v) 0
0
w0 (u)w0 (v) 2
e−tv h(s, t)d(s, t).
λ w0 (v) ∫
=−
(1 − ρ0 )2 0
[0,M ]2

Finally, compounding above covariance structures yields ζ(⋅, ⋅).
The covariance structure ζ(⋅, ⋅) depends on the unknown objects. However, one can use
the suggested estimators and plug them in place of the true ones. This might be helpful to
implement the problem and the provided consistency results ensure accuracy as long as the
sample size is large enough.

31

Funding: This work was supported by the Deutsche Forschungsgemeinschaft (German
Research Foundation) within the programme "Statistical Modeling of Complex Systems and
Processes—Advanced Nonparametric Approaches", grant GRK 1953.
Corresponding Author
Dr. Cornelia Wichelhaus
Technische Universität Darmstadt
Schlossgartenstraße 7
64289 Darmstadt
Germany
wichelhaus@mathematik.tu-darmstadt.de

Moritz von Rohrscheidt
Ruprecht-Karls Universität Heidelberg
Berliner Straße 41-49
69120 Heidelberg
Germany
rohrscheidt@uni-heidelberg.de

32

References
Armero, C. and M. Bayarri
1994. Prior Assessments for Prediction in Queues. The Statistician, Pp. 139–153.
Armero, C. and D. Conesa
1998. Inference and Prediction in Bulk Arrival Queues and Queues with Service in Stages.
Applied Stochastic Models and Data Analysis, 14(1):35–46.
Armero, C. and D. Conesa
2004. Statistical Performance of a Multiclass Bulk Production Queueing System. European
Journal of Operational Research, 158(3):649–661.
Armero, C. and D. Conesa
2006. Bayesian Hierarchical Models in Manufacturing Bulk Service Queues. Journal of
Statistical Planning and Inference, 136(2):335–354.
Asmussen, S.
2008. Applied Probability and Queues, volume 51. Springer Science & Business Media.
Ausín, M. C., R. E. Lillo, and M. P. Wiper
2004. Bayesian estimation for the M/G/1 queue using a phase-type approximation. Journal
of Statistical Planning and Inference, 118(1):83–101.
Ausín, M. C., R. E. Lillo, and M. P. Wiper
2007. Bayesian Control of the Number of Servers in a GI/M/c Queueing System. Journal
of Statistical Planning and Inference, 137(10):3043–3057.
Ausín, M. C., R. E. Lillo, and M. P. Wiper
2008. Bayesian Prediction of the Transient Behaviour and Busy Period in short-and longtailed GI/G/1 Queueing Systems. Computational Statistics & Data Analysis, 52(3):1615–
1635.
Conti, P. L.
1999. Large Sample Bayesian analysis for Geo/G/1 discrete-time Queueing Models. The
Annals of Statistics, Pp. 1785–1807.
De Finetti, B.
1974. Theory of Probability, volume 1. Wiley New York.
de Laplace, P.
1774. Mémoire sur les suites récurro-récurrentes et sur leurs usages dans la théorie des
hasards. Mém. Acad. Roy. Sci. Paris, 6:353–371.
Dey, J., R. Erickson, and R. Ramamoorthi
2003. Some Aspects of neutral-to-right Priors. International Statistical Review, 71(2):383–
401.
Diaconis, P. and D. Freedman
1986. On the Consistency of Bayes Estimates. The Annals of Statistics, Pp. 1–26.

33

Diaconis, P. and D. Ylvisaker
1985. Quantifying Prior Opinion. Bayesian Statistics, 2:133–156.
Doksum, K.
1974. Tailfree and Neutral Random Probabilities and their Posterior Distributions. The
Annals of Probability, Pp. 183–201.
Doob, J. L.
1949. Application of the Theory of Martingales. Le calcul des probabilites et ses applications,
Pp. 23–27.
Epifani, I., A. Lijoi, and I. Prünster
2003. Exponential Functionals and Means of neutral-to-the-right Priors.
90(4):791–808.

Biometrika,

Feigin, P. D. and R. L. Tweedie
1989.
Linear Functionals and Markov Chains associated with Dirichlet Processes.
105(03):579–585.
Ferguson, T. S.
1973. A Bayesian Analysis of some Nonparametric Problems. The Annals of Statistics,
Pp. 209–230.
Freedman, D. A.
1963. On the Asymptotic Behavior of Bayes’ Estimates in the Discrete Case. The Annals
of Mathematical Statistics, Pp. 1386–1403.
Ghosh, J. K. and R. Ramamoorthi
2003. Bayesian Nonparametrics. Springer Series in Statistics.
Gill, R. D. and S. Johansen
1990. A survey of Product-Integration with a View toward Application in Survival Analysis.
The Annals of Statistics, Pp. 1501–1555.
Grübel, R. and S. M. Pitts
1992. A Functional Approach to the Stationary Waiting Time and Idle Period Distributions
of the GI/G/1 Queue. The Annals of Probability, Pp. 1754–1778.
Haigh, J.
2004. Probability Models. Springer Science & Business Media.
Hewitt, E. and L. J. Savage
1955. Symmetric Measures on Cartesian Products. Transactions of the American Mathematical Society, 80(2):470–501.
Hjort, N. L.
1990. Nonparametric Bayes Estimators based on Beta Processes in Models for Life History
Data. The Annals of Statistics, Pp. 1259–1294.
Insua, D. R., M. Wiper, and F. Ruggeri
1998. Bayesian Analysis of M/Er/1 and M/H_k/1 Queues. Queueing Systems, 30(3-4):289–
308.
34

Johnson, R. A.
1967. An Asymptotic Expansion for Posterior Distributions. The Annals of Mathematical
Statistics, Pp. 1899–1906.
Kechris, A. S.
1995. Classical Descriptive Set Theory. Springer.
Kendall, D. G.
1953. Stochastic Processes occurring in the Theory of Queues and their Analysis by the
Method of the imbedded Markov Chain. The Annals of Mathematical Statistics, Pp. 338–
354.
Kim, Y. and J. Lee
2001. On Posterior Consistency of Survival Models. The Annals of Statistics, Pp. 666–686.
Kim, Y. and J. Lee
2004. A Bernstein-von Mises Theorem in the Nonparametric Right-Censoring Model. The
Annals of Statistics, Pp. 1492–1512.
Kingman, J. F.
1978. Uses of Exchangeability. The Annals of Probability, Pp. 183–197.
Kosorok, M.
2008. Introduction to Empirical Processes and Semiparametric Inference. Springer, New
York.
Nelson, R.
2013. Probability, Stochastic Processes, and Queueing Theory: the Mathematics of Computer
Performance Modeling. Springer Science & Business Media.
Phadia, E. G.
2015. Prior Processes and their Applications. Springer.
Regazzini, E., A. Lijoi, and I. Prünster
2003. Distributional Results for Means of Normalized Random Measures with Independent
Increments. The Annals of Statistics, Pp. 560–585.
Sato, K.
1999. Lévy Processes and Infinitely Divisible Distributions.
Schervish, M. J.
1995. Theory of Statistics. Springer Science & Business Media.
Walker, S. and P. Muliere
1997. Beta-Stacy Processes and a Generalization of the Pólya-Urn Scheme. The Annals of
Statistics, Pp. 1762–1780.
Wiper, M.
1998. Bayesian Analysis of Er/M/1 and Er/M/c Queues. Journal of Statistical Planning
and Inference, 69(1):65–79.

35

