Integer Sequences from Queueing Theory

Joseph Abate 900 Hammond Road Ridgewood, NJ 07450-2908 USA

Ward Whitt Department of Industrial Engineering and Operations Research Columbia University New York, NY 10027-6699 USA

                      ww2040@columbia.edu

Abstract Operators on probability distributions can be expressed as operators on the associated moment sequences, and so correspond to operators on integer sequences. Thus, there is an opportunity to apply each theory to the other. Moreover, probability models can be sources of integer sequences, both classical and new, as we show by considering the classical M/G/1 single-server queueing model. We identify moment sequences that are integer sequences. We establish connections between the M/M/1 busy period distribution and the Catalan and Schroeder numbers.

1 Introduction

Given the well established link between integer sequences and combinatorics, it is surprising that there are so few connections between integer sequences and queueing theory, because two prominent scholars, John Riordan (1902–1988) and Lajos Takács (1924– ), were active in both combinatorics and queueing theory [22, 23, 24, 28, 29, 30]. Of course, these authors saw connections; e.g., the index of Riordan’s queueing book [23] contains entries for Bell polynomials, binomial moments, Catalan numbers, cumulant generating function, Lagrange expansion, and rooted trees.

In the OEIS [26], a search on “Riordan” yields 1221 entries, while a search on “queueing” produces 2. We think Riordan would want less disparity. One exception is the entry submitted by A. Harel (\underline{A122525}) making connection to the Erlang delay formula associated with the classical M/M/s queue. A search on “queue” produces 15 entries, but these primarily focus on combinatorial problems [11, 27] or a queue as a tool in computer science rather than queueing theory [12]; queueing theory concerns stochastic process models describing congestion, e.g., the probability distribution of customer waiting times [23, 28].

The purpose of this paper is to point out connections between the theories of probability and integer sequences, and to exhibit some integer sequences that arise in queueing theory. As a branch of probability theory, queueing theory has exploited the classical analytical approach to probability theory using transforms [23, 28], whose series expansions involve the integer moments of the probability distributions and close relatives. The entire probability distribution is characterized by the sequence of moments in great generality [10, 17]. Thus we were motivated to introduce an operational calculus for manipulating probability distributions on the positive halfline by either manipulating the associated Laplace transforms or the associated moment sequences [7]; a quick overview is provided by [7, Tables 1-3]. A more recent paper in the same spirit is [18].

The operational calculus for probability distributions on the positive halfline in the framework of moment sequences is closely related to operators commonly used to analyze integer sequences. Since many operators on integer sequences can be applied to moment sequences arising in probability theory, both when they are integers and when they are not, there is an opportunity for experts on integer sequences to contribute to probability theory through moment sequences. (The connection is also discussed in [13].) The probability connection also provides concrete models where integer sequences arise.

Here is how the present paper is organized. We start in §2 by reviewing moment sequences in probability theory and relating operators on probability distributions to operators on sequences. Then in §3 we review the classical M/G/1 single-server queueing model and identify integer sequences arising in that context. We consider the M/M/1 busy period distribution, its stationary excess distribution and the equilibrium time to emptiness. In that context we identify random quantities whose probability density functions have the Catalan and Large Schroeder numbers as moments. Finally, we consider the moment sequence for the M/G/1 steady-state waiting time. We give the proofs of Theorems 5 and 11 in §4.

2 Moment Sequences in Probability

Let Z be a nonnegative random variable with cumulative distribution function (cdf) F and probability density function (pdf) f, i.e.,

F (t) ≡ P (Z ≤ t) = ∫_{0}^{t} f(u) du, t ≥ 0,

where ≡ means “defined as.” Let f̂(s) be the Laplace transform (LT) of f (and thus Z) and let ϕ(x) be the associated moment generating function (mgf) of f, i.e.,

f̂(s) ≡ ∫_{0}^{∞} e^{−st}f(t) dt ≡ E[e^{−sZ}] and
ϕ(x) ≡ f̂(−x) = ∫_{0}^{∞} e^{xt}f(t) dt ≡ E[e^{xZ}],  (1)

where s in (1) is understood to be a complex number with positive real part, while x in (1) is a positive real number. The LT is always well defined; we assume that there exists x^{∗} > 0 such that ϕ(x) < ∞ for x < x^{∗}.

We obtain sequences by considering series expansions for the mgf ϕ; in particular, we can write

ϕ(x) = ∑_{n=0}^{∞} \frac{m_{n}x^{n}}{n!} = ∑_{n=0}^{∞} μ_{n}x^{n},  (2)

where of course we must have μ_{n} = m_{n}/n!, n ≥ 0. From the probability perspective, the object of primary interest is the cdf F, but the moments m_{n} provide a useful partial characterization. The moment sequence can be calculated and employed to derive other quantities of interest in probability models [1, 7, 16]. There is a developing operational calculus for manipulating probability distributions via their moment sequences [7].

From the sequence perspective, we may instead regard the sequences {m_{n} : n ≥ 1} and {μ_{n} : n ≥ 1} as the objects of primary interest. The two perspectives meet with the mgf ϕ. From the sequence perspective, ϕ(x) arises as the generating function (gf) of the sequence {μ_{n} : n ≥ 1}, and we speak of {μ_{n}} as being its coefficients, while ϕ(x) arises as the exponential generating function (egf) of the sequence {m_{n} : n ≥ 1}. Of course, {m_{n} : n ≥ 1} is always an integer sequence whenever {μ_{n} : n ≥ 1} is, but not necessarily conversely. We will be giving examples where both are integer sequences.

A simple example from probability theory is the exponential distribution with mean m, specified by

F (t) ≡ 1 − e^{−t/m}, f(t) ≡ (1/m)e^{−t/m}, t ≥ 0, and ϕ(x) ≡ (1 − mx)^{−1},

which has associated sequences

m_{n} ≡ n!m^{n} and μ_{n} ≡ m^{n}, n ≥ 0.

It is immediate that {μ_{n}} and {m_{n}} are both elementary integer sequences whenever m is an integer. Since the mean (first moment) is probabilistically only a scale parameter, depending on how the measurement units are defined, it is natural to follow the convention that the first moment is 1; here that gives m_{1} = 1. Then we obtain the fundamental integer sequences m_{n} = n! and μ_{n} = 1, n ≥ 1.

As a consequence of the relations outlined above, we see that results about probability distributions can be translated into results about integer sequences, provided that the sequence {μ_{n}} or {m_{n}} is indeed an integer sequence. Similarly, results about integer sequences can be translated into results about probability distributions, provided that the egf or gf of the integer sequence is indeed the mgf of a bonafide pdf.

