- Email: [email protected]
Precise formulae for Bravo coefficients
Accepted Manuscript Precise formulae for Bravo coefficients Saralees Nadarajah, Tibor K. Pogány
PII: DOI: Reference:
S0167-6377(18)30004-X https://doi.org/10.1016/j.orl.2018.01.001 OPERES 6324
To appear in:
Operations Research Letters
Received date : 27 February 2016 Revised date : 2 January 2018 Accepted date : 2 January 2018 Please cite this article as: S. Nadarajah, T.K. Pogány, Precise formulae for Bravo coefficients, Operations Research Letters (2018), https://doi.org/10.1016/j.orl.2018.01.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to view linked References
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Precise Formulae for Bravo Coefficients (Short Title: Precise Formulae for Bravo Coefficients) by Saralees Nadarajah 1 University of Manchester, Manchester, UK Tibor K. Pog´ any University of Rijeka, Rijeka, CROATIA and ´ Obuda University, Budapest, HUNGARY
Abstract: Daley et al. assessed the BRAVO (balancing reduces asymptotic variance of outputs) effect of a finite-buffer Markovian many-server system in terms of an asymptotic variance ratio. The expressions given for the ratio in just the two cases considered both appear incorrect. Here, we derive the correct expressions. Keywords: Error function; Integral; Normal distribution
1 Corresponding author: Saralees Nadarajah, School of Mathematics, University of Manchester, Manchester M13 9PL, UK, email: [email protected]
1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1
Introduction
In queuing systems, there has been recent phenomena called the BRAVO (balancing reduces asymptotic variance of outputs) effect. It arises whenever the departure variance decreases significantly as the result of the service capacity being set to match the arrival rate. The BRAVO effect was first observed by [8] for M/M/1/K systems. Since then several papers have appeared on the BRAVO effect, see [1], [2], [7] and [6]. According to https: // www.smp.uq.edu.au / people / YoniNazarathy / research / BRAVO Cambridge v1.pdf , some open problems are: How can BRAVO be harnessed in practice? Why does BRAVO occur? The numerical experiments in [1] suggest that the BRAVO effect may occur in critically loaded systems. So, a characterization of the BRAVO effect (if found) could be of practical use. Suppose a queuing system has finite waiting capacity K, arrival rate λ, s servers and service rate µ/s per server. The BRAVO effect has been studied in terms of the following variance-mean ratio for various queuing systems: D=
Var [Ndep (0, t]] , s,K→∞ t→∞ E [Ndep (0, t]] lim
lim
(1)
where Ndep (0, t] denotes the number of serviced customers in (0, t]. For simple Poisson processes, D = 1. For M/M/1/K systems, D = 23 when λ = µ, as shown by [8]. For M/M/1/∞ systems, D = 1 when λ 6= µ and D = 2 (1 − 2/π) when λ = µ, as shown by [1]. [1] went on to derive expressions for D for GI/GI/1 and GI/GI/s systems. The most recent paper on the BRAVO effect is that by [3]. They investigated the BRAVO effect √ in the quality and efficiency driven regime for the M/M/s system such that s (1 − λ/µ) → β and √ K/ s → η as λ → ∞, s → ∞, λ/µ → 1 and K → ∞. The parameter β describes the “scaled shortfall from one in the system capacity” and the parameter η describes the “relative buffer size”. The main result in [3] gave expressions for D := Dβ,η in (1) for the two cases: λ/µ = 1 and √ λ/µ = 1 − β/ s, β 6= 0. Unfortunately, the expressions given in both cases appear incorrect. The aim of this note is to derive the correct expressions for D := Dβ,η in (1) in both cases. Sec√ tion 2 considers the case λ/µ = 1. Section 3 considers the case λ/µ = 1−β/ s, β 6= 0. Throughout, Φ(·) and φ(·) denote, respectively, the cumulative distribution function and the probability density function of a standard normal random variable.
2
The case λ/µ = 1
In this case, β = 0. [3] showed that p 3 2 Z ∞ η 2 π2 + η3 2b2 u √ − Φ(−u) 1 − b exp Φ(−u) du, (2) D0,η = 1 − p 3 π 2 2π 0 + η 2 √ pπ where b = 2π/ 2 + η . The integral in the right hand side of (2) can be calculated as follows. Lemma 2.1. It holds that √ 2 Z ∞ u 1 − b log 2 √ . Φ(−u) 1 − b exp I+ = Φ(−u) du = 2 2π 0 2
(3)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
√ Proof: Note that Φ(−u) = (1/2)erfc u/ 2 , where erfc(·) denotes the complementary error function defined by Z ∞ 2 √ erfc (z) = exp −t2 dt. π z So,
I+ = = =
2 Z u u u 1 ∞ b erfc √ erfc √ 1 − exp du 2 0 2 2 2 2 2 Z Z u u u b ∞ 1 ∞ 2 du − erfc √ du erfc √ exp 2 0 4 0 2 2 2 b log 2 1 √ − √ , 2π 2 2π
where the last equality follows by equations (2.8.2.1) and (2.8.20.11) in [9], volume 2. Hence, the result. 2 Combining (2) and (3), we obtain D0,η
p √ 3 η 2 π2 + η3 2b2 1 − b log 2 √ . =1− p 3 − √ π 2π 2π 2 +η
The corresponding expression given in [3] is √ 2 − π2 η + 2π 1 − log 2 − 2 D0,η = − 3 pπ 3 +η 2
π 12
.
(4)
Furthermore, Lemma 3.1 in [3], showing 2 Z ∞ √ √ u Φ(−u) 1 − exp Φ(−u) du = 2π 1 − log 2 , 2 0 and used to prove (4), is incorrect. In fact, √ 2 Z ∞ u 1 − log 2 √ Φ(−u) 1 − exp , Φ(−u) du = 2 2π 0
(5)
by our Lemma 2.1 above. We have verified (5) numerically using the R software [3]. The following code f=function (u) {exp(pnorm(-u,log.p=T)+log(1-exp(u*u/2+pnorm(-u,log.p=T))))} integrate(f,lower=0,upper=Inf)$value returned 0.2606794, the value of the right hand side of (5) up to seven decimal places. The proof of Lemma 3.1 in [3] appears correct, only the statement of Lemma 3.1 is incorrect. But their proof appears rather long. Our proof based on the integration book bypasses their detailed derivations, our proof is only a few lines long. Furthermore, (4) appears correct and not affected by the incorrect statement in Lemma 3.1 in [3]. 3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
√ The case λ/µ = 1 − β/ s, β 6= 0
3
In this case, [3] established the expression (see their equation (2.1) and also their graphics) Dβ,η = 1 −
2β 2 exp(−βη)h2 (η, β) f (η, β) + g(η, β), φ(β)
(6)
where 1 , 1 − exp(−βη) + βΦ(β)/φ(β) f (η, β) = K−β − β exp(−βη)h(η, β)J−β , h(η, β) =
g(η, β) = 2 exp(−βη)h(η, β) [1 − exp(−βη)h(η, β)] [1 − βη − exp(−βη)]
+2 exp(−βη)h2 (η, β) [1 − exp(−βη)h(η, β)] [1 − 2βη exp(−βη) − exp(−2βη)] , Z ∞ Φ(−u)du, Kβ = Jβ =
Z
β ∞
β
[Φ(−u)]2 du. φ(u)
[3], page 246, equations (4.9), (4.10) went on to show that Z ∞ Φ (− | β |) β β 2 x2 1 Jβ = − √ log 2 + exp − log 1 + 2 dx 2π 1 2 x 2π
(7)
and J−|β|
1 √ 4 log 2 2|β| − exp = J|β| − J0 + du + √ Φ 2 2π 0 Z √2 2 2 β x 1 2|β| exp − log 1 + 2 dx − π 2 x 0 Z
|β|
u2 2
(8)
for β > 0 and β < 0, respectively. Both (7) are (8) are incorrect: to see this note that log 2 lim RHS(7) = √ β↓0 2 2π and lim RHS(8) = 0, β↑0
√ but J0 = log 2 by the proof of Lemma 2.1. Furthermore, (6), (7) are (8) are all non-explicit in that they are expressed in terms of integrals. In the rest of this note, we give formulas for Kβ and Jβ expressed explicitly in terms of the hypergeometric and Kamp´e de F´eriet functions. Using Lemma 2.1, we can write Jβ = log
Z √ 2−
β 0
4
[Φ(−u)]2 du φ(u)
(9)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
and Jβ = log
√
2+
Z
−β 0
[1 − Φ(−u)]2 du φ(u)
(10)
for β > 0 and β < 0, respectively. Using the fact 2z erfc(z) = 1 − √ 1 F1 π
1 3 ; ; −z 2 , 2 2
(11)
where 1 F1 (α; β; z) denotes the confluent hypergeometric function defined by 1 F1 (α; β; z) =
∞ X (α)k z k k=0
(β)k k!
,
where (α)k = α(α + 1) · · · (α + k − 1) denotes the ascending factorial, we can rewrite (9) and (10) as 2 Z β √ 1 1 u 1 3 u2 du (12) − √ 1 F1 ; ;− Jβ = log 2 − 2 2 2 2π 0 φ(u) 2 and Z √ Jβ = log 2 +
0
−β
2 1 3 u2 u 1 1 + √ 1 F1 ; ;− du, φ(u) 2 2 2 2 2π
(13)
respectively, for β > 0 and β < 0, respectively. Lemma 3.1 derives an explicit expression for the integral Ka . Lemmas 3.2 to 3.4 derive explicit expressions for the integrals in (12) and (13). Lemma 3.1. For a > 0, Ka =
Z
∞ a
1 a 1 1 a2 Φ(−u)du = √ 1 F1 − ; ; − − . 2 2 2 2 2π
Proof: By the proof of Lemma 2.1, we can write Z ∞ Φ(−u)du Ka = a Z a Z ∞ Φ(−u)du Φ(−u)du − = 0 0 Z a 1 √ = Φ(−u)du − 2π 0 Z a u 1 1 erfc √ du. = √ − 2π 2 0 2
5
(14)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
By (11), the last term in (14) can be expressed as Z a Z 1 a u u 1 1 3 u2 − √ 1 F1 ; ;− erfc √ du = du 2 0 2 2 2 2 2 2π 0 k 2k+2 ∞ a 1 X (1/2)k a 1 = −√ − 2 2 2k + 2 2π k=0 (3/2)k k! 2 k+1 ∞ (−1/2)k+1 a 1 X a = −√ − 2 (1/2) (k + 1)! 2 2π k=0 k+1 1 a 1 1 a2 −√ −1 . = 1 F1 − ; ; − 2 2 2 2 2π The result follows. 2 Lemma 3.2. For a > 0,
Z
a 0
1 du = πerfi φ(u)
a √ 2
,
where erfi(·) denotes the imaginary error function defined by Z z 2 erfi (z) = √ exp t2 dt. π 0 Proof: Follows from the definition of the imaginary error function. 2 Lemma 3.3. For a > 0, √ Z a a2 1 3 u2 3 2πa2 u ; ;− du = , 1 F1 2 F2 1, 1; , 2; 2 2 2 2 2 2 0 φ(u) where 2 F2 (α, β; γ, δ; z) denotes the hypergeometric function defined by 2 F2 (α, β; γ, δ; z) =
∞ X (α)k (β)k z k . (γ)k (δ)k k! k=0
Proof: By the transformation formula 1 F1 (β − α; β; z) = exp(z) 1 F1 (α; β; −z), we can write 2 Z a 1 3 u2 u ; ;− du u exp 1 F1 2 2 2 2 0 Z a 3 u2 u1 F1 1; ; = du 2 2 0 Z a X ∞ (1)k u2k u = du (3/2)k k! 2k 0 k=0 Z a ∞ X (1)k 1 u2k+1 du = (3/2)k k! 2k 0 k=0
= =
∞ X
(1)k 1 a2k+2 (3/2)k k! 2k 2(k + 1)
k=0 ∞ a2 X
2
k=0
(1)k (1)k a2k . (3/2)k (2)k k! 2k 6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
The result now follows from the definition of the 2 F2 (α, β; γ, δ; z) hypergeometric function. 2 Lemma 3.4. For a > 0, √ 2 Z a 2 1 3 u2 u 2πa3 1,1,1 du = ; ;− F1,1,1 1 F1 2 2 2 3 0 φ(u) where p,r,t Fq,s,u
=
a1 , . . . , ap c1 , . . . , cr b1 , . . . , bq d1 , . . . , ds
p Y
∞ ∞ X X j=1 k=0 ℓ=0
(aj )k+ℓ
q Y
j=1
(bj )k+ℓ
r Y
j=1 s Y
j=1
(cj )k
(dj )ℓ
3/2 1 5/2 3/2
! a2 1/2 a2 , , − 3/2 2 2
e1 , . . . , at f1 , . . . , fu
t Y
j=1 u Y
(ej )ℓ
(fj )ℓ
! x, y
xk y ℓ k! ℓ!
j=1
denotes the Kamp´e de F´eriet function. Proof: By the transformation formula 1 F1 (β − α; β; z) = exp(z) 1 F1 (α; β; −z), we can write 2 2 Z a u 1 3 u2 2 du ; ;− u exp 1 F1 2 2 2 2 0 Z a 3 u2 1 3 u2 u2 1 F1 1; ; = F ; ; − du 1 1 2 2 2 2 2 0 Z a X ∞ X ∞ (−1)ℓ (1)k (1/2)ℓ u2(k+ℓ) u2 = du (3/2)k (3/2)ℓ k!ℓ! 2k+ℓ 0 k=0 ℓ=0 Z a ∞ X ∞ X (−1)ℓ (1)k (1/2)ℓ u2+2(k+ℓ) du = (3/2)k (3/2)ℓ k!ℓ! 2k+ℓ 0 k=0 ℓ=0 k ℓ ∞ ∞ a2 /2 −a2 /2 (1)k (1/2)ℓ a3 X X = 2 (k + ℓ + 3/2) (3/2)k (3/2)ℓ k! ℓ! k=0 ℓ=0 k ℓ ∞ ∞ −a2 /2 a3 X X (3/2)k+ℓ (1)k (1/2)ℓ a2 /2 = . 3 (5/2)k+ℓ (3/2)k (3/2)ℓ k! ℓ! k=0 ℓ=0
The results now follows from the definition of the Kamp´e de F´eriet function. 2
Theorem 3.1 derives explicit expressions for (12) and (13). Its proof is immediate from Lemmas 3.2 to 3.4. Theorem 3.1. For β > 0, √ π β β2 β3 β2 3 1,1,1 F1,1,1 + 2 F2 1, 1; , 2; − √ Jβ = log 2 − erfi √ 4 2 2 2 2 3 2π For β < 0, β β2 π β2 β3 3 1,1,1 F1,1,1 Jβ = log 2 + erfi − √ + + √ 2 F2 1, 1; , 2; 4 2 2 2 2 3 2π √ Also J0 = log 2 (by the proof of Lemma 2.1). √
7
3/2 1 5/2 3/2
! β2 1/2 β 2 . , − 3/2 2 2
3/2 1 5/2 3/2
! β2 1/2 β 2 . , − 3/2 2 2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Combining Lemma 3.1 and Theorem 3.1, we obtain an explicit expression for (6). In-built routines for computing the imaginary error, hypergeometric and Kamp´e de F´eriet functions are available in packages like Maple, Matlab and Mathematica. See also [4], Section 1.3.2 and [5].
Acknowledgments The authors would like to thank the Editor for careful reading and comments which improved the paper.
References [1] Al Hanbali, A., Mandjes, M., Nazarathy, Y. and Whitt, W. (2011). The asymptotic variance of departures in critically loaded queues. Advances in Applied Probability, 43, 243-263. [2] Daley, D. J. (2011). Revisiting queueing output processes: A point process viewpoint. Queueing Systems, 68, 395-405. [3] Daley, D. J., Van Leeuwaarden, J. S. H. and Nazarathy, Y. (2015). BRAVO for many-server QED systems with finite buffers. Advances in Applied Probability, 47, 231-250. [4] Exton, H. (1978). Handbook of Hypergeometric Integrals: Theory, Applications, Tables, Computer Programs. Chichester, England. [5] Exton, H. and Krupnikov, E. D. (1998). A register of computer-oriented reduction identities for the Kamp´e de F´eriet function. Draft manuscript. [6] Hautphenne, S., Kerner, Y., Nazarathy, Y. and Taylor, P. (2015). The intercept term of the asymptotic variance curve for some queueing output processes. European Journal of Operational Research, 242, 455-464. [7] Nazarathy, Y. (2011). The variance of departure processes: Puzzling behavior and open problems. Queueing Systems, 68, 385-394. [8] Nazarathy, Y. and Weiss, G. (2008). The asymptotic variance rate of the output process of finite capacity birth-death queues. Queueing Systems, 59, 135-156. [9] Prudnikov, A. P., Brychkov, Y. A. and Marichev, O. I. (1986). Integrals and Series, volumes 1, 2 and 3. Gordon and Breach Science Publishers, Amsterdam. [10] R Development Core Team (2017). A Language and Environment for Statistical Computing: R Foundation for Statistical Computing. Vienna, Austria.
8
PII: DOI: Reference:
S0167-6377(18)30004-X https://doi.org/10.1016/j.orl.2018.01.001 OPERES 6324
To appear in:
Operations Research Letters
Received date : 27 February 2016 Revised date : 2 January 2018 Accepted date : 2 January 2018 Please cite this article as: S. Nadarajah, T.K. Pogány, Precise formulae for Bravo coefficients, Operations Research Letters (2018), https://doi.org/10.1016/j.orl.2018.01.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to view linked References
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Precise Formulae for Bravo Coefficients (Short Title: Precise Formulae for Bravo Coefficients) by Saralees Nadarajah 1 University of Manchester, Manchester, UK Tibor K. Pog´ any University of Rijeka, Rijeka, CROATIA and ´ Obuda University, Budapest, HUNGARY
Abstract: Daley et al. assessed the BRAVO (balancing reduces asymptotic variance of outputs) effect of a finite-buffer Markovian many-server system in terms of an asymptotic variance ratio. The expressions given for the ratio in just the two cases considered both appear incorrect. Here, we derive the correct expressions. Keywords: Error function; Integral; Normal distribution
1 Corresponding author: Saralees Nadarajah, School of Mathematics, University of Manchester, Manchester M13 9PL, UK, email: [email protected]
1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1
Introduction
In queuing systems, there has been recent phenomena called the BRAVO (balancing reduces asymptotic variance of outputs) effect. It arises whenever the departure variance decreases significantly as the result of the service capacity being set to match the arrival rate. The BRAVO effect was first observed by [8] for M/M/1/K systems. Since then several papers have appeared on the BRAVO effect, see [1], [2], [7] and [6]. According to https: // www.smp.uq.edu.au / people / YoniNazarathy / research / BRAVO Cambridge v1.pdf , some open problems are: How can BRAVO be harnessed in practice? Why does BRAVO occur? The numerical experiments in [1] suggest that the BRAVO effect may occur in critically loaded systems. So, a characterization of the BRAVO effect (if found) could be of practical use. Suppose a queuing system has finite waiting capacity K, arrival rate λ, s servers and service rate µ/s per server. The BRAVO effect has been studied in terms of the following variance-mean ratio for various queuing systems: D=
Var [Ndep (0, t]] , s,K→∞ t→∞ E [Ndep (0, t]] lim
lim
(1)
where Ndep (0, t] denotes the number of serviced customers in (0, t]. For simple Poisson processes, D = 1. For M/M/1/K systems, D = 23 when λ = µ, as shown by [8]. For M/M/1/∞ systems, D = 1 when λ 6= µ and D = 2 (1 − 2/π) when λ = µ, as shown by [1]. [1] went on to derive expressions for D for GI/GI/1 and GI/GI/s systems. The most recent paper on the BRAVO effect is that by [3]. They investigated the BRAVO effect √ in the quality and efficiency driven regime for the M/M/s system such that s (1 − λ/µ) → β and √ K/ s → η as λ → ∞, s → ∞, λ/µ → 1 and K → ∞. The parameter β describes the “scaled shortfall from one in the system capacity” and the parameter η describes the “relative buffer size”. The main result in [3] gave expressions for D := Dβ,η in (1) for the two cases: λ/µ = 1 and √ λ/µ = 1 − β/ s, β 6= 0. Unfortunately, the expressions given in both cases appear incorrect. The aim of this note is to derive the correct expressions for D := Dβ,η in (1) in both cases. Sec√ tion 2 considers the case λ/µ = 1. Section 3 considers the case λ/µ = 1−β/ s, β 6= 0. Throughout, Φ(·) and φ(·) denote, respectively, the cumulative distribution function and the probability density function of a standard normal random variable.
2
The case λ/µ = 1
In this case, β = 0. [3] showed that p 3 2 Z ∞ η 2 π2 + η3 2b2 u √ − Φ(−u) 1 − b exp Φ(−u) du, (2) D0,η = 1 − p 3 π 2 2π 0 + η 2 √ pπ where b = 2π/ 2 + η . The integral in the right hand side of (2) can be calculated as follows. Lemma 2.1. It holds that √ 2 Z ∞ u 1 − b log 2 √ . Φ(−u) 1 − b exp I+ = Φ(−u) du = 2 2π 0 2
(3)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
√ Proof: Note that Φ(−u) = (1/2)erfc u/ 2 , where erfc(·) denotes the complementary error function defined by Z ∞ 2 √ erfc (z) = exp −t2 dt. π z So,
I+ = = =
2 Z u u u 1 ∞ b erfc √ erfc √ 1 − exp du 2 0 2 2 2 2 2 Z Z u u u b ∞ 1 ∞ 2 du − erfc √ du erfc √ exp 2 0 4 0 2 2 2 b log 2 1 √ − √ , 2π 2 2π
where the last equality follows by equations (2.8.2.1) and (2.8.20.11) in [9], volume 2. Hence, the result. 2 Combining (2) and (3), we obtain D0,η
p √ 3 η 2 π2 + η3 2b2 1 − b log 2 √ . =1− p 3 − √ π 2π 2π 2 +η
The corresponding expression given in [3] is √ 2 − π2 η + 2π 1 − log 2 − 2 D0,η = − 3 pπ 3 +η 2
π 12
.
(4)
Furthermore, Lemma 3.1 in [3], showing 2 Z ∞ √ √ u Φ(−u) 1 − exp Φ(−u) du = 2π 1 − log 2 , 2 0 and used to prove (4), is incorrect. In fact, √ 2 Z ∞ u 1 − log 2 √ Φ(−u) 1 − exp , Φ(−u) du = 2 2π 0
(5)
by our Lemma 2.1 above. We have verified (5) numerically using the R software [3]. The following code f=function (u) {exp(pnorm(-u,log.p=T)+log(1-exp(u*u/2+pnorm(-u,log.p=T))))} integrate(f,lower=0,upper=Inf)$value returned 0.2606794, the value of the right hand side of (5) up to seven decimal places. The proof of Lemma 3.1 in [3] appears correct, only the statement of Lemma 3.1 is incorrect. But their proof appears rather long. Our proof based on the integration book bypasses their detailed derivations, our proof is only a few lines long. Furthermore, (4) appears correct and not affected by the incorrect statement in Lemma 3.1 in [3]. 3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
√ The case λ/µ = 1 − β/ s, β 6= 0
3
In this case, [3] established the expression (see their equation (2.1) and also their graphics) Dβ,η = 1 −
2β 2 exp(−βη)h2 (η, β) f (η, β) + g(η, β), φ(β)
(6)
where 1 , 1 − exp(−βη) + βΦ(β)/φ(β) f (η, β) = K−β − β exp(−βη)h(η, β)J−β , h(η, β) =
g(η, β) = 2 exp(−βη)h(η, β) [1 − exp(−βη)h(η, β)] [1 − βη − exp(−βη)]
+2 exp(−βη)h2 (η, β) [1 − exp(−βη)h(η, β)] [1 − 2βη exp(−βη) − exp(−2βη)] , Z ∞ Φ(−u)du, Kβ = Jβ =
Z
β ∞
β
[Φ(−u)]2 du. φ(u)
[3], page 246, equations (4.9), (4.10) went on to show that Z ∞ Φ (− | β |) β β 2 x2 1 Jβ = − √ log 2 + exp − log 1 + 2 dx 2π 1 2 x 2π
(7)
and J−|β|
1 √ 4 log 2 2|β| − exp = J|β| − J0 + du + √ Φ 2 2π 0 Z √2 2 2 β x 1 2|β| exp − log 1 + 2 dx − π 2 x 0 Z
|β|
u2 2
(8)
for β > 0 and β < 0, respectively. Both (7) are (8) are incorrect: to see this note that log 2 lim RHS(7) = √ β↓0 2 2π and lim RHS(8) = 0, β↑0
√ but J0 = log 2 by the proof of Lemma 2.1. Furthermore, (6), (7) are (8) are all non-explicit in that they are expressed in terms of integrals. In the rest of this note, we give formulas for Kβ and Jβ expressed explicitly in terms of the hypergeometric and Kamp´e de F´eriet functions. Using Lemma 2.1, we can write Jβ = log
Z √ 2−
β 0
4
[Φ(−u)]2 du φ(u)
(9)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
and Jβ = log
√
2+
Z
−β 0
[1 − Φ(−u)]2 du φ(u)
(10)
for β > 0 and β < 0, respectively. Using the fact 2z erfc(z) = 1 − √ 1 F1 π
1 3 ; ; −z 2 , 2 2
(11)
where 1 F1 (α; β; z) denotes the confluent hypergeometric function defined by 1 F1 (α; β; z) =
∞ X (α)k z k k=0
(β)k k!
,
where (α)k = α(α + 1) · · · (α + k − 1) denotes the ascending factorial, we can rewrite (9) and (10) as 2 Z β √ 1 1 u 1 3 u2 du (12) − √ 1 F1 ; ;− Jβ = log 2 − 2 2 2 2π 0 φ(u) 2 and Z √ Jβ = log 2 +
0
−β
2 1 3 u2 u 1 1 + √ 1 F1 ; ;− du, φ(u) 2 2 2 2 2π
(13)
respectively, for β > 0 and β < 0, respectively. Lemma 3.1 derives an explicit expression for the integral Ka . Lemmas 3.2 to 3.4 derive explicit expressions for the integrals in (12) and (13). Lemma 3.1. For a > 0, Ka =
Z
∞ a
1 a 1 1 a2 Φ(−u)du = √ 1 F1 − ; ; − − . 2 2 2 2 2π
Proof: By the proof of Lemma 2.1, we can write Z ∞ Φ(−u)du Ka = a Z a Z ∞ Φ(−u)du Φ(−u)du − = 0 0 Z a 1 √ = Φ(−u)du − 2π 0 Z a u 1 1 erfc √ du. = √ − 2π 2 0 2
5
(14)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
By (11), the last term in (14) can be expressed as Z a Z 1 a u u 1 1 3 u2 − √ 1 F1 ; ;− erfc √ du = du 2 0 2 2 2 2 2 2π 0 k 2k+2 ∞ a 1 X (1/2)k a 1 = −√ − 2 2 2k + 2 2π k=0 (3/2)k k! 2 k+1 ∞ (−1/2)k+1 a 1 X a = −√ − 2 (1/2) (k + 1)! 2 2π k=0 k+1 1 a 1 1 a2 −√ −1 . = 1 F1 − ; ; − 2 2 2 2 2π The result follows. 2 Lemma 3.2. For a > 0,
Z
a 0
1 du = πerfi φ(u)
a √ 2
,
where erfi(·) denotes the imaginary error function defined by Z z 2 erfi (z) = √ exp t2 dt. π 0 Proof: Follows from the definition of the imaginary error function. 2 Lemma 3.3. For a > 0, √ Z a a2 1 3 u2 3 2πa2 u ; ;− du = , 1 F1 2 F2 1, 1; , 2; 2 2 2 2 2 2 0 φ(u) where 2 F2 (α, β; γ, δ; z) denotes the hypergeometric function defined by 2 F2 (α, β; γ, δ; z) =
∞ X (α)k (β)k z k . (γ)k (δ)k k! k=0
Proof: By the transformation formula 1 F1 (β − α; β; z) = exp(z) 1 F1 (α; β; −z), we can write 2 Z a 1 3 u2 u ; ;− du u exp 1 F1 2 2 2 2 0 Z a 3 u2 u1 F1 1; ; = du 2 2 0 Z a X ∞ (1)k u2k u = du (3/2)k k! 2k 0 k=0 Z a ∞ X (1)k 1 u2k+1 du = (3/2)k k! 2k 0 k=0
= =
∞ X
(1)k 1 a2k+2 (3/2)k k! 2k 2(k + 1)
k=0 ∞ a2 X
2
k=0
(1)k (1)k a2k . (3/2)k (2)k k! 2k 6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
The result now follows from the definition of the 2 F2 (α, β; γ, δ; z) hypergeometric function. 2 Lemma 3.4. For a > 0, √ 2 Z a 2 1 3 u2 u 2πa3 1,1,1 du = ; ;− F1,1,1 1 F1 2 2 2 3 0 φ(u) where p,r,t Fq,s,u
=
a1 , . . . , ap c1 , . . . , cr b1 , . . . , bq d1 , . . . , ds
p Y
∞ ∞ X X j=1 k=0 ℓ=0
(aj )k+ℓ
q Y
j=1
(bj )k+ℓ
r Y
j=1 s Y
j=1
(cj )k
(dj )ℓ
3/2 1 5/2 3/2
! a2 1/2 a2 , , − 3/2 2 2
e1 , . . . , at f1 , . . . , fu
t Y
j=1 u Y
(ej )ℓ
(fj )ℓ
! x, y
xk y ℓ k! ℓ!
j=1
denotes the Kamp´e de F´eriet function. Proof: By the transformation formula 1 F1 (β − α; β; z) = exp(z) 1 F1 (α; β; −z), we can write 2 2 Z a u 1 3 u2 2 du ; ;− u exp 1 F1 2 2 2 2 0 Z a 3 u2 1 3 u2 u2 1 F1 1; ; = F ; ; − du 1 1 2 2 2 2 2 0 Z a X ∞ X ∞ (−1)ℓ (1)k (1/2)ℓ u2(k+ℓ) u2 = du (3/2)k (3/2)ℓ k!ℓ! 2k+ℓ 0 k=0 ℓ=0 Z a ∞ X ∞ X (−1)ℓ (1)k (1/2)ℓ u2+2(k+ℓ) du = (3/2)k (3/2)ℓ k!ℓ! 2k+ℓ 0 k=0 ℓ=0 k ℓ ∞ ∞ a2 /2 −a2 /2 (1)k (1/2)ℓ a3 X X = 2 (k + ℓ + 3/2) (3/2)k (3/2)ℓ k! ℓ! k=0 ℓ=0 k ℓ ∞ ∞ −a2 /2 a3 X X (3/2)k+ℓ (1)k (1/2)ℓ a2 /2 = . 3 (5/2)k+ℓ (3/2)k (3/2)ℓ k! ℓ! k=0 ℓ=0
The results now follows from the definition of the Kamp´e de F´eriet function. 2
Theorem 3.1 derives explicit expressions for (12) and (13). Its proof is immediate from Lemmas 3.2 to 3.4. Theorem 3.1. For β > 0, √ π β β2 β3 β2 3 1,1,1 F1,1,1 + 2 F2 1, 1; , 2; − √ Jβ = log 2 − erfi √ 4 2 2 2 2 3 2π For β < 0, β β2 π β2 β3 3 1,1,1 F1,1,1 Jβ = log 2 + erfi − √ + + √ 2 F2 1, 1; , 2; 4 2 2 2 2 3 2π √ Also J0 = log 2 (by the proof of Lemma 2.1). √
7
3/2 1 5/2 3/2
! β2 1/2 β 2 . , − 3/2 2 2
3/2 1 5/2 3/2
! β2 1/2 β 2 . , − 3/2 2 2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Combining Lemma 3.1 and Theorem 3.1, we obtain an explicit expression for (6). In-built routines for computing the imaginary error, hypergeometric and Kamp´e de F´eriet functions are available in packages like Maple, Matlab and Mathematica. See also [4], Section 1.3.2 and [5].
Acknowledgments The authors would like to thank the Editor for careful reading and comments which improved the paper.
References [1] Al Hanbali, A., Mandjes, M., Nazarathy, Y. and Whitt, W. (2011). The asymptotic variance of departures in critically loaded queues. Advances in Applied Probability, 43, 243-263. [2] Daley, D. J. (2011). Revisiting queueing output processes: A point process viewpoint. Queueing Systems, 68, 395-405. [3] Daley, D. J., Van Leeuwaarden, J. S. H. and Nazarathy, Y. (2015). BRAVO for many-server QED systems with finite buffers. Advances in Applied Probability, 47, 231-250. [4] Exton, H. (1978). Handbook of Hypergeometric Integrals: Theory, Applications, Tables, Computer Programs. Chichester, England. [5] Exton, H. and Krupnikov, E. D. (1998). A register of computer-oriented reduction identities for the Kamp´e de F´eriet function. Draft manuscript. [6] Hautphenne, S., Kerner, Y., Nazarathy, Y. and Taylor, P. (2015). The intercept term of the asymptotic variance curve for some queueing output processes. European Journal of Operational Research, 242, 455-464. [7] Nazarathy, Y. (2011). The variance of departure processes: Puzzling behavior and open problems. Queueing Systems, 68, 385-394. [8] Nazarathy, Y. and Weiss, G. (2008). The asymptotic variance rate of the output process of finite capacity birth-death queues. Queueing Systems, 59, 135-156. [9] Prudnikov, A. P., Brychkov, Y. A. and Marichev, O. I. (1986). Integrals and Series, volumes 1, 2 and 3. Gordon and Breach Science Publishers, Amsterdam. [10] R Development Core Team (2017). A Language and Environment for Statistical Computing: R Foundation for Statistical Computing. Vienna, Austria.
8