-The moment generating function completely characterizes the random variable distribution. At the same time, the natural desire of the researcher is to obtain the density and the distribution function as more informative characteristics. The article analyzes popular transition methods from the moment generating function to the probability density, in particular, the saddlepoint method, the method based on the Heaviside theorem and Pade´ approximation, a significant number of the inverse Laplace transform numerical methods. Computational complexity and practical feasibility of the methods are investigated for a number of practical problems of the inverse Laplace transform. Recommendations on the optimal method choice are given and the directions for the further development of the topic are suggested as well. Index Terms-probability density function, moment generating function, numerical Laplace inversion methods, Pade´ approximation
I. INTRODUCTION
Moment generating function (MGF) is an important statistical characteristic that uniquely describes the random variable distribution. Since the probability density function (PDF) is related to the MGF through the inverse Laplace transform [1]: f (t) = L 1hM ( s)i; the computational complexity can be significantly reduced by using the MGF instead of PDF. For example, MGF of the two random variables sum is calculated as the product of the MGF of these variables [1].
Flowgraph models operate MGF to obtain the stochastic process total time distribution. They can be applied in different areas. For example, in reliability theory for the uptime estimation, in medical statistics for the disease progression and the treatment effectiveness analysis [2], [3], for the power system survivability analysis and it’s weak points identification [4].
Pade´ approximation, Heaviside’s theorem, Saddlepoint approximation and numerical inverse Laplace transform methods are used to obtain the PDF as the MGF inversion (while the MGF is acquired as the result of the flowgraph specific analysis). The transform method must provide an acceptable approximation error with a relatively low computational complexity and must generate non-negative PDF values.
The paper is organised as follows. Section II introduces the MGF transform via Pade´ approximation, section III presents the Saddlepoint approximation, section IV provides the review of eight numerical inverse Laplace transform methods. In section V we select methods for the software implementation and comparison on the basis of the preliminary analysis. Finally section VI is devoted to the results discussing.
The Heaviside theorem is suggested to obtain the PDF applying the inverse Laplace transform of the MGF [5], [6]: L 1h U (s) i = Xq U ( k)
R(s) k=1 R0( k) e k t; where k are the roots of the polynomial R(s).
MGF Pade´ approximation can be straightforward provided as series:
Pade´ approximation is applicable to MGF because: 1 X ci xi = i=0
Pp
j=0 aj Pqk=0 bk xj xk 1 M (s) = X n=0 n n! sn; where n is n-th moment.
Ren [5] compares three methods: Pade´ approximation, maximum entropy method and saddlepoint approximation. It is noted that the Pade´ approximation is more informative and always provides an analytical approximation for the unknown original PDF.
But this approach has a significant disadvantage: Pade´ approximation requires the high order derivatives computation (because the n-th moment is the n-th derivative of the MGF). This fact leads to great computational complexity often exceeding the hardware possibilities for the high order derivatives calculation in the symbolic form.
f (t) = 2 n where K(s) = log(M (s)), K0(s) = t, n is the number of random variables. For one random variable PDF approximation is equal to [7]:
There are over a hundred different numerical inverse Laplace transform (ILT) methods. Being applied to the various types of functions they differ in computational complexity and approximation accuracy.
Some ILT methods have Gibbs oscillations – an oscillatory error which occurs when a discontinuous function is approximated with Fourier series. Gibbs oscillations do not disappear when more terms are added to the approximation (the ”size” of the excess does not decrease vertically, but only horizontally [9]).
Different ILT numerical methods can be combined within
the same mathematical framework. Abate and Whitt propose
the unified formula (
Abate–Whitt framework [11]:
(
N X k F ( k );
t t k=1 where t > 0, nodes k and weights k (internal and external scaling constants) are real or complex numbers that depend on N , but do not depend on the transform F or the time argument t [12].
Abate–Whitt framework has low computational complexity, but it is characterized by a deterioration of the approximation on periodic functions at large values of t [9]. Since PDF is not a periodic function, this disadvantage can be ignored.
Various ILT methods in this framework differ only by approaches to k and k calculation.
Several methods belonging to the Abate-Whitt framework are considered below.
1) CME: Method is based on the concentrated matrixexponential distribution [11].
Matrix-exponential distributions of order N contain positive random variables with PDF f (t) = AeAt1, at t 0, where is a real row vector of length N , A is an NxN real matrix, and 1 is a column vector of ones of size N [11].
Let A be diagonalizable with spectral decomposition A = PkN=1 uk k k, where k are the eigenvalues, uk are the right eigenvectors and k are the left eigenvectors of A with kuk = 1, then the PDF of the matrix-exponential distribution can be written as:
N X cke kt; k=1 where ck = Aukvk1.
A matrix-exponential distribution is called concentrated if its squared coefficient of variation (SCV) is minimal [13]. The density of the matrix-exponential distribution with the minimal SCV is known analytically only for the order N < 3. Therefore, the matrix-exponential distributions required for the approximation were calculated using numerical optimization by the expression: f (t) = ce
N2 1 t Y cos2(!t j=0
N j ) = X ke k=1 kt; where c; ; !; j are positive real values.
Thus, the values k and k providing the minimum SCV for the matrix-exponential distribution are the coefficients of the Abate-Whitt framework and are obtained from numerical optimization [14].
It is guaranteed that the CME method is free from Gibbs oscillations and preserves monotonicity [11]. Many of the known methods do not have these properties: for example, the fn(t) function for the Euler and Gaver–Stehfest methods can take negative values.
Advantages of the method: it does not cause overshoot, maintains monotonicity, is accurate when using floating point arithmetic, the quality of the approximation improves with increasing order [11].
2) Gaver–Stehfest method: Method is based on the threeparameter exponential PDF properties (Gaver method [15]). Stehfest [16] improves method by taking the weighted average of a Gaver approximations sequence for fixed t.
Gaver–Stehfest algorithm does not use complex numbers; weights and nodes are real numbers and are calculated as follows [10]:
k = k ln(
X j=b k+21 c ( N 2 j)!j!(j
N j 2 2j! 1)!(k j)!(2j k)! ; where N is the number of terms used in the equation (the method is defined only for even N ). Wang in [10] notices that N is recommended to be taken in the range from 10 to 14.
The main disadvantage of the method is Gibbs oscillations [11].
Abate and Whitt [12] show that Euler and Talbot methods are more efficient than show that Euler and Talbot methods are more efficient than those of Gaver–Stehfest. It is also noted that the algorithm implementation requires high computational accuracy due to manipulations with large numbers [12], [17]. However, the Gaver–Stehfest algorithm has the advantage of using only real numbers.
3) Euler method: is an implementation of the Fourier series method using Euler summation to accelerate convergence [18]. Method is applicable only to odd orders and assumes nodes with positive imaginary parts.
Nodes and weights are calculated as follows [11]: k = (N 1) ln(10) 6 + i(k
1) where 1 = 12 ; k = 1, 2 < k < n+21 ; n = n1 1 ; n k =2 2n k+1 + 2 n2 1 nk2 1 , 1 < k < n2 1 Method suffers from Gibbs oscillations [11].
4) Talbot method: Method is based on deforming the contour integral in the Bromwich inversion [19]. It assumes nodes with positive imaginary parts and also applies values of k with a negative real part, which is usually outside the region of convergence, but there can be an analytic continuation of F (s) [9].
Nodes and weights are calculated as follows [11]: k = 2(k
1) 5 (k
1) N
+ i 1 = 2N 5 cot Method is not free from Gibbs oscillations [11].
5) Zakian method: Method is based on Fourier series method with Pade´ approximation [11].
General formula is [10]:
k=1
f (t) = 1t XN RefKkF ( tk )g;
(
Zakian proposes two methods for choosing coefficients Kk and k. In the first he compares the Laplace transform (t; u) (a rational function) with the Laplace transform of n(u t) (an exponential function) and chooses the coefficients so that the rational functions are equal to the classical Pade´ approximations of the exponential function [21]. Another method includes least squares optimization [22].
based on the Laplace transform inverse kernel approximation by expressions containing hyperbolic functions sinh and cosh or their combinations [23].
The general formula [24] is based on the combining two approaches to the Laplace transform inverse kernel approximation est 2sinhe(aa st) and est coshe(aa st) to increase accuracy: where M is the maximum absolute value of the original f (t).
The parameter a is had to be found empirically. It is recommended to have a in the following range 2 a 6 [23], [25]. Brancˇ´ık and Smith [26] also notice that one of the ways to improve the accuracy of this approximation is to increase the parameter a.
f (t) :=
(
N N t
F (N 1)( N ); t where F (n 1)(x) is the n-th derivative of F (x).
The main issue is that method requires high order derivatives computation and is characterized by slow convergence [27].
Post-Widder method does not cause the overshoot and maintains monotonicity [11]. However, according to Horvath’s research [11], the CME method gives better approximation. Davies [28] also notices that the method rarely gives a high accuracy of the approximation.
Method is based on Laguerre polynomials. ILT is calculated by the formula [29]: f (t) 1 X qnln(t); n=0
x ln(x) = e 2 Ln(x); where ln is Laguerre function, qn are Laguerre coefficients, dependent on F .
Laguerre function: where Ln is Laguerre polynomial.
Laguerre polynomial:
n
Ln(x) = X
m=0
n ( x)m
m m!
(
Q(z) = 1 X qnzn = n=0 1 1 z
F
Cohen method is based on the power series F (s) as a function 1s . However, this method is not applicable when F (s) has a pole at 0, and the approximation also gets rapidly worse for larger values of t [11].
Honig-Hirdes and de Hoog methods are the variation of the Fourier series method. Wellekens uses the Fourier series method based on the Pade´ approximation [11]. Schapery method is described in [30], but according to [28], it rarely gives the high accuracy of the approximation.
V. CHOOSING METHODS FOR THE IMPLEMENTATION We will determine the most promising solutions for the considered methods.
Pade´ approximation, Post-Widder and Laguerre methods provide high computational complexity, because of calculating high order derivatives.
Gaver–Stehfest, Euler, Talbot methods suffer from the Gibbs oscillations, and since PDF is restricted to have negative values, these methods are also beyond our consideration.
So we choose four of the most promising methods for the comparison: the saddlepoint approximation, CME method, Zakian method and hyperbolic kernel approximation method.
VI. IMPLEMENTATION AND COMPARISON OF THE
SELECTED METHODS
The selected methods are implemented as Matlab functions (Matlab R2020a version was used), because Matlab [31] has one of the most powerful symbolic toolboxes. For the comparison with other methods the CME method is taken in the format suggested in [32] with some parameters pre-calculated by the authors of the method. The saddlepoint method, Zakian method, hyperbolic kernel approximation method are Matlab implemented according to [33]. Zakian method is realised with pre-calculated parameters Ki and i borrowed from [20]. The a parameter value is selected empirically for the hyperbolic kernel approximation method. Since Matlab has the built-in inverse Laplace transform function, it is reasonable to compare it with the selected methods.
The method’s average running time is stated as the performance criteria. The hardware in use is Intel(R) Core(TM) i5-1035G1 CPU @ 1.00GHz 1.19 GHz, RAM 8,00 GB with Windows 10.
The accuracy of the methods tested is assessed by two parameters: maximum absolute deviation "abs, to check for the sharp exceeds in the approximation, and 1-norm of the numerical error "n calculated by the formula: "n =
Z T 0 jf (t) fappr(t)jdt 1 XM jf (m) M m=1 fappr(m)j
Original function f (t) is compared with approximated PDF fappr(t) at M = 200 equidistant points.
The problem of the PDF mining from the flowgraph is stated for the beforehand known PDFs just for the testing purpose. The methods are checked on the approximation of exponential, normal, triangular and uniform distributions with well known Laplace transforms. The flowgraph with 5 nodes shown on fig. 1 is also investigated on the subject of the transition time distribution from the node 1 to the node 5. The flowgraph parameters are presented in the Table I. It should be noted that for this graph the 10-th MGF derivative calculation takes 2 seconds, the 12-th derivative calculation takes 4.5 seconds. It confirms that methods based on the high order derivatives calculation are not suitable for the implementation.
Built-in Matlab function demonstrates the best performance and the highest approximation accuracy (tables II - III) for the exponential distribution ( = 0:5) among all the methods in question. A visual comparison (Fig. 2) confirms these results.
The Matlab built-in function is unable to provide the inverse Laplace transform of the normal (Gauss) distribution with the expected value = 30 and standard deviation = 3. It returns an unevaluated call to ilaplace function. Saddlepoint approximation provides the set of values that exactly approximates the Gauss distribution only for the bounded segment [0;116]. Zakian and CME methods show sharp spikes in the values of the approximated function. Only the hyperbolic method gives a sufficiently accurate approximation (Table V). A visual comparison shows that CME method is characterized by the sharp exceed of values in small t, but Zakian method cannot provide an accurate approximation for t > 30 (Fig. 3). Zakian method ensures maximum performance, but because of low accuracy, we recommend to use hyperbolic approximation for this distribution. Fig. 2. PDF approximation. Exponential distribution
The saddlepoint approximation for the triangular distribution with parameters a = 10; b = 50; c = 30 is unable to find the answer: method fails with error ”Unable to find explicit solution”. The highest approximation accuracy and performance is demonstrated by the built-in Matlab function (Tables VI, VII). Zakian method does not provide the required approximation accuracy (Fig. 4).
The built-in Matlab function provides the highest performance with the approximation error of the order 10 2, that makes it the best option among the considered methods. The hyperbolic approximation method suffers from Gibbs oscillations. The Zakian method does not provide the required approximation accuracy.
The saddlepoint approximation for the uniform distribution with parameterss a = 20; b = 40 fails with the error ”Unable to find explicit solution”. The error for the built-in Matlab function and CME method is explained by the PDF approximation corners rounding (Fig. 5). The CME method provides the best accuracy for the order of N = 50 (Table IX).
For the flowgraph with different transition time distributions (Fig. 1, Table I) the saddlepoint approximation method fails with the error ”Unable to find explicit solution” and the same result is demonstrated by the built-in Matlab function that returns an unevaluated call to ilaplace function. The Zakian method does not provide the required approximation accuracy (Fig. 6). The best approximation is provided by the CME method (Table XI). For this case the CME method is faster than the hyperbolic approximation and the Zakian method as well (Table X).
In this paper we consider some well-known MGF inversion methods. Next, we have identified a group of methods that, according to the sources, are characterized by the least computational complexity and the absence of Gibbs oscillations. This group, which includes saddlepoint approximation, CME method, Zakian method and hyperbolic kernel approximation method, was tested on inverse Laplace transform problems for a number of distributions (exponential, normal, triangular, uniform). We also tested these methods on the meaningful problem represented by the five-state flowgraph. All methods are compared with the built-in Matlab function according to the criteria of accuracy and performance. The following conclusions were made on the test results: the Zakian method with the recommended in [10], [20] order N = 10 failed in most cases. So the calculation of parameters for higher orders is required; the hyperbolic kernel approximation method is not free from Gibbs oscillations and also requires the empirical selection of the a parameter, which makes it difficult to use method in the real cases when true distribution is unknown; the saddlepoint method gives an exact approximation for exponential and normal distributions, but it unable to provide the acceptable approximation of the PDF for uniform and triangular distributions; the built-in Matlab function provides the most accurate and fast result for exponential, triangular and uniform distributions, but does not cope with the MGF transform for the normal distribution and the flowgraph (fig. 1); thus, we can recommend the built-in Matlab function or the saddlepoint approximation only for a number of distributions; the CME method is the most stable among the considered methods. It is not subject to Gibbs oscillations, uses pre-calculated parameter values, so it has a relatively low computational complexity. In fact the definition of PDF here is reduced to summation and depends only on the complexity of the MGF expression. It is advisable to test the method on flowgraphs with large number of probabilistic transitions. [8] D. H. Collins, R. L. Warr, and A. V. Huzurbazar, “An introduction to statistical flowgraph models for engineering systems,” Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, vol. 227, no. 5, pp. 461–470, 2013. [9] I. Horva´th, Z. Talyiga´s, and M. Telek, “An optimal inverse laplace transform method without positive and negative overshoot–an integral based interpretation,” Electronic Notes in Theoretical Computer Science, vol. 337, pp. 87–104, 2018. [10] Q. Wang and H. Zhan, “On different numerical inverse laplace methods for solute transport problems,” Advances in Water Resources, vol. 75, pp. 80–92, 2015. [11] G. Horva´th, I. Horva´th, S. A.-D. Almousa, and M. Telek, “Numerical inverse laplace transformation using concentrated matrix exponential distributions,” Performance Evaluation, vol. 137, p. 102067, 2020. [12] J. Abate and W. Whitt, “A unified framework for numerically inverting laplace transforms,” INFORMS Journal on Computing, vol. 18, no. 4, pp. 408–421, 2006. [13] G. Horva´th, I. Horva´th, S. A.-D. Almousa, and M. Telek, “High order low variance matrix-exponential distributions,” in The Tenth International Conference on Matrix-Analytic Methods in Stochastic Models (MAM10), vol. 2, 2019. [14] I. Horva´th, O. Sa´fa´r, M. Telek, and B. Za´mbo´ , “Concentrated matrix exponential distributions,” in European Workshop on Performance Engineering. Springer, 2016, pp. 18–31. [15] D. P. Gaver Jr, “Observing stochastic processes, and approximate transform inversion,” Operations Research, vol. 14, no. 3, pp. 444–459, 1966. [16] H. Stehfest, “Algorithm 368: Numerical inversion of laplace transforms [d5],” Communications of the ACM, vol. 13, no. 1, pp. 47–49, 1970. [17] R. Jacquot, J. Steadman, and C. Rhodine, “The gaver-stehfest algorithm for approximate inversion of laplace transforms,” IEEE Circuits & Systems Magazine, vol. 5, no. 1, pp. 4–8, 1983. [18] J. Abate, G. L. Choudhury, and W. Whitt, “An introduction to numerical transform inversion and its application to probability models,” in Computational probability. Springer, 2000, pp. 257–323. [19] A. Talbot, “The accurate numerical inversion of laplace transforms,”
IMA Journal of Applied Mathematics, vol. 23, no. 1, pp. 97–120, 1979. [20] V. Zakian and M. Edwards, “Tabulation of constants for full grade iMN approximants,” Mathematics of Computation, vol. 32, no. 142, pp. 519– 531, 1978. [21] V. Zakian, “Optimisation of numerical inversion of laplace transforms,”
Electronics Letters, vol. 6, no. 21, pp. 677–679, 1970. [22] V. Zakian and R. Littlewood, “Numerical inversion of laplace transforms by weighted least-squares approximation,” The Computer Journal, vol. 16, no. 1, pp. 66–68, 1973. [23] J. Valsa and L. Brancˇik, “Approximate formulae for numerical inversion of laplace transforms,” International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, vol. 11, no. 3, pp. 153–166, 1998. [24] N. Smith and L. Brancik, “Comparative study on one-dimensional numerical inverse laplace transform methods for electrical engineering,” Elektrorevue, vol. 18, no. 1, pp. 1–6, 2016. [25] N. A.-Z. R-Smith and L. Brancˇ´ık, “Convergence acceleration techniques for proposed numerical inverse laplace transform method,” in 2016 24th Telecommunications Forum (TELFOR). IEEE, 2016, pp. 1–4. [26] L. Brancˇ´ık and N. Smith, “Two approaches to derive approximate formulae of nilt method with generalization,” in 2015 38th International Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO). IEEE, 2015, pp. 155–160. [27] A. M. Cohen, “Quadrature methods,” in Numerical methods for Laplace transform inversion. Springer, 2007, pp. 71–101. [28] B. Davies and B. Martin, “Numerical inversion of the laplace transform: a survey and comparison of methods,” Journal of computational physics, vol. 33, no. 1, pp. 1–32, 1979. [29] J. Abate, G. L. Choudhury, and W. Whitt, “On the laguerre method for numerically inverting laplace transforms,” 1995. [30] K. L. Kuhlman, “Review of inverse laplace transform algorithms for laplace-space numerical approaches,” Numerical Algorithms, vol. 63, no. 2, pp. 339–355, 2013. [31] “[electronic source] matlab - mathworks - matlab simulink,” (Access date: 25.04.2021). [Online]. Available: https://uk.mathworks.com/products/matlab.html