We first observe that there is a natural probabilistic setting in which, not only is {m_{n}} a moment sequence, but so is the associated sequence {μ_{n}}. That occurs in spectral representations, where one pdf serves as a mixing pdf for another. In particular, suppose that the pdf f can be represented as a continuous mixture of exponential pdf’s via

f(t) = ∫_{τ_{1}}^{τ_{2}} y^{−1}e^{−t/y}g(y) dy, t ≥ 0,  (3)

in which case we call g the mixing pdf for f [4]. Append a superscript f to {m_{n}} and {μ_{n}} to denote dependence upon f. We observe that the sequence {μ_{n}} is itself the moment sequence of the associated mixing pdf g. Hence we call {μ_{n}} the mixing moments of f. (We omit the elementary proofs in this section.)

Proposition 1. (mixing moments) For pdf’s f and g related via (3), m_{n}^{g} = μ_{n}^{f}, n ≥ 1.

A canonical operation in probability theory is convolution. If two independent nonnegative random variables Z_{1} and Z_{2} with pdf’s f_{Z_{1}} and f_{Z_{2}} are added, then the sum Z_{1} + Z_{2} has a pdf that is the convolution of the two component pdf’s, i.e.,

f_{Z_{1}+Z_{2}}(t) = ∫_{0}^{t} f_{Z_{1}}(y)f_{Z_{2}}(t − y) dy, t ≥ 0.

One reason transforms are so frequently used in probability is that they convert convolution into a simple product; i.e., the associated mgf’s are related by

ϕ_{Z_{1}+Z_{2}}(x) = ϕ_{Z_{1}}(x)ϕ_{Z_{2}}(x).

Because of the assumed stochastic independence, the moments are related by the binomial theorem,

m_{n}^{Z_{1}+Z_{2}} = ∑_{k=1}^{n} \binom{n}{k}m_{k}^{Z_{1}}m_{n−k}^{Z_{2}}.

Thus, if {m_{n}^{Z_{1}}} and {m_{n}^{Z_{2}}} are integer sequences, then so is {m_{n}^{Z_{1}+Z_{2}}}. For example, by above, that occurs when Z_{1} and Z_{2} have exponential or deterministic distributions with integer means.

With this background, we can interpret [7, Tables 1-3], which show various operators mapping one probability distribution into another. There are four columns: The first column contains the name and notation for the operator; the second column shows how the operator acts on LT’s; the third column shows how it acts on pdf’s; and the fourth column shows how it acts on moment sequences. From the perspective of integer sequences, the fourth column shows how it acts on the coefficients of the egf; we could then add a fifth column showing how it acts on the corresponding coefficients of the gf. Those familiar with integer sequences might want to translate the LT into the mgf by replacing s with −x, and then interpret that mgf as the egf of the given k^{th} moment.

In the rest of this section we highlight a few striking connections between operators on probability distributions, as in [7], and operators on integer sequences. First, a standard operator on integer sequences is the simple shift to the left, e.g., converting 1, 1, 2, 6, 22, 90, … to 1, 2, 6, 22, 90, . . .. We now show that, probabilistically, the simple shift applied to the coefficients μ_{n} corresponds to constructing the stationary-excess cdf of a cdf on the positive halfline having mean 1.

Given a nonnegative real-valued random variable Z with cdf F having finite moments m_{k}, k ≥ 1, let Z_{e} be a random variable with the associated stationary-excess cdf F_{e} (a.k.a. the equilibrium excess or stationary residual-life cdf), defined by

F_{e}(t) ≡ P (Z_{e} ≤ t) ≡ \frac{1}{m_{1}} ∫_{0}^{t}(1 − F (u)) du, t ≥ 0.  (4)

The stationary-excess cdf frequently arises in renewal theory; see of [25, Examples 7.16,7.17, 7.23, 7.24] and [31]; it appears in [7, Table 2]. (A search on “renewal theory” in the OEIS gives two unrelated entries.) For us, the important fact is that the random variable Z_{e} has moments

m_{e,k} ≡ E[Z_{e}^{k}] = \frac{m_{k+1}}{(k + 1)m_{1}} , k ≥ 1.  (5)

Hence, the transformation from a cdf of a nonnegative random variable to its associated stationary-excess cdf produces a simple shift on the gf coefficients.

Proposition 2. (the stationary-excess operator) Let F be the cdf of a nonnegative random variable with mean 1, mgf ϕ in (1) and associated sequence of mixing moments {μ_{n}} in (2) (coefficients of ϕ(x) when it is regarded as a gf). The associated stationary-excess cdf F_{e} in (4) has mgf ϕ_{e}(x) = (ϕ(x) − 1)/x. The mixing moments of ϕ_{e}(x) (coefficients of ϕ_{e} when it is regarded as a gf) are μ_{e,k} = μ_{k+1}, k ≥ 1.

A similar relationship holds for the stationary-lifetime operator, mapping a pdf f into the associated pdf

f_{s}(t) ≡ \frac{tf(t)}{m_{1}}, t ≥ 0.  (6)

The stationary-lifetime pdf also frequently arises in renewal theory, e.g., [25, §7.7], and also appears in [7, Table 2]. For us, the important fact is that the moments are related by

m_{s,k} ≡ \frac{m_{k+1}}{m_{1}} , k ≥ 1,

just like (5) without the k + 1 in the denominator. Hence, the transformation from a pdf of a nonnegative random variable to its associated stationary-lifetime pdf produces a simple shift on the moments m_{n} (the egf coefficients).

Proposition 3. (the stationary-lifetime operator) Let f be the pdf of a nonnegative random variable with mean 1, mgf ϕ in (1) and associated sequence of moments {m_{n}} in (2). The associated stationary-lifetime pdf f_{s} in (6) has mgf ϕ_{s}(x) = ϕ′(x). The coefficients of ϕ_{s} when it is regarded as an egf are m_{s,k} = m_{k+1}, k ≥ 1.

We now turn to another basic probability operator, which is conveniently related to a continued fraction representation of the mgf [8, 9, 13, 15]. If f̂ is the LT of a pdf f, then the associated exponential mixture pdf has LT

f̂_{EM}(s) ≡ (1 + sf̂(s))^{−1};  (7)

it appears in [7, Table 3]. The special case of exponential mixtures of inverse Gaussian (EMIG) distributions is discussed in [7, §8] and in [9].

