Joseph Abate

The Fourier-series method for inverting transforms of probability distributions

This paper reviews the Fourier-series method for calculating cumulative distribution functions (cdfs) and probability mass functions (pmfs) by numerically inverting characteristic functions, Laplace transforms and generating functions. Some variants of the Fourier-series method are remarkably easy to use, requiring programs of less than fifty lines. The Fourier-series method can be interpreted as numerically integrating a standard inversion integral by means of the trapezoidal rule. The same formula is obtained by using the Fourier series of an associated periodic function constructed by aliasing; this explains the name of the method. This Fourier analysis applies to the inversion problem because the Fourier coefficients are just values of the transform. The mathematical centerpiece of the Fourier-series method is the Poisson summation formula, which identifies the discretization error associated with the trapezoidal rule and thus helps bound it. The greatest difficulty is approximately calculating the infinite series obtained from the inversion integral. Within this framework, lattice cdfs can be calculated from generating functions by finite sums without truncation. For other cdfs, an appropriate truncation of the infinite series can be determined from the transform based on estimates or bounds. For Laplace transforms, the numerical integration can be made to produce a nearly alternating series, so that the convergence can be accelerated by techniques such as Euler summation. Alternatively, the cdf can be perturbed slightly by convolution smoothing or windowing to produce a truncation error bound independent of the original cdf. Although error bounds can be determined, an effective approach is to use two different methods without elaborate error analysis. For this purpose, we also describe two methods for inverting Laplace transforms based on the Post-Widder inversion formula. The overall procedure is illustrated by several queueing examples.

Numerical Inversion of Laplace Transforms of Probability Distributions

