=Paper=
{{Paper
|id=Vol-2953/SEIM_2021_paper_15
|storemode=property
|title=Moment Generating Function to Probability Density Function Transform Methods
|pdfUrl=https://ceur-ws.org/Vol-2953/SEIM_2021_paper_15.pdf
|volume=Vol-2953
|authors=Daria Gaponenko,Aleksandr Sidnev
}}
==Moment Generating Function to Probability Density Function Transform Methods==
https://ceur-ws.org/Vol-2953/SEIM_2021_paper_15.pdf
Moment Generating Function to Probability Density
Function Transform Methods
1st Darya Gaponenko 2nd Aleksandr Sidnev
Peter the Great St.Petersburg Polytechnic University Peter the Great St.Petersburg Polytechnic University
St.Petersburg, Russia St.Petersburg, Russia
daria.gaponenko@yandex.ru alex2922@gmail.com
Abstract—The moment generating function completely char- of eight numerical inverse Laplace transform methods. In
acterizes the random variable distribution. At the same time, the section V we select methods for the software implementation
natural desire of the researcher is to obtain the density and the and comparison on the basis of the preliminary analysis.
distribution function as more informative characteristics. The
article analyzes popular transition methods from the moment Finally section VI is devoted to the results discussing.
generating function to the probability density, in particular, the
saddlepoint method, the method based on the Heaviside theorem II. PAD É APPROXIMATION
and Padé approximation, a significant number of the inverse
Laplace transform numerical methods. Computational complex-
The Heaviside theorem is suggested to obtain the PDF
ity and practical feasibility of the methods are investigated for a applying the inverse Laplace transform of the MGF [5], [6]:
number of practical problems of the inverse Laplace transform. h U (s) i q
Recommendations on the optimal method choice are given and −1
X U (αk )
L = ∗ eαk ∗t ,
the directions for the further development of the topic are R(s) R0 (αk )
suggested as well. k=1
Index Terms—probability density function, moment generating where αk are the roots of the polynomial R(s).
function, numerical Laplace inversion methods, Padé approxima-
tion
MGF Padé approximation can be straightforward provided
as series:
I. I NTRODUCTION ∞ Pp j
j=0 aj ∗ x
X
i
Moment generating function (MGF) is an important statisti- ci ∗ x = Pq k
cal characteristic that uniquely describes the random variable i=0 k=0 bk ∗ x
distribution. Since the probability density function (PDF) is Padé approximation is applicable to MGF because:
related to the MGF through the inverse Laplace transform [1]:
∞
h i X µn
f (t) = L−1 M (−s) , M (s) = ∗ sn ,
n=0
n!
the computational complexity can be significantly reduced by
where µn is n-th moment.
using the MGF instead of PDF. For example, MGF of the two
Ren [5] compares three methods: Padé approximation, max-
random variables sum is calculated as the product of the MGF
imum entropy method and saddlepoint approximation. It is
of these variables [1].
noted that the Padé approximation is more informative and
Flowgraph models operate MGF to obtain the stochastic
always provides an analytical approximation for the unknown
process total time distribution. They can be applied in different
original PDF.
areas. For example, in reliability theory for the uptime estima-
But this approach has a significant disadvantage: Padé
tion, in medical statistics for the disease progression and the
approximation requires the high order derivatives computa-
treatment effectiveness analysis [2], [3], for the power system
tion (because the n-th moment is the n-th derivative of the
survivability analysis and it’s weak points identification [4].
MGF). This fact leads to great computational complexity
Padé approximation, Heaviside’s theorem, Saddlepoint ap-
often exceeding the hardware possibilities for the high order
proximation and numerical inverse Laplace transform methods
derivatives calculation in the symbolic form.
are used to obtain the PDF as the MGF inversion (while
the MGF is acquired as the result of the flowgraph specific III. S ADDLEPOINT APPROXIMATION
analysis). The transform method must provide an acceptable
approximation error with a relatively low computational com- Ren [5] also considers the technique of saddlepoint approx-
plexity and must generate non-negative PDF values. imation for the MGF-PDF transform:
The paper is organised as follows. Section II introduces the
MGF transform via Padé approximation, section III presents n 12
the Saddlepoint approximation, section IV provides the review f (t) = ∗ exp(n ∗ (K(s) − st)),
2 ∗ π ∗ K 00 (s)
Copyright © 2021 for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
where K(s) = log(M (s)), K 0 (s) = t, n is the number of ran- where ck = −αAuk vk 1.
dom variables. For one random variable PDF approximation A matrix-exponential distribution is called concentrated if
is equal to [7]: its squared coefficient of variation (SCV) is minimal [13].
− 12 The density of the matrix-exponential distribution with the
f (t) = 2 ∗ π ∗ K 00 (s) ∗ exp(K(s) − st) minimal SCV is known analytically only for the order N < 3.
1 Therefore, the matrix-exponential distributions required for the
The error of this approximation is estimated as O(n− 2 ) [5].
approximation were calculated using numerical optimization
Collins [8] notes this method to have less efficiency than the
by the expression:
inverse Laplace transform.
IV. I NVERSE L APLACE TRANSFORM METHODS N −1
2 N
There are over a hundred different numerical inverse Y X
Laplace transform (ILT) methods. Being applied to the various f (t) = ce−λt cos2 (ωt − φj ) = ηk e−βk t ,
j=0 k=1
types of functions they differ in computational complexity and
approximation accuracy. where c, λ, ω, φj are positive real values.
Some ILT methods have Gibbs oscillations – an oscillatory Thus, the values ηk and βk providing the minimum SCV
error which occurs when a discontinuous function is approxi- for the matrix-exponential distribution are the coefficients of
mated with Fourier series. Gibbs oscillations do not disappear the Abate-Whitt framework and are obtained from numerical
when more terms are added to the approximation (the ”size” of optimization [14].
the excess does not decrease vertically, but only horizontally
It is guaranteed that the CME method is free from Gibbs
[9]).
oscillations and preserves monotonicity [11]. Many of the
A. Abate–Whitt framework known methods do not have these properties: for example,
Different ILT numerical methods can be combined within the fn (t) function for the Euler and Gaver–Stehfest methods
the same mathematical framework. Abate and Whitt propose can take negative values.
the unified formula (1) and notice, that such algorithms as Advantages of the method: it does not cause overshoot,
Gaver–Stehfest, Talbot and Fourier series method with Euler maintains monotonicity, is accurate when using floating point
summation can fit into this framework [10]. arithmetic, the quality of the approximation improves with
Abate–Whitt framework [11]: increasing order [11].
N 2) Gaver–Stehfest method: Method is based on the three-
X ηk βk
f (t) ≈ F( ), (1) parameter exponential PDF properties (Gaver method [15]).
t t Stehfest [16] improves method by taking the weighted average
k=1
where t > 0, nodes βk and weights ηk (internal and external of a Gaver approximations sequence for fixed t.
scaling constants) are real or complex numbers that depend on Gaver–Stehfest algorithm does not use complex numbers;
N , but do not depend on the transform F or the time argument weights and nodes are real numbers and are calculated as
t [12]. follows [10]:
Abate–Whitt framework has low computational complexity, βk = k ∗ ln(2)
but it is characterized by a deterioration of the approximation
on periodic functions at large values of t [9]. Since PDF is N
not a periodic function, this disadvantage can be ignored. ηk = (−1)( 2 +k) ln(2)∗
Various ILT methods in this framework differ only by
approaches to ηk and βk calculation. min(k, N
2 ) N
X j 2 2j!
Several methods belonging to the Abate-Whitt framework ∗ ,
( N2 − j)!j!(j − 1)!(k − j)!(2j − k)!
are considered below. j=b k+1
2 c
1) CME: Method is based on the concentrated matrix-
exponential distribution [11]. where N is the number of terms used in the equation (the
Matrix-exponential distributions of order N contain positive method is defined only for even N ). Wang in [10] notices that
random variables with PDF f (t) = −αAeAt 1, at t ≥ 0, where N is recommended to be taken in the range from 10 to 14.
α is a real row vector of length N , A is an NxN real matrix, The main disadvantage of the method is Gibbs oscillations
and 1 is a column vector of ones of size N [11]. [11].
PLet A be diagonalizable with spectral decomposition A = Abate and Whitt [12] show that Euler and Talbot methods
N
k=1 uk λk νk , where λk are the eigenvalues, uk are the right are more efficient than show that Euler and Talbot methods are
eigenvectors and νk are the left eigenvectors of A with νk uk = more efficient than those of Gaver–Stehfest. It is also noted
1, then the PDF of the matrix-exponential distribution can be that the algorithm implementation requires high computational
written as: accuracy due to manipulations with large numbers [12], [17].
N
f (t) =
X
ck eλk t , However, the Gaver–Stehfest algorithm has the advantage of
k=1
using only real numbers.
3) Euler method: is an implementation of the Fourier series 6) Hyperbolic kernel approximation method: Method is
method using Euler summation to accelerate convergence [18]. based on the Laplace transform inverse kernel approximation
Method is applicable only to odd orders and assumes nodes by expressions containing hyperbolic functions sinh and cosh
with positive imaginary parts. or their combinations [23].
Nodes and weights are calculated as follows [11]: The general formula [24] is based on the combining two
approaches to the Laplace transform inverse kernel approxi-
(N − 1) ln(10) ea
mation est ≈ 2sinh(a−st) ea
and est ≈ cosh(a−st) to increase
βk = + πi(k − 1)
6 accuracy:
N −1
ηk = 10 6 (−1)k ξk , ea 1 a XN
a π
f (t) = F( ) + (−1)n Re F + in +
where ξ1 = 12 ; 2t 2 t n=1
t t
ξk = 1, 2 < k < n+1
2 ; a 1 π
1 +Im F + i(n − ) (3)
ξn = n−1 ;
2 2 t 2 t
n−1
− n−1
, 1 < k < n−1
ξn−k = ξn−k+1 + 2 2
k
2
2 The absolute error depends on the a as follows [24]:
Method suffers from Gibbs oscillations [11].
4) Talbot method: Method is based on deforming the εaM ≈ M e−4a ,
contour integral in the Bromwich inversion [19]. It assumes where M is the maximum absolute value of the original f (t).
nodes with positive imaginary parts and also applies values of
The parameter a is had to be found empirically. It is
βk with a negative real part, which is usually outside the region
recommended to have a in the following range 2 ≤ a ≤ 6
of convergence, but there can be an analytic continuation of
[23], [25]. Brančı́k and Smith [26] also notice that one of
F (s) [9].
the ways to improve the accuracy of this approximation is to
Nodes and weights are calculated as follows [11]: increase the parameter a.
2N
β1 = B. Post–Widder method
5
ILT is calculated by the formula [27]:
2(k − 1)π (k − 1)π
βk = cot +i (−1)N −1 N N (N −1) N
5 N f (t) := F ( ), (4)
(N − 1)! t t
1
η1 = eβ1
5 where F (n−1) (x) is the n-th derivative of F (x).
The main issue is that method requires high order derivatives
2h (k − 1)π (k − 1)π 2 computation and is characterized by slow convergence [27].
ηk = 1 + i 1 + cot −
5 N N Post-Widder method does not cause the overshoot and
(k − 1)π i βk maintains monotonicity [11]. However, according to Horvath’s
−i cot e research [11], the CME method gives better approximation.
N
Davies [28] also notices that the method rarely gives a high
Method is not free from Gibbs oscillations [11]. accuracy of the approximation.
5) Zakian method: Method is based on Fourier series
method with Padé approximation [11]. C. Laguerre method
General formula is [10]: Method is based on Laguerre polynomials. ILT is calculated
N by the formula [29]:
1 X αk
f (t) = Re{Kk F ( )}, (2) ∞
t t f (t) ≈
X
qn ln (t), (5)
k=1
n=0
where Kk and αk are real or complex constants. N represents
the number of terms used in the summation. According to where ln is Laguerre function, qn are Laguerre coefficients,
[10], [20] a sufficient value of N in most cases is N = 10. dependent on F .
Coefficients for 1 ≤ N ≤ 10 are calculated in [20]. Laguerre function:
Zakian proposes two methods for choosing coefficients Kk x
ln (x) = e− 2 Ln (x), (6)
and αk . In the first he compares the Laplace transform δ(t, u)
(a rational function) with the Laplace transform of δn (u − t) where Ln is Laguerre polynomial.
(an exponential function) and chooses the coefficients so that Laguerre polynomial:
the rational functions are equal to the classical Padé approx- n
imations of the exponential function [21]. Another method X n (−x)m
Ln (x) = (7)
includes least squares optimization [22]. m=0
m m!
Laguerre coefficients generating function: i5-1035G1 CPU @ 1.00GHz 1.19 GHz, RAM 8,00 GB with
∞ Windows 10.
X 1 1+z
Q(z) = qn z n = F (8) The accuracy of the methods tested is assessed by two
n=0
1−z 2(1 − z) parameters: maximum absolute deviation εabs , to check for
Formulas 5 – 8 can be rewritten as follows [11]: the sharp exceeds in the approximation, and 1-norm of the
N −1 N −1
numerical error εn calculated by the formula:
X 1 X k!
f (t) = F (j) ( ) l (t),
2 (k − j)! k
(9)
2 (j!) Z T M
j=0 k=j 1 X
εn = |f (t) − fappr (t)|dt ≈ |f (m) − fappr (m)|
where F (j) (s) is the j-th derivative of F (s). 0 M m=1
Davies [28] notices that Laguerre polynomials give an
Original function f (t) is compared with approximated PDF
accurate approximation over a wide range of functions.
fappr (t) at M = 200 equidistant points.
The disadvantage of this method is that computation re-
The problem of the PDF mining from the flowgraph is stated
quires high order differentiation. Laguerre method also suffers
for the beforehand known PDFs just for the testing purpose.
from Gibbs oscillations [11].
The methods are checked on the approximation of exponential,
D. Other methods normal, triangular and uniform distributions with well known
Cohen method is based on the power series F (s) as a Laplace transforms. The flowgraph with 5 nodes shown on
function 1s . However, this method is not applicable when F (s) fig. 1 is also investigated on the subject of the transition time
has a pole at 0, and the approximation also gets rapidly worse distribution from the node 1 to the node 5. The flowgraph
for larger values of t [11]. parameters are presented in the Table I. It should be noted
Honig-Hirdes and de Hoog methods are the variation of that for this graph the 10-th MGF derivative calculation takes
the Fourier series method. Wellekens uses the Fourier series 2 seconds, the 12-th derivative calculation takes 4.5 seconds.
method based on the Padé approximation [11]. Schapery It confirms that methods based on the high order derivatives
method is described in [30], but according to [28], it rarely calculation are not suitable for the implementation.
gives the high accuracy of the approximation.
V. C HOOSING METHODS FOR THE IMPLEMENTATION
We will determine the most promising solutions for the
considered methods.
Padé 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. Fig. 1. Flowgraph
VI. I MPLEMENTATION AND COMPARISON OF THE
SELECTED METHODS
TABLE I
The selected methods are implemented as Matlab functions F LOWGRAPH PARAMETERS
(Matlab R2020a version was used), because Matlab [31] has Edge Probability Distribution Parameters
one of the most powerful symbolic toolboxes. For the com- q12 0.2 Exponential λ = 0.2
parison with other methods the CME method is taken in the q14 0.8 Uniform a = 2.4, b = 9.6
format suggested in [32] with some parameters pre-calculated q23 1 Triangular a = 1.6, b = 6.4, c = 4
q34 0.9 Exponential λ = 0.5
by the authors of the method. The saddlepoint method, Zakian q35 0.1 Uniform a = 6, b = 24
method, hyperbolic kernel approximation method are Matlab q43 0.5 Triangular a = 4, b = 16, c = 10
implemented according to [33]. Zakian method is realised with q45 0.5 Exponential λ = 0.3
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 A. Exponential distribution
inverse Laplace transform function, it is reasonable to compare Built-in Matlab function demonstrates the best performance
it with the selected methods. and the highest approximation accuracy (tables II - III) for the
The method’s average running time is stated as the per- exponential distribution (λ = 0.5) among all the methods in
formance criteria. The hardware in use is Intel(R) Core(TM) question. A visual comparison (Fig. 2) confirms these results.
TABLE II sharp exceed of values in small t, but Zakian method cannot
E XPONENTIAL DISTRIBUTION APPROXIMATION TIME provide an accurate approximation for t > 30 (Fig. 3). Zakian
Matlab function, s Zakian, s Saddlepoint, s method ensures maximum performance, but because of low
0.0110 ± 0.0005 0.055 ± 0.005 0.111 ± 0.002 accuracy, we recommend to use hyperbolic approximation for
Order CME, s Hyperbolic, s this distribution.
5 0.319 ± 0.009 0.098 ± 0.004
10 0.356 ± 0.009 0.182 ± 0.012
TABLE IV
20 0.406 ± 0.007 0.353 ± 0.009
N ORMAL DISTRIBUTION APPROXIMATION TIME
50 0.5915 ± 0.0104 0.860 ± 0.012
100 0.888 ± 0.019 1.714 ± 0.021 Matlab function, s Zakian, s Saddlepoint, s
- 0.072 ± 0.004 0.112 ± 0.004
Order CME, s Hyperbolic, s
TABLE III 5 0.319 ± 0.005 0.096 ± 0.003
E XPONENTIAL DISTRIBUTION APPROXIMATION ERRORS 10 0.376 ± 0.014 0.219 ± 0.009
20 0.44 ± 0.02 0.420 ± 0.018
Matlab function Zakian Saddlepoint 50 2.92 ± 0.03 2.955 ± 0.009
εn εabs εn εabs εn εabs 100 4.55 ± 0.06 5.89 ± 0.02
0 0 2.51E-03 6.81E-03 3.25E-04 2.56E-02
CME Hyperbolic
Order εn εabs εn εabs
5 3.40E-04 3.70E-03 2.68E-03 1.79E-02
TABLE V
10 5.62E-05 8.02E-04 1.62E-03 6.11E-03
N ORMAL DISTRIBUTION APPROXIMATION ERRORS
20 1.03E-05 1.74E-04 1.02E-03 3.15E-03
50 9.60E-07 2.30E-05 8.29E-04 4.36E-03 Matlab function Zakian Saddlepoint, t ∈ [0; 116]
100 3.60E-07 9.00E-06 5.54E-04 5.00E-03 εn εabs εn εabs εn εabs
- - p. inf. p. inf. 0 0
CME Hyperbolic
Order εn εabs εn εabs
B. Normal distribution 5 2.94E-03 6.03E-02 5.03E-03 1.95E-02
10 1.75E+00 3.49E+02 2.82E-03 6.01E-03
The Matlab built-in function is unable to provide the inverse 20 p. inf. p. inf. 8.91E-04 2.92E-03
Laplace transform of the normal (Gauss) distribution with the 50 p. inf. p. inf. 4.66E-05 2.44E-03
100 p. inf. p. inf. 1.85E-05 2.44E-03
expected value µ = 30 and standard deviation σ = 3. It returns
*p. inf. (practically infinite) stands for values larger than 1000
an unevaluated call to ilaplace function. Saddlepoint approx-
imation 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
Fig. 3. PDF approximation. Normal distribution
C. Triangular distribution
The saddlepoint approximation for the triangular distribu-
tion with parameters a = 10; b = 50; c = 30 is unable
Fig. 2. PDF approximation. Exponential distribution to find the answer: method fails with error ”Unable to find
explicit solution”. The highest approximation accuracy and The built-in Matlab function provides the highest performance
performance is demonstrated by the built-in Matlab function with the approximation error of the order 10−2 , that makes it
(Tables VI, VII). Zakian method does not provide the required the best option among the considered methods. The hyperbolic
approximation accuracy (Fig. 4). approximation method suffers from Gibbs oscillations. The
Zakian method does not provide the required approximation
TABLE VI accuracy.
T RIANGULAR DISTRIBUTION APPROXIMATION TIME
Matlab function, s Zakian, s Saddlepoint, s TABLE VIII
0.0108 ± 0.0005 0.16 ± 0.03 - U NIFORM DISTRIBUTION APPROXIMATION TIME
Order CME, s Hyperbolic, s
Matlab function, s Zakian, s Saddlepoint, s
5 0.46 ± 0.08 0.144 ± 0.005
0.0114 ± 0.0016 0.18 ± 0.08 -
10 0.68 ± 0.02 2.56 ± 0.03
20 1.068 ± 0.018 5.03 ± 0.04 Order CME, s Hyperbolic, s
50 2.34 ± 0.08 13 ± 1 5 0.331 ± 0.004 0.21 ± 0.08
100 3.39 ± 0.04 25.56 ± 0.17 10 0.6 ± 0.10 0.42 ± 0.16
20 0.83 ± 0.03 3.55 ± 0.03
50 1.70 ± 0.07 10.2 ± 0.2
100 2.98 ± 0.11 21.6 ± 0.8
TABLE VII
T RIANGULAR DISTRIBUTION APPROXIMATION ERRORS
Matlab function Zakian Saddlepoint TABLE IX
εn εabs εn εabs εn εabs U NIFORM DISTRIBUTION APPROXIMATION ERRORS
2.50E-04 5.00E-02 1.17E-03 5.28E-02 - -
CME Hyperbolic Matlab function Zakian Saddlepoint
Order εn εabs εn εabs εn εabs εn εabs εn εabs
5 9.24E-04 5.92E-02 3.86E-03 5.23E-02 2.50E-04 2.50E-02 3.14E-03 2.72E-02 - -
10 4.08E-04 5.43E-02 1.25E-03 5.15E-02 CME Hyperbolic
20 2.88E-04 5.20E-02 4.13E-04 5.07E-02 Order εn εabs εn εabs
50 2.58E-04 5.07E-02 2.84E-04 5.03E-02 5 1.93E-03 2.56E-02 5.59E-03 2.68E-02
100 2.55E-04 5.05E-02 2.59E-04 5.02E-02 10 8.93E-04 2.55E-02 2.50E-03 2.60E-02
20 4.45E-04 2.52E-02 1.55E-03 2.58E-02
50 2.61E-04 2.51E-02 8.91E-04 2.52E-02
100 2.52E-04 2.51E-02 6.91E-04 2.52E-02
Fig. 4. PDF approximation. Triangular distribution
D. Uniform distribution Fig. 5. PDF approximation. Uniform distribution
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 E. Flowgraph
Matlab function and CME method is explained by the PDF For the flowgraph with different transition time distributions
approximation corners rounding (Fig. 5). The CME method (Fig. 1, Table I) the saddlepoint approximation method fails
provides the best accuracy for the order of N = 50 (Table IX). with the error ”Unable to find explicit solution” and the same
result is demonstrated by the built-in Matlab function that according to the sources, are characterized by the least com-
returns an unevaluated call to ilaplace function. The Zakian putational complexity and the absence of Gibbs oscillations.
method does not provide the required approximation accuracy This group, which includes saddlepoint approximation, CME
(Fig. 6). The best approximation is provided by the CME method, Zakian method and hyperbolic kernel approximation
method (Table XI). For this case the CME method is faster method, was tested on inverse Laplace transform problems
than the hyperbolic approximation and the Zakian method as for a number of distributions (exponential, normal, triangular,
well (Table X). uniform). We also tested these methods on the meaningful
problem represented by the five-state flowgraph. All methods
TABLE X are compared with the built-in Matlab function according
F LOWGRAPH APPROXIMATION TIME to the criteria of accuracy and performance. The following
Matlab function, s Zakian, s Saddlepoint, s conclusions were made on the test results:
- 5.09 ± 0, 14 - • the Zakian method with the recommended in [10], [20]
Order CME, s Hyperbolic, s
5 2.48 ± 0.05 10.20 ± 0.12
order N = 10 failed in most cases. So the calculation of
10 4.83 ± 0.11 20.25 ± 0.19 parameters for higher orders is required;
20 8.77 ± 0.12 32.40 ± 0.18 • the hyperbolic kernel approximation method is not free
50 21.8 ± 0.6 81.50 ± 0.82 from Gibbs oscillations and also requires the empirical
100 34.5 ± 1.6 182 ± 8
selection of the a parameter, which makes it difficult to
use method in the real cases when true distribution is
unknown;
TABLE XI • the saddlepoint method gives an exact approximation for
F LOWGRAPH APPROXIMATION ERRORS
exponential and normal distributions, but it unable to
Matlab function Zakian Saddlepoint provide the acceptable approximation of the PDF for
εn εabs εn εabs εn εabs uniform and triangular distributions;
- - 2.26E-03 9.85E-03 - -
• the built-in Matlab function provides the most accurate
CME Hyperbolic
Order εn εabs εn εabs
and fast result for exponential, triangular and uniform
5 5.52E-04 7.08E-03 3.87E-03 9.15E-03 distributions, but does not cope with the MGF transform
10 5.34E-04 7.95E-03 2.26E-03 1.06E-02 for the normal distribution and the flowgraph (fig. 1);
20 5.67E-04 9.83E-03 1.04E-03 1.08E-02 • thus, we can recommend the built-in Matlab function
50 5.80E-04 1.05E-02 6.21E-04 1.04E-02
100 5.81E-04 1.04E-02 5.90E-04 1.05E-02 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.
R EFERENCES
[1] M. G. Bulmer, Principles of statistics. Courier Corporation, 1979.
[2] B. Garcı́a-Mora, C. Santamarı́a, G. Rubio, and J. L. Pontones, “Bayesian
prediction for flowgraph models with covariates. an application to blad-
der carcinoma,” Journal of Computational and Applied Mathematics,
vol. 291, pp. 85–93, 2016.
[3] B. Garcı́a-Mora, C. Santamarı́a, and G. Rubio, “Markovian modeling
for dependent interrecurrence times in bladder cancer,” Mathematical
Methods in the Applied Sciences, vol. 43, no. 14, pp. 8302–8310, 2020.
[4] A. Fetanat, “A decision support tool for modeling time-to-event data
in power systems.” Journal of Theoretical & Applied Information
Technology, vol. 60, no. 2, 2014.
[5] Y. Ren, “The methodology of flowgraph models,” Ph.D. dissertation,
Fig. 6. PDF approximation of the transition time between the first and the The London School of Economics and Political Science (LSE), 2011.
last nodes of the graph [6] A. Sidnev, “Business process analytical modeling — business analysis
body of knowledge required component,” in The European Proceedings
of Social Behavioural Sciences. Edited by: Prof. Valeria Chernyavskaya
VII. C ONCLUSION and Prof. Holger Kuße. Future Academy, 2018, pp. 1337–1348.
[7] A. V. Huzurbazar, “Flowgraph models: a bayesian case study in con-
In this paper we consider some well-known MGF inversion struction engineering,” Journal of statistical planning and inference, vol.
methods. Next, we have identified a group of methods that, 129, no. 1-2, pp. 181–193, 2005.
[8] D. H. Collins, R. L. Warr, and A. V. Huzurbazar, “An introduction to [32] “[electronic source] github - ghorvath78/iltcme: Inverse laplace trans-
statistical flowgraph models for engineering systems,” Proceedings of form based on concentrated matrix-exponential functions,” (Access date:
the Institution of Mechanical Engineers, Part O: Journal of Risk and 15.07.2020). [Online]. Available: https://github.com/ghorvath78/iltcme
Reliability, vol. 227, no. 5, pp. 461–470, 2013. [33] “[electronic source] github - dariagap/mgf to pdf transform methods,”
[9] I. Horváth, Z. Talyigás, and M. Telek, “An optimal inverse laplace (Access date: 14.02.2021). [Online]. Available:
transform method without positive and negative overshoot–an integral https://github.com/dariaGap/mgf to pdf transform methods
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. Horváth, I. Horvá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. Horváth, I. Horváth, S. A.-D. Almousa, and M. Telek, “High order
low variance matrix-exponential distributions,” in The Tenth Interna-
tional Conference on Matrix-Analytic Methods in Stochastic Models
(MAM10), vol. 2, 2019.
[14] I. Horváth, O. Sáfár, M. Telek, and B. Zámbó, “Concentrated matrix
exponential distributions,” in European Workshop on Performance En-
gineering. 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 numer-
ical 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 iM N
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 trans-
forms by weighted least-squares approximation,” The Computer Journal,
vol. 16, no. 1, pp. 66–68, 1973.
[23] J. Valsa and L. Branč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. Brančı́k, “Convergence acceleration techniques
for proposed numerical inverse laplace transform method,” in 2016 24th
Telecommunications Forum (TELFOR). IEEE, 2016, pp. 1–4.
[26] L. Brančı́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