Laplace Transforms of Probability Density

Joseph Abate, Ward Whitt

چکیده

In order to numerically invert Laplace transforms to calculate probability distributions in queueing and related models, we need to be able to calculate the Laplace transform values. In many cases the desired Laplace transform values (e.g., of a waiting-time distribution) can be computed when the Laplace transform values of component probability density functions (pdf's) (e.g., of a service-time pdf) can be computed. However, in some cases explicit expressions for Laplace transforms of component pdf's are not available. Hence, we propose the construction of in nite-series representations for Laplace transforms of pdf's and show how they can be used to calculate transform values. We use the Laplace transforms of Laguerre functions, Erlang pdf's and exponential pdf's as basis elements in the series representation. We develop several speci c parametric families of pdf's in this in nite series framework. We show how to determine the asymptotic form of the pdf from the series representation and how to truncate so as to preserve the asymptotic form for a time of interest. Keywords: Laplace transforms, numerical inversion of Laplace transforms, moments, Laguerre polynomials, Laguerre-series representation, Laguerre transform, in nite series, long-tail probability distributions, Erlang transform, Bell numbers, asymptotics for tail probabilities. 1. Introduction Many descriptive quantities of interest in queueing models and other probability models arising in operations research can be e ectively computed by numerically inverting Laplace transforms; e.g., by using the Fourier-series method or the Laguerre-series method; see Hosono (1984), Abate and Whitt (1992, 1995), Choudhury, Lucantoni and Whitt (1994), Abate, Choudhury and Whitt (1996, 1997, 1998) and references therein. The most di cult step in performing the numerical inversion, if there is any di culty at all, is usually computing the Laplace transform values; e.g., see the application to polling models in Choudhury and Whitt (1996). In many cases, numerical inversion is straightforward provided that Laplace transforms are available for component probability density functions (pdf's). A familiar example is the steady-state waiting-time distribution in the M/G/1 queue. Numerical inversion can be applied directly to the classical Pollaczek-Khintchine transform provided that the Laplace transform of the service-time pdf is available. More generally, the steady-state waiting-time distribution in the GI/G/1 queue can be computed by numerical transform inversion provided that the Laplace transforms of both the interarrival-time pdf and the service-time pdf are available; see Abate, Choudhury and Whitt (1993, 1994). In the GI/G/1 case, an extra numerical integration is required to calculate the required waiting-time transform values. For these inversion algorithms to be e ective, the pdf's also need to be suitably smooth; otherwise preliminary smoothing may need to be performed; see Section 6 of Abate and Whitt (1992). In this paper, all the pdf's considered will be continuous. Another example is the renewal function. Since the Laplace transform of the renewal function can be expressed directly in terms of the Laplace transform of the interrenewal-time pdf, the renewal function can be computed by numerical inversion provided that the Laplace transform of the interrenewal-time pdf is available; see Section 13 of Abate and Whitt (1992). A remaining di culty, however, is that convenient closed-form Laplace transforms are not readily available for some pdf's. This is particularly true for long-tail (also called heavytail) pdf's, which are currently of interest to model features of communication networks. To address this problem, Abate, Choudhury and Whitt (1994) introduced a family of long-tail distributions with convenient Laplace transforms, the Pareto mixtures of exponential (PME) distributions. To provide additional possibilities, Abate and Whitt (1996) initiated the development of an operational calculus for manipulating Laplace transforms, and Abate and Whitt 1 (1998) introduced two classes of beta mixtures of exponential (BME and B2ME) distributions. Feldmann and Whitt (1997) also provided an algorithm for approximating probabilities density functions with decreasing failure rate, which include the Pareto, Weibull and other long-tail distributions, by hyperexponential distributions ( nite mixtures of exponential distributions), whose transforms are readily available. Asmussen, Nerman and Olsson (1996) and Turin (1996) showed that the expectation-maximization (EM) algorithm can often be used to t phase-type distributions to general probability distributions. The phase-type distributions also have relatively convenient Laplace transforms. In this paper we make further contributions in this direction. We suggest a general approach for constructing probability density functions for which the Laplace transform values can be computed. In particular, we suggest representing the Laplace transform as an in nite series and then numerically calculating the sum, using acceleration methods if necessary. Perhaps the main point is to recognize that for the inversion we do not actually need a convenient closed-form expression for the transform. It su ces to have an algorithm to compute the transform value f̂(s) for the required arguments s. Thus an in nite series can be a satisfactory representation. The rst method we propose is the Laguerre-series representation. In Abate, Choudhury and Whitt (1996) we studied the Laguerre-series representation as a tool to numerically invert Laplace transforms, given that the Laplace transform values are available. Our purpose here is di erent. Here we want to compute Laplace transform values of a Laplace transform of a component pdf, such as a service-time pdf in a queueing model. We do so as an intermediate step to compute the nal Laplace transform values of some other pdf or cdf, perhaps of a waiting-time distribution, which will be used to perform a numerical inversion. The nal inversion may be done by any method, possibly but not necessarily by the Laguerre-series method. Given that we decide to seek a Laguerre-series representation, there are two problems: (1) calculating the Laguerre coe cients, and (2) calculating the sum for the Laplace transform. It turns out that both problems can be addressed very e ectively when the pdf is shorttailed, in particular, when the rightmost singularity of its Laplace transform is strictly less than 0. For this we use two approaches. The rst way is to compute the coe cients by numerical integration using the standard integral representation of the coe cients, assuming that the pdf is known. The second way is to compute the coe cients from the moments, assuming that the moments are known. In exploiting the moments, we follow Keilson and 2 Nunn (1979), who showed that the Laguerre coe cients can be calculated from the moments under a regularity condition. However, we obtain a more elementary connection between the Laguerre coe cients and moments by altering the way the Laguerre coe cients are de ned. The resulting Laguerre-series representation for the Laplace transform is convenient because the series tends to converge geometrically fast, even when the corresponding Laguerre-series in the time domain converges slowly. Moreover, it is sometimes possible to analytically determine the Laguerre coe cients, as we illustrate here. In this work, as in Abate, Choudhury and Whitt (1996), we draw on previous studies of Laguerre-series representation by Keilson and Nunn (1979), Keilson, Nunn and Sumita (1981), Sumita (1984) and Sumita and Kijima (1988), but we do not focus on manipulations that can be performed on the Laguerre coe cients, i.e., the Laguerre transform calculus. Applications of the Laguerre transform calculus are described by Keilson and Sumita (1982) and Litko (1989). These authors show that the Laguerre transform calculus can be e ective, but much can also be done with Laplace transforms. We also propose two other series representations for Laplace transforms of pdf's that are useful for long-tail pdf's. The rst is the Erlang-series representation, which was introduced by Keilson and Nunn (1979) as the Erlang transform. Then the nth basis function is an Erlang pdf of order n (En) or its Laplace transform (1+s) n. The coe cients of the Erlang transform are obtained from the coe cients of a power series representation of the pdf, assuming that it exists. The Erlang-series representation is a convenient alternative to the Laguerre-series representation because it often applies when the Laguerre-series representation does not and because the Erlang series for the Laplace transform tends to converge even more rapidly. The nal series representation we consider is the exponential-series representations, in which the nth basis element is an exponential pdf with mean an or its Laplace transform (1+ans) 1, where an ! 1 as n ! 1. Since the Laplace transform has poles at s = 1=an for all n, the Laplace transform must have its rightmost singularity at 0, so these are long-tail distributions. We show how the exponential-series representation enables us to construct pdf's with given moment sequences. To illustrate, we construct pdf's with the famous Bell numbers and the ordered Bell numbers as moment sequences; e.g., see pp. 20, 42, 175 of Wilf (1994). We also show how to determine the asymptotic form of a pdf represented by an in nite series and how to truncate the in nite series so that the asymptotic form is preserved for times of interest. In particular, we apply the Euler-Maclaurin formula and Laplace method to determine the asymptotics and the appropriate truncation point (as a function of time) to 3 capture this asymptotic form. Here is how the rest of this paper is organized. In Section 2 we present the Laguerre-series representation that we will use and derive our new algorithm for computing the Laguerre coe cients from the moments. In Section 3 we discuss the application of this algorithm to develop Laguerre-series representations of Laplace transforms of mixtures of exponential pdf's. An important class for which explicit results hold is the beta mixture of exponential (BME) pdf's. Properties of BME distributions, including this representation, were established in Abate and Whitt (1998). The BME example is very nice because we are able to apply the moment algorithm analytically to obtain simple expressions for the Laguerre coe cients. The moment algorithm is genuinely useful for BME pdf's because neither the pdf nor its transform has a convenient closed form. In Section 4 we discuss the Erlang-series representation for pdf's. We illustrate in Section 5 by establishing results for mixtures of exponential distributions, including an explicit Erlangseries expansion for a third class of beta mixtures of exponentials. This class includes a class of Pareto mixtures of exponentials pdf's with a pure power tail, i.e., for which the next term in the asymptotic expansion decays exponentially. We consider the exponential-series representations in Section 6. In Section 7 we show how to determine the asymptotic form of the pdf and how to truncate the in nite series to preserve the asymptotic form at times of interest. Finally, we draw conclusions in Section 8. We close this introduction by mentioning that multivariate extensions are possible. For work on multidimensional Laguerre-series representations, see Sumita (1984), Sumita and Kijima (1985) and Abate, Choudhury and Whitt (1997). 2. Laguerre-Series Representations Let f be a pdf and f̂ its Laplace transform, i.e., f̂(s) = Z 1 0 e stf(t)dt ; (2.1) which is well de ned and analytic for all complex s with Re(s) > 0; see Doetsch (1974). If f̂ is analytic for Re(s) > s for s > 0, then f̂ is analytic at 0 and has the power series (about the origin) f̂ (s) = 1 Xn=0 sn f (n)(0) n! = 1 X n=0( s)nmn n! ; (2.2) where mn is the nth moment of f . Unfortunately, though, even if (2.2) is valid and the moments mn are known, (2.2) is usually not a good way to compute f̂(s). Indeed, we may 4 want to compute f̂(s) for values of s outside the radius of convergence of (2.2). For all s with Re(s) > s , f̂ (s) will be analytic, so that a power-series representation for f̂ (s) about each of those points s is possible, but the domain of convergence for the power series about the origin need not include these points. A Laguerre-series representation may be a much better series representation for computation of the Laplace transform. The Laguerre-series representation of f is f(t) = 1 X n=0 qnLn(t); t 0 ; (2.3) where qn are real numbers andLn(t) n Xk=0 nk ( t)k k! ; t 0 ; (2.4) are the Laguerre polynomials, e.g., see 22.3.9 of Abramowitz and Stegan (1972). The Laguerre polynomials form a complete orthonormal basis for the space L2[0;1) of real-valued functions on the nonnegative real line [0;1) that are square integrable with respect to the weight function e t, using the inner product hf1; f2i = Z 1 0 e tf1(t)f2(t)dt (2.5) and norm kfk = hf; fi1=2. Then (2.3) is valid in the sense of convergence in L2[0;1) with qn = hf; Lni Z 1 0 e tf(t)Ln(t)dt (2.6) being the Laguerre coe cient in (2.3); e.g., see Chapter 22 of Abramowitz and Stegun (1972) and Szeg o (1975). Assuming that f is continuous, the series (2.3) converges pointwise for all t as well. The Laguerre polynomial Ln(t) can be calculated in a numerically stable way via the recursion Ln(t) = 2n 1 t n Ln 1(t) n 1 n Ln 2(t) : (2.7) Our interest, however, is in the associated series representation for the Laplace transform f̂ . Since the Laplace transform of Ln(t) has the simple form L̂n(s) = Z 1 0 e stLn(t)dt = (s 1)n sn+1 ; (2.8) the associated Laguerre-series representation for f̂ is f̂(s) = 1 Xn=0 qn Z 1 0 e stLn(t)dt = 1 X n=0 qn (s 1)n sn+1 ; (2.9) 5 where again qn is the Laguerre coe cient in (2.3) and (2.6). Above we have used the Laguerre polynomials and the weight e t. An essentially equivalent approach is to use the associated orthonormal Laguerre functions e t=2Ln(t) without a weight. Although that approach is essentially the same, the resulting coe cients qn are di erent. In order to exploit the Laguerre-series representation for computing the Laplace transform, we primarily work with the Laguerre series representation of etf(t). This requires that etf(t) be square integrable with respect to the weight e t or, equivalently, that et=2f(t) be square integrable directly. That, in turn, clearly requires that f(t) be asymptotically dominated by the exponential e t=2, which means that its Laplace transform f̂(s) should have its rightmost singularity s satisfy s 0 because the Laguerre coe cients are bounded and js=(s+1)j 0. The Laguerre coe cients are bounded because jLn(t)j et=2 for t 0, see 22.14.12 on p. 786 of Abramowitz and Stegun (1972), and jqnj Z 1 0 jf(t)jjLn(t)jdt Z 1 0 et=2jf(t)jdt M 0 and each N 1, jf̂(s) f̂N (s)j M 1 s+ 1 N 1 ; (2.21) where fN (s) = N X n=0 1 s + 1qn s s+ 1 n (2.22) with qn begin the Laguerre coe cients in (2.11) and (2.12), and = 4 + 2(1 + 2) + 2 4 + 2 2 + 1 1=2 1=2, so that scaling is needed in this case. As indicated earlier, the nal decay rate has to be at least 1/2. From numerical experience, we found that the scaling becomes more e ective in this example as the nal decay rate approaches 1/2. Hence, we use the scale 11 parameter = e=2, making the ccdf Gc(t) F c( t) with momentsmn(G) = (n+1)n+1(2=e)n. Then Gc(t) 2e p t 3=2e t=2 as t!1 (2.38) and mn(G) 2e p (2n) 3=22n as n!1 : (2.39) For the ccdf, we obtain the Laguerre-series expansion Gc(t) = 1 Xn=0(qn qn+1)Ln(t) (2.40) where qn can be obtained from the moments using (2.12). We now compute Gc(t) for several values of t using the method of moments. As in Abate, Choudhury and Whitt (1996), we also use Wynn's (1956, 1966) -algorithm to accelerate convergence. In particular, we use the fourth-order version; i.e., we use en2m for m = 4 and n = 80. This requires that we compute qn for 0 n 89. To determine the correct digits, we redo the computation for n = 100, and keep only the digits that agree. The numerical results are displayed in Table 1. Table 1 shows that we obtain four correct digits except for the largest values of t considered, t = 30. The asymptotic values from (2.38) are also displayed in Table 1. As should be anticipated, the asymptotics do not perform too well for the displayed times because the next term is smaller only by a factor t 1; see Abate and Whitt (1997) for further discussion. However, the asymptotics would be useful for larger times. time asymptotic t Gc(t) values in (2.38) :1 :7563 :5 :4110 1 :2346 2 :9334 e 1 3 :4178 e 1 5 :9843 e 2 7 :2599 e 2 :501 e 2 10 :3957 e 3 :654 e 3 15 :2029 e 4 :292 e 4 20 :1169 e 5 :156 e 5 25 :7224 e 7 :914 e 7 30 :466 e 8 :571 e 8 Table 1: Numerical values of the CEP ccdf computed from the Laguerre-series representation using the method of moments. 12 Because of the combinatorial sum, the calculation of qn from the moments requires about (n 225n=2 + 15) digits of precision. Hence, for n = 80, we need about 70 digits precision. For this purpose, we used UBASIC; see Kida (1990). 3. Mixtures of Exponential Distributions In this section we obtain some results about Laguerre-series representations for Laplace transforms of completely monotone pdf's, i.e., for pdf's that are mixtures of exponential pdf's. We rst obtain some general results. Then we describe an application of the method of moments in Theorem 2.1 to obtain a Laguerre-series representation for a class of pdf's called a beta mixture of exponential (BME) pdf's. These BME pdf's are investigated further in Abate and Whitt (1998). The important point here is that the method of moments can be applied analytically to determine simple expressions for the Laguerre coe cients. We rst establish integral representations for the Laguerre coe cients for general mixtures of exponentials. Theorem 3.1. If a pdf f can be represented as f(t) = Z b a xe xtdG(x) (3.1) for a cdf G, where 1=2 0, q > 0 and a third positive scale parameter, which has been omitted from (3.6). (The parameter q should not be confused with the Laguerre coe cients qn, which have the subscript.) In general, a more explicit form for the BME pdf is not known, but special cases are described in Abate and Whitt (1998). Especially tractable are the BME pdf's where both p and q are integer multiples of 1/2.The Laplace transform v̂(p; q; s) has rightmost singularity s = 1. The BME pdf is class II with asymptotic form v(p; q; t) (p+ q)e t (p)tq as t!1 (3.8) and nth moment mn(p; q) = (p)nn! (p+ q)n (3.9) where (x)n = x(x+ 1) : : :(x+ n 1) is the Pochammer symbol with (x)0 = 1. Theorem 2.1, (3.9) and a combinatorial identity, (7.1) on p. 58 of Gould (1972), yields the Laguerre coe cients in (2.11), as shown in Theorem 2.2 and 2.3 of Abate and Whitt (1998); in particular, qn = (q)n (p+ q)n ; n 0 : (3.10) Since 1 Y has a beta pdf with parameter pair (q; p) when Y has a beta pdf with parameter pair (p; q), we could also apply Theorem 3.2 to obtain (3.10). 15 Given the explicit expression for the Laguerre coe cients qn in (3.10), we can apply the series representation (2.17) to compute the BME transform values, but in this instance it turns out to be more e ective to use continued fractions; see Abate and Whitt (1998a). We can use the Laguerre-series representation for the BME pdf's to obtain Laguerre-series representations for some related pdf's. In Abate and Whitt (1998) a second class of beta mixtures of exponential (B2ME) is de ned with pdf v2(p; q; t) = (p+ q) (p) (q) Z 1 0 y 1e t=yyp 1(1 + y) (p+q)dy : (3.11) In (3.11) the mixing pdf is a beta pdf of the second kind. However, by (1.21) of Abate and Whitt (1998), the B2ME pdf is an exponentially undamped version of a BME pdf, i.e.,v2(p; q; t) = q p+ q etv(p; q+ 1; t); t 0 ; (3.12) so that v2(p; q; t) q (p+ q) (p)tq+1 as t!1 : (3.13) From (3.13), we see that the B2ME pdf v2(p; q; t) is a long-tail (class III) pdf. From (3.13) we can obtain Laguerre-series representations for v2(p; q; t) and its Laplace transform v̂2(p; q; s); in particular, from (3.10) and (3.12), we obtain v2(p; q; t) = 1 X n=0 q p+ q (q + 1)n (p+ q + 1)nLn(t) (3.14) and v̂2(p; q; s) = 1 Xn=0 q p+ q (q + 1)n (p+ q + 1)n (s 1)n sn+1 : (3.15) Unfortunately, however, as noted in Section 2, the Laguerre-series representation (3.15) is not as useful as (2.17) with (3.10). Because of the undamping, s in (2.17) has been shifted to s 1. 4. The Erlang-Series Representation In this section, following Section 4 of Keilson and Nunn (1979), we consider another series representation for the Laplace transform f̂ of a pdf f based on a Taylor series expansion of g(t) etf(t), i.e., g(t) etf(t) = 1 Xn=0 g(n)(0) tn n! ; t 0 ; (4.1) 16 where g(n)(0) is the nth (right) derivative of g evaluated at 0 and the coe cients f (n)(0) and g(n)(0) are related by g(n)(0) = n Xk=0 nk f (n k)(0) : (4.2) Then f(t) can be rewritten asf(t) = 1 X n=1 g(n 1)(0)e(n; t); t 0 ; (4.3) where, for n 1, en(t) is the Erlang (En) pdf e(n; t) = e ttn 1 (n 1)! ; t 0 ; (4.4) which has mean n, variance n and Laplace transform ê(n; s) = (1+s) n . As n increases e(n; t) approaches a normal cdf with mean and variance n. The squared coe cient of variation (SCV, variance divided by the square of the mean) of e(n; t) is 1=n, showing that the variability is asymptotically negligible compared to the mean. From (4.3) we immediately obtain an associated Erlang-series representation for the Laplace transform, i.e., f̂(s) = 1 Xn=1 g(n 1)(0)(1 + s) n : (4.5) Closely paralleling Theorem 2.3, we can easily relate the Erlang-series representation of pdf's and cdf's. Theorem 4.1. Let g be a pdf with ccdf Gc. Then the following are equivalent: (i) Gc(t) = 1 P n=1 cnen(n; t) (ii) g(t) = 1 P n=1(cn cn+1)e(n; t). Proof. Note that ĝ(s) = 1 +Gc(s) (s + 1)Gc(s). Hence, Gc(s) = 1 Xn=1 cn(1 + s) n if and only if ĝc(s) = 1 X n=1 dn(1 + s) n with P1n=1 dn = c1 for the coe cients as indicated. Remarkably, many exponential pdf's can be represented as a mixtures of Erlangs. Theorem 4.2. The Erlang coe cients are cn = (1 )n 1 for > 0, if and only if the pdf f is exponential with mean 1= . 17 Proof. Starting with the exponential pdf e t with 6= 1 we get the expansion for the Laplace transform f̂ (s) = 1 1 X n=1 (1 )n (1 + s)n = + s ; (4.6) and vice versa. For = 1, we get c1 = 1 and cn = 0 for n 1. For > 1, the coe cients in (4.6) alternate in sign; otherwise they are probabilities. Corollary. A pdf f is a geometric probabilistic mixture of Erlang pdf 's if and only if f is an exponential pdf with mean m1 > 1, i.e., f(t) = e t for 0 0 (5.12) so that f(t) must be completely monotone at 0. Example 5.1. In this example, the conditions of Theorem 5.1 and Lemma 5.1 are not satis ed, so we run into di culties. Consider Feller's second Bessel pdf for the case r = 0, as discussed in Section 11 of Abate and Whitt (1996), f(t) = e (1+t)I0(2pt); t 0 ; (5.13) 21 which is not completely monotone, because f 00(t) (1 4t=3)=2 as t ! 0. We can see that f(t) has a generalized Laguerre-series representation from the Laplace transform f̂(s) = 1 1 + s exp( s=(1 + s)) = 1 1 + s 1 X n=0 ( 1)n n! s 1 + s n : (5.14) We can construct the transform Ĝc(s) = s 1f̂(s 1) = 1 1 + s exp( 1=(1 + s)) (5.15) which has inverse Gc(t) = e tJ0(2pt); t 0 ; (5.16) which is not monotone and thus not a bona de ccdf. Example 5.2. As a quick illustration of Theorem 5.1, we draw on Example 3.1 and let the mixing density be w(x) = log(1 x), 0 x 1. Then (5.1) and (5.3) are valid with qn = 1=(n+ 1)2. Example 5.3. We now illustrate how we can apply Theorem 5.1 to characterize the dual ccdf Gc(t) de ned in (5.2). Suppose that f(t) = p2e 3t=2I0(t=2); t 0 ; (5.17) where I0(t) is the Bessel function. We can then show that f(t) = Z 1 1=2 y 1e t=y dy p1 ypy 1=2 = e t 1 Xn=0 8 n 2n n Ln(t) : (5.18) Thus the dual ccdf is Gc(t) = Z 1 1=2 e tx dx p1 xpx 1=2 = Z 2 1 e t=y dy ypy 1p1 y=2 : (5.19) From Theorem 5.1, Gc(t) = 1 Xk=1 8 (k 1) 2k 2 k 1 e(k; t) : (5.20) By Lemma 5.1, Ĝc(s) = s 1f̂(s 1) = 1 p1 + sp(1=2) + s (5.21) from which we can deduce that Gc(t) = e 3t=4I0(t=4); t 0 (5.22) see 29.3.49 of Abramowitz and Stegun (1972).22 Next we apply Theorem 5.1 to develop an Erlang-series representation for the third beta mixture of exponential (B3ME) pdf v3(p; q; t) = Z 1 0 xe txb(p; q; x)dx ; (5.23) for b(p; q; x) in (3.7), which is dual to the BME pdf in (3.6); i.e., it is dual in that the exponential means are averaged in (3.6) while the rates are averaged in (5.23), both with respect to the same standard beta pdf. The BME pdf (3.6) has a Laguerre-series representation, but no Erlang-series representation. Paralleling the connection between the BME pdf and the Tricomi con uent hypergeometric function U(a; b; z) established in Theorem 1.7 of Abate and Whitt (1998), we now establish a connection between the B3ME pdf and the Kummer con uent hypergeometric function M(a; b; z); see Chapter 13 of Abramowitz and Stegun (1972). From the integral representation 13.2.1 there, we obtain the following result. Theorem 5.2. For all p > 0 and q > 0, v3(p; q; t) = p p+ qM(p+ 1; p+ q + 1; t); t 0 : (5.24) Proof. Starting from the de nition (5.23), we obtain v3(p; q; t) = Z 1 0 xe txb(p; q; x)dx = p p+ q e tM(q; p+ q + 1; t) = p p+ qM(p+ 1; p+ q + 1; t) : Remark 6.1. The B3ME pdf in the form (5.24) was introduced by Mathai and Saxena (1966); see (67) on p. 32 of Johnson and Kotz (1970). We can apply Theorem 5.2 to describe the asymptotic form of the B3ME tail probabilities. We see that a B3ME pdf has a long tail. Corollary. For each p > 0 and q > 0, v3(p; q; t) (p+ q)p (q)tp+1 as t!1 : (5.25) As before, it is signi cant that we can evaluate the coe cients analytically. The transform series also tend to converge quite rapidly. 23 Theorem 5.3. For each p > 0 and q > 0, v3(p; q; t) has an Erlang-series representation, in particular, v3(p; q; t) = p 1 Xn=1 (q)n 1 (p+ q)n e(n; t) (5.26) and v̂3(p; q; s) = p 1 Xn=1 (q)n 1 (p+ q)n (1 + s) n : (5.27) First Proof. Make the change of variables x = 1 y in (5.23) to obtain v3(p; q; t) = (p+ q) (p) (q)e t Z 1 0 etyyq 1(1 y)pdy : (5.28) Now expand ety in a power series and integrate term by term. Second Proof. Apply (5.4) in Theorem 5.1 with (3.10). We get the Erlang-series representation with qn 1 qn = (q)n 1 (p+ 1)n 1 (q)n (pq)n = p(q)n 1 (p+ q)n : (5.29) Let v3e(p; q; t) m(p; q) 1V c 3 (p; q; t) denote the stationary-excess pdf associated with v3(p; q; t). Just as for the BME pdf, we can integrate to see that the B3ME stationary-excess pdf is again a B3ME pdf with di erent parameters. Theorem 5.4. For each p > 1 and q > 0, v3e(p; q; t) = v3(p 1; q; t); t 0 : (5.30) Example 5.4. Consider the BME pdf v(2; 1; t) = Z 1 0 y 1e t=y(2y)dy = e t 1 X n=0 2 (n+ 1)(n+ 2)Ln(t) : The dual cdf is V c 3 (2; 1; t) = Z 1 0 e tx(2x)dx = Z 1 1 e t=y2y 3dy : From Theorem 5.1, V c 3 (2; 1; t) = 1 Xk=1 2 k(k + 1)e(k; t) : Also, since v̂(2; 1; s) = 2s 1 2s 2 log(1 + s), V̂ c 3 (2; 1; s) = 2 2s log(1 + s 1) ; so that V3(2; 1; t) = 2(1 (1 + t)e t); t 0 : 24 We now observe that the special case of a B3ME pdf with q = 1 corresponds to a Pareto mixture of exponentials (PME); the PME is a minor modi cation of PME considered in Abate, Choudhury and Whitt (1994); a scale factor there made the mean 1. For each r > 0, let the Pareto ccdf be F c(r; y) = y r; y > 1 : (5.31) The associated Pareto pdf is f(r; y) = ry (r+1); y > 1 : (5.32) The nth moment is mn(F ) = r=(r n) for n 0 v3(p; 1; t) p (p+ 1) tp+1 pe t t 1 + pt + p(p 1) t2 + as t!1 : (5.36) 25 Proof. Apply 6.5.3 and 6.5.32 of Abramowitz and Stegun (1972). We now relate the ccdf to the pdf. Theorem 5.6. For each p > 1, V c 3 (p; 1; t) = p p 1v3(p 1; 1; t) = t pv3(p; 1; t) + e t; t 0 : (5.37) Proof. Start with the second relation in (5.33). Then use the rst relation in (5.33) with 6.5.22 of Abramowitz and Stegun (1972). 6. Exponential-Series Representations We now consider a simple alternative to the Erlang-series representation in Section 5: an exponential-series representation, which is just a countably in nite mixture of exponential pdf's, i.e., f(t) = 1 Xk=1 pk e t=ak ak ; t 0 ; (6.1) and f̂ (s) = 1 Xk=1 pk(1 + sak) 1 ; (6.2) where fpk : k 1g is a probability mass function (pmf) and fak : k 1g is the sequence of means of the component exponential pdf's. The standard case is ak 0. The idea is that, for suitably large t, the integral will be dominated by the contribution in the neighborhood of x (t). Then Laplace's method yields Z 1 0 e (x;t)dx s 2 00(x (t); t)e (x (t);t) as t!1 : (7.5) In summary, the asymptotic form is given by (7.5) and the truncation point as a function of t is x (t). In particular, if we are interested in time t0, then the truncation point should be at least x (t0). Example 7.1. Suppose that the weights in (7.1) are pn = (1 ) n (a geometric pmf with an atom (1 ) at 0) and Fn is exponential with mean an = bnc. For the special case = 1=2, b = 1 and c = 1, we obtain the second Bell pdf in Example 6.2. In general, by (7.3), (n; t) = log(1 ) n log + t bnc : (7.6) As indicated above, x (t) is the root of 0(x; t) = log ct bxc+1 = 0 ; (7.7) so that x (t) = ct b( log ) 1=(c+1) ; (7.8) 00(x; t) = c(c+ 1)t bx(c+2) ; (7.9) 28 00(x (t); t) = (1 + c) b( log )2+c ct 1=(1+c) ; (7.10) and F c(t) At1=2(1+c)e Bt1=(1+c) as t!1 ; (7.11) where A = (1 ) 2 c1=(1+c) (1 + c)b1=(1+c)( log )(2+c)=(1+c)!1=2 (7.12) and B = 1 + 1c ( log )cc b 1=(1+c) : (7.13) In the special case = 1=2, b = c = 1 corresponding to Example 6.2, F c(t) At1=4e Bt1=2 as t!1 ; (7.14) where A = p 2(log 2)3=4 and B = 2plog 2 : (7.15) For the case of Example 6.2 we can verify the asymptotic by applying Section 6 of Abate, Choudhury, Lucantoni and Whitt (1995), but we must correct errors in formulas (6.8) and (6.9) there. Assuming that F c(t) t e t as t ! 1, we obtain estimators for the parameters , , and from the moments. First, by p. 176 of Wilf (1994), the ordered Bell numbers satisfy ~b(n) n!=2an+1 as n! 1 for a = log 2. Therefore, by (6.8), mn 1 2 log 2 (n!)2 (log 2)n as n!1 ; (7.16) rn mn mn 1 n2 a as n!1 ; (7.17) n rn r(rn+1 rn) ! 1 2 as n!1 ; (7.18) n n r n ! 2pa as n!1 ; (7.19) n ( =n)1= nrn n+ ( 1 1)=2 ! 1=4 as n!1 ; (7.20) n mn (n+ ) ((n + )= + 1) ! p 2a3=4 as n!1 : (7.21) Note that (7.18){(7.21) agree with (7.15). Formulas (7.20) and (7.21) here correct (6.8) and (6.9) in Abate, Choudhury, Lucantoni and Whitt (1995). 29 Example 7.2. We now produce a pdf with a power tail. Suppose that an = n for > 0 and pn = n for > 0 and > 1. Then (x; t) = log + log x+ t x ; (7.22) so that x (t) = ( t= )1= ; (7.23) 00(x; t) = x2 + ( + 1)t x +2 0 (7.25) and F c(t) At ( 1)= as t!1 (7.26) for A = e = ( = ) ( 1)= : (7.27) If = 1, then F c(t) e t 1 as t!1 (7.28) and mn = n! 1 Xk=1 n = n! ( n) ; (7.29) where (z) is the Riemann zeta function; e.g., see p. 19 of Magnus et al. (1966). The moment mn is nite for n the Laplace transform values in applications. At the same time, other properties of the pdf's can be deduced, such as moments and the asymptotic form, so that the pdf can be chosen to satisfy desired properties in applications. For example, exponential series representations always yield completely monotone pdf's whose moments can be computed. We showed how to determine the asymptotic form in Section 7. By imposing additional structure, we obtain relatively simple parametric families of pdf's such as the BME pdf's in (3.6), the B3ME pdf's in (5.23), the PME pdf's in (5.33), the generalized Bell pdf's in Example 6.1, and the geometric mixture of exponential pdf's in Example 7.1. Since the Laplace transform values for these pdf's can be computed, these pdf's are natural candidates to use in applications. Then the parameters can be chosen to match the moments, the asymptotic form and other features of interest. Toward that end, we studied the BME pdf's further in Abate and Whitt (1998). The connection to the con uent hypergeometric function M(a; b; z) established in Theorem 5.2 makes it possible to establish analogous results for B3ME pdf's. Over all, our analysis help extend the analyst's toolkit. 32 References Abate, J., Choudhury, G. L., Lucantoni, D. M. and Whitt, W. (1995). Asymptotic analysis of tail probabilities based on the computation of moments. Ann. Appl. Prob. 5, 983{1007. Abate, J., Choudhury, G. L. and Whitt, W. (1993). Calculation of the GI/G/1 waiting- time distribution and its cumulants from Pollaczek's formulas. Archiv f ur Elcktronik und Ubertragungstechnik 47, 311{321. Abate, J., Choudhury, G. L. and Whitt, W. (1994). Waiting-time tail probabilities in queues with long-tail service-time distributions. Queueing Systems 16, 311{338. Abate, J., Choudhury, G. L. and Whitt, W. (1996). On the Laguerre method for numerically inverting Laplace transforms. INFORMS J. Computing 8, 413{427. Abate, J., Choudhury, G. L. and Whitt, W. (1997). Numerical inversion of multidimensional Laplace transforms by the Laguerre method. Performance Evaluation 31, 229{243. Abate, J., Choudhury, G. L. and Whitt, W. (1998). An introduction to numerical trans- form inversion and its application to probability models. Advances in Computational Probability, W. Grassman (ed.), Kluwer, Boston, to appear. Abate, J. and Whitt, W. (1992). The Fourier-series method for inverting transforms of probability distributions. Queueing Systems 10, 5{88. Abate, J. and Whitt, W. (1995). Numerical inversion of Laplace transforms of probability distributions. ORSA J. Computing 7, 36{43. Abate, J. and Whitt, W. (1996). An operational calculus for probability distributions via Laplace transforms. Adv. Appl. Prob. 28, 75{113. Abate, J. and Whitt, W. (1997). Asymptotics for M/G/1 low-priority waiting-time tail probabilities. Queueing Systems 25, 173{233. Abate, J. and Whitt, W. (1998). Modeling service times with densities having non-exponential tails, submitted. Abate, J. andWhitt, W. (1998a). Computing transforms for numerical inversion via continued fractions. AT&T Labs, Florham Park, NJ. 33 Abramowitz, M. and Stegun, I. A. (1972). Handbook of Mathematical Functions, National Bureau of Standards, Washington, D. C. Asmussen, S., Nerman, O. and Olsson, M. (1996). Fitting phase type distributions via the EM algorithm. Scand. J. Statist. 23, 419{441. Choudhury, G. L. Lucantoni, D. M. and Whitt, W. (1994). Multidimensional transform inversion with applications to the transient M/G/1 queue. Ann. Appl. Prob. 4, 719{ 740. Choudhury, G. L. and W. Whitt (1996). Computing distributions and moments in polling models by numerical transform inversion. Performance Evaluation 25, 267{292. Daley, D. J. and Jacobs, D. R., Jr. (1969). The total waiting time in a busy period of a stable single-server queue, II, J. Appl. Prob. 6, 565{572. Doetsch, G. (1974). Introduction to the Theory and Application of the Laplace Transforma- tion, Springer, New York. Feldmann, A. and Whitt, W. (1997). Fitting mixtures of exponentials to long-tail distribu- tions to analyze network performance models. Performance Evaluation 31, 245{279. Gould, H. W. (1972). Combinatorial Identities, Morgantown, WV. Guillemin, F. and Pinchon, D. (1998). On the area swept under the occupation process of an M/M/1 queue in a busy period, France Telecom/CNET. Hardy, G. H. (1991). Divergent Series, Chelsea, New York (reprinting of 1949 Oxford edition). Hosono, T. (1984). Fast Inversion of Laplace Transform by BASIC, Kyoritsu Publishers, Japan (in Japanese). Johnson, N. L. and Kotz, S. (1970). Continuous Univariate Distributions 1, Wiley, New York. Keilson, J. and Nunn, W. R. (1979). Laguerre transformation as a tool for the numerical solution of integral equations of convolution type. Applied Math. Computation 5, 313{ 359. Keilson, J., Nunn, W. and Sumita, U. (1981). The bilateral Laguerre transform. Applied Math. Computation 8, 137{174. 34 Keilson, J. and Sumita, U. (1982). Waiting time distribution response to tra c surges via the Laguerre transform. In Applied Probability | Computer Science: The Interface, R. L. Disney and T. J. Ott (eds.), Birkh auser, Boston, 109{130. Kida, Y. (1990). UBASIC Version 8.12, Faculty of Science, Kanazawa University, 1-1 Marunouchi, Kanazawa 920, Japan. Litko, J. R. (1989). GI/G/1 interdeparture time and queue-length distributions via the Laguerre transform. Queueing Systems 4, 367{382. Magnus, W., Oberhettinger, F. and Soni, R. P. (1966). Formulas and Theorems for the Special Functions of Mathematical Physics, Springer, New York. Mathai, A. M. and Saxena, R. K. (1966). On a generalized hypergeometric distribution. Metrika 11, 127{132. Olver, F. W. J. (1974). Asymptotics and Special Functions, Academic Press, New York. Sumita, U. (1984). The Laguerre transform and a family of functions with nonnegative Laguerre coe cients. Math. Oper. Res. 9, 510{521. Sumita, U. (1984a). The matrix Laguerre transform. Applied Math. Computation 15, 1{28. Sumita, U. and Kijima, M. (1985). The bivariate Laguerre transform and its applications: numerical exploration of bivariate processes. Adv. Appl. Prob. 17, 683{708. Sumita, U. and Kijima, M. (1988). Theory and algorithms of the Laguerre transform, I. Theory. J. Opns. Res. Soc. Japan 31, 467{494. Szeg o, G. (1975). Orthogonal Polynomials, 4th ed., American Math. Soc. Colloq. Pub. 23, Providence, RI. Turin, W. (1996). Fitting probabilistic automata via the EM algorithm. Stochastic Models 12, 405{424. Wilf, H. W. (1994). Generatingfunctionology, second ed., Academic Press, San Diego. Wynn, P. (1956). On a device for computing the em(Sn) transformation. Math. Tables Aids Comput. 10, 91{96. Wynn, P. (1966). On the convergence and stability of the epsilon algorithm. SIAM J. Numer. Anal. 3, 91{122. 35