Asymmetrical uncertainties

In several disciplines, measurement results occasionally are expressed using coverage intervals that are asymmetric relative to the measured value. The conventional treatment of such results, when there is the need to propagate their uncertainties to derivative quantities, is to replace the asymmetric uncertainties by ‘symmetrized’ versions thereof. We show that such simplification is unnecessary, illustrate how asymmetry may be modeled and recognized explicitly, and propagated using standard Monte Carlo methods. We present three distributions (Fechner, skew-normal, and generalized extreme value), among many available alternatives, that can be used as models for asymmetric uncertainties associated with scalar input quantities, in the context of the measurement model considered in the GUM. We provide an example where such uncertainties are propagated to the uncertainty of a ratio of mass fractions. We also show how a similar, model-based approach can be used in the context of data reductions from interlaboratory studies and other consensus building exercises where the reported uncertainties are expressed asymmetrically, illustrating the approach to obtain consensus estimates of the absorption cross-section of ozone, and of the distance to galaxy M83 in the Virgo cluster.


Introduction
The following are but a few illustrative measurement results that comprise uncertainty evaluations expressed in the form of coverage regions that are not centered at, hence are asymmetrical relative to the measured value: • Gavrilyuk et al [19] report a measurement of the half-life of the double-β decay of 78 Kr, as ( may then be regarded as expressing the belief that, with approximately 68% probability, it includes the true value of the half-life: that is, as a ±1-sigma interval assuming that the uncertainty is expressed using a Gaussian distribution. • The Event Horizon Telescope Collaboration [41] estimate the distance from Earth to the black hole at the center of the M87 galaxy in the Virgo cluster as 16.8 +0.8 −0.7 Mpc. • Crowther et al [16] report an asymmetric uncertainty associated with their estimate of the age of the R136 star cluster in the Tarantula nebula as 1.5 +0.3 −0.7 × 10 6 year. • The measurement result that Yoshino et al [44] obtained for the absorption cross-section of ozone at 253.65 nm was reported as 1145 +7.1 −14.4 × 10 −20 cm 2 molecule −1 . • Nelson et al [32, page 8657] use an asymmetric probability distribution to characterize the uncertainty associated with the purity, 95.97%, of 3-epi-25(OH)D 3 In several disciplines, measurement results occasionally are expressed using coverage intervals that are asymmetric relative to the measured value. The conventional treatment of such results, when there is the need to propagate their uncertainties to derivative quantities, is to replace the asymmetric uncertainties by 'symmetrized' versions thereof. We show that such simplification is unnecessary, illustrate how asymmetry may be modeled and recognized explicitly, and propagated using standard Monte Carlo methods. We present three distributions (Fechner, skew-normal, and generalized extreme value), among many available alternatives, that can be used as models for asymmetric uncertainties associated with scalar input quantities, in the context of the measurement model considered in the GUM. We provide an example where such uncertainties are propagated to the uncertainty of a ratio of mass fractions. We also show how a similar, model-based approach can be used in the context of data reductions from interlaboratory studies and other consensus building exercises where the reported uncertainties are expressed asymmetrically, illustrating the approach to obtain consensus estimates of the absorption cross-section of ozone, and of the distance to galaxy M83 in the Virgo cluster.
Keywords: uncertainty, ozone, black hole, skew-normal, interlaboratory study, symmetrization, consensus S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. determined using quantitative 1 H-nuclear magnetic resonance spectroscopy, and report a 95% coverage interval that is asymmetric relative to the measured value, [94.73%, 97.59%]. Nelson et al [33, table 6] report an asymmetric, 95% coverage interval for the mass fraction of benzoic acid in a primary standard for quantitative NMR, as 0.999 92 g/g − 0.000 06 g/g + 0.000 04 g/g.
• The Bank of England uses the Fechner distributiondiscussed in section 2.1-to characterize the uncertainty associated with forecasts of future inflation, thus capturing perceived asymmetries between upside and downside risks [15]. , which together define a box that includes the true source with 90% probability [22].
This contribution addresses the problem of propagating measurement uncertainties expressed asymmetrically, when the corresponding measured values are combined using a suitable function to produce an estimate of a derivative quantity, or when comparable, independent measurements of the same quantity are combined to form a consensus value.
Rather than 'symmetrizing' the uncertainties reported asymmetrically before combining the measurement results, for example as proposed by Barlow [7,8] and Audi et al [4, appendix A], we address the problem directly by preserving and modeling such uncertainties, first (section 2) in the context of the measurement model considered in the Guide to the expression of uncertainty in measurement (GUM) [26], and second (section 3) in the context of statistical models used in consensus building exercises, for example in interlaboratory studies or meta-analyses [27]. Section 4 highlights some conclusions and recommendations.
In section 2 we present three simple models, selected from among the plethora of asymmetric distributions that may be used for the same purpose: the Fechner distribution, the skew-normal distribution, and the generalized extreme value distribution. In section 2.4, we apply all three to the estimation and evaluation of the uncertainty associated with a ratio of mass fractions of two rare earth elements, and compare the results.
In section 3 we demonstrate how the skew-normal distribution may be used to model measurement uncertainties reported asymmetrically, in a meta-analysis of results for the absorption cross-section of ozone, first when evaluating the uncertainty associated with the DerSimonian-Laird estimator of the consensus value, and second when a Bayesian statistical model is used to compute both an estimate of the consensus value and an evaluation of the associated uncertainty. We also demonstrate how the same techniques can be used to blend three determinations of the distance to galaxy M38 in the Virgo cluster, whose associated uncertainties are expressed asymmetrically.

Asymmetric uncertainties in measurement equations
The measurement model considered in the GUM expresses an output quantity y as a known function of several input quantities, y = f (x 1 , . . . , x n ), and handles the uncertainties associated with the input quantities by modeling them as random variables.
The approach to uncertainty propagation described in the GUM uses only the standard deviations of these random variables (and any correlations between them) to produce an approximate evaluation of the standard uncertainty u(y), by application of what is variously known as the delta method [ Alternatively, a Monte Carlo method may be used that propagates the full distributions of those random variables, as described by Morgan and Henrion [30] and by the GUM Supplement 1 (GUM-S1) [24]. Both approaches are implemented in the NIST Uncertainty Machine (uncertainty.nist. gov) [28].
The conventional technique described in the GUM requires that the random variables used as models for the input quantities should have finite standard deviations, but not necessarily that their probability distributions be symmetrical. However, even when the input quantities have symmetrical distributions, the distribution of the output quantity need not be either Gaussian or even symmetric.
The Monte Carlo method will recognize and express any asymmetry that may be present, because it uses the full distributions, not just their standard deviations, and its results can be used to produce valid coverage intervals regardless of whether the distribution of the output quantity is symmetric or asymmetric [36].
We assume that each scalar input x is qualified with a left uncertainty u L (x) and a right uncertainty u R (x), both positive, but the underlying probability distributions are left unspecified. Very often, these are interpreted so that the interval covers the true value of x with probability approximately 68%; in other words, it is a 1-sigma interval.
In our general treatment, we assume that probabilities where X denotes the random variable modeling the uncertainty associated with the input x. In the examples below, p M = 0.5 (that is, x is regarded as an estimate of the median of the underlying distribution), p L is either 0.025 or 0.16, and p R is either 0.84 or 0.975.
In sections 2.1-2.3 we review three different models for asymmetric probability distributions that reproduce the aforementioned, limited information provided about each input quantity. In the online supplementary materials (stacks.iop. org/MET/56/045009/mmedia) we present R code [37] that finds the values of the parameters of the model that best reproduce x − u L (x), x, and x + u R (x) as specified percentiles of the corresponding distribution, using a non-linear, weighted least-squares method. In section 2.4 we present an illustrative example and compare the accuracy that the different models achieve.

Fechner distribution
The Fechner distribution, also known as the split normal or two-piece normal distribution [43], consists of two halfnormal distributions with the same mode, one concentrated to the left of the mode, the other to the right of the mode, and with their respective densities suitably rescaled so that the resulting probability density is continuous. Nelson et al [32] use it to report uncertainties in studies of chemical purity of vitamin D metabolites.
The distribution is determined by three parameters: the mode, µ, and the left and right uncertainties, σ L > 0 and σ R > 0. The mean is µ + 2/π(σ R − σ L ), and the variance is ( This distribution can approximate left or right skewed, asymmetric distributions whose Pearson's moment coefficient of skewness has absolute value no larger than 0.9953, approximately. Figure 1 shows two examples of the Fechner probability density, and algorithm F describes how a Fechner distribution may be chosen to reproduce specified percentiles. R package fanplot [1] provides computational facilities supporting the Fechner distribution. The Fechner distribution may be used for the 'symmetrization' of asymmetric uncertainties very much along the same lines as in Method 2 of Audi et al [4, appendix A]. The idea is to fit an asymmetric distribution to the measurement result (measured value and associated uncertainty, expressed using an asymmetric interval), and then use the mean and standard deviation of the fitted distribution as estimate of the measurand and as standard measurement uncertainty, respectively. Table 1 displays the symmetrization produced by Audi et al [4], and the symmetrization based on the mean and standard deviation of the Fechner distribution, fitted by non-linear weighted least squares to the original measurement results for the half-life T1 /2 , of each of four different radioactive isotopes.
The choice of values for the weights corresponding to the left, central, and right percentiles depends on where the greatest modeling accuracy is required. In most cases where the left and right percentiles are not far out in the tails, equal weights for all three percentiles produce reasonable results. But when one wishes to ensure that the tails are modeled accurately because they have the greatest impact on the uncertainty, possibly at the price of a less accurate reproduction of the central percentile, then left and right weights larger than the central weight may prove helpful.   Table 2. Mean (R), standard deviation (u(R)), and lower ( L 95% (R)) and upper (U 95% (R)) endpoints of a 95% coverage interval for the true value of the ratio w S (La)/w C (Yb Algorithm F has coverage probability about 68%, and 0.645 < u R (x)/u L (x) < 1.55, then the Fechner distribution may be a suitable model for the reported asymmetry. If the coverage probability is approximately 95%, then the model may be suitable when 0.410 < u R (x)/u L (x) < 2.44. F-2: Use R function parFechner defined in Listing 1 of the online supplementary materials to find the values of µ, σ L , and σ R of the best fitting Fechner distribution.

Skew-normal distribution
The skew-normal distribution [6], with computational facilities implemented in R package sn [5], generalizes the Gaussian distribution to accommodate a modicum of either left or right skewness (figure 1). It is specified by three realvalued parameters: location ξ, scale ω > 0, and shape α. The Gaussian distribution corresponds to α = 0. For distributions in this family, Pearson's moment coefficient of skewness does not exceed 0.995. Given a measured value, x, its associated left and right uncertainties, u L (x) and u R (x), and the probabilities p L , p X and p R such that x − u L (x), x, and x + u R (x) are the 100p L , 100p X , and 100p R percentiles of the fitted distribution, the best fitting skew-normal distribution is determined by the following procedure.

SN-1: If the interval
has coverage probability about 68%, and 0.645 < u R (x)/u L (x) < 1.55, then the skew-normal distribution may be a suitable model for the reported asymmetry. If the coverage probability is approximately 95%, then the model may be suitable when 0.410 < u R (x)/u L (x) < 2.44. SN-2: Use R function parSkewNormal defined in Listing 2 of the online supplementary materials to find the values of ξ, ω, and α of the best fitting skew-normal distribution.

Generalized extreme value distribution
The family of generalized extreme value (GEV) distributions [23, chapter 22], with computational facilities implemented in R package evd [39], includes the Gumbel, Fréchet, and reverse Weibull distributions as special cases. The members of this family are specified by real-valued location (µ), scale (σ > 0), and shape (ξ) parameters. When ξ > 0, the support of the distribution is the semi-axis extending from µ − σ/ξ to plus infinity, and when ξ < 0 it is the semi-axis extending from minus infinity to µ − σ/ξ. The Gumbel distribution corresponds to ξ = 0, and its support extends from minus infinity to plus infinity. Given a measured value, x, its associated left and right uncertainties, u L (x) and u R (x), and the probabilities p L , p X and p R such that x − u L (x), x, and x + u R (x) are the 100p L , 100p X , and 100p R percentiles of the fitted distribution, the best fitting GEV distribution is determined by the following procedure. For the GEV densities depicted in figure 1, the fitted values of ξ both are negative, −0.1465 for w La and −0.4967 for w Yb .

Algorithm GEV
GEV-1: The GEV may be an appropriate model for asymmetric distributions with a very wide range of values of the ratio u R (x)/u L (x), regardless of whether the coverage probability of [x − u L (x), x + u R (x)] is around 68% or 95%. GEV-2: Use R function parGEV defined in Listing 3 of the supplementary materials to find the values of µ, σ, and ξ of the best fitting GEV distribution.
Suppose also that the measured mass fraction of lanthanum in samples of a marine sediment has median w S (La) = 2277 µg kg −1 , and that, with 95% probability, its true value lies between L 95% (w S (La)) = 1041 µg kg −1 and U 95% (w S (La)) = 3988 µg kg −1 , suggesting an asymmetric probability distribution skewed to the right.
Since mass fractions of rare earths in terrestrial samples are often reported in the form of ratios to the mass fraction of ytterbium in carbonaceous chondrites [3], this example focuses on the ratio w S (La)/w C (Yb) and evaluates the associated uncertainty, which involves propagating the two asymmetrical uncertainties given above.
Application of each of the three models described in section 2 involved the same steps, and produced the results listed in table 2 and depicted in figure 2: (1) Fit the model distribution to the measurement result for w S (La), matching its 2.5th, 50th, and 97.5th percentiles to L 95% (w S (La)), w S (La), and U 95% (w S (La)); (2) Fit the model distribution to the measurement result for w C (Yb), matching its 2.5th, 50th, and 97.5th percentiles to L 95% (w C (Yb)), w C (Yb), and U 95% (w C (Yb)); (3) Draw a large sample w S,1 (La), …, w S,K (La) from the distribution fitted in step (1); (4) Draw a large sample w C,1 (Yb), …, w C,K (Yb) from the distribution fitted in step (2); (5) Use the sample of ratios w S,1 (La)/w C,1 (Yb), …, w S,K (La)/w C,K (Yb) to produce an estimate of the true value of w S (La)/w C (Yb), and an evaluation of the associated uncertainty.
Listing 4 of the online supplementary materials shows R code used to apply the steps above using the skew-normal model. Similar codes produce the corresponding results for the Fechner and GEV models.
To validate the proposed methods, some underlying probability distributions have to be assigned to w S (La) and to w C (Yb). Since we are using the Fechner, skew-normal, or GEV as models for the asymmetric uncertainties for these mass fractions, to probe their effectiveness, albeit only in one particular case, we chose to model w S (La) and w C (Yb) as random variables with asymmetric distributions different from the three that we use as approximants: w S (La) as value of a random variable with a gamma probability distribution with shape 9 and scale 253 µg kg −1 , and w C (Yb) as value of a random variable with a Weibull distribution with shape 64 and scale 163.94 µg kg −1 .
Under these assumptions, we know the 'right' answer for the distribution of the ratio, and can then compare the approximations obtained using the Fechner, skew-normal, and GEV, against it. The probability density of the corresponding ratio is depicted using a dotted line in figure 2, and its mean, standard deviation, and 2.5th and 97.5th percentiles are listed under gam/wei in table 2. Table 3 lists the measurement results for the ozone absorption cross-section σ at 253.65 nm that Hodges et al [20, table 2] used to produce an estimate of σ as a consensus value. Application of algorithm SN from section 2.2 to the cases whose model is SN, under the assumption that each interval ranging from σ − u L (σ) to σ + u R (σ) has approximately 68% coverage and that σ is the median of the corresponding skew-normal distribution, produced the location ξ, scale ω, and shape α of this distribution, hence u(σ) = ω 1 − 2α 2 /(π(1 + α 2 )) [5].

Absorption cross-section of ozone
Hodges et al [20] computed the consensus value σ DL = 1.1329 × 10 −17 cm 2 molecule −1 using the DerSimonian-Laird procedure for consensus building [17,27]. Because this is a method-of-moments estimate involving only the estimates of σ and the standard measurement uncertainties u(σ), it is indifferent to any asymmetries of the reported uncertainties and is expected to perform well when the asymmetries are mild, as they are in this dataset.
The asymmetries were taken into account in the evaluation of u(σ DL ) = 0.0035 cm 2 molecule −1 by applying a Monte Carlo method similar to the method described by Koepke et al [27, section 5.2], except that the measurement 'errors' were sampled from skew-normal distributions for the laboratories whose probabilistic model (in table 3) was SN, and from Gaussian distributions for the laboratories whose probabilistic model was G. The output of this Monte Carlo method produced a 95% coverage interval for σ DL , ranging from 1.1260 × 10 −17 cm 2 molecule −1 to 1.1397 × 10 −17 cm 2 molecule −1 .
An alternative approach, patterned after the hierarchical Bayesian procedure implemented in the NIST Consensus Builder [27, section 5.3], involves the following steps: (a) Use a random-effects model for the measured values of the absorption cross-section, σ j = θ j + ε j for j = 1, . . . , n, where n = 14 is the number of results in table 3, the {θ j } are a sample from a Gaussian distribution with mean σ 0 (the true value of the absorption cross-section) and standard deviation τ (dark uncertainty [27,42]), and the {ε j } are measurement 'errors' with either skew-normal or Gaussian probability distributions, whose parameter values are listed in table 3. (b) Assign a vague Gaussian prior distribution to σ 0 , with mean 0 × 10 −17 cm 2 molecule −1 and standard deviation 1000 × 10 −17 cm 2 molecule −1 , and a mildly informative prior distribution to τ , in the spirit of empirical Bayes: a half-Cauchy distribution with median equal to one half of the MAD of the {σ j }) (MAD, as defined in R function mad, is the median of the absolute deviations from the median, rescaled so that it reproduces the standard Table 3. Determinations of the ozone absorption cross-section σ at 253.65 nm, and evaluations of the associated uncertainty produced by Hodges et al [20], expressed asymmetrically in the form of left and right uncertainties u L (σ) and u R (σ). The probability distribution used to model the reported uncertainty is either the Gaussian distribution (G) when u L (σ) = u R (σ), or the skew-normal distribution (SN) in the other cases, with location ξ, scale ω, and shape α. In all cases, u(σ) denotes the standard deviation of the model used to describe the reported uncertainty. The labels under lab are as described by Hodges et al [20], who fully reference the corresponding publications.  For MCMC we employed a Hamiltonian or Hybrid Monte Carlo (HMC) sampler [9,18,31]. HMC samples the parameters in the model by simulating a discretized trajectory of a Hamiltonian dynamical system whose potential function is the negative of the logarithm of the probability density of the posterior distribution of interest. HMC samplers often produce samples with autocorrelations of smaller absolute value than a simple Metropolis sampler, and the Markov chain also tends to reach equilibrium faster than for Metropolis samplers. We implemented HMC sampling using facilities available in R package rstan [38]. In particular, we used the so-called No-U-Turn sampler [21] that chooses the number of steps taken in the simulated trajectory adaptively.
Since the standardized difference, |σ DL − σ BAYES | / u 2 (σ DL ) + u 2 (σ BAYES ), between the DerSimonian-Laird and Bayes estimates of the absorption cross-section equals 0.24, the two estimates, σ DL and σ BAYES , are judged not to be significantly different.
The posterior distribution of the dark uncertainty τ is not only concentrated near 0, but the posterior probability is 99.2% that τ is smaller than the median of the standard deviations of the measurement errors for the different laboratories. This is consistent with the fact that the DerSimonian-Laird procedure estimates τ to be 0 (that is, the measured values are not overdispersed by comparison with their reported measurement uncertainties), which motivated the choice of value assigned to the median of the half-Cauchy prior distribution for τ .

Distance to black-hole in galaxy M87
The Event Horizon Telescope Collaboration [41, table 9] lists three measurements of the distance from Earth to galaxy M87 in the Virgo cluster, whose combination produced 16.8 +0. 8 −0.7 Mpc (approximately 55 million light-years) as estimate of that distance. At the center of this galaxy there lies the first black hole ever to be delineated in images, which were acquired by the Event Horizon Telescope [40]. The measurement of that distance is an input to the measurement of the mass of the black hole.
The reported distance measurements are 16 [13], considered to be independent. Each of these results was modeled using a skew-normal distribution. For example, for the first, (16.75 − 1.04) Mpc, 16.75 Mpc, and (16.75 + 1.11) Mpc were treated as the 16th, 50th, and 84th percentiles of the best fitting skew-normal distribution, and similarly for the other two. The 16th and 84th percentiles are chosen because we interpret these asymmetrically reported uncertainties as nominally equivalent to 1-sigma uncertainties. Figure 3 shows that these percentiles reproduce the measurement results exactly.
The same modeling choices and techniques described in section 3.1 produced a sample of size 16 000, resulting from thinning and merging the samples drawn from the posterior distribution of the distance, by 4 MCMC chains of length 250 000 each. The median of this sample is 16.88 Mpc. The corresponding, complete measurement result may then be expressed as 16.88 +0.55 −0.59 Mpc, whose difference from the result reported in the original study is statistically insignificant. Probability density estimates, based on samples of size K = 10 6 , of the probability distribution of the ratio w S (La)/w C (Yb), when the uncertainties associated with the numerator and denominator are modeled using the Fechner, skew-normal, or GEV distributions. The dotted curve, labeled GAM/WEI, reflects the hypothetical true distribution of the ratio assuming that the uncertainty associated with w S (La) is best described using a gamma distribution, and the uncertainty associated with w S (Yb) is best described using a Weibull distribution.

Conclusions and suggestions
The problem of propagating uncertainties expressed asymmetrically can be solved, without replacing these uncertainties with symmetrized proxies, by modeling them using appropriate asymmetric probability distributions and sampling from these using a suitable Monte Carlo procedure.
We have described how the Fechner, the skew-normal, and the generalized extreme value distributions may be used as models for asymmetric uncertainties associated with mutually independent scalar input quantities of a measurement model of the sort considered in the GUM. Many other asymmetric distributions may be used for the same purpose.
When such input quantities are correlated, then either multivariate versions of these or other distributions may be used similarly, or a copula may be employed that imposes the required dependence structure as Possolo [34] illustrates.
The Fechner and skew-normal distributions can model only modest asymmetry while remaining close to the conventional, Gaussian model for measurement error. However, the generalized extreme value distribution can accommodate more extreme skewness, either to the left or to the right. All three performed comparably well in the example involving a ratio of two mass fractions, where associated uncertainties were expressed asymmetrically and whose 'true' models were different from those three.
We have also provided two alternative approaches to the combination of measurement results obtained in an interlaboratory study, using a collection of determinations of the absorption cross-section of ozone as an example. The model for the measured values comprised Gaussian random effects, and either skew-normal or Gaussian laboratory-specific, measurement error models. Both choices can easily be modified, for example, by replacing the Gaussian distribution with a Laplace distribution for the random effects, or by adopting still other asymmetric distributions for the measurement error models.
In one of the alternatives, based on the DerSimonian-Laird procedure, the asymmetries are recognized only during the evaluation of the uncertainty associated with the consensus value. The other alternative was a full-fledged hierarchical Bayesian model, hence was likelihood based, and the asymmetries were recognized fully both in the estimation of the consensus value and in the evaluation of the associated uncertainty. For the data set we chose as illustration, the two alternatives produced statistically indistinguishable results.
We applied the same models and techniques to blend three independent estimates of the distance from Earth to galaxy M87, as an alternative to the procedure that The Event Horizon Telescope Collaboration [41] used for the same purpose, having obtained a result that is statistically indistinguishable from its counterpart in the original study.
The same techniques may be used to model the uncertainties associated with the standard atomic weights that are expressed as intervals whose endpoints are not equidistant from the conventional atomic weight: for example, the standard atomic weight of lithium is [6.938, 6.997], and its conventional atomic weight is 6.94 (www.ciaaw.org, retrieved June 26, 2019).
In all cases, great caution is recommended in verifying that the model adopted to describe the asymmetries can indeed reproduce the measured value, and the endpoints of the reported coverage interval, closely enough for the intended purpose. This verification should be done both numerically and graphically, similarly to what we have done in figure 1. also was improved in consequence of comments and suggestions generously offered by Juris Meija (NRC Canada) and by two anonymous referees.