The probabilistic exponential mixing operator has a simple manifestation in the continued fraction representation of the mgf, regarding that mgf as a gf. Starting with the LT f̂(s) of a pdf f, if we represent the associated mgf as a formal power series by

ϕ(x) ≡ f̂(−x) = 1 + μ_{1}x + μ_{2}x^{2} + μ_{3}x^{3} + μ_{4}x^{4} + μ_{5}x^{5} + …,  (8)

then the corresponding continued fraction (CF) is

ϕ(x) = \frac{1}{1−} \frac{h_{1}x}{1−} \frac{h_{2}x}{1−} \frac{h_{3}x}{1−} \frac{h_{4}x}{1−} …  (9)

where h_{1} ≡ μ_{1}, h_{2} ≡ (μ_{2} − μ_{1}^{2})/μ_{1}, etc. When h_{n} > 0 for all n we have an S-fraction and the underlying pdf f is completely monotone (CM) [8].

Proposition 4. (the exponential-mixture operator) Let f be the pdf of a nonnegative random variable with mean 1, mgf ϕ in (1) with associated sequence of CF coefficients {h_{n}} in (9). The associated exponential-mixture pdf f_{EM} with LT in (7) has CF coefficients

h_{EM,1} = 1, h_{EM,n} = h_{n−1}, n ≥ 2;

i.e., it produces a shift to the right. The inverse exponential mixture operator [7, (7.3), p. 94] gives the corresponding shift to the left.

Proof. From the transform expression in (7), the conclusion is immediate: Given the CF representation for f̂(s), it is immediate that the corresponding CF for (1 + sf̂(s))^{−1} shifts the coefficients one to the right; i.e., we write the CF for (1 + sf̂(s))^{−1} as 1/(1 + sf̂(s)), inserting the CF for f̂(s).

3 Queueing Examples

3.1 The M/G/1 Model

The M/G/1 queue is a basic model in queueing theory, usually discussed in queueing textbooks; e.g., [25, §8.5]. There is a single server with unlimited waiting room. Customers arrive according to a Poisson process (the first M, for Markov) with rate λ, 0 < λ < ∞. If the system is empty, then the customer goes immediately into service; otherwise the customer waits in queue. The successive service times come from a sequence of independent and identically distributed (i.i.d.) random variables with cdf G having mean 1/μ, 0 < μ < ∞. We will mostly consider the easiest case, in which the service-time cdf G is exponential; then the model is denoted by M/M/1.

A waiting customer enters service immediately upon service completion. Let Q(t) be the number of customers in the system at time t for t ≥ 0. In the M/M/1 model, the stochastic process Q ≡ {Q(t) : t ≥ 0} is a birth-and-death stochastic process, with constant birth rate λ and constant death rate μ. Let ρ ≡ λ/μ be the traffic intensity. If ρ < 1, then P (Q(t) = j|Q(0) = i) → (1 − ρ)ρ^{j} as t → ∞ for each i and j; i.e., Q(t) converges in distribution to a geometric distribution on the nonnegative integers, having mean ρ/(1 − ρ). We assume that ρ < 1, under which the system is said to be stable. (If ρ ≥ 1, then P (Q(t) ≤ j|Q(0) = i) → 0 as t → ∞ for each i and j.) For the M/G/1 model, the steadystate distribution of Q(t) is characterized by the Pollaczek-Khintchine transform [23, 28].

3.2 The M/M/1 Busy Period

A busy period is the time from the arrival of a customer finding an empty system until the system is empty again [23, §4.8]. The first passage time of Q from any state j > 0 to j − 1 is distributed as a busy period. Without loss of generality, we can measure time in units of mean service times, so that we let μ ≡ 1. Then the model has the single parameter ρ (= λ). Let X_{ρ} denote a random variable with the busy period distribution, as a function of the traffic intensity ρ. Let B_{ρ} be the cdf of X_{ρ}, i.e., B_{ρ}(t) ≡ P (X_{ρ} ≤ t), t ≥ 0, and let b_{ρ} be the associated pdf. It turns out that

b_{ρ}(t) ≡ \frac{1}{t\sqrt{ρ}}e^{−(1+ρ)t}I_{1}(2t\sqrt{ρ}), t ≥ 0,

where I_{1}(t) is the Bessel function of the first kind [23, (39), p. 63]. The LT of b_{ρ} (and thus X_{ρ}) is

b̂_{ρ}(s) ≡ ∫_{0}^{∞} e^{−st}b_{ρ}(t) dt ≡ E[e^{−sX_{ρ}}] = \frac{1 + ρ + s − \sqrt{(1 + ρ + s)^{2} − 4ρ}}{2ρ},

[23, (38), p. 63]. As usual, the moments can be obtained by differentiating the Laplace transform. The first two moments are E[X_{ρ}] = (1 − ρ)^{−1} and E[X_{ρ}^{2}] = 2(1 − ρ)^{−3}.

3.3 The Busy-Period Moment Sequence After Scaling

Our goal here is to obtain interesting integer sequences from the sequence of successive integer moments of a busy period X_{ρ}. To obtain integer sequences, we first perform a change of variables, introducing σ ≡ ρ/(1 − ρ) (the mean steady state number in system) or, equivalently, ρ ≡ σ/(1 + σ). Notice that the first two moments become 1 + σ and 2(1 + σ)^{3}. It turns out that the entire moment sequence becomes an integer sequence whenever σ is an integer. Hence, from this first step, we obtain an entire sequence of integer sequences.

To seek simple integer sequences, we further scale these moment sequences all to have mean (first moment) 1. We remark that this final spatial scaling also plays a role in understanding the performance of the M/M/1 queue as the traffic intensity ρ increases toward its critical value 1. With appropriate scaling of both time and space, the stochastic process Q approaches reflected Brownian motion with negative drift, while the busy period distribution has interesting behavior, in which both small values and large values play a role [3, 6, 32].

In terms of random variables, let

Y_{σ} ≡ \frac{X_{σ/(1+σ)}}{1 + σ} = (1 − ρ)X_{ρ}, where ρ ≡ \frac{σ}{1 + σ} and σ = \frac{ρ}{1 − ρ}.  (10)

Let b(x; σ) be the mgf of Y_{σ}. From above,

b(x; σ) ≡ E[e^{xY_{σ}}] = \frac{1 + 2σ − x − Ψ(x)}{2σ}, where Ψ(x) ≡ \sqrt{1 − 2(1 + 2σ)x + x^{2}}. (11)