We present a simple algorithm for numerically inverting Laplace transforms. The algorithm is designed especially for probability cumulative distribution functions, but it applies to other functions as well. Since it does not seem possible to provide effective methods with simple general error bounds, we simultaneously use two different methods to confirm the accuracy. Both methods are variants of the Fourier-series method. The first, building on Dubner and Abate (Dubner, H., J. Abate. 1968. Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform. JACM 15 115–123.) and Simon, Stroot, and Weiss (Simon, R. M., M. T. Stroot, G. H. Weiss. 1972. Numerical inversion of Laplace transforms with application to percentage labeled experiments. Comput. Biomed. Res. 6 596–607.), uses the Bromwich integral, the Poisson summation formula and Euler summation; the second, building on Jagerman (Jagerman, D. L. 1978. An inversion technique for the Laplace transform with applications....

Numerical Inversion of Laplace Transforms by Relating Them to the Finite Fourier Cosine Transform

In this paper the problem of readily determining the inverse Laplace transform numerically by a method which meets the efficiency requirements of automatic digital computation is discussed. Because the result inverse function is given as a Fourier cosine series, the procedure requires only about ten FORTRAN statements. Furthermore, it does not require the use of involved algorithms for the generation of any special functions, but uses only cosines and exponentials. The basis of the method hinges on the fact that in evaluating the inverse Laplace transform integral there exists a freedom in choosing the contour of integration. Given certain restrictions, the contour may be any vertical line in the right-half plane. Specifying a line, the integral to be evaluated is essentially a Fourier integral. However, the method is concerned with determining the proper line, so that when the integral (along this line) is approximated, the error is as small as desired by virtue of having chosen the particular contour.

Numerical inversion of probability generating functions

Random quantities of interest in operations research models can often be determined conveniently in the form of transforms. Hence, numerical transform inversion can be an effective way to obtain desired numerical values of cumulative distribution functions, probability density functions and probability mass functions. However, numerical transform inversion has not been widely used. This lack of use seems to be due, at least in part, to good simple numerical inversion algorithms not being well known. To help remedy this situation, in this paper we present a version of the Fourier-series method for numerically inverting probability generating functions. We obtain a simple algorithm with a convenient error bound from the discrete Poision summation formula. The same general approach applies to other transforms.

A Unified Framework for Numerically Inverting Laplace Transforms

We introduce and investigate a framework for constructing algorithms to invert Laplace transforms numerically. Given a Laplace transform \hat{f} of a complex-valued function of a nonnegative real-variable, f, the function f is approximated by a finite linear combination of the transform values; i.e., we use the inversion formula f(t) \approx f_n (t) \equiv \frac{1}{t} \sum_{k = 0}^{n}\omega_{k}\hat{f}\biggl(\frac{\alpha_{k}}{t}\biggr),\quad 0 where the weights ωk and nodes αk are complex numbers, which depend on n, but do not depend on the transform \hat{f} or the time argument t. Many different algorithms can be put into this framework, because it remains to specify the weights and nodes. We examine three one-dimensional inversion routines in this framework: the Gaver-Stehfest algorithm, a version of the Fourier-series method with Euler summation, and a version of the Talbot algorithm, which is based on deforming the contour in the Bromwich inversion integral. We show that these three building blocks can be combined to produce different algorithms for numerically inverting two-dimensional Laplace transforms, again all depending on the single parameter n. We show that it can be advantageous to use different one-dimensional algorithms in the inner and outer loops.

An Introduction to Numerical Transform Inversion and Its Application to Probability Models

Numerical transform inversion has an odd place in computational probability. Historically, transforms were exploited extensively for solving queueing and related probability models, but only rarely was numerical inversion attempted. The model descriptions were usually left in the form of transforms. Vivid examples are the queueing books by Takacs [Takacs, 1962] and Cohen [Cohen, 1982]. When possible, probability distributions were calculated analytically by inverting transforms, e.g., by using tables of transform pairs. Also, moments of probability distributions were computed analytically by differentiating the transforms and, occasionally, approximations were developed by applying asymptotic methods to transforms, but only rarely did anyone try to compute probability distributions by numerically inverting the available transforms. However, there were exceptions, such as the early paper by Gaver [Gaver, 1966]. (For more on the history of numerical transform inversion, see our earlier survey [Abate and Whitt, 1992a].) Hence, in the application of probability models to engineering, transforms became regarded more as mathematical toys than practical tools. Indeed, the conventional wisdom was that numerical transform inversion was very difficult. Even numerical analysts were often doubtful of the numerical stability of inversion algorithms. In queueing, both theorists and practitioners lamented about the “Laplace curtain” obscuring our understanding of system behavior.

Exponential Approximations for Tail Probabilities in Queues, I: Waiting Times

This paper focuses on simple exponential approximations for tail probabilities of the steady-state waiting time in infinite-capacity multiserver queues based on small-tail asymptotics. For the GI/GI/s model, we develop a heavy-traffic asymptotic expansion in powers of one minus the traffic intensity for the waiting-time asymptotic decay rate. We propose a two-term approximation for the asymptotic decay rate based on the first three moments of the interarrival-time and service-time distributions. We also suggest approximating the asymptotic constant by the product of the mean and the asymptotic decay rate. We evaluate the exponential approximations based on the exact asymptotic parameters and their approximations by making comparisons with exact results obtained numerically for the BMAP/GI/1 queue, which has a batch Markovian arrival process, and the GI/GI/s queue. Numerical examples show that the exponential approximations are remarkably accurate, especially for higher percentiles, such as the 90th percentile and beyond.

Asymptotics for M/G/1 low-priority waiting-time tail probabilities

We consider the classical M/G/1 queue with two priority classes and the nonpreemptive and preemptive-resume disciplines. We show that the low-priority steady-state waiting-time can be expressed as a geometric random sum of i.i.d. random variables, just like the M/G/1 FIFO waiting-time distribution. We exploit this structures to determine the asymptotic behavior of the tail probabilities. Unlike the FIFO case, there is routinely a region of the parameters such that the tail probabilities have non-exponential asymptotics. This phenomenon even occurs when both service-time distributions are exponential. When non-exponential asymptotics holds, the asymptotic form tends to be determined by the non-exponential asymptotics for the high-priority busy-period distribution. We obtain asymptotic expansions for the low-priority waiting-time distribution by obtaining an asymptotic expansion for the busy-period transform from Kendalls functional equation. We identify the boundary between the exponential and non-exponential asymptotic regions. For the special cases of an exponential high-priority service-time distribution and of common general service-time distributions, we obtain convenient explicit forms for the low-priority waiting-time transform. We also establish asymptotic results for cases with long-tail service-time distributions. As with FIFO, the exponential asymptotics tend to provide excellent approximations, while the non-exponential asymptotics do not, but the asymptotic relations indicate the general form. In all cases, exact results can be obtained by numerically inverting the waiting-time transform.

Waiting-time tail probabilities in queues with long-tail service-time distributions

We consider the standardGI/G/1 queue with unlimited waiting room and the first-in first-out service discipline. We investigate the steady-state waiting-time tail probabilitiesP(W>x) when the service-time distribution has a long-tail distribution, i.e., when the service-time distribution fails to have a finite moment generating function. We have developed algorithms for computing the waiting-time distribution by Laplace transform inversion when the Laplace transforms of the interarrival-time and service-time distributions are known. One algorithm, exploiting Pollaczeks classical contourintegral representation of the Laplace transform, does not require that either of these transforms be rational. To facilitate such calculations, we introduce a convenient two-parameter family of long-tail distributions on the positive half line with explicit Laplace transforms. This family is a Pareto mixture of exponential (PME) distributions. These PME distributions have monotone densities and Pareto-like tails, i.e., are of orderx−r forr>1. We use this family of long-tail distributions to investigate the quality of approximations based on asymptotics forP(W>x) asx→∞. We show that the asymptotic approximations with these long-tail service-time distributions can be remarkably inaccurate for typicalx values of interest. We also derive multi-term asymptotic expansions for the waiting-time tail probabilities in theM/G/1 queue. Even three terms of this expansion can be remarkably inaccurate for typicalx values of interest. Thus, we evidently must rely on numerical algorithms for determining the waiting-time tail probabilities in this case. When working with service-time data, we suggest using empirical Laplace transforms.

ASYMPTOTICS FOR STEADY-STATE TAIL PROBABILITIES IN STRUCTURED MARKOV QUEUEING MODELS

We apply Tauberian theorems with known transforms to establish asymptotics for the basic steady-state distributions in the BMAP/G/l queue. The batch Markovian arrival process (BMAP)is equivalent to the versatile Markovian point process or Neuts (N) process; it generalizes the Markovian arrival process (MAP) by allowing batch arrivals. We consider the waiting time, the workload (virtual waiting time) and the queue lengths at an arbitrary time, just before an arrival and just after a departure. We begin by establishing asymptotics for steady-state distributions of M/G/1-type Markov chains. Then we treat steady-state distributions in the BMAP/G/l and MAP/MSP/l queues. The MSPis a MAPindependent of the arrival process generating service completions during