Template for an Acta IMEKO event paper ACTA IMEKO ISSN: 2221-870X November 2016, Volume 5, Number 3, 32-44 ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 32 Numerical inversion of a characteristic function: An alternative tool to form the probability distribution of output quantity in linear measurement models Viktor Witkovský Institute of Measurement Science, Slovak Academy of Sciences, Dúbravská cesta 9, SK-841 04 Bratislava, Slovakia Section: RESEARCH PAPER Keywords: GUM; linear measurement model; probability density function; characteristic function; numerical inversion; Gil-Pelaez inversion formulae; FFT algorithm Citation: Viktor Witkovský, Numerical inversion of a characteristic function: An alternative tool to form the probability distribution of output quantity in linear measurement models, Acta IMEKO, vol. 5, no. 3, article 6, November 2016, identifier: IMEKO-ACTA-05 (2016)-03-06 Section Editor: Franco Pavese, Italy Received April 27, 2016; In final form April 27, 2016; Published November 2016 Copyright: © 2016 IMEKO. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited Funding: The work was supported by the Slovak Research and Development Agency, project APVV-15-0295, and by the Scientific Grant Agency VEGA of the Ministry of Education of the Slovak Republic and the Slovak Academy of Sciences, by the projects VEGA 2/0047/15 and VEGA 2/0011/16 Corresponding author: Viktor Witkovský, e-mail: witkovsky@savba.sk 1. INTRODUCTION The basic working tool in measurement uncertainty analysis, as advocated in the current revision (under preparation) of the Guide to the expression of uncertainty in measurement (GUM) [1], and consistent with its Supplement 1 – Propagation of distributions using a Monte Carlo method [2], is the state-of-knowledge PDF about the quantity (true value of measurand), based on the currently available information. The state-of-knowledge PDF quantifies the degree of belief about the values that can be assigned to the quantity based on the available information. The expectation and the standard deviation of this PDF (if they exist) are used to report the measurement result and the associated (standard) measurement uncertainty. Although the latest GUM development emphasizes the Bayesian view of probability in the evaluation of measurement uncertainty, it should be clearly stated and understood that this approach is not based on the strict Bayesian principles of statistical inference (i.e. straightforward application of the ABSTRACT Measurement uncertainty analysis based on combining the state-of-knowledge distributions requires evaluation of the probability density function (PDF), the cumulative distribution function (CDF), and/or the quantile function (QF) of a random variable reasonably associated with the measurand. This can be derived from the characteristic function (CF), which is defined as a Fourier transform of its probability distribution function. Working with CFs provides an alternative and frequently much simpler route than working directly with PDFs and/or CDFs. In particular, derivation of the CF of a weighted sum of independent random variables is a simple and trivial task. However, the analytical derivation of the PDF and/or CDF by using the inverse Fourier transform is available only in special cases. Thus, in most practical situations, a numerical derivation of the PDF/CDF from the CF is an indispensable tool. In metrological applications, such approach can be used to form the probability distribution for the output quantity of a measurement model of additive, linear or generalized linear form. In this paper we propose new original algorithmic implementations of methods for numerical inversion of the characteristic function which are especially suitable for typical metrological applications. The suggested numerical approaches are based on the Gil-Pelaez inverse formulae and on using the approximation by discrete Fourier transform and the fast Fourier transform (FFT) algorithm for computing PDF/CDF of the univariate continuous random variables. As illustrated here, for typical metrological applications based on linear measurement models the suggested methods are an efficient alternative to the standard Monte Carlo methods. mailto:witkovsky@savba.sk ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 33 Bayes' theorem). For more details and further discussion see, e.g., [3]-[5], or [6], and also [7] and [8]. In fact, the GUM approach is based on using a well-defined functional relationship between the mutually inter-related quantities for propagating the state-of-knowledge PDFs of the input quantities, represented by the random variables (RVs), into the state-of-knowledge PDF of the output quantity – which is believed to be a RV reasonably associated with the measurand. Frequently, it is suggested to use a well-known functional relationship (based, e.g., on the physical and/or geometrical laws) between the true value of the measurand and the true values of the other influencing input variables, which is typically expressed by the measurement equation of the measurement model. Obviously, such PDF of the output quantity represents currently available knowledge (limited, but hopefully the best to date) about the measurand, i.e. it expresses probability distribution of the values being attributed to a quantity (the measurand), based on information used (which could be rather limited and/or heavily biased). This interpretation is consistent to the original GUM definition of the uncertainty in measurement (see [1], clause 2.2.3), which is defined as a parameter, associated with the result of a measurement, that characterizes the dispersion of the values that could be reasonably attributed to the measurand. However, the derived term coverage interval is inconsistent with this interpretation, for more details and discussion see [8]. In fact, without imposing further (well and clearly defined) model assumptions and optimality criteria for selecting and combining the information, it can be only hardly expected that the presented result shall represent the best (in what sense?) estimate of the true measurand value. On the other hand, the proposed GUM approach could be well accepted as a (simple) method for combining experimental results with the expert judgment in order to get comprehensive characterization of our knowledge about the true value of measurand, based on all currently available information, albeit without the possibility of guaranteeing the (otherwise naturally) required statistical properties and/or optimality criteria. If this is the goal, other means and/or subsequent analysis should be applied and properly used. As already mentioned, the term coverage interval (introduced in [2], and defined as the interval containing the (true) value of a quantity with a stated probability (say 95 %), based on the information available) is not properly used in this context. Hence, as an alternative to the 95 % coverage interval, here we shall use a more appropriate term – the 95 % state-of-knowledge interval. This should read as the interval of 95 % values that could be reasonably attributed to the unknown value of measurand based on the current state-of-knowledge (i.e. based on the measurement model, the currently available information, and method used for combining the information). Of course, further study is necessary for characterizing the optimality properties of the used method, e.g., under repeatability conditions. A standard approach to derive the state-of-knowledge PDF is based on the propagation of distributions using a Monte Carlo method, as suggested in Supplement 1 of the GUM, [2]. For more details and discussion on applicability of the uncertainty evaluation methods based on the GUM and its Supplement 1 see, e.g., [9]-[12]. A principal advantage of the Monte Carlo methods is their simplicity and asymptotic consistency (convergence to the true values with growing number of the Monte Carlo simulations). A disadvantage of the Monte Carlo methods is their principal ambiguity (i.e. non-uniqueness/variability of the results under independently repeated computations with given number of simulations) and typically large demand on the computational resources. Often, to achieve the pre-selected accuracy level, a need of a very large number of simulations is required. Alternatively, in linear measurement models, the state-of- knowledge PDF of the output quantity can be derived by inverting its characteristic function, which can be easily derived from the known characteristic functions of the input quantities. A principal advantage of the methods based on inversion of the characteristic function is their principal exactness (as there is a theoretical one-to-one mapping between the CF and the CDF) and efficiency for all values of the computed CDF/PDF/QF. However, in real situations, the precision of the used computational methods is also subject to possible numerical errors. A typical numerical error of such results is due to particular algorithmic implementation (as e.g., the truncation error and/or the integration error), which can be and should be properly controlled. As we shall illustrate below, in typical metrological applications with linear measurement models and the well behaved probability distributions of the input quantities, the numerical precision and efficiency of the here presented simple methods and algorithms is superior to the standard Monte Carlo methods. We notice that among other possible alternative approaches to evaluate the propagated probability distribution of the output quantity we can include the advanced methods for arithmetic computations with random variables and their distributions, see e.g. [13], [14], and also [15], [16]. However, applicability of these methods is still limited to a relatively small number of the input random variables. 2. LINEAR MEASUREMENT MODEL AND THE CHARACTERISTIC FUNCTIONS As mentioned above, an alternative tool to form the state- of-knowledge probability distribution of the output quantity in a linear measurement model is based on the numerical inversion of its characteristic function (CF), which is defined as a Fourier transform of its PDF, see (2). Computing the (inverse) Fourier transform numerically is a well-known problem, frequently connected with the problem of computing integrals of highly oscillatory (complex) functions. The problem was studied for a long time in general, but also with focus on specific applications, see, e.g., [17]-[22], to show just a few. In particular, the methods suggested for inverting the characteristic function for obtaining the probability distribution function include [23]-[27]. Approximations of the continuous Fourier transform by the discrete Fourier transform and by using the FFT algorithm are widely used in different fields of engineering. However, using the FFT for evaluation of the PDF/CDF from the characteristic function is not widespread in statistical applications (one important exception is the field of financial mathematics and econometrics), and in general, not well implemented in relevant software packages. In [28], Korczynski, Cox, and Harris suggested and illustrated the use of convolution principles in metrology applications. Their suggested approach was based on consecutive replacing of the convolution integrals by the convolution sums evaluated by using the fast Fourier transform (i.e. without direct using the exact characteristic functions), to ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 34 form the probability distribution for the output quantity in the measurement model of additive, linear or generalized linear form. If compared with the here proposed approach based on combining and inverting the exact characteristic functions, the numerical precision of their suggested approach can quickly deteriorate with the growing number of the required convolution integrals. In fact, in metrological applications a number of measurement models used in uncertainty evaluation are, at least approximately (up to reasonable level), of the additive linear form 𝑌 = 𝑐1𝑋1 + ⋯ + 𝑐𝑛𝑋𝑛, (1) where the input quantities 𝑋1, … , 𝑋𝑛 are independent random variables with known probability distributions, 𝑋𝑗 ∼ 𝐹𝑋𝑗, for 𝑗 = 1, … , 𝑛, possibly parametrized by 𝜃𝑗. Here, 𝑐1, … , 𝑐𝑛 denote the known constants and 𝑌 represents the univariate output quantity (a random variable with an unknown distribution to be determined). The characteristic function of a continuous univariate random variable 𝑋 ∼ 𝐹𝑋, with its probability density function pdf𝑋(𝑥), is defined as a Fourier transform of its PDF, cf𝑋(𝑡) = ∫ 𝑒𝑖𝑖𝑖 ∞ −∞ pdf𝑋(𝑥) 𝑑𝑥, 𝑡 ∈ 𝑹 (2) Analytical expressions of the characteristic functions are known for many standard probability distributions, see e.g. [29]-[31], or other available sources. Otherwise, CF could be derived either analytically, expressed by using computer-based tools, e.g. by using MATHEMATICA, or evaluated numerically. In Table 1 we present some selected characteristic functions of the univariate distributions, frequently used in metrological applications. Compare the presented distributions with those in Table 1 in [2]. Notice that the characteristic functions of the symmetric zero-mean distributions are purely real functions of the argument 𝑡 ∈ 𝑹. Deriving CF of a weighted sum of independent random variable is a simple and trivial task. Let cf𝑋𝑗(𝑡) denote the characteristic function of 𝑋𝑗. The characteristic function of 𝑌 defined by (1) is cf𝑌(𝑡) = cf𝑋1(𝑐1𝑡) ⋯ cf𝑋𝑛(𝑐𝑛𝑡). (3) For illustration, in Figure 1 we plotted the CF of a linear combination of two independent chi-squared random variables with 𝜈1 = 1 and 𝜈10 = 10 degrees of freedom, evaluated for 𝑡 ∈ (−1,1). Here we shall assume that the considered characteristic functions of the input and/or output quantities, the random variables 𝑋1, … , 𝑋𝑛 and 𝑌 are known or can be easily derived. Then, by the Fourier inversion theorem, the PDF of the random variable 𝑌 is given by pdf𝑌(𝑦) = 1 2𝜋 ∫ 𝑒−𝑖𝑖𝑖 ∞ −∞ cf𝑌(𝑡) 𝑑𝑡, 𝑦 ∈ 𝑹. (4) Analytical derivation of the PDF by using the (inverse) Fourier transform (4) is available only in special cases. Thus, in most practical situations, a numerical derivation of the PDF/CDF from the CF is an indispensable tool. In the next section we present two simple but frequently very efficient approaches (approximate numerical methods) for inversion of the characteristic function together with a detailed description of their algorithmic implementations, which are suitable for typical metrological applications. 3. NUMERICAL INVERSION OF THE CHARACTERISTIC FUNCTION The exact inverse Fourier transform (4) can be naturally approximated by the truncated integral form, pdf𝑌(𝑦) = 1 2𝜋 ∫ 𝑒−𝑖𝑖𝑖 𝑇 −𝑇 cf𝑌(𝑡) 𝑑𝑡, (5) Table 1. Characteristic functions of continuous univariate distributions used in metrological applications (selected symmetric zero-mean distributions and non-negative distributions). Here, Kν(z) denotes the modified Bessel function of the second kind, Jν(z) is the Bessel function of the first kind, and U(a, b, z) is the confluent hypergeometric function of the second kind. Probability distribution Characteristic function (CF) Gaussian 𝑁(0,1) cf(𝑡) = 𝑒− 1 2 𝑖2 Student’s t 𝑡𝜈 cf(𝑡) = 1 2 𝜈 2−1Γ�𝜈 2 � �𝜈 1 2 |𝑡|� 𝜈 2 𝐾𝜈 2 �𝜈 1 2 |𝑡|� Rectangular 𝑅(−1,1) cf(𝑡) = sin (𝑖) 𝑖 Triangular 𝑇(−1,1) cf(𝑡) = 2−2cos (𝑖) 𝑖2 Arcsine 𝑈(−1,1) cf(𝑡) = 𝐽0(𝑡) Exponential 𝐸𝑥𝐸(𝜆) cf(𝑡) = 𝜆 𝜆−𝑖𝑖 𝜆 > 0 rate Gamma Γ(𝛼, 𝛽) cf(𝑡) = �1 − 𝑖𝑖 𝛽 � −𝛼 𝛼 > 0 shape, 𝛽 > 0 rate Chi-squared 𝜒𝜈2 cf(𝑡) = (1 − 2𝑖𝑡)− 𝜈 2 𝜈 > 0 degrees of freedom Fisher-Snedecor’s F 𝐹𝜈1,𝜈2 cf(𝑡) = Γ�ν1 2 +ν2 2 � Γ�𝜈2 2 � 𝑈 �𝜈1 2 , 1 − 𝜈2 2 , − 𝜈2 𝜈1 𝑖𝑡 � 𝜈1 > 0, 𝜈2 > 0 degrees of freedom Figure 1. Real (blue) and imaginary (red) part of the characteristic function of 𝑌 = 10 𝑋𝜒12 + 𝑋𝜒10 2 — the linear combination of two independent chi- squared random variables with 𝜈1 = 1 and 𝜈1 = 10 degrees of freedom, evaluated for 𝑡 ∈ (−1,1). ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 35 where 𝑇 is a sufficiently large (real) value, and the integrand is a complex (oscillatory) function, which is assumed to decay to zero reasonably fast for |𝑡| → ∞ (as is typical for continuous distributions used in metrological applications). Note that the selection of the integration limit 𝑇 and the speed of integrand function decay contribute significantly to the truncation error of this approximation. In general, the required integral can be evaluated by any suitable numerical quadrature method. Frequently, a simple trapezoidal rule gives fast and satisfactory results. However, caution is necessary if the integrand is a highly oscillatory function, what is the case if abs(𝑦) is a large value from the tail area of the distribution, as e.g. the 99 % quantile. In such situations a more advanced quadrature method should be used, typically in combination with efficient root-finding algorithms and algorithms for accelerated computing of the limit values of the sums of alternating series (not considered here), see e.g., [21], [32]. Note that the selection of the integration method significantly contributes to the integration error of this approximation. Figure 2 illustrates the problem of integrating the highly oscillatory integrand function. For typical metrological applications we suggest to consider algorithms based on the Gil-Pelaez inversion formulae (6)-(7) and the approximate discrete Fourier transform (21) resp. (22), based on the simple trapezoidal rule for computing the required sub-integrals. The possible numerical error of the results (i.e. the truncation error and the integration error) can be controlled by selecting a proper truncation limit 𝑇 (sufficiently large, leading to a small truncation error), and by dividing the selected integration interval (0, 𝑇) to a large number, say 𝑁, of small sub-intervals, where the trapezoidal rule provides a satisfactorily good numerical approximation (with small integration error) of the true integral value. For theoretical results on how to control the truncation and integration errors in such and similar situations see, e.g., [35]. 3.1. The Gil-Pelaez inversion formulae In [33], Gil-Pelaez derived the inversion formulae, suitable for numerical evaluation of the PDF and/or the CDF, which require integration of a real-valued function, only. The PDF is given by pdf𝑌(𝑦) = 1 𝜋 ∫ ℜ�𝑒−𝑖𝑖𝑖cf𝑌(𝑡)� ∞ 0 𝑑𝑡. (6) Further, if 𝑦 is a continuity point of the cumulative distribution function of 𝑌, the CDF is given by cdf𝑌(𝑦) = 1 2 − 1 𝜋 ∫ ℑ� 𝑒−𝑖𝑖𝑖cf𝑌(𝑖) 𝑖 � ∞0 𝑑𝑡. (7) By ℜ(𝑓(𝑡)) and ℑ(𝑓(𝑡)) we denote the real and imaginary part of the complex function 𝑓(𝑡), respectively. Numerical inversion of the characteristic function based on (6) and (7) have been successfully implemented for evaluation of the distribution function of a linear combination of independent chi-squared RVs by Imhof in [34] and by Davies in [35]. Further, Gil-Pelaez's method has been implemented in the algorithm tdist, see [36] and [37], for computing the distribution of a linear combination of independent Student's t random variables and/or other symmetric zero-mean random variables, and also for computing the distribution of a linear combination of independent inverted gamma random variables suggested in [38], and the distribution of a linear combination of independent log-Lambert 𝑊 × 𝜒𝜈2 RVs, [39]. In [40], the algorithm tdist has been suggested and applied for computing the 95 % state-of-knowledge interval (considered as the approximate 95 % confidence interval) for the common mean value in the inter-laboratory comparisons with systematic effects (biases). In general, the integrals in (6) and (7) can be computed by any numerical quadrature method, possibly in combination with efficient root-finding algorithms and accelerated computing of limits of the alternating series, as considered, e.g., in [20] and [25]. Frequently, (6) and (7) can be efficiently approximated by a simple trapezoidal quadrature: pdf𝑌(𝑦) ≈ 𝛿𝑖 𝜋 �𝑤𝑗 𝑁 𝑗=0 ℜ�𝑒−𝑖𝑖𝑗𝑖cf𝑌�𝑡𝑗�� ≈ 𝛿𝑖 𝜋 � 𝑤0 + ∑ 𝑤𝑗𝑁𝑗=1 cos�𝑡𝑗𝑦�ℜ�cf𝑌�𝑡𝑗�� + ∑ 𝑤𝑗𝑁𝑗=1 sin�𝑡𝑗𝑦�ℑ�cf𝑌�𝑡𝑗�� � (8) cdf𝑌(𝑦) ≈ 1 2 − 𝛿𝑖 𝜋 �𝑤𝑗 𝑁 𝑗=0 ℑ� 𝑒−𝑖𝑖𝑗𝑖cf𝑌�𝑡𝑗� 𝑡𝑗 � ≈ 1 2 − 𝛿𝑖 𝜋 ⎝ ⎜ ⎛ 𝑤0(mean(𝑌) − 𝑦) + + ∑ 𝑤𝑗𝑁𝑗=1 cos�𝑡𝑗𝑦�ℑ� cf𝑌�𝑖𝑗� 𝑖𝑗 � + ∑ 𝑤𝑗𝑁𝑗=1 sin�𝑡𝑗𝑦�ℜ� cf𝑌�𝑖𝑗� 𝑖𝑗 � ⎠ ⎟ ⎞ , (9) where 𝑁 is a sufficiently large number of (equidistant) sub- intervals of (0, 𝑇), 𝑤𝑗 are the appropriate quadrature weights, and 𝑡𝑗 denote the appropriate (equidistant) nodes from the interval (0, 𝑇), for sufficiently large 𝑇. In particular, for the trapezoidal quadrature rule we set • 𝛿𝑖 = 𝑇 𝑁 or alternatively 𝛿𝑖 = 2𝜋 𝐵−𝐴 , which gives 𝑇 = 𝑁 2𝜋 𝐵−𝐴 , • 𝑤0 = 𝑤𝑁 = 1 2 , and 𝑤𝑗 = 1 for 𝑗 = 1, … , 𝑁 − 1, • 𝑡𝑗 = 𝑗𝛿𝑖 for 𝑗 = 0, … , 𝑁, with 𝑇 = 𝑡𝑁 = 𝑁𝛿𝑖. Here, the interval (𝐴, 𝐵) specifies the range of typical values 𝑦, i.e. a large part of the distribution support of the random variable 𝑌. Figure 2. Integrand functions for computing the PDF/CDF of the chi-squared distributed random variable, 𝑌 ∼ 𝜒12, at 𝑦 = 15, computed by the Gil-Pelaez inversion formulae of its characteristic function cfY(𝑡) = (1 − 2𝑖𝑡) −1 2. The plotted integrand functions of the integrals in (6) and (7) are evaluated for 𝑡 ∈ (0, 𝑇), 𝑇 = 5. The red circles depict the zeros (roots) of the integrand functions on (0, 𝑇). The total integral can be computed as an infinite sum of sub-integrals - the limit value of the sum of alternating series. http://www.mathworks.com/matlabcentral/fileexchange/4199-tdist http://www.mathworks.com/matlabcentral/fileexchange/4199-tdist ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 36 As a simple rule of thumb, if the (optimum) value of 𝑇 is unknown we suggest to start with the application of the six- sigma-rule, i.e. set the typical range (𝐴, 𝐵) as an intersection of the natural parametric space of 𝑌 with the interval (𝐿, 𝑈) (e.g., (𝐴, 𝐵) = (𝐿, 𝑈) ∩ 𝑹 or (𝐴, 𝐵) = (𝐿, 𝑈) ∩ 𝑹+, with • 𝐿 = mean(𝑌) − 6 std(𝑌), • 𝑈 = mean(𝑌) + 6 std(𝑌), where mean(𝑌) and std(𝑌) represent the expectation and the standard deviation of the probability distribution of 𝑌, and define 𝑇 = 𝑁 2𝜋 𝐵−𝐴 for some pre-selected fixed 𝑁. Note that the value of 𝑇 should be kept sufficiently large, such that the characteristic function is sufficiently small for all 𝑡 > 𝑇, i.e. 𝑎𝑎𝑎(cf𝑌(𝑡)) < 𝜖 for small 𝜖, say 𝜖 = 10−12. This can be effectively controlled by using a proper value (sufficiently large) of 𝑁. Further, for computing the leading term in (9), we use the result from [36]: If the mean (expectation) of 𝑌 exists, then lim𝑖→0 ℑ� 𝑒−𝑖𝑖𝑖cf𝑌(𝑖) 𝑖 � = mean(𝑌) − 𝑦 . (10) The required mean(𝑌) and std(𝑌) can be evaluated analytically, from the moments of the input variables, or approximately, by using numerical differentiation of the characteristic function of 𝑌, cf𝑌(𝑡). In particular, mean(𝑌) ≈ 1 12𝑖ℎ � cf𝑌(−2ℎ) − 8cf𝑌(−ℎ) +8cf𝑌(−ℎ) − cf𝑌(2ℎ) � , (11) std(𝑌) ≈ �m2(𝑌) − mean2(𝑌)� , (12) where m2(𝑌) ≈ 1 144ℎ2 ⎝ ⎜ ⎛ cf𝑌(−4ℎ) − 16cf𝑌(−3ℎ) +64cf𝑌(−2ℎ) + 16cf𝑌(−ℎ) −130 +16cf𝑌(ℎ) + 64cf𝑌(2ℎ) −16cf𝑌(3ℎ) + cf𝑌(4ℎ) ⎠ ⎟ ⎞ , (13) for any small ℎ > 0, e.g., ℎ = 10−4. This numerical approach is applicable even in the cases when the theoretical moments of the considered distribution, defined by cf𝑌(𝑡), do not exist (e.g., for the Student’s t distribution with 1 or 2 degrees of freedom). The truncation limit 𝑇 = 𝑁𝛿𝑖 = 𝑁 2𝜋 𝐵−𝐴 depends on 𝑁 and (𝐴, 𝐵). The trade-off between the values of 𝐵 − 𝐴, 𝑁, 𝛿𝑖 and 𝑇, strongly depends on the particular distribution of 𝑌 and its CF. For example, for a highly precise numerical inversion of the standard normal CF, cfY(𝑡) = 𝑒 −1 2 𝑖2, computed in double precision arithmetic with precision better than 𝜖 = 10−14, it is sufficient to set (𝐴, 𝐵) = (−8,8) with 𝐵 − 𝐴 = 16, 𝑁 = 25 = 32, leading to 𝛿𝑖 = 𝜋 8 and 𝑇 = 4 𝜋 with cfY(𝑇) = 𝑒 −1 2 𝑇2 ≈ 5 × 10−35. On the other hand, the numerical inversion of the rectangular CF given by cfY(𝑡) = sin (𝑖) 𝑖 requires in double precision arithmetic the value 𝑇 = 1014 to get 𝑎𝑎𝑎�cf(𝑡)� ≤ 10−14 for 𝑡 > 𝑇. For the natural choice of (𝐴, 𝐵) = (−1,1) with 𝐵 − 𝐴 = 2, this suggests to set 𝑁 ≈ 1013, what is an unacceptable value, and thus it reveals that the simple trapezoidal rule is not a suitable integration method for the highly precise numerical inversion of the rectangular distribution CF. Fortunately, CF of the output quantity 𝑌 in typical metrological situations based on linear measurement models is a well behaved function, and the methods based on simple trapezoidal quadrature are efficient to provide good numerical approximations of the true values. In general, the pre-selected values of 𝑇 and 𝑁 should be checked and/or properly corrected. A simple diagnostic check is based on evaluating abs(cf𝑌(𝑇)). For large 𝑇 this should be a small value (smaller that the accepted truncation error). This check allows to control also the level of 𝑁, given the already fixed and sufficiently large interval (𝐴, 𝐵). We note that the presented quadrature method requires only one evaluation of the characteristic function cf𝑌�𝑡𝑗� at 𝑡𝑗, 𝑗 = 1, … , 𝑁, for any 𝑦 ∈ (𝐴, 𝐵) in pdf𝑌(𝑦) and cdf𝑌(𝑦), respectively. Moreover, the computation is further simplified if 𝑌 is a continuous random variable with a symmetric zero-mean distribution, i.e. with purely real CF, pdf𝑌(𝑦) ≈ 𝛿𝑖 𝜋 � 1 2 + ∑ cos�𝑡𝑗𝑦�cf𝑌�𝑡𝑗�𝑁−1𝑗=1 + 1 2 cos(𝑡𝑁𝑦)cf𝑌(𝑡𝑁) � , (14) cdf𝑌(𝑦) ≈ 1 2 − 𝛿𝑖 𝜋 � − 𝑖 2 + ∑ sin�𝑡𝑗𝑦� cf𝑌�𝑖𝑗� 𝑖𝑗 𝑁−1𝑗=1 + 1 2 sin(𝑡𝑁𝑦) cf𝑌(𝑖𝑁) 𝑖𝑁 � . (15) Finally, the quantile function (QF) can be evaluated by using the iterative Newton-Raphson scheme, based on repeated evaluations of the PDF/CDF, see (8)-(9) and/or (14)-(15). In particular, for a fixed probability level 𝐸 ∈ (0,1), the 𝐸- quantile of the (continuous) distribution of 𝑌, say 𝑞 = qf𝑌(𝐸), is given as a solution (fixed point) of the following iterative scheme, qf𝑌𝑘+1(𝐸) = qf𝑌𝑘(𝐸) − �cdf𝑌�qf𝑌 𝑘(𝑝)�−𝑝� pdf𝑌�qf𝑌 𝑘(𝑝)� , (16) where 𝑘 = 0,1, …, and the starting value qf𝑌0(𝐸) is set as, e.g., qf𝑌0(𝐸) = mean(𝑌), defined by (11). %% EXAMPLE (MATLAB ALGORITHM CF2DISTGP) % % PDF and CDF of a linear combination of RVs: % Y = c1*X1 + c2*X2 + c3*X3 + c4*X4 + c5*X5, % where, % X1 ~ Normal(0,1) with c1=1, % X2 ~ Student's t with 1 df and c2=1, % X3 ~ Rectangular on (-1,1) with c3=5, % X4 ~ Triangular on (-1,1) with c4=1, % X5 ~ U-distribution on (-1,1) with c5=10 cfN = @(t) exp(-t.^2/2); cft = @(t,nu) min(1,besselk(nu/2, ... abs(t).*sqrt(nu),1) .* ... exp(-abs(t).*sqrt(nu)) .* ... (sqrt(nu).*abs(t)).^(nu/2) / ... 2^(nu/2-1)/gamma(nu/2)); cfR = @(t) min(1,sin(t)./t); cfT = @(t) min(1,(2-2*cos(t))./t.^2); cfU = @(t) besselj(0,t); c = [1 1 5 1 10]; nu = 1; cfY = @(t) ... cfN(c(1)*t) .* ... cft(c(2)*t,nu) .* ... cfR(c(3)*t) .* ... cfT(c(4)*t) .* ... cfU(c(5)*t); y = linspace(-50,50,201)'; [result,cdf,pdf] = cf2DistGP(cfY,y); ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 37 A working version of the MATLAB algorithm cf2DistGP for computing the PDF/CDF by numerical inversion of the characteristic function, based on the Gil-Pelaez inversion formulae, is presented in Appendix A. For illustration, the MATLAB code presented above evaluates the PDF and CDF of the output variable 𝑌, which is a linear combination of the independent random variables with a normal, Student's t, rectangular, triangular and arcsine distributions, i.e. 𝑌 = 𝑋𝑁 + 𝑋𝑖𝜈 + 5𝑋𝑅 + 𝑋𝑇 + 10𝑋𝑈, by using the algorithm cf2DistGP with its default settings of 𝑇 and 𝑁. Similarly, the PDF/CDF of the random variable 𝑌 = 𝑋𝑁 + 𝑋𝑖𝜈 + 5𝑋𝑅 + 𝑋𝑇 + 10𝑋𝑈 can be evaluated by using the MATLAB algorithm tdist, see also Figure 3: 3.2. Numerical inversion of the characteristic function by using the FFT algorithm This approach for computing the PDF by numerical inversion of the characteristic function by using the FFT algorithm is based on the results by Hürlimann in [41]. Alternatively, for other applications based on using the fractional fast Fourier transform (FRFT), see [42] and also [43]- [46]. Here we shall approximate the continuous Fourier transform (CFT), say 𝐹(𝑦) = ∫ 𝑒−𝑖2𝜋𝜋𝑖 ∞ −∞ 𝑓(𝑢) 𝑑𝑢, (17) by a discrete Fourier transform (DFT). The DFT can be efficiently evaluated by using the FFT algorithm that computes the same result as the DFT, but much faster. For complex numbers 𝑓0, … , 𝑓𝑁−1 the DFT is defined as 𝐹𝑘 = ∑ 𝑒 −𝑖2𝜋𝑘𝑗𝑁 𝑓𝑗𝑁−1𝑗=0 , 𝑘 = 0, … , 𝑁 − 1. (18) Formally, here we shall use the following notation, 𝐅𝑁 = 𝐹𝐹𝑇(𝐟𝑁), (19) where 𝐟𝑁 = (𝑓0, … , 𝑓𝑁−1) and 𝐅𝑁 = (𝐹0, … , 𝐹𝑁−1). The relationship between the CF and the PDF is given by the (inverse) continuous Fourier transform defined by (4). For a sufficiently large interval (−𝑇, 𝑇), it is possible to approximate a PDF by (5). Here we shall consider the integral approximation, based on the mid-point integration rule, ∫ 𝑓(𝑥)𝑑𝑥 ≈ 𝑓(𝑎)+𝑓(𝑏) 2 (𝑎 − 𝑎)𝑏𝑎 , which corresponds to the trapezoidal quadrature. For more alternative approaches based on more sophisticated integration rules see, e.g., [39]. Similarly as before, let (𝐴, 𝐵) denote a sufficiently large interval, where the distribution of 𝑌 is concentrated. A reasonable rule for determining (𝐴, 𝐵) can be based, for example, on the six-sigma-rule, or its modifications, by using an altered multiplication coefficient (e.g., use 10 instead of 6). Let 𝑦𝑘 = 𝐴 + 𝑘𝛿𝑖, with 𝛿𝑖 = 𝐵−𝐴 𝑁 , and 𝑘 = 0, … , 𝑁 − 1. For 𝑁 large, 𝑇 = 𝜋/𝛿𝑖 is also large, and from (5), by using 𝑡 = 2𝜋𝑢, 𝑑𝑡 = 2𝜋𝑑𝑢, and 𝑑𝑢 = 1 𝐵−𝐴 , we get pdf𝑌(𝑦𝑘) ≈ 1 2𝜋 ∫ 𝑒−𝑖2𝜋𝜋𝑖𝑘 1 2𝛿𝑖 − 12𝛿𝑖 cf𝑌(2𝜋𝑢) 𝑑𝑢 . (20) Now, we shall approximate the integral (19) by using the approximate trapezoidal quadrature rule, for each of the 𝑁 sub- intervals. Thus, pdf𝑌(𝑦𝑘) ≈ 1 𝐵−𝐴 ∑ 𝑒−𝑖2𝜋𝜋𝑗𝑖𝑘𝑁−1𝑗=0 cf𝑌�2𝜋𝑢𝑗� , (21) Figure 3. The probability density function (PDF) and the cumulative distribution function (CDF) of a random variable 𝑌 = ∑ 𝑐𝑗𝑋𝑗 5 𝑗=1 , with 𝑋1 ∼ 𝑁(0,1), 𝑋2 ∼ 𝑡𝜈=1, 𝑋3 ∼ 𝑅(−1,1), 𝑋4 ∼ 𝑇(−1,1), 𝑋5 ∼ 𝑈(−1,1), and coefficients 𝑐 = (𝑐1, … , 𝑐5) = (1,1,5,1,10), evaluated by numerical inversion of its characteristic function by the MATLAB algorithm tdist, see also the Examples. %% EXAMPLE (MATLAB ALGORITHM TDIST) % % TDIST at Matlab Central File Exchange: % http://www.mathworks.com/matlabcentral/ % /fileexchange/4199-tdist % % PDF and CDF of a linear combination of RVs % Y = c1*X1 + c2*X2 + c3*X3 + c4*X4 + c5*X5 % with: % X1 ~ Normal(0,1) [we set df1=Inf] with c1=1, % X2 ~ Student's t with 1 df [set df2=1], c2=1, % X3 ~ Rectangular on (-1,1) [set df3=-1],c3=5, % X4 ~ Triangular on (-1,1) [set df4=-2], c4=1, % X5 ~ U-distribution on (-1,1) [df5=-3], c5=10 df = [Inf 1 -1 -2 -3]; coefs = [1 1 5 1 10]; [pdf,y] = tdist([],df,coefs,'PDF'); cdf = tdist(y,df,coefs,'CDF'); figure; plot(y,pdf); grid figure; plot(y,cdf); grid http://www.mathworks.com/matlabcentral/fileexchange/4199-tdist http://www.mathworks.com/matlabcentral/fileexchange/4199-tdist ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 38 where 𝑢𝑗 = 1 2+𝑗− 𝑁 2 𝐵−𝐴 , 𝑗 = 0, … , 𝑁 − 1. From that, by using 𝑒𝑖𝜋 = −1, the expressions for 𝑢𝑗 and 𝑦𝑘, and the DFT defined by (17), we finally get the formal relationship 𝐩𝐩𝐟 = 𝐂 ⊙ 𝐹𝐹𝑇(𝐃 ⊙ 𝐜𝐟), (22) where ⊙ denotes the dot product (element wise multiplication), • 𝐩𝐩𝐟 = �pdf𝑌(𝑦0), … , pdf𝑌(𝑦𝑁−1)�, • 𝐂 = (𝐶0, … , 𝐶𝑁−1), with • 𝐶𝑘 = 1 𝐵−𝐴 (−1)�� 1−1 𝑁�� 𝑁𝑁 𝐵−𝑁+𝑘�� , 𝑘 = 0, … , 𝑁 − 1, • 𝐃 = (𝐷0, … , 𝐷𝑁−1), with • 𝐷𝑘 = (−1) − 2𝑁𝐵−𝑁𝑘, 𝑘 = 0, … , 𝑁 − 1, • 𝐜𝐟 = �cf𝑌(𝑡0), … , cf𝑌(𝑡𝑁−1)�, with • 𝑡𝑘 = 2𝜋 𝐵−𝐴 �1 2 + 𝑘 − 𝑁 2 �, 𝑘 = 0, … , 𝑁 − 1. Further, CDF is evaluated here by simple cumulative sum from the evaluated PDF values, and QF is evaluated by interpolation from the CDF. A working version of the MATLAB algorithm cf2DistFFT for computing the PDF/CDF/QF by numerical inversion of the characteristic function, based on the FFT algorithm, is presented in Appendix B. For illustration, the MATLAB code presented below evaluates the PDF/CDF of the output variable 𝑌 = 10 𝑋𝜒12 + 𝑋𝜒102 , by using the algorithm cf2DistFFT with its default settings of (𝐴, 𝐵) and 𝑁. Other specific versions of the algorithm for computing the PDF/CDF/QF of a linear combination of independent random variables with the Fisher-Snedecor’s F-distributions and the log-normal distributions, by numerical inversion of the characteristic function by using the FFT algorithm, are available at the MATLAB CENTRAL FILE EXCHANGE as the algorithm Fdist, file ID: 56262, and the algorithm logNdist, file ID: 56512, respectively. 4. COMPARISON OF THE CF APPROACH AND THE MONTE CARLO METHOD The main advantage of the Monte Carlo methods lies in their simplicity and asymptotic consistency. A disadvantage of the Monte Carlo methods lies in the large computational demands, typically required in order to achieve the pre-specified accuracy level. On the other hand, application of the CF approach offers principal theoretical exactness which could be, however, influenced by unacceptable numerical errors (i.e. the truncation and the integration errors), if not properly used. For illustration, here we consider the linear measurement model for calibration of a coaxial step attenuator as considered in [47], a model typical for metrological applications. We consider this example in order to compare results based on the CF approach and the Monte Carlo method and to illustrate the applicability and/or advantage of the suggested methods for potential users in metrology applications. For more details about this specific example see [47], example S7, and/or [48]. The linear measurement model of the attenuation 𝐿𝑋 of the attenuator to be calibrated is given by 𝐿𝑋 = 𝑐𝑐𝑛𝑎𝑡 + 𝐿𝑆 + 𝛿𝐿𝑆 + 𝛿𝐿𝐷 + 𝛿𝐿𝑀 + 𝛿𝐿𝐾 + 𝛿𝐿𝑖𝑏 − 𝛿𝐿𝑖𝑎 + 𝛿𝐿0𝑏 − 𝛿𝐿0𝑎, (23) with the following information about the distributions of the input quantities (available from the given uncertainty budget, see S7.12 in [47] with the correction as presented in [48]): • 𝑐𝑐𝑛𝑎𝑡 = 30.04 + 0.003 = 30.043, • 𝐿𝑆 ∼ 0.0090 × 𝑁(0,1), • 𝛿𝐿𝑆 ∼ (0.0025/� 1 3 ) × 𝑅(−1,1), • 𝛿𝐿𝐷 ∼ (0.0011/� 1 2 ) × 𝑈(−1,1), • 𝛿𝐿𝑀 ∼ (0.0200/� 1 2 ) × 𝑈(−1,1), • 𝛿𝐿𝐾 ∼ (0.0017/� 1 2 ) × 𝑈(−1,1), • 𝛿𝐿𝑖𝑏 ∼ (0.0003/� 1 3 ) × 𝑅(−1,1) • 𝛿𝐿𝑖𝑎 ∼ (0.0003/� 1 3 ) × 𝑅(−1,1), • 𝛿𝐿0𝑏 ∼ 0.0020 × 𝑁(0,1), • 𝛿𝐿0𝑎 ∼ 0.0020 × 𝑁(0,1). The algorithms cf2DistGP, with default option parameters based on the six-sigma-rule and pre-selected 𝑁 = 210 = 1024, evaluated the support interval (𝐴, 𝐵) = (−0.1341,0.1341), and hence 𝑇 = 23989, with 𝛿𝑖 = 𝑇 𝑁 = 23.4. Then, the 97.5 %-quantile of 𝑌 = 𝐿𝑋 − 𝑐𝑐𝑛𝑎𝑡 was calculated as 𝑞 = 0.03900448275179. Thus, the calculated 95 % state-of-knowledge interval for the attenuation 𝐿𝑋 is given as 30.043 ∓ 0.039. The basic diagnostic check ensures the highest precision of the presented calculations in double precision arithmetic, as abs(cf𝑌(𝑇))=0 (i.e. the value equals to the machine zero, which is here defined as 4.94 × 10−324). The used computer time was 𝑡 = 3.9 × 10−4 s. Further, here we also present the estimated values of the required 97.5 % quantile computed by the Monte Carlo method for sample sizes 𝑀, with 𝑀 = 104, 105, 106, 107, and 108. %% EXAMPLE (MATLAB ALGORITHM CF2DISTFFT) % % Distribution of a linear combination of RVs % (chi-squared RVs with 1 and 10 DFs) % Y = 10*X_{\chi^2_1} + X_{\chi^2_10} df1 = 1; df2 = 10; cfChi2_1 = @(t) (1-2i*t).^(-df1/2); cfChi2_10 = @(t) (1-2i*t).^(-df2/2); cfY = @(t) cfChi2_1(10*t) .* cfChi2_10(t); clear options options.isPositiveSupport = true; result = cf2DistFFT(cfY,[],[],options); % PLOT THE CF of Y t = linspace(-1,1,501); figure plot(t,real(cfY(t)),t,imag(cfY(t)));grid xlabel('t'); ylabel('Characteristic function'); title('Y = 10*X_{\chi^2_1}+X_{\chi^2_{10}}') http://www.mathworks.com/matlabcentral/fileexchange/56262-fdist http://www.mathworks.com/matlabcentral/fileexchange/56512-logndist ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 39 The estimated values of the 97.5 % quantile 𝑞, and the computer time 𝑡 (in seconds), used for the sample size 𝑀: • 𝑀 = 104, 𝑞 = 0.039343812021804, 𝑡 = 0.07, • 𝑀 = 105, 𝑞 = 0.038985244853524, 𝑡 = 0.08, • 𝑀 = 106, 𝑞 = 0.038952503603421, 𝑡 = 0.75, • 𝑀 = 107, 𝑞 = 0.038998144396413, 𝑡 = 7.49, • 𝑀 = 108, 𝑞 = 0.039003626106614, 𝑡 = 153.2. This clearly illustrates the advantages and computational efficiency of the proposed approach over the standard Monte Carlo methods, especially if high precision of the computed (estimated) quantiles is required. 5. CONCLUSIONS We suggest to consider numerical methods for derivation of the state-of-knowledge PDF/CDF of the output quantity in linear measurement models from its characteristic function. Such approach can be used to form the probability distribution for the output quantity of a measurement model of additive, linear or generalized linear form, and can be considered as an alternative tool to the uncertainty evaluation based on the Monte Carlo methods. Here we have presented two simple but efficient approaches for numerical inversion of the characteristic function, which are especially suitable for metrological applications. The suggested numerical approaches are based on the Gil- Pelaez inverse formula and on the approximation by discrete Fourier transform (DFT) and the FFT algorithm for computing the PDF/CDF of (univariate) continuous random variables. As we have explained and illustrated, the suggested CF approach should be considered as an alternative to the standard approach based on the Monte Carlo methods (as considered e.g. in Supplement 1 – Propagation of distributions using a Monte Carlo method [2]) in specific situations, i.e. for evaluating the state-of- knowledge probability distributions of the output quantity in a linear measurement model, especially if the highest precision of the reported distribution quantiles is required. However, numerical errors of such results should be properly controlled. For illustration purposes, here we present working versions of the MATLAB codes of the suggested algorithms (cf2DistGP and cf2DistFFT), as well as some simple examples in order to illustrate the applicability of the suggested methods. Although the methods for inverting characteristic functions for obtaining the probability distribution functions have been studied for long time, especially in statistical literature, and the possible applications are much more general than those motivated by metrology, surprisingly, such methods are still not widespread and used in engineering applications and only rarely used among statisticians. One possible reason might be that the characteristic functions and the algorithms for numerical inversions are not directly available in standard (statistical) software packages, like e.g., R or MATLAB. Systematic development of the methods, algorithms and the software toolbox (developed for R and/or MATLAB) for computing, combining and inverting the characteristic functions is our highly desirable goal for the next future. ACKNOWLEDGEMENT The work was supported by the Slovak Research and Development Agency, project APVV-15-0295, and by the Scientific Grant Agency VEGA of the Ministry of Education of the Slovak Republic and the Slovak Academy of Sciences, by the projects VEGA 2/0047/15 and VEGA 2/0011/16. %% COMPARISON / CF APPROACH & MONTE CARLO METHOD % CF APPROACH: % Distribution and quantile of the output % quantity in linear measurement model % Y = L_S + δL_S + δL_D + δL_M + δL_K + % δL_ib - δL_ia + δL_0b - δL_0a % const = 30.043; prob = 0.975; c = [0.009 0.0025/sqrt(1/3) ... 0.0011/sqrt(1/2) 0.0200/sqrt(1/2) ... 0.0017/sqrt(1/2) 0.0003/sqrt(1/3) ... -0.0003/sqrt(1/3) 0.0020 -0.0020 ]; cfN = @(t) exp(-t.^2/2); cfR = @(t) min(1,sin(t)./t); cfU = @(t) besselj(0,t); cfY = @(t) cfN(c(1)*t) .* cfR(c(2)*t) .* ... cfU(c(3)*t) .* cfU(c(4)*t) .* ... cfU(c(5)*t) .* cfR(c(6)*t) .* ... cfR(c(7)*t) .* cfN(c(8)*t) .* ... cfN(c(9)*t); [result,cdf,pdf,q_CF] = cf2DistGP(cfY,[],prob) % 95% state-of-knowledge interval / CF APPROACH: I_CF = const + [-q_CF,q_CF] % MONTE CARLO METHOD: M = [10^4 10^5 10^6 10^7 10^8]; q_MC = zeros(5,1); time = zeros(5,1); rng; for m = 1:5 tic; N = randn(M(m),3); R = 2*rand(M(m),3)-1; U = 2*betarnd(0.5,0.5,M(m),3)-1; Y = c(1)*N(:,1) + c(2)*R(:,1) ... + c(3)*U(:,1) + ... + c(4)*U(:,2) + c(5)*U(:,3) ... + c(6)*R(:,2) + c(7)*R(:,3) ... + c(8)*N(:,2) + c(9)*N(:,3); Y = sort(Y); q_MC(m) = Y(ceil(M(m)*prob)); time(m) = toc; end % 95% state-of-knowledge interval / MC METHOD: q = q_MC(5); I_MC = const + [-q,q] ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 40 APPENDIX A. MATLAB ALGORITHM CF2DISTGP FOR NUMERICAL INVERSION OF THE CHARACTERISTIC FUNCTION BASED ON THE GIL-PELAEZ INVERSION FORMULAE function [result,cdf,pdf,qf] = ... cf2DistGP(cf,x,prob,options) % Calculates CDF/PDF/QF from the characteristic % function by the Gil-Pelaez inversion formulae, % integrated by a simple trapezoidal quadrature. % % SYNTAX: % [result,cdf,pdf,qf] = ... % cf2DistGP(cf,x,prob,options) %% CHECK THE INPUT PARAMETERS tic; narginchk(1, 4); if nargin < 4, options = []; end if nargin < 3, prob = []; end if nargin < 2, x = []; end if ~isfield(options, 'N') options.N = 2^10; end if ~isfield(options, 'T') options.T = []; end if ~isfield(options, 'SixSigmaRule') options.SixSigmaRule = 6; end if ~isfield(options, 'xMean') options.xMean = []; end if ~isfield(options, 'xStd') options.xStd = []; end if ~isfield(options, 'tolDiff') options.tolDiff = 1e-4; end if ~isfield(options, 'qf0') options.qf0 = (cf(1e-4)-cf(-1e-4))/(2e-4*1i); end if ~isfield(options, 'crit') options.crit = 1e-13; end if ~isfield(options, 'maxiter') options.maxiter = 1000; end if ~isfield(options, 'isPlot') options.isPlot = true; end if ~isfield(options, 'DIST') options.DIST = []; end if ~isempty(options.DIST) xMean = options.DIST.xMean; cft = options.DIST.cft; dt = options.DIST.dt; N = length(cft); t = (1:N)' * dt; range = 2*pi / dt; xMin = xMean - range/2; xMax = xMean + range/2; xStd = []; else N = options.N; T = options.T; SixSigmaRule = options.SixSigmaRule; xMean = options.xMean; xStd = options.xStd; h = options.tolDiff; cft = cf(h*(1:4)); if isempty(xMean) xMean = real((-cft(2) ... + 8*cft(1)-8*conj(cft(1)) ... + conj(cft(2)))/(1i*12*h)); end if isempty(xStd) xM2 = real(-(conj(cft(4)) ... - 16*conj(cft(3)) ... + 64*conj(cft(2)) ... + 16*conj(cft(1)) ... - 130 + 16*cft(1) ... + 64*cft(2) ... - 16*cft(3)+cft(4))/(144*h^2)); xStd = sqrt(xM2 - xMean^2); end if ~isempty(T) dt = T / N; t = (1:N)' * dt; cft = cf(t); cft(N) = cft(N)/2; range = 2*pi / dt; xMin = xMean - range/2; xMax = xMean + range/2; xStd = []; else xMin = xMean - SixSigmaRule * xStd; xMax = xMean + SixSigmaRule * xStd; range = xMax - xMin; dt = 2*pi / range; t = (1:N)' * dt; cft = cf(t); cft(N) = cft(N)/2; end options.DIST.xMean = xMean; options.DIST.cft = cft; options.DIST.dt = dt; end %% ALGORITHM if isempty(x) x = linspace(xMin,xMax,101)'; end if any(x < xMin) || any(x > xMax) warning(['x out-of-range (the support): ', ... '[xMin, xMax] = [',num2str(xMin),... ', ',num2str(xMax),'] !']); end [n,m] = size(x); x = x(:); E = exp(-1i*x*t'); % CDF cdf = (xMean - x)/2 + imag(E * (cft ./ t)); cdf = 0.5 - (cdf * dt) / pi; cdf = reshape(max(0,min(1,cdf)),n,m); % PDF pdf = 0.5 + real(E * cft); pdf = (pdf * dt) / pi; pdf = reshape(max(0,pdf),n,m); ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 41 APPENDIX B. MATLAB ALGORITHM CF2DISTFFT FOR NUMERICAL INVERSION OF THE CHARACTERISTIC FUNCTION BASED ON THE FFT ALGORITHM % QF if ~isempty(prob) isPlot = options.isPlot; options.isPlot = false; [n,m] = size(prob); prob = prob(:); maxiter = options.maxiter; crit = options.crit; qf = options.qf0; criterion = true; count = 0; [res,cdfQ,pdfQ] = ... cf2DistGP([],qf,[],options); options = res.options; while criterion count = count + 1; correction = (cdfQ - prob) ./ pdfQ; qf = qf - correction; [~,cdfQ,pdfQ] = ... cf2DistGP([],qf,[],options); criterion = any(abs(correction) ... > crit * abs(qf)) ... && max(abs(correction)) ... > crit && count < maxiter; end qf = reshape(qf,n,m); prob = reshape(prob,n,m); options.isPlot = isPlot; else qf = []; count = []; correction =[]; end %% RESULT result.cdf = cdf; result.pdf = pdf; result.qf = qf; result.x = x; result.xMean = xMean; result.xStd = xStd; result.xMin = xMin; result.xMax = xMax; result.prob = prob; result.SixSigmaRule = options.SixSigmaRule; result.t = t; result.T = t(end); result.dt = dt; result.cf = cf; result.N = N; result.count = count; result.correction = correction; result.options = options; result.tictoc = toc; %% PLOT the PDF / CDF if length(x)==1, options.isPlot = false; end if options.isPlot figure plot(x,pdf,'.-') grid title('PDF Specified by the CF') xlabel('x') ylabel('pdf') % figure plot(x,cdf,'.-') grid title('CDF Specified by the CF') xlabel('x') ylabel('cdf') end end function [result,cdf,pdf,qf] = ... cf2DistFFT(cfFun,y,prob,options) %cf2DistFFT calculates the approximate values % of CDF, PDF, and QF by numerical inversion of % the characteristic function CF by using the % FFT algorithm. % % SYNTAX: % [result,cdf,pdf,qf] = ... % cf2DistFFT2(cfFun,y,prob,options) % Viktor Witkovsky (witkovsky@savba.sk) % Ver.: 24-Apr-2016 17:12:15 %% CHECK THE INPUT PARAMETERS if nargin < 1 error('Too few inputs'); end if nargin < 4, options = []; end if nargin < 3, prob = []; end if nargin < 2, y = []; end if ~isfield(options,'N') options.N = 2^10; end N = options.N; if ~isfield(options,'SixSigmaRule') options.SixSigmaRule = 6; end if ~isfield(options,'minY') options.minY = []; end if ~isfield(options,'maxY') options.maxY = []; end if ~isfield(options,'isForcedSymmetric') options.isForcedSymmetric = []; end if isempty(options.isForcedSymmetric) isForcedSymmetric = false; else isForcedSymmetric = ... options.isForcedSymmetric; end if ~isfield(options,'isZeroSymmetric') options.isZeroSymmetric = []; end if isempty(options.isZeroSymmetric) if isForcedSymmetric isZeroSymmetric = true; else isZeroSymmetric = false; end else isZeroSymmetric = options.isZeroSymmetric; end if ~isfield(options,'isPositiveSupport') options.isPositiveSupport = []; end ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 42 if isempty(options.isPositiveSupport) isPositiveSupport = false; else isPositiveSupport = ... options.isPositiveSupport; end if ~isfield(options,'isPlot') options.isPlot = true; end if ~isfield(options,'delta') options.tolDiff = 1e-4; end %% MOMENTS AND SUPPORT h = options.tolDiff; t = h*(1:4); cf = cfFun(t); meanY = real((-cf(2) + 8*cf(1) – ... 8*conj(cf(1))+ conj(cf(2)))/(1i*12*h)); m2 = real(-(conj(cf(4)) - 16*conj(cf(3)) + ... 64*conj(cf(2)) + 16*conj(cf(1)) - 130 + ... 16*cf(1) + 64*cf(2) - 16*cf(3) + ... cf(4))/(144*h^2)); stdY = sqrt(m2 - meanY^2); A = meanY - options.SixSigmaRule * stdY; B = meanY + options.SixSigmaRule * stdY; if isPositiveSupport if A <= 0 && ... isempty(options.isForcedSymmetric) A = max(0,A); isForcedSymmetric = true; elseif A > 0 && ... isempty(options.isForcedSymmetric) isForcedSymmetric = false; end end % Use the specified values (if available) if ~isempty(options.minY), A = options.minY; end if ~isempty(options.maxY), B = options.maxY; end % Symmetric support [-B,B] ? if isForcedSymmetric || isZeroSymmetric B = options.SixSigmaRule * ... sqrt(stdY^2 + meanY^2); if ~isempty(options.maxY) B = options.maxY; end A = -B; end %% CHARACTERISTIC FUNCTION CF k = (0:(N-1))'; t = 2*pi * (0.5-N/2+k) / (B-A); cf = cfFun(t(N/2+1:end)); cf = [conj(cf(end:-1:1));cf]; % CF of the 'SYMETRIZED' distribution if isForcedSymmetric cf = real(cf); end %% PDF BY the FFT algorithm dy = (B-A)/N; C = (-1).^((1-1/N)*(A/dy+k))/(B-A); D = (-1).^(-2*(A/(B-A))*k); pdfFFT = real(C.*fft(D.*cf)); cdfFFT = cumsum(pdfFFT*dy); yFFT = A + k * dy; if options.isZeroSymmetric cdfFFT = cdfFFT + 0.5 - ... (cdfFFT(N/2+1)+cdfFFT(N/2))/2; end % SPECIAL TREATMENT for symmetrized distribution if isForcedSymmetric pdfFFT = max(0,2*pdfFFT(N/2+1:end)); cdfFFT = cdfFFT + 0.5 - ... (cdfFFT(N/2+1)+cdfFFT(N/2))/2; cdfFFT = min(1, ... max(0,2*cdfFFT(N/2+1:end)-1)); yFFT = yFFT(N/2+1:end); else pdfFFT = max(0,pdfFFT); cdfFFT = min(1,max(0,cdfFFT)); end yMin = min(yFFT); yMax = max(yFFT); %% INTERPOLATE QUANTILE FUNCTION : QF(prob) if isempty(prob) prob = [0.9,0.95,0.975,0.99,0.995,0.999]; end [cdfU,id] = unique(cdfFFT); yyU = yFFT(id); szp = size(prob); qfFun = @(prob) interp1([-eps;cdfU],... [-eps;yyU+dy/2],prob); qf = reshape(qfFun(prob),szp); % INTERPOLATE CDF/QF/PDF if isempty(y) y = linspace(A,A+(N-1)*dy,100); end szy = size(y); cdfFun = @(x) interp1([-eps;yyU+dy/2],... [-eps;cdfU],x(:)); cdf = reshape(cdfFun(y),szy); try pdfFun = @(x) interp1(yFFT,pdfFFT,y(:)); pdf = reshape(pdfFun(y),szy); catch warning('Unable to interpolate') pdf = NaN*y; pdfFun = []; end %% RESULT result.y = y; result.cdf = cdf; result.pdf = pdf; result.prob = prob; result.quant = qf; result.cdfFun = cdfFun; result.pdfFun = pdfFun; result.qfFun = qfFun; result.yMin = yMin; result.yMax = yMax; result.cdfMin = min(cdfFFT); result.cdfMax = max(cdfFFT); result.N = N; result.Details.yFTT = yFFT; result.Details.pdfFFT = pdfFFT; result.Details.cdfFFT = cdfFFT; result.Details.meanY = meanY; result.Details.stdY = stdY; result.Details.A = A; result.Details.B = B; result.Details.dy = dy; result.Details.dt = 2*pi/(B-A); result.Details.t = t; result.Details.cf = cf; result.Details.cfFun = cfFun; result.options = options; ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 43 REFERENCES [1] JCGM100:2008, Evaluation of measurement data – Guide to the expression of uncertainty in measurement (GUM 1995 with minor corrections), in JCGM - Joint Committee for Guides in Metrology, (ISO, BIPM, IEC, IFCC, ILAC, IUPAC, IUPAP and OIML, 2008). [2] JCGM101:2008, Evaluation of measurement data – Supplement 1 to the Guide to the expression of uncertainty in measurement – Propagation of distributions using a Monte Carlo method, in JCGM - Joint Committee for Guides in Metrology, (ISO, BIPM, IEC, IFCC, ILAC, IUPAC, IUPAP and OIML, 2008). [3] L.J. Gleser, Assessing uncertainty in measurement, Statistical Science, 277 (1998). [4] A. Possolo and B. Toman, Assessment of measurement uncertainty via observation equations, Metrologia 44, p. 464 (2007). [5] A. Forbes and J. Sousa, The GUM, Bayesian inference and the observation and measurement equations, Measurement 44, 1422 (2011). [6] C. Elster, Bayesian uncertainty analysis compared with the application of the GUM and its supplements, Metrologia 51, p. S159 (2014). [7] R.Willink, What can we learn from the GUM of 1995?, Measurement (2016). [8] R.N. Kacker, R. Kessel and J.F. Lawrence, Removing divergence of JCGM documents from the GUM (1993) and repairing other defects, Measurement 88, 194 (2016). [9] M.G. Cox and P.M. Harris, Validating the applicability of the GUM procedure, Metrologia 51, p. S167S175 (2014). [10] P.M. Harris and M.G. Cox, On a Monte Carlo method for measurement uncertainty evaluation and its implementation, Metrologia 51, p. S176-S182 (2014). [11] R. Willink, Measurement Uncertainty and Probability (Cambridge University Press, 2013). [12] R. Willink, An improved procedure for combining Type A and Type B components of measurement uncertainty. International Journal of Metrology and Quality Engineering 4, p. 55-62 (2013). [13] M. Korzen and S. Jaroszewicz, Pacal: A Python package for arithmetic computations with random variables, Journal of Statistical Software 57, 1 (2014). [14] R.C. Williamson and T. Downs, Probabilistic arithmetic. I. Numerical methods for calculating convolutions and dependency bounds, International Journal of Approximate Reasoning 4, 89 (1990). [15] M. Sheppard, Apply a function to a set of probability distributions. Matlab Central File Exchange, file 34108 (Dec 2011 (updated Apr 2012), 2012). [16] T.A. Driscoll, N. Hale and L. N. Trefethen, Chebfun guide, Pafnuty Publications (2014). [17] A. Asheim and D. Huybrechs, Complex Gaussian quadrature for oscillatory integral transforms, IMA Journal of Numerical Analysis , p. drs060 (2013). [18] D. Levin, Fast integration of rapidly oscillatory functions, Journal of Computational and Applied Mathematics 67, 95 (1996). [19] G. Milovanović, Numerical calculation of integrals involving oscillatory and singular kernels and some applications of quadratures, Computers & Mathematics with Applications 36, 19 (1998). [20] A. Sidi, The numerical evaluation of very oscillatory infinite integrals by extrapolation, Mathematics of Computation 38, 517 (1982). [21] A. Sidi, A user-friendly extrapolation method for oscillatory infinite integrals, Mathematics of Computation 51, 249 (1988). [22] A. Sidi, A user-friendly extrapolation method for computing infinite range integrals of products of oscillatory functions, IMA Journal of Numerical Analysis 32, 602 (2012). [23] J. Abate and W. Whitt, The Fourier-series method for inverting transforms of probability distributions, Queueing Systems 10, 5 (1992). [24] N.G. Shephard, From characteristic function to distribution function: A simple framework for the theory, Econometric Theory 7, 519 (1991). [25] L.A. Waller, B.W. Turnbull and J.M. Hardin, Obtaining distribution functions by numerical inversion of characteristic functions with applications, The American Statistician 49,346 (1995). [26] R. Zieliński, High-accuracy evaluation of the cumulative distribution function of 𝛼-stable symmetric distributions, Journal of Mathematical Sciences 105, 2630 (2001). [27] R.L. Strawderman, Computing tail probabilities by numerical Fourier inversion: The absolutely continuous case, Statistica Sinica, 175 (2004). [28] M.J. Korczynski, M. Cox and P. Harris, Convolution and uncertainty evaluation, Series on Advances in Mathematics for Applied Sciences 72, p. 188 (2006). [29] E. Lukacs, Characteristics functions, Griffin, London (1970). [30] N.L. Johnson, S. Kotz, and N. Balakrishnan, Continuous univariate distributions, vol. 1, Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. (1994). [31] N.L. Johnson, S. Kotz, and N. Balakrishnan, Continuous univariate distributions, vol. 2, Wiley Series in Probability and Mathematical Statistics: Applied Probability and Statistics. (1995). [32] H. Cohen, F. R. Villegas and D. Zagier, Convergence acceleration of alternating series, Experimental Mathematics 9, 3 (2000). [33] J. Gil-Pelaez, Note on the inversion theorem, Biometrika 38, 481 (1951). [34] J. Imhof, Computing the distribution of quadratic forms in normal variables, Biometrika 48, 419 (1961). [35] R. Davies, Algorithm AS 155: The distribution of a linear combinations of 𝜒2 random variables, Applied Statistics 29, 232 (1980). [36] V. Witkovský, On the exact computation of the density and of the quantiles of linear combinations of t and F random variables, Journal of Statistical Planning and Inference 94, 1 (2001). [37] V. Witkovský, Matlab algorithm TDIST: The distribution of a linear combination of Student’s t random variables, in COMPSTAT 2004 Symposium, ed. J. Antoch (Physica- Verlag/Springer, Heidelberg, Germany, 2004) pp. 1995–2002. [38] V. Witkovský, Computing the distribution of a linear combination of inverted gamma variables, Kybernetika 37, 79 (2001). [39] V.Witkovský, G.Wimmer and T. Duby, Logarithmic Lambert 𝑊 × 𝐹 random variables for the family of chi-squared distributions and their applications, Statistics & Probability Letters 96, 223 (2015). [40] V. Witkovský, G. Wimmer and S. Ďuriš, On statistical methods for common mean and reference confidence intervals in interlaboratory comparisons for temperature, International Journal of Thermophysics 36, 2150 (2015). [41] W. Hürlimann, Improved FFT approximations of probability functions based on modified quadrature rules, International Mathematical Forum 8, 829 (2013). %% PLOT THE PDF/CDF, if required if options.isPlot figure plot(yFFT,pdfFFT,'-','LineWidth',2) grid title('PDF Specified by the CF') xlabel('y') ylabel('pdf') figure plot(yFFT,cdfFFT,'-','LineWidth',2) grid title('CDF Specified by the CF') xlabel('y') ylabel('cdf') end end ACTA IMEKO | www.imeko.org November 2016 | Volume 5 | Number 3 | 44 [42] D.H. Bailey and P.N. Swarztrauber, The fractional Fourier transform and applications, SIAM Review 33, 389 (1991). [43] P. Carr and D. Madan, Option valuation using the fast Fourier transform, Journal of Computational Finance 2, 61 (1999). [44] K. Chourdakis, Option pricing using the fractional FFT, Journal of Computational Finance 8, 1 (2004). [45] M. Held, CFH Toolbox (Characteristic Function Option Pricing). Matlab Central File Exchange, file 46489 (May 5 (updated Aug 12), 2014). [46] Y.S. Kim, S. Rachev, M.L. Bianchi and F.J. Fabozzi, Computing VaR and AVaR in infinitely divisible distributions, Probability and Mathematical Statistics 30, 223 (2010). [47] EA-4/02, Evaluation of the uncertainty of measurement in calibration, M: 2013 (2013). [48] A. Possolo, Statistical models and computation to evaluate measurement uncertainty, Metrologia 51, p. S228 (2014). Numerical inversion of a characteristic function: An alternative tool to form the probability distribution of output quantity in linear measurement models