Let the associated moments of Y_{σ} be m_{n}(σ) ≡ E[Y_{σ}^{n}], n ≥ 1, where m_{1}(σ) = 1 for all σ > 0. These moments, divided by n!, are the coefficients of the series expansion of b(x; σ)

b(x; σ) = ∑_{n=0}^{∞} \frac{m_{n}(σ)x^{n}}{n!} = ∑_{n=0}^{∞} b_{n}(σ)x^{n}, where b_{n}(σ) ≡ \frac{m_{n}(σ)}{n!}.  (12)

The coefficients b_{n}(σ) (and thus also the moments m_{n}(σ)) are polynomials in σ (and thus integers when σ is an integer); the first few are: b_{0}(σ) = 1, b_{1}(σ) = 1, b_{2}(σ) = 1 + σ and b_{3}(σ) = 1 + 3σ + 2σ^{2}.

More specifically, we now relate the coefficients b_{n}(σ) in (12) to the Catalan numbers, denoted by C_{k}, (\underline{A000108}) starting with 1, 1, 2, 5, 14, 42, 132, 429, 1430; in particular,

C_{n} ≡ \frac{1}{n + 1}\binom{2n}{n} and C_{n+1} = ∑_{i=0}^{n} C_{i}C_{n−i}, n ≥ 1.  (13)

The Catalan numbers can be characterized in terms of their generating function

c(y) ≡ ∑_{k=0}^{∞} C_{k}y^{k} ≡ \frac{1 − \sqrt{1 − 4y}}{2y} .  (14)

Theorem 5. (mixing moments of the busy period pdf) For n ≥ 1,

b_{n+1}(σ) = ∑_{k=0}^{n} \binom{n + k}{n − k}C_{k}σ^{k},  (15)

where b_{n}(σ) is defined in (11) and (12), and C_{k} is the k^{th} Catalan number in (13). In addition, there is a convenient recurrence relation, with b_{0}(σ) ≡ b_{1}(σ) ≡ 1 and

(n + 1)b_{n+1}(σ) = (2n − 1)(1 + σ)b_{n}(σ) − (n − 2)b_{n−1}(σ).  (16)

We provide a proof in §4. For σ = 1, 2, 3, the coefficient sequences are:

\{b_{n}(1)\} = 1, 1, 2, 6, 22, 90, 394, (\underline{A155069})
\{b_{n}(2)\} = 1, 1, 3, 15, 93, 645, 4791, (\underline{A103210})
\{b_{n}(3)\} = 1, 1, 4, 28, 244, 2380, 24868, (\underline{A103211})

Sequence \{b_{n}(1)\} (\underline{A155069}) is a relatively recent addition to the OEIS [26], having the title “Expansion of (3 − x − \sqrt{1 − 6x + x^{2}})/2 in powers of x.” We provide a model context. However, sequence \{b_{n}(1)\} shifted one to the left, i.e., 1, 2, 6, 22, 90, … is (\underline{A006318}), corresponding to the famous “Large Schroeder numbers.” The “little Schroeder numbers” in (\underline{A001003}) are obtained by dividing \underline{A006318} by 2, i.e., the sequence 1, 1, 3, 11, 45, . . ..

The moment sequences themselves, m_{n}(σ) ≡ n!b_{n}(σ) are of course also integer sequences, but they are evidently not in the OEIS [26]. For example, the sequences with terms n!b_{n}(1) and n!b_{n}(1/2), n ≥ 0, yielding 1, 1, 4, 36, 528, 10800 and 1, 1, 3, 18, 171, 2250, respectively, are not found.

The busy period is the first passage time of Q from 1 to 0. Since the first passage time pdf from each state k to 0 is the k-fold convolution of the busy period pdf [3, Theorem 3.1]), the moments and mixing moments of these pdf’s also generate integer sequences for integer σ.

Since we find regularity by multiplying {b_{n}(σ)} by 1/2 = σ/(1 + σ) when σ = 1, we are motivated to consider the associated sequences {σb_{n+1}(σ)/(1 + σ) : n ≥ 1} for integers σ > 1.

Corollary 6. For n ≥ 1,

\frac{σb_{n+1}(σ)}{1 + σ} = ∑_{k=0}^{n}(−1)^{n+k}\binom{n + k}{n − k}C_{k}(1 + σ)^{k}.

We have already observed that Corollary 6 gives the little Schroeder numbers for σ = 1; it is also discussed in [24, Problem 15, p. 168].

3.4 The Busy-Period Stationary-Excess Distribution

Let B_{e} denote the busy period stationary-excess cdf associated with the busy-period cdf B_{ρ}, defined as in (4) after applying the scaling in (10). Let b_{e} be the associated stationary-excess pdf. We can apply Proposition 2 to characterize the mixing moments b_{e,n}(σ) of b_{e}.

Corollary 7. (mixing moments of the busy-period stationary-excess pdf) The M/M/1 busyperiod stationary excess pdf b_{e} has mgf

b_{e}(x, σ) ≡ ∫_{0}^{∞} e^{xt}b_{e}(t) dt ≡ ∑_{n=0}^{∞} b_{e,n}(σ)x^{n} = \frac{b(x; σ) − 1}{x}
= \frac{2}{1 − x + \sqrt{1 − 2(1 + 2σ)x + x^{2}}}
= \frac{c(σx/(1 − x)^{2})}{1 − x} = \frac{2}{1 − x + Ψ(x)},  (17)

where c(y) is the generating function of the Catalan numbers in (14) and Ψ is defined in (11). For n ≥ 1,

b_{e,n}(σ) = b_{n+1}(σ), n ≥ 1,  (18)

for which an expression is given in Theorem 5. The mean of b_{e} is σ.

3.5 Large Schroeder and Catalan Numbers as Moments of PDF’s

In this section we identify the pdf’s whose moments are (i) the Large Schroeder numbers and (ii) the Catalan numbers. By Corollary 7 and our observation after Theorem 5, the sequence \{b_{e,n}(1)\} coincides with the Large Schroeder numbers (\underline{A006318}). From [3, Theorem 5.1 and Corollary 5.2.1] and [4, Theorem 4.1], we can obtain the spectral representation for the pdf b_{e}. (We need to account for the different scaling of time used there.) That identifies the desired pdf and, by Proposition 1, the mixing moments.

Corollary 8. (Large Schroeder numbers as moments) The large Schroeder numbers arise as the moments of the mixing pdf associated with the stationary-excess busy-period pdf b_{e} when σ = 1. The mixing pdf for b_{e} as a function of σ is

