Laplace Transforms of Probability Density

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 < 1=2. This method thus applies to short-tail pdf's, but not longtail pdf's. In the terminology of Abate, Choudhury and Whitt (1994) and Abate and Whitt (1997), the algorithm applies to pdf's in classes I and II but not class III. Class I contains pdf's with a pure-exponential tail, in particular, pdf's f for which the rightmost singularity of its Laplace transform f̂ is at s < 0 and f̂( s ) = 1. For class II, again s < 0, but f̂( s ) < 1. The rightmost singularity of f̂(s) is 0 for class III. Class II routinely arises for busy-period distributions and low-priority waiting-time distributions in queueing; e.g., see Abate and Whitt (1997). Remark 2.1. Some pdf's f such as the gamma pdf with shape parameter 1/2, i.e., (1=2; t) e t=p t, are not square integrable with respect to the exponential weight because f(t) tp as t ! 0 for 1=2 p < 1. As a consequence, such pdf's do not t directly into the L2 theory. Indeed, then P1n=1 q2 n = 1, as can be seen from Examples 2.1{2.4 of Abate, Choudhury and Whitt (1996). however, the L2 theory can be applied by using the generalized Laguerre polynomials L( ) n (t) Pnk=0( 1)k n+ n k tk=k! with respect to the weight t e t for = 1. The expansion for etf(t) can then be re-expressed in terms of the standard Laguerre polynomials and coe cients, i.e., etf(t) = 1 Xn=1 q(1) n L(1) n (t) = 1 Xn=1 qnLn(t) ; (2.10) using relationships among the Laguerre coe cients and polynomials, i.e., q(1) n = qn qn+1, L(1) n = L0n+1(t) and L0n+1(t) = L0n(t) Ln(t), see p. 241 of Magnus et al. (1966). As a consequence, the direct Laguerre-series representations are still valid, with the understanding that the L2 properties require the generalized Laguerre polynomials. Even for class I and II pdf's, we nd it necessary to scale appropriately so that the rightmost singularity s of the Laplace transform f(s) satis es s < 1=2. Such scaling can easily 6 be done for classes I and II, but cannot be done for class III. Class III can be transformed to satisfy this condition too, e.g., by exponential damping as discussed in Abate, Choudhury and Whitt (1994) and Abate and Whitt (1996), but then the moments and Laguerre coe cients are altered, so we are unable to exploit the transformation. (This point is made in Section 7 of Keilson and Nunn (1979) too.) In contrast, the moments are simply scaled under our linear scaling of class I and II pdf's. Assuming that we are indeed considering the Laguerre-series representation for etf(t), then the Laguerre coe cient qn becomesqn = Z 1 0 f(t)Ln(t)dt : (2.11) Since Ln(t) is a polynomial of degree n, we see that qn is a linear combination of the rst n moments of f . Indeed, we can combine (2.4) and (2.11) to obtain the following result. Theorem 2.1. If the function etf(t) has a Laguerre-series representation with respect to the weight e t then all moments of f are nite and the nth Laguerre coe cient qn of etf(t) can be expressed as qn = n Xk=0 nk ( 1)kmk k! ; (2.12) where mk is the kth moment of f with m0 = 1. Since the rightmost singularity is less than 1=2, all moments are nite. Remark 2.2. In establishing a connection between moments and Laguerre coe cients, we follow Section 7 of Keilson and Nunn (1979). However, (2.12) is a more elementary relation than their in nite series. We obtain the simplicity by letting qn be the Laguerre coe cient of etf(t) with respect to the weight e t. Both approaches require that the rightmost singularity of the Laplace transform f̂(s) satisfy s < 1=2. Example 2.1. To better understand the condition requiring that the rightmost singularity of f̂(s) be less than 1=2, consider the case of an exponential pdf with mean 3. The Laplace transform is f̂(s) = (1 + 3s) 1 = 1 Xn=0 3n( s)n (2.13) the rightmost singularity is 1=3 and the series has radius of convergence 1/3. If we were to apply (2.12), we would getqn = n Xk=0( 3)k nk = (1 3)n = ( 2)n ; (2.14) 7 which is growing geometrically in n. Moreover, the Laguerre-series representation for f̂(s) would be f̂(s) = 1 1 + s 1 Xn=0( 2)n s 1 + s n ; (2.15) which is absolutely convergent for jsj < 1=3, but fails to converge for real s with s 1. With the Laguerre coe cients as in (2.11) and (2.12), we obtain f(t) = e t 1 X n=0 qnLn(t) (2.16) and f̂(s) = 1 Xn=0 qn sn (s + 1)n+1 : (2.17) It is noteworthy that the convergence of the transform series (2.17) tends to be faster than the convergence of the basic Laguerre series (2.16), because Ln(t) converges to 0 slowly as n!1. Indeed, Ln(t) et=2 ( 2nt)1=4 cos(2pnt =4) as n! 1 ; (2.18) see p. 245 of Magnus, Oberhettinger and Soni (1966). The Laguerre series in (2.17) tends to converge geometrically for all s with Re(s) > 0 because the Laguerre coe cients are bounded and js=(s+1)j < 1 for Re(s) > 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 <1 : (2.19) Using the Hilbert space theory of the function space L2[0;1), we know that, with g(t) = etf(t), kgk22 Z 1 0 etf(t)2dt = 1 Xn=0 q2 n <1 ; (2.20) which implies that qn ! 0. Hence the series (2.17) indeed converges geometrically fast. It is important to note, however, that when we do numerical inversion, we need to compute the transform f̂ (s) not for a single value of s but for several values of s. If we use the Fourierseries method in Abate and Whitt (1992, 1995), then we need to consider values of s = u+ iv with u = A=2t and v = k =t for A = 30, say, and k 1. A serious di culty occurs because js=(1 + s)j ! 1 as v ! 1 and, thus, as k ! 1. Thus, the Euler summation used in the Fourier-series method helps greatly, because it enables us to restrict attention to k 40, say. Then u = 15=t and v 125=t. We now give a bound on the error from truncation. 8 Theorem 2.2. Let M be the bound on jqnj in (2.19). For s = u + iv with u > 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.23) with = v=(1 + u) and = u=(1 + u). Proof. If s = u+ vi for real u and v, then z s s+ 1 = A+ Bi C +Bi (2.24) for A = u, B = v and C = u+ 1. If (2.24) holds with 0 jAj < C, then jzj2 = 4 + 2(1 + 2) + 2 4 + 2 2 + 1 < 1 (2.25) for = B=C and = jAj=C. Now let us consider again the Fourier-series method of numerical transform inversion with Euler summation. We have observed that we need to consider s = u + iv with u = 15=t and v 125=t. From (2.23) we see that s s + 1 1 u v2 when u << v ; (2.26) so that here s s + 1 1 t10 3 (2.27) in the worst case. Hence, to have truncation error jf̂N(s) f̂(s)j M , we need N log( t10 3) log(1 t10 3) log( t10 3) t10 3 : (2.28) For example, for = 10 15 and M = 1, we require N 4:4t 1 104. Thus, for t = 1, we require N = 44; 000. The resulting computation is feasible, but not easy. Hence, if alternative methods are available, then they might well be perferred. To directly apply (2.17) for computing transform values given a pdf f , we can compute the Laguerre coe cients qn associated with etf(t) for the given pdf f by numerically integrating 9 the integral in (2.11) and then approximate the sum in (2.17) by f̂N(s) in (2.22) to achieve desired error in (2.21). However, it is often possible to obtain the coe cients more easily from the moments of f using (2.12). To apply Theorem 2.1, we can rst scale the transform so that it has rightmost singularity 1. This is achieved by considering the scaled Laplace transform f̂s (s) f̂(s s). In terms of a random variable X with pdf f and moments mn, this corresponds to replacing X by X=s , which has pdf s f(s t), t 0, and moments mn=(s )n. For the exponential pdf with mean 3 in Example 2.1, the rightmost singularity is s = 1=3. The scaling converts the pdf to an exponential pdf with mean 1. However, for nonexponential pdf's, we usually do not have s = 1=m1. To carry out the scaling, we need to determine the rightmost singularity s . Fortunately, algorithms to compute s from the moments are contained in Abate, Choudhury, Lucantoni and Whitt (1995), where s is denoted . Under considerable generality, s = lim n!1 nmn 1 mn ; (2.29) but the proposed algorithms often do much better by exploiting Richardson extrapolation; see p. 987 of Abate, Choudhury, Lucantoni and Whitt (1995). It is convenient that we can relate the Laguerre coe cients of a pdf f(t) to its complementary cdf (ccdf) F c(t) 1 F (t) directly. Theorem 2.3. If etf(t) has a Laguerre-series representation for a pdf f(t), then so does etF c(t) for the associated ccdf F c(t), and the expansions are related by f(t) = e t 1 X n=0 qnLn(t) (2.30) and F c(t) = e t 1 Xn=0(qn qn+1)Ln(t) : (2.31) First Proof. Using Laplace transforms, F̂ c(s) = 1 f̂ (s) s = 1s + f̂(s) s+ 1 s f̂ (s) = 1 1 + s 1 X n=0 qn s s + 1 n 1 1 + s 1 X n=0 qn+1 s s + 1 n using the fact that q0 = Z 1 0 f(t)L0(t)dt = Z 1 0 f(t)dt = 1 : 10 Second Proof. Use one of the basic di erentiation formulas for Laguerre polynomials L0n(t) = nLn(t) nLn 1(t) t (2.32) or L0n(t) L0n+1(t) = Ln(t) (2.33) e.g., see p. 241 of Magnus, Oberhettinger and Soni (1966). We now give two examples illustrating how the Laguerre-series representation of the transform obtained from the moments can be helpful. Example 2.2. In Theorem 5.4 of Abate and Whitt (1998) we indicated how new pdf's can be constructed from others using their moment sequences. In particular, suppose that we have two pdf's f(t) and g(t) where their moments are related by mn(G) = mn(F ) 2n+ 1 ; n 1 : (2.34) Then their Laplace transform are related by ĝ(s) = 1 2ps Z s 0 f̂(z) pz dz : (2.35) If both the moments and the Laplace transform of f are known, it may be easier to compute ĝ(s) using the Laguerre-series representation based on Theorem 2.1 and (2.34) than to perform the integration in (2.35). Example 2.3. An interesting example from p. 86 of Abate and Whitt (1996) and p. 986 of Abate, Choudhury, Lucantoni and Whitt (1995) is the Cayley-Einstein-Polya (CEP) pdf with mean 1 and Laplace transform satisfying the functional equation f̂(s) = exp( sf̂ (s)) : (2.36) Evidently, the associated pdf is unknown, but the moments are mn = (n+ 1)n 1, n 1. The ccdf has the asymptotic formF c(t) r e5 2 t 3=2e t as t!1 (2.37) for = 1=e. Even though the mean is m1 = 1, the rightmost singularity of f̂(s) is s = e 1 > 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 < a < b 1, then etf(t) has a Laguerre-series representation with Laguerre coe cients qn = Z b a x 1 x n dG(x) : (3.2) Proof. First, because of the exponentials in (3.1), the rightmost singularity of f̂(s) satis es s < 1=2. We exploit the fact that the Laplace transform of Ln(t) is (s 1)n=sn+1. We change the order of integration to obtain qn = Z 1 0 e t(etf(t))Ln(t)dt = Z b a Z 1 0 xe xtLn(t)dt dG(x) = Z b a x 1 x n dG(x) : Corollary. If a pdf f can be represented as f(t) = Z d c y 1e t=yw(y)dy (3.3) 13 for 0 c < d < 2, where w is a pdf, then etf(t) has a Laguerre-series representation with Laguerre coe cients qn = Z c 1 d 1 x 1 x n x 2w(1=x)dx : (3.4) Proof. Make the change of variables x = 1=y to go from the mixture in (3.3) to f(t) = Z c 1 d 1 xe xtx 2w(1=x)dx : Then apply Theorem 3.1. We next characterize Laguerre coe cients of mixtures of exponentials where the mixing pdf has support [0; 1]. Recall that mixtures can be represented as products of independent random variables. Theorem 3.2. Let X and Y be independent random variables with X exponentially distributed having mean 1 and Y having support on [0; 1]. If f is the pdf of XY , then etf(t) has a Laguerreseries representation with the nth Laguerre coe cient qn equal to the nth moment of (1 Y ). Hence 1 = q0 qn qn+1 for all n and qn ! 0. Proof. The rightmost singularity of the Laplace transform is less than or equal to 1, so that the Laguerre-series representation exists. By Theorem 2.1, the Laguerre coe cient is qn = n Xk=0 nk ( 1)kmk k! = n Xk=0 nk ( 1)kyk where mk (yk) is the kth moment of XY (Y ). The nth moment of 1 Y clearly has the nal expression too. Example 3.1. Consider the pdff(t) = Z 1 0 y 1e t=yw(y)dy ; (3.5) where w(y) = log(1 y), 0 y 1. The nth moment of the pdf w(1 y) = log y has nth moment (n+1) 2 by 4.1.50 of Abramowitz and Stegan (1972), so that the Laguerre coe cient of etf(t) is qn = (n+ 1) 2 by Theorem 3.2. 14 Corollary. If, in addition to the assumptions of Theorem 3.2, the mixing pdf w(y) of Y is symmetric on [0; 1], i.e., if w(y) = w(1 y), 0 y 1, then the nth Laguerre coe cient qn of etf(t) equals the nth moment of Y . Example 3.2. Consider the pdf (3.5), where the mixing pdf w(y) is uniform on [0; 1]. By the Corollary to Theorem 3.2, qn = (n+ 1) 1. We now obtain even more explicit representations for the Laguerre coe cients of etf(t) for the class of BME pdf's. A BME pdf can be expressed as v(p; q; t) = Z 1 0 y 1e t=yb(p; q; y)dy; t 0 ; (3.6) where b(p; q; y) is the standard beta pdf, i.e., b(p; q; y) = (p+ q) (p) (q)yp 1(1 y)q 1; 0 y 1 ; (3.7) and (x) is the gamma function. A BME pdf has three parameters p > 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 < < 1. In general it is necessary to check that the Erlang-series representation for f̂(s) in (4.5) converges. Example 4.1. Consider the pdf f(t) = (log 2) 1e t=(2 e t); t 0 : (4.7) Using the Taylor series expansion of (2 e t) 1, we would obtain the associated Erlang-series representation f(t) = (log 2) 1 1 Xk=1( 1)k 1~b(k 1)e(k; t) ; (4.8) where ~b(k) are the ordered Bell numbers; see 5.27 on p. 175 of Wilf (1994). The ordered Bell number ~b(n) is the number of ordered partitions of a set of n elements, where the order of the classes is counted but not the order of the elements in the classes. The associated term-by-term Erlang-series representation for the Laplace transform would be f̂(s) = (log 2) 1 1 Xk=1( 1)k 1~b(k 1)(1 + s) n : (4.9) However, neither (4.8) nor (4.9) is convergent. Thus, we see that the Erlang-series representation for the Laplace transform is not e ective in this example. Example 4.2. As a second more favorable example, consider the pdf f(t) = 2e t(1 cos t); t 0 ; (4.10) as on p. 990 of Abate, Choudhury, Lucantoni and Whitt (1995). Using the Taylor series representation of 1 cos t, we obtain f(t) = 2 1 Xk=1( 1)k+1e(2k + 1; t); t 0 ; (4.11) 18 and f̂(s) = 2 1 Xk=1( 1)k+1(1 + s) (2k+1) : (4.12) Unlike in Example 4.1, the Erlang-series representation in (4.12) is (conditionally) convergent for all s with Re(s) 0. 5. Mixtures of Exponentials Again In this section we obtain results about Erlang-series representations when the pdf is a mixture of exponential pdf's. We start by showing that every mixture of exponential pdf's (spectral representation) where the mixing pdf has support [0; b] for b < 1 has an Erlangseries representation and identify the coe cients. First, we note that it su ces to consider mixing pdf's with support [0; 1], because we can always rescale; i.e., if Z = Y X where Y has support on [0; b], we consider (Y=b) with support [0; 1] and let Z = (Y=b)X . We then consider Z = bZ = Y X . Theorem 5.1. Consider a completely monotone pdf f(t) that has a Laguerre-series representation as in Theorem 3.1, i.e., f(t) = Z b a y 1e t=ydH(y) = e t 1 Xn=0 qnLn(t) ; (5.1) where H is a cdf with support in [a; b] with 0 a < b < 2. If Gc(t) = Z b a e txdH(x); t 0 ; (5.2) Then Gc(t) has the Erlang-series representation Gc(t) = 1 Xk=1 qk 1e(k; t); t 0 ; (5.3) and Gc(t) has a pdf g(t) with g(t) = 1 Xk=1(qk 1 qk)e(k; t) : (5.4) To prove Theorem 5.1, we use the following lemma. Lemma 5.1. If (5.2) and the rst relation (5.1) hold, where 0 a < b 1, then Ĝc(s) = s 1f̂ (s 1). 19 Proof. Note that f̂(s) = Z b a (1 + sy) 1dH(y) and Gc(s) = Z b a (x+ s) 1dH(y) = 1s Z b a (1 + xs 1) 1dH(x) : Proof of Theorem 5.1. From Lemma 5.1, Ĝc(s) = s 1f̂(s 1) = Z 1 0 (s 1e x=s)f(x)dx ; so that Gc(t) = Z 1 0 L 1(s 1e x=s)f(x)dx = Z 1 0 J0(2pxt)f(x)dx (5.5) where J0(2pxt) = e t 1 X n=0 tn n!Ln(t) ; (5.6) see 29.3.75 and 22.9.16 of Abramowitz and Stegun (1972). Therefore, Gc(t) = e t 1 Xk=0 qktk k! = 1 Xk=1 qk 1e(k; t) for qk in (2.11). Remark 5.1. We use the condition b < 2 to apply Theorem 3.1 to have the Laguerre-series representation in (5.1). As remarked before Theorem 5.1, we can rescale to have b < 2 or b = 1 if b <1. By Theorem 3.2, for the special case in which b 1, we have 1 = q0 q1 qn qn+1 for all n. Second Proof of Theorem 5.1 for b = 1. Working with the density g(t) of Gc(t) in (5.2), we have g(t) = e t Z 1 0 xe (1 x)th(x)dx = e t Z 1 0 1 Xn=0 x(1 x)ntn n! h(x)dx = 1 Xn=1 e ttn 1 (n 1)! Z 1 0 x(1 x)nh(x)dx = 1 Xn=1 e(n; t) Z 1 0 (1 x)xnh(1 x)dx : (5.7) 20 On the other hand, the Laplace transform f̂(s) of f in (5.1) can be expressed as f̂(s) = Z 1 0 e st Z 1 0 y 1e t=yh(y)dy = Z 1 0 1 1 + syh(y)dy = Z 1 0 1 1 + s(1 z)h(1 z)dz = 1 1+ s Z 1 0 1 1 + z s 1+s h(1 z)dz = 1 X n=0 1 1 + s s 1 + s n Z 1 0 znh(1 z)dz : Hence, f(t) = 1 Xn=0 qne tLn(t) for qn = Z 1 0 (1 z)nh(z)dz : (5.8) The result follows by comparing the coe cients. Remark 5.2. We need the condition that the pdf f(t) be completely monotone in order for Gc(t) to be a bona de ccdf. To see this, suppose that f(t) = 1 Xn=0 f (n)(0)tn n! (5.9) and that its transform can be computed term by term, i.e., f̂ (s) = 1 Xn=0 f (n)(0)s (n+1) ; (5.10) see Section 30 of Doetsch in (1974) for conditions. Also suppose that Ĝc(s) is analytic at s = 0, so that Gc(s) = 1 Xn=0 mn+1 (n+ 1)!( s)n : (5.11) The representations (5.10) and (5.11) imply that ( 1)nf (n)(0) = mn+1(G) (n+ 1)! > 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 < r. Theorem 5.5. For q = 1, the B3ME pdf is the PME pdf with p = r, i.e., v3(p; 1; t) = p (p+ 1; t) tp+1 = Z 1 1 y 1e t=yfp(y)dy ; t 0 : (5.33) where (p; t) R t 0 p 1e d is the incomplete gamma function, and its moments aremn(p; 1) = p(n!)=(p n). Proof. To prove the rst relation in (5.33), compare the power series for (p+ 1; t) with v3(p; 1; t) = p 1 Xn=1 (n 1)! (p+ 1)n e(n; t) = pe t 1 Xn=1 tn 1 (p+ 1)n ; (5.34) see 6.5.4 and 6.5.29 of Abramowitz and Stegun (1972). For the second relation, apply 6.5.2 of Abramowitz and Stegun (1972) to get v3(p; 1; t) = p tp+1 Z t 0 pe d = p Z 1 0 xpe txdx = p Z 1 1 e t=yy (r+2)dy : (5.35) We now give the full asymptotic form as t ! 1. It is signi cant that the second term is exponentially damped. Corollary. For all p > 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 < ak+1 and ak ! 1 as k !1, so that f̂(s) has poles at 1=ak for all k, implying that f̂(s) has a singularity at 0, so that it is a long-tail pdf. If we truncate the in nite series in (6.1) or (6.2), we obtain a nite mixture of exponential distributions, i.e., a hyperexponential distribution. We consider how to truncate in the next section with the objective of capturing the tail asymptotics at times of interest. Thus this section together with the next constitutes an alternative approach to Feldmann and Whitt (1997). It is signi cant that the nth moment of f in (6.1) can be simply expressed as mn = n! 1 Xk=1 pkank : (6.3) By choosing appropriate sequences fpkg and fakg, we can produce desired moment sequences. In many cases we can describe the asymptotic behavior of the pdf and/or its associated ccdf distribution function F c(t) 1 F (t) from the asymptotic behavior of the moments, e.g., by invoking Abate, Choudhury, Lucantoni and Whitt (1995). 26 Example 6.1. Let pk = e 1=k! and ak = k. Then mn n! = e 1 1 Xk=1 kn k! = b(n) ; (6.4) where b(n) = f1; 2; 5; 15; 52; 203; : : :g are the Bell numbers; see p. 20 of Wilf (1994). The Bell number b(n) is the number of ways a set of n elements can be partitioned. The Bell numbers themselves have the relatively simple exponential generating function B(z) 1 Xn=0 b(n) n! zn = eez 1 ; (6.5) from which we can deduce the recurrence relation b(n) = n 1 Xk=0 n 1 k b(k) (6.6) with b(0) = 1; see p. 23 of Wilf (1994). The discrete distribution with mass pk = e 1=k! at k has Laplace transform b̂(s) = B( s); see (3.9) on p. 86 of Abate and Whitt (1996). We can easily generalize to obtain a two-parameter family; e.g., keep pk = e 1=k! but now let ak = c(k 1 + a). The parameter c is a scale parameter; m1 = 1 if c = 1=(1 + a). In general, for this \generalized Bell pdf", the nth moment is mn = n!cn e 1 Xk=1 (k + 1 a)n k! : (6.7) Example 6.2. Let pk = 2 (k+1) and ak = k, k 0. Then the nth moment is mn = n! 1 Xk=0 kn 2k+1 = n!~b(n) ; (6.8) where f~b(n)g is the sequence of ordered Bell numbers as in Example 4.1. 7. Truncating In nite Series The Erlang-series and exponential-series representations in Sections 4{6 produce in nite series of pdf's, which in applications we will want to truncate. In this section we consider how to truncate. Appropriate truncations depend on the asymptotic form of the pdf or ccdf and the times of interest. We now show how the asymptotic form and the appropriate truncation point (as a function of the time of interest) can be identi ed by applying the Euler-Maclaurin summation formula to approximate the sum by an integral and then Laplace's method to determine the asymptotic behavior of the integral; see pp. 80 and 285 of Olver (1974). 27 We start with the ccdf F c(t) = 1 X n=1 pnF c n(t); t 0 ; (7.1) where F c n(t) is a ccdf for each n, which we rewrite as F c(t) = 1 X n=1 e (n;t) ; (7.2) where (n; t) = log pn logF c n(t) : (7.3) We now regard (x; t) as a continuous function of x and assume that we can approximate the sum by an integral (invoking the Euler-Maclaurin formula to quantify the error), obtaining F c(t) Z 1 0 e (x;t)dx : (7.4) We now use Laplace's method on the integral in (7.4). We assume that there is a unique x (t) that minimizes (x; t), i.e., for which 0(x (t); t) = 0 and 00(x (t); t) > 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 < (7.24) 00(x (t); t) = t 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 < 1 and otherwise in nite. Example 7.3. An exponential series representation for a cdf associated with the M/M/1 queue was recently derived and applied to determine the asymptotic form by Guillemin and Pinchon (1998). They focused on the total time spent in the system by all customers in a busy period (the area under the queue length process) in the M/M/1 queue. Daley and Jacobs (1969) had shown that the transform can be expressed as the ratio of two Bessel functions, i.e., ĝ(s) = 1 p J +1(2p =s) J (2p =s) for = (1 + )=s ; (7.30) where is the tra c intensity. Guillemin and Pinchon showed that the distribution has an exponential series representation, where the exponentials satisfy ak k=(1 + ) as k ! 1. Moreover, they applied the Laplace method directly to the discrete sums to derive the asymptotic form of the ccdf, Gc(t) 1 p2 b 4 t 1=4 exp( p4 t=b) as t!1 ; (7.31) 30 where b (1 + ) 1 and log 2(1 )=(1+ ). We now develop a relatively simple approximation for this ccdf. We propose the threeparameter exponential-series ccdf Gca(t) = A 1 Xk=1 rk k exp( t=bk) ; (7.32) which would be proper if A = 1= log(1 r). Using the results of Example 7.1, we have Gca(t) p2 A b 4!t 1=4 exp(1 p4!t=b) as t!1 (7.33) where ! log r. Hence, we can make the asymptotic values agree by letting A = (1 )=2 and w = , which gives r = exp(2(1 )=(1 + )). Moreover, if the total probability mass is less than 1, we can make the approximating ccdf proper by adding an atom at 0. Numerical results for the case = 0:6 are displayed in Table 2. The exact values are taken from Table 3 of Guillemin and Pinchon (1998). From Table 2, we see that the complicated exact ccdf can be reasonably approximated by the relatively simple exponential-series ccdf Gca(t) in (7.32). time t exact approximation asymptote 0 1:000 1:000 1 0+ 1:000 0:481 1 1 0:482 0:313 0:398 10 0:117 0:113 0:127 20 0:656 e 1 0:688 e 1 0:755 e 1 40 0:324 e 1 0:362 e 1 0:390 e 1 60 0:199 e 1 0:229 e 1 0:242 e 1 80 0:134 e 1 0:157 e 1 0:165 e 1 200 0:274 e 2 0:322 e 2 0:333 e 2 300 0:109 e 2 0:127 e 2 0:130 e 2 400 0:507 e 3 0:586 e 3 0:600 e 3 2000 0:555 e 7 0:598 e 7 0:598 e 7 4000 0:365 e 9 0:384 e 9 0:384 e 9 Table 2: A comparison of the approximate exponential-series ccdf Gca(t) in (7.32) with the exact ccdf and the asymptotic in (7.31) for Example 7.3. 8. Conclusions We have shown how pdf's and their transforms can be constructed via in nite-series representations, in particular, as Laguerre series (Sections 2 and 3), Erlang series (Sections 4 and 31 5) and exponential series (Section 6). These series representations make it possible to computethe Laplace transform values in applications. At the same time, other properties of the pdf'scan be deduced, such as moments and the asymptotic form, so that the pdf can be chosento satisfy desired properties in applications. For example, exponential series representationsalways yield completely monotone pdf's whose moments can be computed. We showed how todetermine the asymptotic form in Section 7.By imposing additional structure, we obtain relatively simple parametric families of pdf'ssuch as the BME pdf's in (3.6), the B3ME pdf's in (5.23), the PME pdf's in (5.33), thegeneralized Bell pdf's in Example 6.1, and the geometric mixture of exponential pdf's inExample 7.1. Since the Laplace transform values for these pdf's can be computed, these pdf'sare natural candidates to use in applications. Then the parameters can be chosen to match themoments, the asymptotic form and other features of interest. Toward that end, we studied theBME pdf's further in Abate and Whitt (1998). The connection to the con uent hypergeometricfunction M(a; b; z) established in Theorem 5.2 makes it possible to establish analogous resultsfor B3ME pdf's.Over all, our analysis help extend the analyst's toolkit.32 ReferencesAbate, J., Choudhury, G. L., Lucantoni, D. M. and Whitt, W. (1995). Asymptotic analysis oftail 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 fur Elcktronik undUbertragungstechnik 47, 311{321.Abate, J., Choudhury, G. L. and Whitt, W. (1994). Waiting-time tail probabilities in queueswith long-tail service-time distributions. Queueing Systems 16, 311{338.Abate, J., Choudhury, G. L. and Whitt, W. (1996). On the Laguerre method for numericallyinverting Laplace transforms. INFORMS J. Computing 8, 413{427.Abate, J., Choudhury, G. L. and Whitt, W. (1997). Numerical inversion of multidimensionalLaplace 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 ComputationalProbability, W. Grassman (ed.), Kluwer, Boston, to appear.Abate, J. and Whitt, W. (1992). The Fourier-series method for inverting transforms ofprobability distributions. Queueing Systems 10, 5{88.Abate, J. and Whitt, W. (1995). Numerical inversion of Laplace transforms of probabilitydistributions. ORSA J. Computing 7, 36{43.Abate, J. and Whitt, W. (1996). An operational calculus for probability distributions viaLaplace transforms. Adv. Appl. Prob. 28, 75{113.Abate, J. and Whitt, W. (1997). Asymptotics for M/G/1 low-priority waiting-time tailprobabilities. Queueing Systems 25, 173{233.Abate, J. and Whitt, W. (1998). Modeling service times with densities having non-exponentialtails, submitted.Abate, J. andWhitt, W. (1998a). Computing transforms for numerical inversion via continuedfractions. AT&T; Labs, Florham Park, NJ.33 Abramowitz, M. and Stegun, I. A. (1972). Handbook of Mathematical Functions, NationalBureau of Standards, Washington, D. C.Asmussen, S., Nerman, O. and Olsson, M. (1996). Fitting phase type distributions via theEM algorithm. Scand. J. Statist. 23, 419{441.Choudhury, G. L. Lucantoni, D. M. and Whitt, W. (1994). Multidimensional transforminversion 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 pollingmodels 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 astable 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 anM/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, NewYork.Keilson, J. and Nunn, W. R. (1979). Laguerre transformation as a tool for the numericalsolution of integral equations of convolution type. Applied Math. Computation 5, 313{359.Keilson, J., Nunn, W. and Sumita, U. (1981). The bilateral Laguerre transform. AppliedMath. Computation 8, 137{174.34 Keilson, J. and Sumita, U. (1982). Waiting time distribution response to tra c surges viathe Laguerre transform. In Applied Probability | Computer Science: The Interface, R.L. Disney and T. J. Ott (eds.), Birkhauser, 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 theLaguerre transform. Queueing Systems 4, 367{382.Magnus, W., Oberhettinger, F. and Soni, R. P. (1966). Formulas and Theorems for theSpecial 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 nonnegativeLaguerre 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.Szego, 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 Models12, 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 AidsComput. 10, 91{96.Wynn, P. (1966). On the convergence and stability of the epsilon algorithm. SIAM J. Numer.Anal. 3, 91{122.35