Measurement Science and Technology ACCEPTED MANUSCRIPT Performance evaluation of Monte Carlo simulation: Case study of Monte Carlo approximation vs. analytical solution for a chi-squared distribution To cite this article before publication: Hanno Ohvril et al 2019 Meas. Sci. Technol. in press https://doi.org/10.1088/1361-6501/ab57a9 Manuscript version: Accepted Manuscript Accepted Manuscript is “the version of the article accepted for publication including all changes made as a result of the peer review process, and which may also include the addition to the article by IOP Publishing of a header, an article ID, a cover sheet and/or an ‘Accepted Manuscript’ watermark, but excluding any other editing, typesetting or other changes made by IOP Publishing and/or its licensors” This Accepted Manuscript is © 2019 IOP Publishing Ltd. During the embargo period (the 12 month period from the publication of the Version of Record of this article), the Accepted Manuscript is fully protected by copyright and cannot be reused or reposted elsewhere. As the Version of Record of this article is going to be / has been published on a subscription basis, this Accepted Manuscript is available for reuse under a CC BY-NC-ND 3.0 licence after the 12 month embargo period. After the embargo period, everyone is permitted to use copy and redistribute this article for non-commercial purposes only, provided that they adhere to all the terms of the licence https://creativecommons.org/licences/by-nc-nd/3.0 Although reasonable endeavours have been taken to obtain all necessary permissions from third parties to include their copyrighted content within this article, their full citation and copyright line may not be present in this Accepted Manuscript version. Before using any content from this article, please refer to the Version of Record on IOPscience once published for full citation and copyright details, as permissions will likely be required. All third party content is fully copyright protected, unless specifically stated otherwise in the figure caption in the Version of Record. View the article online for updates and enhancements. This content was downloaded from IP address 129.137.5.42 on 18/11/2019 at 11:14 Page 1 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 Performance Evaluation of Monte Carlo Simulation: Case study of Monte Carlo 4 5 6 approximation vs. analytical solution for a chi-squared distribution 7 8 Hanno Ohvril1, Alan H. Tkaczyk1,2*, Peeter Saari1, Tõnu Kollo3, Koit Mauring1, Piia 9 10 Post1, Martin Vilbaste4, Jüri Vedru1, Cagatay Ipbüker1 11 12 1Institute of Physics, University of Tartu, Estonia 13 14 2 15 Institute of Technology, University of Tartu, Estonia 16 17 3Institute of Mathematics and Statistics, University of Tartu, Estonia 18 19 4Institute of Chemistry, University of Tartu, Estonia 20 21 22 * Corresponding Author: alan@ut.ee 23 24 25 26 ABSTRACT 27 28 29 The guide to the expression of uncertainty in measurement (GUM) describes the law of propagation 30 31 of uncertainty for linear models based on the first-order Taylor series approximation of Y = f(X1, X2, 32 33 …, XN). However, for non-linear models this framework leads to unreliable results while estimating 34 the combined standard uncertainty of the model output [u(y)]. In such instances, it is possible to 35 36 implement the method(s) described in Supplement 1 to GUM – Propagation of distributions using a 37 38 Monte Carlo Method. As such, a numerical solution is essential to overcome the complexity of the 39 analytical approach to derive the probability density functions of the output. In this paper, Monte 40 41 Carlo simulations are performed with the aim of providing an insight into the analytical 42 43 transformation of the probability density function (PDF) for Y = X2 where X is normally distributed 44 45 and a detailed comparison of analytical and Monte Carlo approach results are provided. This paper 46 displays how the used approach enables to find PDF of Y = X2 without the use of special functions. 47 48 In addition, the singularity of the PDF and the nonsymmetric coverage interval are also discussed. 49 50 51 52 53 Keywords: GUM; Uncertainty estimation; Monte Carlo method; Non-central non- 54 55 standard chi-squared distribution 56 57 58 59 60 1. Introduction to uncertainty and chi-squared distribution 1 Acc pted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 2 of 26 1 2 3 Any obtained quantity as a result of observations, measurements, modelling or prediction 4 5 6 is associated with an uncertainty that emerges from the followed procedure. The concept 7 8 of errors in measurement was established in the beginning of the 20th century (Wallis & 9 10 Roberts, 1956; Traub, 1997; Lane, 2011). The misconception of attributing the term 11 12 ‘error’ on ‘measurement uncertainty’ was resolved with the publication of the Guide to 13 14 15 the expression of uncertainty in measurement (GUM). GUM established a standard 16 17 procedure for assessing uncertainty (GUM-1993; GUM-1995; JCGM:2008). In addition, 18 19 GUM states that evaluation of uncertainty is not a routine task, but it depends on the 20 21 22 understanding and analysis of the performed method as well as the evaluation of the 23 24 practitioner itself. GUM also accepts approaches to uncertainty evaluation, including 25 26 analytical methods used to derive an exact algebraic form for the probability distribution 27 28 29 for the output Y, or a Monte Carlo method (MCM) with controlled accuracy, etc. 30 31 GUM is mainly concerned about the expression of uncertainty of the measurable 32 33 34 quantity, called the ‘measurand’ – Y. The measurand is determined from N other input 35 36 quantities, X1, X2, ..., XN, through a multivariate functional relationship, Y = f (X1, X2, X3, 37 38 …, XN), where xi denote possible values of corresponding random variable Xi, 39 40 41 respectively. Each input quantity in this relationship has its own uncertainty, expressed as 42 43 u(x1), u(x2), u(x3), …, u(xN), whereby x1, x2,…xN are the best estimates of input quantities 44 45 X1, X2,…XN. The standard uncertainties of input quantities are either evaluated as 46 47 standard deviations of repeated measurement values (type-A uncertainties) or by standard 48 49 50 deviations of the assumed probability density functions (type-B uncertainties). GUM 51 52 defines the standard uncertainty of the measurand as follows: 53 54 𝑁 𝑁 55 𝜕𝑦 𝜕𝑦2 (1) ( 56 𝑢 𝑦 ) = ∑ ∑ ( ) ( ) 𝑢(𝑥 , 𝑥 ) 𝜕𝑥𝑖 𝜕𝑥 𝑖 𝑗 𝑗 57 𝑖=1 𝑗=1 58 59 60 2 Acce ted Manuscript Page 3 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 𝑁 𝑁 4 𝜕𝑦 2 𝜕𝑦 𝜕𝑦 (2) = ∑ ( 𝑢(𝑥𝑖)) + ∑ ( ) ( ) 𝑢(𝑥𝑖, 𝑥5 𝜕𝑥 𝑗 ) 𝑖 𝜕𝑥 𝜕𝑥 6 𝑖=1 𝑖,𝑗=1 𝑖 𝑗 7 𝑖≠𝑗 8 𝜕𝑦 9 where the partial derivatives stand as the sensitivity coefficients, 𝑢(𝑥 𝜕𝑥 𝑖 , 𝑥𝑖) = 10 𝑖 11 12 𝑢(𝑥𝑖) × 𝑢(𝑥𝑖) = 𝑢 2(𝑥𝑖) as the estimated variance of xi, and 𝑢(𝑥𝑖, 𝑥𝑗); 𝑖 ≠ 𝑗 as the 13 14 estimated covariance associated with xi and xj. The input quantities are often assumed to 15 16 17 be mutually uncorrelated (𝑢(𝑥𝑖, 𝑥𝑗) = 0, when 𝑖 ≠ 𝑗), which helps simplify Eq. (2), 18 19 considering only 𝑖 = 𝑗: 20 21 𝑁 22 𝜕𝑦 2 𝑢2(𝑦) = ∑ ( ) 𝑢2 (3) 23 (𝑥𝜕𝑥 𝑖 ) 24 𝑖𝑖=1 25 26 Equations (2) and (3), often called as the law of propagation of uncertainty, are based on 27 28 29 a first-order Taylor series approximation of Y = f (X1, X2, ..., XN) and they express the 30 31 basic GUM framework recommendation for evaluation of uncertainty of a multivariate 32 33 system. 34 35 36 However, there are situations where the application of the present GUM-framework leads 37 38 to unreliable results. If the model is non-linear coupled with high relative uncertainties of 39 40 41 input quantities the present GUM framework provides unreliable estimate for the 42 43 combined standard uncertainty of model output u(y). Also, if the distribution of the 44 45 output Y cannot be assumed to be a Gaussian or a Student’s t-distribution it is not correct 46 47 48 to use the coverage factor k = 2 or corresponding Student’s t-coefficient to calculate the 49 50 expanded uncertainty at P=95 % coverage probability. In these cases it can be 51 52 recommended to use the Monte Carlo method (MCM) based on the Supplement 1 to 53 54 GUM (GUM-S1, 2008). The shortcomings of the GUM are currently being dealt with a 55 56 57 new revision of the GUM which is expected to be consistent with GUM supplements 58 59 (Bich et al., 2012; Bich, 2014; Bich et al., 2016). All the distributions of input quantities 60 3 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 4 of 26 1 2 3 will be estimated from a Bayesian point of view eliminating the need to distinguish the 4 5 6 uncertainties as type A and type B uncertainties (Bich et al., 2016). Besides some other 7 8 important changes, the new GUM will also recommend using MCM if one has little 9 10 knowledge about the distribution of the model output (Bich et al., 2016). One good 11 12 example about non-linear model is Y = X2: 13 14 2 15 𝜕𝑦2 16 𝑢 (𝑦) = ( ) 𝑢 2(𝑥) (4) 𝜕𝑥 17 18 19 𝜕𝑦 (5) 20 𝑢(𝑦) = 𝑢(𝑥) = 2𝑥𝑢(𝑥) 𝜕𝑥 21 22 23 Y = X2 stands as the simplest nonlinear model with widespread applications, for example 24 25 its use in kinematics and in fluid mechanics with velocity profiles. However, it is also 26 27 28 implemented in different measurement systems, e.g. in remote sensing, especially in the 29 30 evaluation of coverage of areas with certain specification or with cloud cover, as well as 31 32 measurements of irradiation from large territories. In addition, in electrical engineering, 33 34 power meters actually detect X2 (Carobbi, 2014). A different perspective to problems of 35 36 37 metrology in measurement systems was also presented in a very systematical manner by 38 39 Danilov (2016). Moreover, another important application is for evaluating the Word 40 41 Error Rate (WER) in automated measurement systems, e.g. speech recognition and 42 43 44 Analog to Digital Converters, etc. (Catelani et al., 2010). 45 46 The probability density function (PDF) of Y = X2 is asymmetric, as it cannot be negative 47 48 49 and this leads to problems in constructing its coverage intervals. For a symmetric PDF 50 51 output, symmetric coverage intervals are usually used, but in the case of an asymmetric 52 53 PDF, the user must have an insight into the shape and properties of the PDF to proceed 54 55 56 with design of coverage interval (Bich, 2014; Bich et al., 2012, 2016; Lira, 2019; Willink, 57 58 2016). 59 60 4 Acc pted Manu cript Page 5 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 For example, when X, a priori, follows a normal distribution with non-zero mean X and 4 5 6 non-unit standard deviation x, X ~ N(X;  ), then the X 2 X has non-central non-standard 7 8 chi-squared distribution with one degree of freedom. 9 10 11 The chi-squared distribution, derived from a set of n-independent standard normal 12 13 variables, Xi ~ N(0; 1) with i = 1, 2, … , n, is widely used for evaluating the goodness of 14 15 16 fit of an observed distribution to a theoretical one. Unfortunately, a general case of Y = 17 18 X2, for a single non-central non-standard input X ~ N(X  0; X  1) is not studied in the 19 20 books nor in the literature (Ventsel, 1969; Papoulis, 1990, 1991; Traub, 1997; Rice, 2007; 21 22 23 Fornasini, 2008; Veerarajan, 2009; Lane, 2011; Suhov and Kelbert, 2014; Sahoo, 2015; 24 25 Thomopoulos, 2017; Kelbert and Suhov, 2018). There is however a comparative study in 26 27 the literature focusing on the random measurement errors and indirect measurement 28 29 30 errors by Monte-Carlo method by Labutin and Pugin (2000). 31 32 The study of Kent & Hainsworth (1995) concludes that in the absence of any clear 33 34 35 optimality criteria for choosing confidence intervals of the χ 2-distribution, a ‘symmetric 36 37 range’ interval is the best choice. Furthermore, in the study of Attivissimo et al. (2012) 38 39 where the use of frequentist and the Bayesian approach to measurement uncertainty is 40 41 42 discussed, the authors consider an electric circuit consisting only of a battery of voltage X 43 44 and a noisy unit-value resistor consuming the power W = X2. A thorough uncertainty 45 46 analysis of the circuit requires the computation of a number of PDFs along with 47 48 49 expectations and variances. This study involves uniform and normal distributions for 50 51 input X, Bayesian, frequentist and Monte Carlo approaches and the results are compared. 52 53 An interesting result from this study is the plot from which the PDF of X2 for normal 54 55 noncentral input X can be depicted, although the analytical expression of PDF for W is 56 57 58 given only in a general form of the Bayes’ formula. 59 60 5 Accepted Manus ript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 6 of 26 1 2 3 The aim of this study is to present a comprehensive analysis of the PDF in the output of Y 4 5 6 = X 2 for non-central non-standard normal input X, supported with Monte Carlo simulation 7 8 for extracting the PDF of the normal squared. 9 10 11 A literature review on the PDF of Y = X 2 reveals a series of publications by Cox & Harris 12 13 (1999; 2003; 2006) and Cox & Siebert (2006), which show that numerical solution using 14 15 MCM is an effective tool to approximate the output PDF by a frequency distribution. In 16 17 18 addition, Cox & Siebert (2006) demonstrated that in the case of the simplest (i.e. 19 20 uniform) input PDF, use of the Heaviside step function and the Markov formula allows to 21 22 derive an analytic expression of the PDF of Y = X2. 23 24 25 The aim of this study is to present a comprehensive analysis of the PDF in the output of Y 26 27 = X2 for a priori known non-central non-standard normal input X ~ N(x, x) determined 28 29 30 with Monte Carlo simulation for extracting the PDF of the squared normally distributed 31 32 variable. This paper improves upon the findings of the previous studies presented above. 33 34 It also introduces important information to help avoid potential mistakes while obtaining 35 36 37 results, displays the key aspects of implementation, as well as a special emphasis on a 38 39 simpler PDF equation for practical calculations. This paper also covers the singularity of 40 41 the PDF which is the main problem of the Monte Carlo simulation and the nonsymmetric 42 43 44 coverage interval. 45 46 47 48 49 2. On the two Monte Carlo experiments 50 51 The series of works carried out by Cox et al. (1999; 2003; 2006) on the univariate model 52 53 54 Y = X 2 can be considered as the basis of this problem with significant contributions to the 55 56 area. Particularly, Cox & Harris (1999; 2003) considered Gaussian input X ~ N(0.5; 0.2), 57 58 with mean X = 0.5, and standard deviation σX = 0.2. The PDF of this input, although 59 60 6 Acc pted Manuscript Page 7 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 symmetric, is right-shifted and narrowed compared to the standard normal distribution X 4 5 6 ~ N(0; 1). The authors ran M = 10 000 MCS trials to draw a rough PDF histogram of 20 7 8 columns for output Y. Obtained set of output quantities yi enables a quick evaluation of 9 10 the two main statistical output quantities, mathematical expectation Y and standard 11 12 13 deviation σY, as well as the median, mode and coverage intervals. The PDF for output Y = 14 15 X2 is asymmetric and includes only nonnegative values which reflect its non-Gaussian 16 17 origin. It is not correct to represent a coverage interval for Y in a symmetric form,  18 Y  k 19 20 σY, where k is a coverage factor (usually, k = 1, 2 or 3). When the GUM framework 21 22 cannot be applied, Monte Carlo simulation (MCS) presents itself as a good alternative 23 24 (Cox & Harris, 2006; Cox & Siebert, 2006; GUM-S1; GUM-Introduction; GUM-S2). 25 26 27 Rearrangement of yi into a non-decreasing order enables determining the quantiles to 28 29 define the required expanded uncertainties through possible coverage intervals [ylow, 30 31 yhigh], where the endpoints depend on the particular output PDF. For example, for the 32 33 34 output 0.025 and 0.975 quantiles define a 95% coverage interval. Evidently, the set yi 35 36 allows depicting the same value coverage interval using any another appropriate pair of 37 38 39 quantiles such as 0.015 and 0.965, or 0.040 and 0.990, etc. 40 41 Using MCS for the statistical description of output quantities of nonlinear models 42 43 presents two possible problems: 44 45 46 - MCS easily overlooks sharp peaks at possible singular points of an output PDF 47 48 while using MCS-derived histograms for visualization of the shape of the PDF of 49 50 51 an output quantity; 52 53 - superficial analysis of empirical cumulative distribution functions (CDF) for Y 54 55 can lead to erroneous generalizations about the median value for the set of yi. 56 57 58 For instance, a histogram of 20 columns for Y = X2 with X ~ N(0.5; 0.2), obtained using 59 60 10 000 MCS trials, did not reveal singularity in vicinity of y = 0 (Cox & Harris, 1999; 7 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 8 of 26 1 2 3 2003). However, this was solved when Cox & Harris (2006) readdressed the problem of 4 5 6 Y = X 2 by inserting X ~ N(1.2; 0.5), performing M = 50 000 MCS trials and drawing a Y 7 8 histogram of 70 columns. The enhanced resolution displayed discernible features in the 9 10 histogram which was previously not evident. The more detailed PDF appeared to be 11 12 bimodal with a sharp peak depicted by the first column. The authors also calculated for 13 14 15 the output quantity the estimates of the mathematical expectation Y and the associated 16 17 standard uncertainty 𝜎𝑌 as provided by the law of propagation of uncertainty: 18 19 2 20 𝜇𝑌 = 𝜇𝑋 = 1. 2 2 = 1.44 (6) 21 22 23 𝜎𝑌 = 2𝜇𝑋𝜎𝑋 = 2 × 1.2 × 0.5 = 1.20 (7) 24 25 26 The Monte Carlo experiment resulted in different values y = 1.70, Y = 1.26. On the 27 28 other hand, considered Y = X2 stands out as a relatively simple non-linear model and 29 30 31 enables analytical explanation of discrepancies in calculation of Y and Y as well as the 32 33 appearance of a sharp peak of the PDF. 34 35 3. Expectation and variance of Y = X2 36 37 38 In order to get the expectation of the output: 39 40 41 Y = EX2 (8) 42 43 44 the relationship for variance DX of a random variable X, 45 46 𝜎2 = 𝜇 − 𝜇2 (9) 47 𝑋 𝑌 𝑋 48 49 gives the expectation for output Y: 50 51 52 2 2 𝜇𝑌 = 𝜇𝑋 + 𝜎𝑋 (10) 53 54 55 The obtained formula is universal, regardless of which PDFs are assigned to X, 56 57 unconstrained by the requirement of the Gaussian distribution as input. When (10) is 58 59 applied to X ~ N(1.2; 0.5): 60 8 Accepted Manuscript Page 9 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 4 𝜇𝑌 = 1. 2 2 + 0. 52 = 1.44 + 0.25 = 1.69 (11) 5 6 which is in agreement with the Monte Carlo experiment cited above, Y = 1.70, but 7 8 9 contradicts the result of Eq. (6), Y = 1.44. Derivation of a formula for variance (σ 2 𝑌) of 10 11 the output requires the calculation of the 4th noncentral moment for input X: 12 13 14 𝜎2𝑌 = 𝐷𝑌 = 𝐸(𝑌 − 𝜇 ) 2 𝑌 = 𝐸𝑌 2 − 𝜇2𝑌 (12) 15 16 considering Eq. (10) and since: 17 18 19 𝐸𝑌2 = 𝐸𝑋4 (13) 20 21 equation (12) can be rewritten as: 22 23 24 𝜎2𝑌 = 𝐸𝑋 4 − (𝜇2 2 2𝑋 + 𝜎𝑋 ) (14) 25 26 th 27 here the 4 non-central moment of the normal distribution: 28 29 ∞ 30 4 (15) 𝐸𝑋 = ∫ 𝑥4𝑝𝑋(𝑥)𝑑𝑥 = 𝜇 4 𝑋 + 6𝜇 2 𝑋𝜎 2 4 31 𝑋 + 3𝜎𝑋 −∞ 32 33 where pX(x) is the PDF for a Gaussian input, X. Combining of Eqs. (14) and (15) gives 34 35 36 for the variance of output, Y: 37 38 𝜎2 = 4𝜇2 𝜎2𝑌 𝑋 𝑋 + 2𝜎 4 𝑋 (16) 39 40 41 which is only valid for the Gaussian input. Applying Eq. (16) applied to a normal 42 43 44 quantity discussed above, X ~ N(1.2; 0.5): 45 46 𝜎2𝑌 = 4  1.2 2  0.52 + 20.54 = 1.565 (17) 47 48 49 50 𝜎𝑌 = √1.565 = 1.251 (18) 51 52 matches the MCM result of Cox & Harris (2006), Y = 1.26, and proves the use of Eq. (7) 53 54 55 less accurate. 56 57 4. The PDF for Y = X2 58 59 60 9 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 10 of 26 1 2 3 There are two methods to derive the PDF for Y = X2. One is the PDF transformation 4 5 6 technique which converts the input probability density function, pX(x), into output, pY(y), 7 8 and the other one is the CDF differentiation technique that starts from the cumulative 9 10 distribution functions, GX(y) and GY(y). 11 12 13 The equality of probability elements at both sides of the model must consider both 14 15 inverse functions (𝑥 = ±√𝑦) and are written as: 16 17 18 𝑝𝑋(−𝑥)𝑑𝑥 + 𝑝𝑋(𝑥)𝑑𝑥 = 𝑝𝑌(𝑦)𝑑𝑦 (19) 19 20 21 𝑑𝑥 𝑑𝑥 (20) 22 𝑝𝑌(𝑦) = 𝑝𝑋(−√𝑦) + 𝑝𝑋(√𝑦) 𝑑𝑦 𝑑𝑦 23 24 25 1 1 26 𝑝𝑌(𝑦) = 𝑝 (− 𝑦) + 𝑝 ( 𝑦) (21) 𝑋 √ 𝑋 √ 2√𝑦 2√𝑦 27 28 29 (21) is the generic form and is valid for any PDF in the input (Papoulis, 1990; 1991). In 30 31 32 the case of normal non-standard input, X ~ N(X; X): 33 34 2 1 1 𝑥−𝜇− ( 𝑋) (22) 35 𝑝𝑋(𝑥) = 𝑒 2 𝜎𝑋 36 𝜎𝑋√2𝜋 37 38 the general formula (22) for the output PDF transforms into: 39 40 41 2 21 1 √𝑦+𝜇𝑋 1 1 √𝑦−𝜇𝑋 (23) − ( ) − ( ) 42 𝑝𝑌(𝑦) = 𝑒 2 𝜎𝑋 + 𝑒 2 𝜎𝑋 43 2𝜎𝑋√2𝜋√𝑦 2𝜎𝑋√2𝜋√𝑦 44 45 46 In (23) the presence of √𝑦 in the denominators means that the obtained PDF has a 47 48 singularity at y = 0: 49 50 51 lim 𝑝𝑌(𝑦) = ∞ (24) 52 𝑦→0 53 54 For a central but non-standard normal distribution (X = 0, X  1), (23) equals to: 55 56 57 (25) 58 59 60 and for the standard normal distribution (X = 0, X = 1): 10 Accepted Manuscript Page 11 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 4 (26) 5 6 which is the chi-squared or χ2-distribution with one degree of freedom. 7 8 9 Alternatively, the PDF for Y = X2 (for a non-centered X) is by definition: 10 11 12 𝑑𝐺𝑌(𝑦) (27) 𝑝𝑌(𝑦) = 13 𝑑𝑦 14 15 where GY(y) is the CDF for output, defined as: 16 17 18 𝐺𝑌(𝑦) = 𝑃(𝑌 ≤ 𝑦) = 𝑃(𝑋 2 ≤ 𝑦) = 𝑃(−√𝑦 ≤ 𝑋 ≤ √𝑦) (28) 19 20 √𝑦 21 𝐺𝑌(𝑦) = ∫ 𝑝𝑋(𝑥)𝑑𝑥 (29) 22 − 𝑦 23 √ 24 25 Splitting the area of integration at x = 0: 26 27 0 𝑦28 √ (30) 29 𝐺𝑌(𝑦) = ∫ 𝑝𝑋(𝑥)𝑑𝑥 + ∫ 𝑝𝑋(𝑥)𝑑𝑥 30 −√𝑦 0 31 32 −√𝑦 √𝑦 (31) 33 𝐺𝑌(𝑦) = − ∫ 𝑝𝑋(𝑥)𝑑𝑥 + ∫ 𝑝𝑋(𝑥)𝑑𝑥 34 0 0 35 36 Differentiating with respect to y gives an interim result of: 37 38 39 𝑑𝐺𝑌(𝑦) 𝑑 −√𝑦 𝑑 √𝑦 40 𝑝𝑌(𝑦) = = − ∫ 𝑝 (𝑥)𝑑𝑥 + ∫ 𝑝 (𝑥)𝑑𝑥 (32) 𝑑𝑦 𝑑𝑦 𝑋0 𝑑𝑦 𝑋 41 0 42 43 here we denote the first and second integral as g1(y) and g2(y), respectively. Before using 44 45 the rule of differentiating with respect to the upper limit of integration, changes in 46 47 48 variables should be performed. For the first integral: 49 50 51 𝑑𝜏(𝑦) −1 𝜏 = −√𝑦; = (33) 𝑑𝑦 2√𝑦 52 53 54 𝑑 −√𝑦 𝑑 𝜏 −1 (34) 55 𝑔1(𝑦) = − ∫ 𝑝𝑋(𝑥)𝑑𝑥 = −( ∫ 𝑝𝑑𝑦 𝑑𝜏 𝑋 (𝑥)𝑑𝑥) × 56 0 0 2√𝑦 57 58 1 1 (35) 59 𝑔1(𝑦) = 𝑝𝑋(𝜏) × = 𝑝𝑋(−√𝑦) 60 2√𝑦 2√𝑦 11 Accept d Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 12 of 26 1 2 3 using  = √𝑦 for the second integral in (32): 4 5 6 1 (36) 7 𝑔2(𝑦) = 𝑝𝑋(√𝑦) 8 2√𝑦 9 10 11 both PDF parts, g1(y) and g2(y), given by (35) and (36), respectively, together equal again 12 13 to equation (21), confirming the steps taken. 14 15 16 By involving modified Bessel functions of the first kind, general formula (23) for PDF of 17 18 Y = X2 can be represented in a form usually applied in theory of chi-squared distribution: 19 20 21 1 1𝐴 = and 𝑎 = (37) 22 2𝜎 2𝑋√2𝜋√𝑦 2𝜎𝑋 23 24 25 Eq. (23) can be written as: 26 27 2𝑝 (𝑦) = 𝐴(𝑒−𝑎(√𝑦+𝜇𝑋) + 𝑒−𝑎(√𝑦−𝜇𝑋) 2 28 𝑌 ) (38) 29 30 −𝑎(𝑦+𝜇231 𝑝 (𝑦) = 𝐴𝑒 𝑋) (𝑒2𝑎√𝑦𝜇 𝑦𝜇 𝑋 𝑌 + 𝑒 −2𝑎√ 𝑋) (39) 32 33 34 using (Andras, Baricz, 2008): 35 36 37 𝑒 𝑧 + 𝑒−𝑧 𝜋𝑧 (40) 38 cosh 𝑧 = = √ 𝐼−1/2(𝑧) 2 2 39 40 41 where I−1/2 is modified Bessel function with −1/2 degrees of freedom, Eq. (39) can be 42 43 rewritten in a more desired form: 44 45 46 1 √𝜇 2 𝑋 𝑦 + 𝜇𝑋 𝜇𝑋 (41) 47 𝑝𝑌(𝑦) = 2 𝑒𝑥𝑝 (−2𝜎 4 2𝜎2 ) 𝐼−1/2 (√𝑦 ) 48 𝑋 √𝑦 𝑋 𝜎 2 𝑋 49 50 51 which, compared to (23), is not so convenient for practical calculations. 52 53 54 The interpretation of the PDF for Y = X 2 can be better understood by analyzing the two 55 56 examples presented in the previous section. For both cases, it is assumed the input X is 57 58 normally distributed, according to the two Monte Carlo experiments, X ~ N(0.5; 0.2) and 59 60 12 Acc pted Manuscript Page 13 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 X ~ N(1.2; 0.5), respectively (Cox & Harris, 1999; 2003; 2006). The general formula for 4 5 6 output PDF (24) transforms respectively into: 7 8 2.5 1 2 2.5 1 2 (42) 9 𝑝𝑌(𝑦) = 𝑒 −12.5(√𝑦+0.5) + 𝑒−12.5(√𝑦−0.5) 10 √2𝜋 √𝑦 √2𝜋 √𝑦 11 12 13 1 1 2 1 1 2 (43) 𝑝𝑌(𝑦) = 𝑒 −2(√𝑦+1.2) + 𝑒−2(√𝑦−1.2) 14 √2𝜋 √𝑦 √2𝜋 √𝑦 15 16 17 Plots of the two input PDFs, and their outputs, Y = X2, according to (42) and (43), are 18 19 presented in Fig. 1. Symmetric coverage intervals,   σ and 20  2σ, can be seen for the 21 22 input and the output, but since Y is neither symmetrically distributed nor have any 23 24 negative values, symmetric coverage intervals with respect to the expectation of the 25 26 output do not present any valuable information. However, the output coverage intervals 27 28 29 can be calculated by appropriate integration of a particular pY(y) both for the analytical 30 31 approach and the MCS. 32 33 34 35 36 Output: Y = X2, 37 38 Input: normal, X~N(0.5; 0.2) 𝜇𝑌 = 0.29, 𝜎𝑌 = 0.2078 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 Input: normal, X~N(1.2; 0.5) Output: Y = X 2, 𝜇𝑌 = 1.69, 𝜎𝑌 = 1.251 59 60 13 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 14 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 Fig. 1. On the left: probability density functions for two Gaussian input quantities, X ~ 20 21 N(0.5; 0.2), and X ~ N(1.2; 0.5) respectively. On the right: probability density functions 22 23 24 for corresponding output quantities, Y = X 2. The dashed horizontal lines indicate the 25 26 coverage intervals ( and  2). For Y the  2 intervals include infeasible negative X2 27 28 values. The filled triangles indicate locations of mathematical expectations at the 29 30 31 horizontal axes and the empty triangles of the medians for outputs. The continuous 32 33 horizontal lines seen on the graphs on the right show the 95% coverage intervals, plotted 34 35 according to 2.5% and 97.5% quantiles. 36 37 38 5. The CDF for Y = X 2 39 40 41 The cumulative distribution functions (CDF) are essential for evaluating the 42 43 normalization condition of the PDF for a random quantity, as well as for a rapid 44 45 assessment of median and quantile values, peak-event probabilities and the coverage 46 47 2 48 intervals. For Y = X , in the case of a normal non-standard input, X ~ N(X; X), the CDF 49 50 is defined as: 51 52 𝑧 𝑧 21 1 √𝑦+𝜇 − ( 𝑋) 𝑑𝑦53 2 𝜎 (44) 𝐺𝑌(𝑧) = ∫ 𝑝𝑌(𝑦) 𝑑𝑦 == ∫ 𝑒 𝑋 54 0 2𝜎𝑋 √2𝜋 0 √𝑦 55 56 57 where the first and second parts of the equation are denoted as Part I, G1(z), and Part II, 58 59 60 G2(z), respectively. Applying the change in variables, from y to t, for Part I and Part II: 14 A cept d Manuscript Page 15 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 1 4 𝑡 = (√𝑦 + 𝜇𝑋) (45) 5 √2𝜎𝑋 6 7 1 8 𝑡 = (√𝑦 − 𝜇𝑋) (46) 9 √2𝜎𝑋 10 11 G1(z) and G2(z) are obtained as: 12 13 14 1 𝐵1(𝑧) 2 1 1 (47) 15 𝐺 (𝑧) = ∫ 𝑒−𝑡1 𝑑𝑡 = erf(𝐵1) - erf(𝐴) 16 √𝜋 𝐴 2 2 17 18 19 20 21 1 𝐵2(𝑧) 1 1 𝐺 (𝑧) = ∫ 𝑒−𝑡 2 2 𝑑𝑡 = erf(𝐵2) + erf(𝐴) (48) 22 √𝜋 −𝐴 2 2 23 24 25 where the new limits of integration: 26 27 𝜇𝑋 𝐴 = (49) 28 √2𝜎𝑋 29 30 31 √𝑧 + 𝜇𝑋𝐵 (𝑧) = (50) 32 1 √2𝜎𝑋 33 34 35 √𝑧 − 𝜇𝑋𝐵 (𝑧) = (51) 36 2 √2𝜎𝑋 37 38 39 and the error function: 40 𝑥 41 2 2erf(𝑥) = ∫ 𝑒−𝑡 𝑑𝑡 (52) 42 𝜋 43 0 44 45 The sum of (47) and (48) gives (after returning from z to y): 46 47 1 1 (53) 48 𝐺𝑌(𝑦) = erf(𝐵1) + erf (𝐵 ) 49 2 2 2 50 51 52 where the coefficients B1 and B2 depend on y, X and X. If the first particular input is 53 54 considered again, X ~ N(0.5; 0.2), from Cox & Harris (1999; 2003), B1 and B2 become: 55 56 57 √𝑦 + 𝜇𝑋𝐵1(𝑦) = = √12.5𝑦 + √3.125 (54) 58 √2𝜎𝑋 59 60 15 Accept d Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 16 of 26 1 2 3 √𝑦 − 𝜇4 𝑋𝐵2(𝑦) = = √12.5𝑦 − √3.125 (55) 5 √2𝜎𝑋 6 7 8 For the second input, X ~ N(1.2; 0.5), following Cox & Harris (2006), the CDF is: 9 10 1 1 11 𝐺𝑌(𝑦) = erf(𝐵3) + erf(𝐵4) (56) 12 2 2 13 14 where: 15 16 √𝑦 + 𝜇17 𝑋𝐵3(𝑦) = = √2𝑦 + 1.2√2 (57) 18 √2𝜎𝑋 19 20 21 √𝑦 − 𝜇 (58) 𝑋 22 𝐵4(𝑦) = = √2𝑦 − 1.2√2 23 √2𝜎𝑋 24 25 The two obtained CDFs GY(y) are plotted in Fig. 2 and denoted as “Cumulative 1” and 26 27 “Cumulative 2”, respectively. 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 Fig. 2. Two examples of the cumulative distribution function of Y = X2: 45 46 47 Cumulative 1 for X ~ N(0.5; 0.2); Cumulative 2 for X ~ N(1.2; 0.5). 48 49 50 51 52 6. Width of narrow peaks in the Y = X2 PDF curves 53 54 55 In this section, the obtained CDFs are used for evaluating the contribution of a peak near 56 57 the origin of the PDF Y = X2 (as seen on the right of the Fig. 1). The area under the peak 58 59 60 is calculated for both cases. For the first input, X ~ N(0.5; 0.2) the local minimum is 16 Accepted Man script Page 17 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 located at y = 0.0145, which defines the range of the narrow peak as an interval, 0 < y  4 5 6 0.0145. Contribution of this peak to the whole PDF is equal to GY(0.0145). Before 7 8 performing (53), limits of integration, B1 and B2 are calculated: 9 10 11 √0.0145 + 0.5 (59) 12 𝐵1(𝑦 = 0.0145) = = 2.19350 √2 × 0.2 13 14 15 √0.0145 − 0.5𝐵2(𝑦 = 0.0145) = = −1.34203 (60) 16 √2 × 0.2 17 18 19 1 1 (61) 20 𝐺𝑌(0.0145) = erf(2.19350) - erf(1.34203) = 0.02789 21 2 2 22 23 The result indicates that a narrow peak near the origin covers approximately 2.8% of the 24 25 2 26 area under the PDF curve of Y = X for the input condition of X ~ N(0.5; 0.2). For the 27 28 second input, X ~ N(1.2; 0.5), the extension of the narrow peak is, 0 < y  0.100. The 29 30 limits of integration, B3 and B4: 31 32 33 √0.100 + 1.2 (62) 34 𝐵3(𝑦 = 0.100) = = 2.14427 35 √2 × 0.5 36 37 √0.100 − 1.238 𝐵4(𝑦 = 0.100) = = −1.24984 (63) 39 √2 × 0.5 40 41 substitution of these into (49): 42 43 44 1 1𝐺𝑌(0.100) = erf(2.14427) - erf(1.24984) = 0.03736 (64) 45 2 2 46 47 48 indicating that the contribution of the narrow peak covers approximately 3.7% of the area 49 50 under the PDF curve of Y = X2 for X ~ N(1.2; 0.5). 51 52 53 54 55 7. Quantiles and coverage intervals for Y = X2 56 57 58 Analytical presentation of the cumulative distribution function enables calculation for Y 59 60 of quantiles and coverage intervals. Often the 95% coverage interval is defined between 17 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 18 of 26 1 2 3 the 0.025 and 0.975 quantiles of the PDF for Y. These quantiles for the first example, X ~ 4 5 6 N(0.5; 0.2), using (49), (50) and (51) are found as: 7 8 𝐺𝑌(0.0125) = 0.025 (65) 9 10 11 𝐺𝑌(0.7957) = 0.975 (66) 12 13 14 which means that the lower and upper quantiles are, y2.5% = 0.0125 and y97.5% = 0.7957, 15 16 respectively. 17 18 19 Quantiles for the second input, X ~ N(1.2; 0.5): 20 21 22 𝐺𝑌(0.0561) = 0.025 (67) 23 24 𝐺𝑌(4.7524) = 0.975 25 (68) 26 27 which means that y2.5% = 0.0561 and y97.5% = 4.7524. These 95% of coverage intervals are 28 29 displayed as continuous horizontal lines in Fig. 1 (on the right). The same coverage 30 31 32 interval of 95% may be also given by any other appropriate pair of quantiles such as y1.5% 33 34 and y96.5%, etc. 35 36 37 38 39 8. About the median of Y = X2 40 41 42 The cumulative distribution function is a useful visual aid to understanding the nature of 43 44 a random quantity. The most important characteristic of a CDF itself is the median. A 45 46 detailed examination of tabulated values for the cumulative distribution functions GY(y) 47 48 49 gives the following results: 50 51 GY(y = 0.25) = 0.499999713, for the input X ~ N(0.5; 0.2) (69) 52 53 54 GY(y = 1.44) = 0.499999207, for the input X ~ N(1.2; 0.5) (70) 55 56 57 from here, it appears that with a very high rate of accuracy the following statement can be 58 59 written: 60 18 Accepted M nuscript Page 19 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 median (𝑌 = 𝑋 2) = 𝜇2 4 𝑋 (71) 5 6 However, this statement cannot present the entirety of the variability of Y = X2 as the 7 8 9 median should depend on both of the input quantities, X and X, not solely on X. 10 11 Equation (53) clarifies this situation. With the substitution of: 12 13 14 𝑦 = 𝜇 2 𝑋 (72) 15 16 into (53), the coefficients B1, B2 and the CDF, 𝐺𝑌(𝑦 = μ 2 𝑥), take the following values: 17 18 19 𝑦+𝜇 √2𝜇𝐵 (𝑦) = √ 𝑋 = 𝑋 √ 𝑦−𝜇 𝐵 (𝑦) = 𝑋20 1 2 = 0 (73) √2𝜎𝑋 𝜎𝑋 √2𝜎𝑋 21 22 23 1 1 2𝜇𝐺 (𝑦 = 𝜇2 √ 𝑋 𝑌 𝑋) = erf(𝐵1) = erf( ) (74) 24 2 2 𝜎𝑋 25 26 27 which means that the CDF, 𝐺 (μ 2 𝑌 𝑥), actually depends on the ratio, X /X, and tends to the 28 29 value 0.5 with X over X (Fig. 3). For X /X = 1.6, the CDF already becomes equal to 30 31 0.499. 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 Fig. 3. Cumulative distribution function GY (y) at 𝑦 = μ2𝑥, as a function of the ratio of two 48 49 50 input quantities, X /X. 51 52 53 The validity of (71) can approximated from X /X = 1.6 towards larger values and the 54 55 median Y = X2 as μ2𝑥 can be calculated. The results of medians obtained with (71) can be 56 57 explained by considering the same inputs, X ~ N(0.5; 0.2), and X ~ N(1.2; 0.5), for which 58 59 60 the ratios X /X are 2.5 and 2.4 respectively. 19 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 20 of 26 1 2 3 4 5 6 9. Readdressing a Monte Carlo experiment 7 8 9 As mentioned previously, MCS method is relatively easy to implement and therefore 10 11 presents itself as a good alternative to the analytical derivation of output statistical 12 13 quantities. However, MCS easily overlooks possible sharp narrow peaks in the shape of a 14 15 16 PDF. The proper comprehension of PDF shape obtained with MCS can be achieved by 17 18 considerably increasing the number of Monte Carlo trials and the number of histogram 19 20 columns. In this study, for the first input, X ~ N(0.5; 0.2), a series of M = 1 000 000 21 22 23 Monte Carlo trials has been performed. 24 25 The Monte Carlo run for Y = X2 resulted in the output mean of Y = 0.2902, and a 26 27 28 standard deviation of Y = 0.2079. Both results stand close to the ones obtained from 29 30 analytical method Y = 0.29 and Y = 0.2078. Additionally, the probabilistically 31 32 33 symmetric coverage interval corresponding to 95 % coverage probability was found to lie 34 35 between 0.0126 and 0.7962 while the corresponding analytically calculated coverage 36 37 interval lies between 0.0125 and 0.7957 indicating a good match between the coverage 38 39 40 intervals found by the two methods. The said-one million sets of xi and yi were 41 42 reassembled into histograms of 50, 100, 300 and 1000 columns, and displayed in Fig. 4, 43 44 respectively. It is obvious that 50 columns are not enough for proper visualization of the 45 46 47 output PDF whereas 100 could be used but the output PDF is still not clearly evident. 48 49 However, the histogram of 300 columns reveals the singularity near the origin while the 50 51 histogram of 1000 columns only sharpens it. 52 53 54 Output: Y = X 2, Y = 0.2902, Y = 55 56 Input: normal, X ~ N(0.5; 0.2) 0.2079 57 58 59 60 20 Accepted Manuscript Page 21 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 Fig. 4. Results of a series of 1000000 MC trials as histograms of 50, 100, 300 and 1000 59 60 columns. 21 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 22 of 26 1 2 3 4 5 6 10. Conclusions 7 8 9 For reliable estimation of the model output uncertainty, the model input quantities should 10 11 be specified in terms of probability density functions (PDFs). In order to determine the 12 13 output PDF, the practitioner must choose between analytical and numerical methods. 14 15 16 Analytical methods ask the user to have calculus and probability and statistics knowledge 17 18 as prerequisites (Rice, 2007; Fornasini, 2008; Thomopoulos, 2017), but the user can 19 20 always use certain softwares to help derive the analytical solutions. 21 22 23 In this study, an analytic approach is described for a simple univariate model Y = X2 24 25 where X is the Gaussian input with non-zero expectation and non-unit standard deviation. 26 27 28 The analytic approach enabled a detailed description of singularity in Y near the origin as 29 30 well as to reveal a peculiarity in calculation of the median for Y. However, for example in 31 32 some cases of environmental modelling, either the set of input quantities (e.g. photons of 33 34 35 solar rays, ionizing radiation from polluted territories, etc.) or the model itself (e.g. 36 37 processes of light scattering on a single plant or in the entire plant cover, propagation of 38 39 pollution in an environment, etc.) may be too complicated for an analytical 40 41 representation. In these situations, Monte Carlo simulations (MCS) appear to be the only 42 43 44 alternative method. 45 46 47 This study also demonstrated that for a relatively simple model, Y = X 2, there can be 48 49 unexpected results, such as overlooking narrow peaks, and recommends that a 50 51 sufficiently large number of trials should be chosen to obtain an adequate plot of the 52 53 output histogram. For Y = X2, a plot with the sufficient details that enables detection of 54 55 56 singularity in the output, was achieved using 300 histogram columns on the basis of 1 57 58 million MC trials. 59 60 22 Accepted Manuscript Page 23 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 4 5 6 Acknowledgements 7 8 9 This work has been supported by the TERRITORIES project, which is part of the 10 11 CONCERT EJP, funded from the Euratom research and training program 2014-2018 12 13 under Grant Agreement No 662287.This publication reflects only the authors’ view. 14 15 16 Responsibility for the information and views expressed therein lies entirely with on the 17 18 authors. The European Commission is not responsible for any use that may be made of 19 20 the information it contains. The authors also acknowledge the institutional research 21 22 23 funding IUT20-11 and IUT34-5 of the Estonian Ministry of Education and Research. 24 25 26 27 28 References 29 30 31 32 Andras, S., Baricz, A. (2008). Properties of the probability density function of the non- 33 34 35 central chi-squared distribution. J. Math. Anal. Appl. 346, 395–402. 36 37 Attivissimo, F., Giaquinto, N., and Savino, M. (2012). A Bayesian paradox and its impact 38 39 on the GUM approach to uncertainty. Measurement 45, 2194–2202. 40 41 42 Bich, W., Cox, M.G., Dybkaer, R., Elster, C., Estler, W.T., Hibbert, B., Imai, H., Kool, 43 44 W., Michotte, C., Nielsen, L. and Pendrill, L. (2012). Revision of the ‘Guide to the 45 46 expression of uncertainty in measurement’. Metrologia 49(6), 702-705. 47 48 Bich, W. (2014). Revision of the ‘guide to the expression of uncertainty in measurement’. 49 50 51 Why and how. Metrologia 51(4), S155-S158. 52 53 Bich, W., Cox, M. and Michotte, C. (2016). Towards a new GUM—an update. 54 55 Metrologia 53(5), S149-S159. Metrol. Meas. Syst., 8(2), 195-204. 56 57 58 Carobbi, C. (2014). Bayesian inference on a squared quantity. Measurement 48, 13–20. 59 60 23 Accepted Manu cript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 24 of 26 1 2 3 Catelani, M., Zanobini, A., Ciani, L. (2010). Uncertainty Interval Evaluation Using the 4 5 6 Chi-Square and Fisher Distributions in the Measurement Process. 7 8 Cox, M. and Harris, P. (1999). Up a GUM tree? Counting on IT. Information System 9 10 Engineering. Newsletter of the National Physical Laboratory, 8. Teddington, 11 12 13 Middlesex, UK, 4−5. 14 15 Cox, M. and Harris, P. (2003). Up a GUM tree? Try the Full Monte! Centre for 16 17 Mathematics & Scientific Computing, National Physical Laboratory. Teddington, 18 19 20 Middlesex, UK, 3 p. 21 22 Cox, M. and Harris, P. (2006). Software Support for Metrology. Best Practice Guide No. 23 24 6, Uncertainty Evaluation. Tech. Report DEM-ES-011. National Physical Laboratory. 25 26 Teddington, UK, 167 p. 27 28 29 Cox, M. and Siebert, B. R. L. (2006). The use of a Monte Carlo method for evaluating 30 31 uncertainty and expanded uncertainty. Metrologia 43 S249–59. 32 33 Danilov, A.A. (2016). General Problems of Metrology and Measurement Technique The 34 35 36 Development of Measurement Systems and Their Metrological Support. Measurement 37 38 Techniques, 59-9, 899-903. 39 40 Fornasini, P. (2008). The Uncertainty in Physical Measurements. Springer, 289 p. 41 42 43 GUM-1993. Guide to the expression of uncertainty in measurement. ISO/TMBG 44 45 Technical Management Board, 105 p. 46 47 GUM-1995, JCGM 100:2008. Guide to the Expression of Uncertainty in Measurement 48 49 (GUM-1995 with minor corrections). Joint Committee for Guides in Metrology. Paris, 50 51 52 Geneva, 134 p. 53 54 GUM-S1 (JCGM 101:2008). Evaluation of measurement data. Supplement 1 to the GUM 55 56 − Propagation of distributions using a Monte Carlo method. Joint Committee for 57 58 59 Guides in Metrology, Paris, Geneva, 90 p. 60 24 Accepted Manuscript Page 25 of 26 AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 1 2 3 GUM-S2, JCGM 101:2011. Evaluation of measurement data. Supplement 2 to the GUM 4 5 6 − Extension to any number of output quantities. Joint Committee for Guides in 7 8 Metrology. Paris, Geneva, 80 p. 9 10 GUM-Introduction (JCGM 104:2009). An introduction to the “Guide to the expression 11 12 13 of uncertainty in measurement". Joint Committee for Guides in Metrology. Paris, 14 15 Geneva, 28 p. 16 17 Kelbert, M., and Suhov, Y. (2018). Probability and Statistics by Example: Volume 1, 18 19 20 Basic Probability and Statistics. MCCME Publishing, Moscow, 518 p. (3 rd ed., 21 22 translation from English, in Russian). 23 24 Kent, J. T., and J. Hainsworth, T. J. (1995). Confidence intervals for the noncentral chi- 25 26 27 squared distribution. Journal of Statistical Planning and Inference. 46, 2, 147−159. 28 29 Labutin, S.A., Pugin, M. V. (2000). Summation of Random Measurement Errors and 30 31 Analysis of Indirect-Measurement Errors by Monte-Carlo Method. Measurement 32 33 Techniques, 43-11, 918-922 (translated from Russian). DOI: 34 35 36 10.1023/A:1010968307890 37 38 Lane, D. M. (2011). History of the Normal Distribution. In: Introduction to Statistics. 39 40 Online Edition, 692 p. http://onlinestatbook.com/Online_Statistics_Education.pdf 41 42 43 Lira, I. (2019). The GUM Revision: Where do We Stand? A Personal View. 44 45 Measurement, Proc. 12th Int. Conf., Smolenice, Slovakia, 39−46. 46 47 Papoulis, A., (1990). Probability & Statistics. Prentice-Hall, 454 p. 48 49 50 Papoulis, A., (1991). Probability, Random Variables, and Stochastic Processes. 51 52 McGraw-Hill Series in Electrical Engineering, 666 p. 53 54 Rice, J. A. (2007). Mathematical Statistics and Data Analysis. Brooks/Cole, Cengage 55 56 57 Learning, 603 p. 58 59 60 25 Accepted Manuscript AUTHOR SUBMITTED MANUSCRIPT - MST-109347.R1 Page 26 of 26 1 2 3 Sahoo, P. (2015). Probability and Mathematical Statistics. University of Louisville, 713 4 5 6 p. 7 8 Suhov, Y., and Kelbert, M. (2014). Probability and Statistics by Example. Cambridge 9 10 University Press, 461 p. 11 12 Thomopoulos, N. T. (2017). Statistical Distributions. Applications and Parameter 13 14 15 Estimates. Springer, 176 p. 16 17 Traub, R. E. (1997). Classical Test Theory in Historical Perspective. Educational 18 19 Measurement: Issues and Practice. 16, 4, 8−14. 20 21 22 Veerarajan, T. S. (2009)., Probability, Statistics and Random Processes. McGraw-Hill 23 24 Publ. 595 p. 25 26 Ventsel, E. S. (1969) The Probability Theory (4th ed.). Nauka, Moscow, 576 p. (in 27 28 29 Russian). 30 31 Wallis, W. A. and Roberts, H. V. (1956). Statistics, a new approach. Free Press, 646 p. 32 33 Willink, R. (2016). What can we learn from the GUM of 1995? Measurement, 91, 692- 34 35 36 698. 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 26 Accepted Manuscript