f(y; σ) = \frac{\sqrt{(τ − y)(y − τ^{−1})}}{2σπy}, \frac{1}{τ} < y < τ,

where τ ≡ 1 + 2σ + 2\sqrt{σ(1 + σ)}.

We introduced the Catalan numbers before Theorem 5. The Catalan numbers are related to reflected Brownian motion with drift −1 and diffusion coefficient 1, denoted by {R(t) : t ≥ 0}. It is the limit of the M/M/1 queue length process as ρ ↑ 1 with appropriate scaling of time and space [2, 32]. Since E[R(t)|R(0) = 0] is nondecreasing, H_{1}(t) ≡ E[R(t)|R(0) = 0]/E[R(∞)], t ≥ 0, is a cdf. By [2, Corollary 1.5.2], [7, §7] and [8, (8.7)], the pdf h_{1} of H_{1} has LT ĥ_{1}(s), which can be characterized as the unique fixed point of the exponential mixture operator, i.e.,

ĥ_{1}(s) = \frac{1}{1 + sĥ_{1}(s)}.  (19)

Together with the spectral representation given in [4, Theorem 4.1], that implies the following result. Without the probability model context, the result already appears in the commentary on (\underline{A000108}) [21].

Theorem 9. (Catalan numbers as moments) The generating function c(x) of the Catalan numbers in (13) and (14) coincides with the mgf of h_{1}, ĥ_{1}(−x), associated with the LT in (19). As a consequence, the Catalan numbers arise as the moments of the mixing density of h_{1},

f(y) = \frac{\sqrt{4 − y}}{2π\sqrt{y}} , 0 < y < 4.  (20)

If we take the Laplace transform of f in (20) and replace s by −x, then we obtain the mgf, which is the egf c̃(x) of the Catalan numbers.

Corollary 10. The egf of the Catalan numbers is

c̃(x) ≡ ∑_{n=0}^{∞} \frac{C_{n}x^{n}}{n!} = e^{2x}(I_{0}(2x) − I_{1}(2x)),

where I_{0} and I_{1} are Bessel functions.

3.6 The M/M/1 Equilibrium Time to Emptiness

From a sequence perspective, the shifting by 1 and multiplying by σ/(1 + σ) following Theorem 5 and Corollary 6 lead us to consider the mgf

p(x; σ) ≡ \frac{1}{1 + σ} + \frac{σ}{1 + σ} (\frac{b(x; σ) − 1}{x}) .  (21)

The function p(x; σ) in (21) turns out to be the moment generating function of the equilibrium time to emptiness in the M/M/1 queue, i.e., the first passage time to state 0 by the stochastic process Q, assuming that Q starts at time 0 according to its (geometric) steadystate distribution. Hence, with probability (1 − ρ) = 1/(1 + σ) the process starts at 0, so the first passage time is 0. As a consequence, there is an atom at 0 with mass 1/(1 + σ). Conditional on starting at a positive value, the equilibrium time to emptiness coincides with the stationary-excess pdf b_{e} of the busy-period pdf b [5, Theorem 3]. Additional characterizations appear in [3, Theorem 3.3]. We make a connection to the exponential mixture operator in (7). Paralleling (12), we let p_{n}(σ) be the coefficient of p(x; σ) as a gf, i.e., writing p(x; σ) ≡ ∑_{n=0}^{∞} p_{n}(σ)x^{n}. We provide a proof in §4.

Theorem 11. (equilibrium time to emptiness) The mgf p(x; σ) in (21) can be expressed as

p(x; σ) = \frac{c((1 + σ)x/(1 + x)^{2})}{1 + x} = \frac{2}{1 + x + Ψ(x)},  (22)

where c(y) is the generating function of the Catalan numbers in (14) and Ψ is defined in (11). Hence,

p_{n}(σ) = \frac{σb_{n+1}(σ)}{1 + σ} = \frac{σ}{1 + σ} ∑_{k=0}^{n} \binom{n + k}{n − k}C_{k}σ^{k}, n ≥ 0,  (23)

and the following recursion can be applied with p_{0}(σ) ≡ 1 and p_{1}(σ) ≡ σ

(n + 2)p_{n+1}(σ) = (2n + 1)(1 + 2σ)p_{n}(σ) − (n − 1)p_{n−1}(σ), n ≥ 1.  (24)

In addition, b(x; σ) can be expressed as an exponential mixture of p(x; σ), i.e.,

b(x; σ) = \frac{1}{1 − xp(x; σ)}.  (25)

From (23) and (24), the sequence {p_{n}(σ) : n ≥ 1} is an integer sequence for each positive integer σ. For σ = 1, 2, 3, the coefficient sequences are:

\{p_{n}(1)\} = 1, 1, 3, 11, 45, 197, (\underline{A001003})
\{p_{n}(2)\} = 1, 2, 10, 62, 430, 3194, (\underline{A107841})
\{p_{n}(3)\} = 1, 3, 21, 183, 1785, 18651, (\underline{A131763}).

We conclude this section by remarking that in the theory of integer sequences there also is a convenient invert operator, which can be expressed for LT’s via

ĝ(s) ≡ I(f̂(s)) ≡ Excess(exp mixture(f̂(s)) = f̂_{EM,e}(s)
= Excess (\frac{1}{1 + sf̂(s)}) = \frac{1}{s} (1 − \frac{1}{1 + sf̂(s)})

= \frac{f̂(s)}{1 + sf̂(s)}.  (26)

The following result links the three mgf’s b(x; σ), b_{e}(x; σ) and (p(x; σ); see Corollary 5.1 of [3]. We also consider the excess of the excess, denoted by b_{e.e} [31].

Corollary 12. The mgf’s b_{e}(x; σ), p(x; σ) and b(x; σ) are related by

b_{e}(x; σ) = I(p(x; σ)) = p(x; σ)b(x; σ)

for the invert operator in (26) and

b_{e,e}(x; σ) = p_{e}(x; σ).

Directly from the final expression in (17), we get a heavy-traffic limit for the scaled mgf with our scaling in (10), with the limit being the gf of the Catalan numbers. Let X_{e}(σ) be a random variable with mgf b_{e}(x; σ); i.e., b_{e}(x; σ) ≡ E[e^{xX_{e}(σ)}]. Let ⇒ denote convergence in distribution [32].

Corollary 13. As σ → ∞,

b_{e}(x/σ; σ) → \frac{2}{1 − \sqrt{1 − 4x}} = \frac{1 − \sqrt{1 − 4x}}{2x} ≡ c(x)

for c in (14), so that

\frac{X_{e}(σ)}{σ} ⇒ X as σ → ∞,

where c(x) = E[e^{xX}].

A related heavy-traffic limit for the time-scaled busy-period stationary-excess pdf was obtained in [6, Theorem 1]. (They are consistent.)

3.7 Continued Fractions and Hankel Transforms

We now look at the three probability mgf’s b(x; σ), b_{e}(x,; σ) ≡ (b(x; σ) − 1)/x and p(x; σ) in (11), (17) and (21) from the perspective of continued fractions and Hankel transforms. From the previous subsections, we know that the associated three random variables are closely related. That is shown again through this new perspective. First, the three mgf’s can be represented as S (Stieltjes) fractions. (In [8] we found that S fractions frequently arise in

\begin{matrix}mgf & h_{1} & h_{2} & h_{3} & h_{4} & h_{5} & h_{6} \\ b(x; σ) & 1 & σ & 1 + σ & σ & 1 + σ & σ \\ b_{e}(x; σ) & 1 + σ & σ & 1 + σ & σ & 1 + σ & σ \\ p(x; σ) & σ & 1 + σ & σ & 1 + σ & σ & 1 + σ \end{matrix}

Table 1: The coefficients of the continued fractions representing the three M/M/1 mgf’s. The pattern repeats in each row.

CF’s associated with birth-and-death stochastic processes.) The CF coefficients are given in Table 1.

For the generating functions of our M/M/1 queueing examples, the determination of the CF coefficients is simple because we can apply the algebraic identity

\frac{α}{1+} \frac{γ}{1+} \frac{α}{1+} \frac{γ}{1+} … = \frac{1}{2} (\sqrt{1 + 2(γ + α) + (γ − α)^{2}} − 1 − (γ − α))

for constants α and γ; i.e., we can solve the equation

ζ = \frac{α}{1+} \frac{γ}{1 + ζ}

for ζ.

The Hankel transform of an integer sequence provides a useful partial characterization; it is a many-to-one function mapping an integer sequence into another integer sequence. (For an example of non-uniqueness, see Corollary 14 below.) For example, the Hankel transform of the Catalan numbers is the sequence 1, 1, 1, … [13, 20]. Following [15, 17], starting from a sequence {ω_{n} : n ≥ 0} ≡ ω_{0}, ω_{1}, ω_{2}, ω_{2}, … with ω_{0} ≡ 1, let the Hankel matrix M^{(n)} be the (n + 1) × (n + 1) symmetric matrix with elements M_{i,j}^{(n)} ≡ ω_{i+j−2}, 0 ≤ i ≤ n, 0 ≤ j ≤ n. (The first row contains the first n + 1 elements and M_{n+1,n+1} ≡ ω_{2n}. Let H_{2n} ≡ det(M^{(n)}), the even Hankel determinant. Let the Hankel transform of the sequence {ω_{n} : n ≥ 0} above be the sequence {H_{2n} : n ≥ 0}; it starts with H_{0} = 1.

One way to compute the Hankel transform of a sequence {ω_{n} : n ≥ 0} is to first determine the corresponding continued fraction from the formal power series; i.e., starting from (8), we determine (9) above. The coefficients h_{n} appearing in (9) can be determined by applying the normalized Viskovatov algorithm, as given on [15, p. 112]. Then we apply the iteration

H_{2n} = (∏_{i=1}^{2n} h_{i}) H_{2n−2}, n ≥ 1;

see Theorem 1.4.10 on of [17, p. 23]. Another way to arrive at this result is via the even contraction of the CF in (9), as given in [13, (12.3), p. 270]; note that β_{n} in [13] is equal to h_{2n−1}h_{2n}, n ≥ 1, in our notation. Then (12.2) of [13, (12.2)] and our iteration yield the same result.

Corollary 14. The Hankel transforms of the sequences {b_{n}(σ)} in Theorem 5 , {b_{e,n}(σ)} in Corollary 7 and {p_{n}(σ} in Theorem 11 are

H_{2n}(b) = σ^{ν(n)}(1 + σ)^{ν(n)−n},
H_{2n}(p) = H_{2n}(b_{e}) = (σ + σ^{2})^{ν(n)},

where ν(n) ≡ n(n + 1)/2.

3.8 The Stationary Waiting Time in the M/G/1 Queue

The mgf w(x; σ) ≡ E[e^{xW}] of the steady-state waiting time W in the M/G/1 queue can be expressed in terms of the mgf g_{e}(x; σ) of the stationary-excess of the general service-time distribution using the Pollazcek-Khintchine transform (again assuming that the mean service time is 1) by applying the random sum operator from [7, Table 1], yielding

w(x; σ) = \frac{1 − ρ}{1 − ρg_{e}(x; σ)} = \frac{1}{1 + σ − g_{e}(x; σ)}  (27)

[23, §4.3] and [5]. Riordan [23, (18a), p. 49] and Takács [28] develop a nice recursion for the moments, which as a function of σ becomes

E[W^{n−1}(σ)] = \frac{σ}{n} ∑_{k=2}^{n} \binom{n}{2}g_{k}E[W^{n−k}(σ)] and E[W^{0}(σ)] ≡ 1,  (28)

where g_{k} is the k^{th} moment of the service-time cdf. Clearly, (28) can be the source of many integer sequences.

To illustrate, we now consider the Catalan numbers as service-time moments, which is legitimate by Theorem 9. We exploit the following result, obtained by combining [2, Corollaries 1.3.2 and 1.5.1] with Theorem 9.

Lemma 15. Let h_{1} be the density of the RBM first moment function, which has the Catalan numbers as its mixing moments. Then

ĥ_{1,e}s = ĥ_{1}(s)^{2}, so that c_{e}(x) = c(x)^{2},  (29)

where c(x) is the generating function of the Catalan numbers.

Theorem 16. If the service time pdf is the RBM first moment pdf h_{1} characterized by its LT in (19), which has the Catalan numbers as its moments, then the waiting time mgf in (27) becomes

w(x; σ) = \frac{1}{1 + σ − σc(x)^{2}}  (30)

and the moments of the stationary waiting time pdf are

E[W^{n−1}(σ)] = σ ∑_{k=2}^{n} \frac{(n − 1)!}{(n − k)!}C_{k}E[W^{n−k}(σ)].  (31)

Hence, if σ is an integer, then {E[W^{n}(σ)]} is an integer sequence.

For the case σ = 1, we find {E[W^{n}(1)]} = 1, 2, 18, 252, 4776, …, which is not in the OEIS. Riordan [23, (19), p. 50] shows how the moments E[W^{n}(σ)] can be expressed in terms of the multivariate Bell polynomials.

Let w_{n}(σ) ≡ E[W^{n}(σ)]/n! be the associated mixing moments. From (31), we obtain a recursion for w_{n}(σ).

Corollary 17. A recursion for the mixing moments w_{n}(σ) defined above is

w_{n−1}(σ) = σ ∑_{k=2}^{n} C_{k}w_{n−k}(σ), n ≥ 1.  (32)

For σ = 1, we get {w_{n}(1)} = 1, 2, 9, 42, 199, …, which also is not yet in OEIS.

3.9 General Birth-and-Death Processes

The most relevant previous work connecting queueing theory to integer sequences and providing a suitable framework for generalization evidently is the previous work connecting birth-and-death processes to continued fractions; see [8, 18] and references therein; continued fractions are known to be intimately connected to integer sequences. We have already mentioned that the stochastic process Q ≡ {Q(t) : t ≥ 0}, representing the number of customers in an M/M/1 queueing model at time t, is a birth-and-death stochastic process. Indeed, it is a special birth-and-death process with constant birth rate λ and constant death rate μ. More generally, these rates are functions of the state [25, §6.3]; λ_{k} is the rate up and μ_{k} is the rate down when Q(t) = k. The general birth-and-death process represents a more general queueing model.

For a general birth-and-death process X ≡ {X(t) : t ≥ 0}, an important quantity is the probability

P_{0,0}(t) ≡ P (X(t) = 0|X(0) = 0), t ≥ 0;  (33)

see [3, §11], [18, §3.1] and [19]. From [18, (4.1), p. 771], the LT of P_{0,0}(t), P̂_{0,0}(s), has a representation as an S fraction

P̂_{0,0}(s) = \frac{1}{s+} \frac{λ_{0}}{1+} \frac{μ_{1}}{s+} \frac{λ_{1}}{1+} \frac{μ_{2}}{s+} …
= \frac{z}{1+} \frac{λ_{0}z}{1+} \frac{μ_{1}z}{1+} \frac{λ_{1}z}{1+} \frac{μ_{2}z}{1+} …  (34)

where z ≡ 1/s. (The last line follows from [8, (1.5)] using the sequence {c_{n} : n ≥ 0} with c_{2n} ≡ 1 and c_{2n+1} ≡ z, n ≥ 1.) We can find the corresponding power series

P̂_{0,0}(s) = z(1 − p_{1}z + p_{2}z^{2} − p_{3}z^{3} + …);

see [8, (3.9), (3.10), p. 397]. The sequence {1, p_{1}, p_{2}, p_{3}, …} may be an integer sequence.

For the special case of the M/M/1 queue, with our scaling we have λ_{i} = σ and μ_{i+1} = 1+σ for all i, i ≥ 0, so that

P̂_{0,0}(s) = z (\frac{1}{1+} \frac{σz}{1+} \frac{(1 + σ)z}{1+} \frac{σz}{1+} \frac{(1 + σ)z}{1+} …) .

Then from Table 1 we have the relation

P̂_{0,0}(s) = zp̂(z),  (35)

where p̂(s) is the Laplace transform of the equilibrium time to emptiness.

4 Proofs

In this concluding section we provide proofs for Theorems 5 and 11.

4.1 Proof of Theorem 5.

We give two proofs. The first is direct; the second applies [14].

Direct proof of (15). From [29, pp. 232–233], after a change of scale and notation, we have

b_{n+1} = (\frac{1 + σ}{n + 1}) ∑_{k=1}^{n} \binom{n + k}{k}\binom{n − 1}{k − 1}σ^{k−1}.

(This equation is also given by [24, (21), p. 151] after we make the identification that H_{n}(z) = p_{n}(σ) = σb_{n+1}(σ)/(1 + σ).) On the right side, move (1 + σ) inside the sum and identify the coefficients of σ^{k}, obtaining

(\frac{1}{n + 1}) (\binom{n + k + 1}{k + 1}\binom{n − 1}{k} + \binom{n + k}{k}\binom{n − 1}{k − 1}) = \binom{n + k}{n − k}\binom{2k}{k} (\frac{1}{k + 1}) .

Direct proof of (16). We follow the technique of [23, p. 107] to establish a three-term recursion. By the successive differentiation with respect to x in (11), we find that

Ψ(x)^{2}b′′(x; σ) = (1 + σ)(1 + 2σ − x)b′(x; σ) + (1 + σ)^{2}b(x; σ).

In this equation, make the following substitutions:

b(x; σ) = ∑_{n=0}^{∞} b_{n}(σ)x^{n},

b′(x; σ) = ∑_{n=0}^{∞}(n + 1)b_{n+1}(σ)x^{n},

b′′(x; σ) = ∑_{n=0}^{∞}(n + 1)(n + 2)b_{n+2}(σ)x^{n}.

Then collect and equate the coefficients of x^{n−1} and the result follows.

Application of [14]. The idea in this second proof is to directly make connection to the Catalan numbers via their generating function c in (14). Starting from (11), we apply (17),

(21) and the proof of (22) in Theorem 11 below to conclude that

b_{e}(x; σ) = \frac{1}{1 − x}c (\frac{σx}{(1 − x)^{2}} ) and p(x; σ) = \frac{1}{1 + x}c (\frac{(1 + σ)x}{(1 + x)^{2}}) .  (36)

We then apply [14, Proposition 15, p.12] (with change of variables y ≡ σx and y ≡ (1 + σ)x, respectively) to immediately deduce the conclusions

b_{e,n}(σ) = ∑_{k=0}^{n} \binom{n + k}{n − k}C_{k}σ^{k} and p_{n}(σ) = ∑_{k=0}^{n} \binom{n + k}{n − k}(−1)^{n+k}C_{k}(1 + σ)^{k}  (37)

given in Corollary 7 and Theorem 11. We then can apply (18) and (17) to get the conclusion in Theorem 5 for b from b_{e}. By this line of reasoning, we establish Theorems 5 and 11 and Corollaries 6 and 7 all at once, and have a new perspective on their relationship.

4.2 Proof of Theorem 11.

Proof of (22). From the representation of p(x; σ) in (21), we have

p(x; σ) = \frac{1}{1 + σ} + \frac{σ}{1 + σ} (\frac{2}{1 − x + Ψ(x)})
= \frac{1}{1 + σ} + \frac{σ}{1 + σ} (\frac{1 − x − Ψ(x)}{2σx})
= \frac{1 + x − Ψ(x)}{2(1 + σ)x} = \frac{2}{1 + x + Ψ(x)}.

On the other hand, by (14),

\frac{c((1 + σ)x/(1 + x)^{2})}{1 + x} = \frac{2}{1 + x + (1 + x)\sqrt{1 − 4(1 + σ)x/(1 + x)^{2}}}
= \frac{2}{1 + x + Ψ(x)}.

Proof of (23). Combine (15) and (21).

Proof of (24). Combine (16) and (21).

Proof of (25). From (11),

b(x; σ) = \frac{1 + 2σ − x − Ψ(x)}{2σ} = \frac{2(1 + σ)}{1 + 2σ − x + Ψ(x)}
= \frac{1}{1 − x ([1 + x − Ψ(x)]/[2(1 + σ)x])} = \frac{1}{1 − xp(x; σ)}.

5 Acknowledgment

The second author was supported by NSF Grant CMMI 0948190.

References

[1] J. Abate, G. L. Choudhury, D. M. Lucantoni and W. Whitt, Asymptotic analysis of tail probabilities based on the computation of moments, Ann. Appl. Prob. 5 (1995), 983–1007.

[2] J. Abate and W. Whitt, Transient behavior of regulated Brownian motion, Adv. Appl. Prob. 19 (1987), 560–598.

[3] J. Abate and W. Whitt, Transient behavior of the M/M/1 queue via Laplace transforms, Adv. Appl. Prob. 20 (1988), 145–178.

[4] J. Abate and W. Whitt, Simple spectral representations for the M/M/1 queue, Queueing Systems 3 (1988), 321–346.

[5] J. Abate and W. Whitt, Transient behavior of the M/G/1 workload process, Operations Research 42 (1994), 750–764.

[6] J. Abate and W. Whitt, Limits and approximations for the busy-period distribution in single-server queues, Prob. Engr. Inform. Sci. 9 (1995), 581–602.

[7] J. Abate and W. Whitt, An operational calculus for probability distributions via Laplace transforms, Adv. Appl. Prob. 28 (1996), 75–113.

[8] J. Abate and W. Whitt, Computing Laplace transforms for numerical inversion via continued fractions, INFORMS J. Computing 11 (1999), 394–405.

[9] J. Abate and W. Whitt, Explicit M/G/1 waiting-time distributions for a class of long-tail service-time distributions.. Operations Research Letters 25 (1999), 25–31.

[10] N. I. Akhiezer, The Classical Moment Problem and Some Related Questions in Analysis, Oliver and Boyd, London, 1965.

[11] M. H. Albert, R. E. L. Aldred, M. D. Atkinson, H. P. van Ditmarsch, C. C. Handley and D. A. Holton, Restricted permutations and queue junping, Discrete Math. 287 (2004), 129–133.

[12] M. D. Atkinson, R. Beals, Priority queues and permutations, SIAM J. Comput. 23 (1994), 1225–1230.

[13] P. Barry, A Study of Integer Sequences, Riordan Arrays, Pascal-like Arrays and Hankel Transforms, PhD dissertation Department of Mathematics, University College Cork, Ireland, 2009.

[14] P. Barry, Continued fractions and transformations of integer sequences, J. Integer Sequences 12 (2009), Article 09.7.6.

[15] A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland and W. B. Jones, Handbook of Continued fractions for Special Functions, Springer, 2008.

[16] G. L. Choudhury and D. M. Lucantoni, Numerical computation of the moments of a probability distribution from its transform, Oper. Res. 44 (1996), 368–381.

[17] H. Dette and W. J. Studden, The Theory of Canonical Moments with Applications to Statistics, Probability and Analysis, Wiley, 1997.

[18] P. Flajolet and F. Guillemin, The formal theory of birth-and-death processes, lattice path combinatorics and continued fractions, Adv. Appl. Prob. 32 (2000), 750–778.

[19] J. F. C. Kingman, Regenerative Phonomena, Wiley, 1972.

[20] J. W. Layman, The Hankel transform and some of its properties, J. Integer Sequences 4 (2001), Article 01.1.5.

[21] K. A. Penson, An integral representation of C_{n}, commentary on (\underline{A000108}), in [26].

[22] J. Riordan, An Introduction to Combinatorial Analysis, Wiley, 1958.

[23] J. Riordan, Stochastic Service Systems, Wiley, 1962.

[24] J. Riordan, Combinatorial Identities, Wiley, 1968.

[25] S. M. Ross, Introduction to Probability Models, ninth edition, Academic Press, 2007.

[26] N. J. A. Sloane, The On-Line Encylopedia of Integer Sequences (OEIS), 2010.

[27] R. P. Stanley, Queue problems revisited, Suomen Tehtavaniekat (Proceedings of the Finnish Chess Problem Society) 59 (2005), 193–203.

[28] L. Takács, Introduction to the Theory of Queues, Oxford University Press, 1962.

[29] L. Takács, Combinatorial Methods in the Theory of Stochastic Processes, Wiley, 1967.

[30] L. Takács, On the number of distinct forests, SIAM J. Discrete Math. 3 (1990), 574–581.

[31] W. Whitt, The renewal-process stationary-excess operator. J. Appl. Prob. 22 (1985), 156–167.

[32] W. Whitt, Stochastic-Process Limits, Springer, 2002.

2000 Mathematics Subject Classification: Primary 11B25; Secondary 60K25. Keywords: integer sequence, moment sequence, probability, probability density function, probability moment sequence, the moment problem, queue, queueing, Laplace transforms, moment generating functions, generating functions, continued fractions, Hankel transform, renewal theory, birth-and-death process, stationary excess distribution, stationary lifetime distribution, exponential mixture, single-server queue, M/M/1 queue, M/G/1 queue, busy period, equilibrium time to emptiness, waiting time.

(Concerned with sequences \underline{A000108}, \underline{A001003}, \underline{A006318}, \underline{A093252}, \underline{A093526}, \underline{A103210}, \underline{A103211}, \underline{A107841}, \underline{A122525}, \underline{A131763}, \underline{A137216}, and \underline{A155069}.)

Received February 19 2010; revised version received May 2 2010. Published in Journal of Integer Sequences, May 3 2010.

Return to Journal of Integer Sequences home page.