Provides a printed summary of the results obtained by propagate, with statistics relevant to JCGM 100:2008 (GUM), JCGM 101:2008 (GUM Supplement 1) and JCGM 102:2011 (GUM Supplement 2), as outlined in 'Details'. This includes a validation of the GUM uncertainty framework by the Monte Carlo method (JCGM 101:2008, 8) and checks of the Monte Carlo input against the assigned input distribution.
# S3 method for propagate
summary(object, normality = FALSE, validate = TRUE, ndig = 2,
interval = c("symmetric", "shortest"), sort = FALSE, width = 25, ...)A printed output with the items listed under 'Details'. The propagate object is returned invisibly.
an object returned from propagate.
logical. If TRUE, the three repeated subsampling normality tests are performed. See 'Details'.
logical. If TRUE, the GUM uncertainty framework is validated by the Monte Carlo method according to JCGM 101:2008, 8. See 'Details'.
the Monte Carlo coverage interval against which the GUF interval is validated: the probabilistically symmetric interval (default) or the shortest interval, as in JCGM 101:2008, 9.3. For the stability check, the blocks use the same type of interval.
the number of significant decimal digits of the standard uncertainty \(u(y)\) that defines the numerical tolerance \(\delta\) for validate.
logical. If TRUE, the rows of the first- and second-order uncertainty budgets are sorted by decreasing index.
the maximal width (characters) of the bars of the uncertainty budgets; in the second-order budget, the bars are scaled to the largest absolute share if one exceeds 100 %.
other parameters for future methods.
Andrej-Nikolai Spiess
summary.propagate prints the following items, in this order. Items that are not available for a given object (e.g. the Hessian if no symbolic or numerical derivatives could be obtained, or the correlation checks for a single input) are omitted.
Monte Carlo input model:
The input distribution assigned for the Monte Carlo simulation, mc.margins ("S1" or "mvt"), with the corresponding JCGM clause, or "user-supplied samples" for type = "sim". See 'Details' in propagate.
The degrees of freedom df.MC actually used for sampling. These are the per-input degrees of freedom \(\nu_i\), or the common \(\nu = n - N\) for "mvt" (\(n\) observations, \(N\) input quantities).
Results (JCGM 100:2008 and JCGM 101:2008):
Results of the first- (mean.1, u.1) and second-order (mean.2, u.2; u.2m with the input moments, see 'Details' of propagate) uncertainty propagation, with the coverage interval \(y \pm U\) of the GUM uncertainty framework (GUF).
Results of the Monte Carlo simulation (mean.MC, u.MC, median.MC, mad.MC) with the probabilistically symmetric coverage interval, formed from the order statistics of the \(M\) = nsim model evaluations as in JCGM 101:2008, 7.7.1.
The shortest coverage interval from the same order statistics (JCGM 101:2008, 7.7.2), preferable for asymmetric output distributions.
Moments extension: the standard uncertainties of the Taylor expansion next to the Monte Carlo value, u.1 (first order), u.2 (second order, Gaussian inputs), u.2m (second order with the skewness and excess kurtosis of each input) and u.MC, with their differences to u.MC in percent. For quadratic models with independent inputs, u.2m equals u.MC within the Monte Carlo noise (e.g. 1.622, 3.768, 2.694 and 2.696 W for the power of a heater with a rectangular phase angle, see the vignette); the remaining difference u.2m - u.MC shows higher-order terms of the model (and, for correlated inputs, their joint moments). Printed if u.2m is available.
The effective degrees of freedom \(\nu_{\mathrm{eff}}\) (Welch-Satterthwaite formula or its generalization to correlated components), as obtained from WelchSatter, or df.tot if supplied.
The coverage factor \(k = t_{\max(1, \lfloor\nu_{\mathrm{eff}}\rfloor)}(1-\alpha/2)\), or the supplied k.
The expanded uncertainty \(U = k \cdot u\), where \(u\) is either the first-order u.1 or the second-order u.2, see exp.uc in propagate.
Nonlinearity of the model:
Two measures of the nonlinearity of \(f\) at the input estimates, from the second-order Taylor expansion (see 'Details' of propagate):
The share of the second-order terms in the second-order variance with the input moments (moments extension),
$$\frac{u_{2m}^2(y) - u_1^2(y)}{u_{2m}^2(y)},$$
i.e. the part of the variance that the first-order law of propagation of uncertainty misses because of the curvature of \(f\). For Gaussian inputs \(u_{2m} = u_2\), and the share is \(\frac{1}{2}\mathrm{tr}(\mathbf{H\Sigma H\Sigma})/u_2^2(y)\), which is never negative because \(\mathrm{tr}(\mathbf{H\Sigma H\Sigma}) = \mathrm{tr}((\mathbf{\Sigma}^{1/2}\mathbf{H}\mathbf{\Sigma}^{1/2})^2) \ge 0\); with skewed inputs, the term \(c_i h_{ii} \gamma_{1,i} u_i^3\) can make it negative. It is zero for models that are linear at the estimates. Without u.2m (objects of earlier versions), u.2 is used.
The shift of the estimate by the second-order term, \((\mathrm{mean.2} - \mathrm{mean.1})/u_1(y) = \frac{1}{2}\mathrm{tr}(\mathbf{H\Sigma})/u_1(y)\), in units of the first-order standard uncertainty.
Both measures are exact for a quadratic model with independent inputs (the shift for any input distribution, the share with the input moments) and approximations otherwise; they do not include the third-derivative terms of the Note to JCGM 100:2008, 5.1.2, nor nonlinearity beyond second order (e.g. kinks of abs or min). Large values indicate that the first-order GUF may be inadequate; the validation by the Monte Carlo method (below) is the decisive criterion.
Uncertainty budget (first order):
For each input quantity the estimate \(x_i\), standard uncertainty \(u(x_i)\), the distribution used for the Monte Carlo simulation ("t", "normal" for degrees of freedom \(\ge 1000\), "mvt", "sample" for type = "sim"), the degrees of freedom, the sensitivity coefficient \(c_i\), the contribution \(u_i(y) = |c_i|\,u(x_i)\) (JCGM 100:2008, 5.1.3) and the index
$$I_i = \frac{c_i^2\,u^2(x_i)}{u_1^2(y)}, \qquad u_1^2(y) = \sum_i \sum_j c_i c_j \Sigma_{ij},$$
the share of input \(i\) in the first-order variance, with a bar of proportional length (block characters in UTF-8 locales, "#" otherwise). For correlated inputs, the covariance terms \(2\sum_{i<j} c_i c_j \Sigma_{ij}/u_1^2(y)\) are given in an additional row correlations, so that all rows sum to 100 %; negative covariance terms (cancellation) are shown with a bar of "-" and can make single indices exceed 100 %. The indices are shares of \(u_1^2(y)\) itself, as in the budget tables of EA-4/02; the relative contributions rel.contr of the propagate object, \(\mathrm{RC}_{ij} = |c_i c_j \Sigma_{ij}| / \sum_{k,\ell} |c_k c_\ell \Sigma_{k\ell}|\), are shares of the sum of absolute terms instead. If all first-order terms are zero (\(u_1(y) = 0\)), the budget is not defined, and a note is printed instead. Inputs with \(c_i = 0\) whose second derivatives do not vanish (e.g. \(x_2 = 0\) in \(Y = X_1^2 + X_2^2\)) contribute only through the second-order terms; a note points this out, as their index is 0. Variance-based shares from the Monte Carlo simulation are given by sobol.
Uncertainty budget (second order):
Printed for models with a non-zero Hessian matrix. It decomposes the second-order variance with the input moments (moments extension, u.2m of propagate),
$$u_{2m}^2(y) = \mathbf{c}^T\mathbf{\Sigma}\mathbf{c} + \frac{1}{2}\mathrm{tr}(\mathbf{H\Sigma H\Sigma}) + \sum_i \left(\frac{1}{4} h_{ii}^2 \gamma_{2,i} u^4(x_i) + c_i h_{ii} \gamma_{1,i} u^3(x_i)\right),$$
into the terms of the input quantities and of their interactions, with the skewness \(\gamma_{1,i}\) and excess kurtosis \(\gamma_{2,i}\) of each input (component moments of the propagate object; both 0 for Gaussian inputs, where \(u_{2m} = u_2\)). For independent inputs,
$$u_{2m}^2(y) = \sum_i \left[c_i^2 u^2(x_i) + \frac{1}{4}h_{ii}^2 (\gamma_{2,i} + 2) u^4(x_i) + c_i h_{ii} \gamma_{1,i} u^3(x_i)\right] + \sum_{i<j} h_{ij}^2 u^2(x_i) u^2(x_j),$$
which are the first-order (main effect) and interaction terms of the variance-based (Sobol') decomposition of the second-order Taylor model (JCGM 100:2008, 5.1.2, note). The table gives for each input the first-order term \(c_i^2 u^2(x_i)\), the second-order term (\(\frac{1}{2}h_{ii}^2 u^4(x_i)\) for Gaussian inputs), their sum and its share Index of \(u_{2m}^2(y)\) with a bar, next to the first-order index (Index 1st) for comparison; each non-zero interaction \(h_{ij}^2 u^2(x_i) u^2(x_j)\) is given in a row x_i x x_j, and interactions with shares below 0.05 % are pooled in the row other interactions. For correlated inputs, all covariance terms of first and second order are given in the row correlations, as they cannot be uniquely allocated to the inputs. The table thus shows inputs that the first-order budget misses: in \(Y = X_1^2 + X_2^2\) with \(x_1 = 0.001\), \(x_2 = 0\), \(u(x_1) = u(x_2) = 0.005\) and independent inputs, \(X_2\) has a first-order index of 0 % but a second-order index of 48 %. If the second-order terms are less than 0.1 % of \(u_{2m}^2(y)\), only a note is printed. With the input moments, the shares agree with the Sobol' indices of sobol for quadratic models with independent inputs (heater example of the vignette: 64 % for the rectangular phase angle; with Gaussian moments, as before version 1.3-0, it would be 81 %). Correlated pairs keep Gaussian joint moments, and the third-derivative term of JCGM 100:2008, 5.1.2, is not included; for strongly nonlinear models, the Monte Carlo results are decisive.
Validation of the GUM uncertainty framework by the Monte Carlo method (validate = TRUE, JCGM 101:2008, 8):
The endpoints of the GUF coverage interval \(y \pm U_p\) are compared with those of the Monte Carlo coverage interval \([y_{\mathrm{low}}, y_{\mathrm{high}}]\) selected by interval (probabilistically symmetric by default):
$$d_{\mathrm{low}} = |y - U_p - y_{\mathrm{low}}|, \quad d_{\mathrm{high}} = |y + U_p - y_{\mathrm{high}}|$$
The numerical tolerance \(\delta\) is derived from the GUF standard uncertainty \(u(y) = U_p/k\): expressing \(u(y) = c \times 10^l\), with \(c\) an integer of ndig digits, gives \(\delta = \frac{1}{2} \times 10^l\) (JCGM 101:2008, 7.9.2). The GUF is validated if both \(d_{\mathrm{low}}\) and \(d_{\mathrm{high}}\) are not larger than \(\delta\); otherwise, the Monte Carlo results should be used. The comparison is only conclusive if the Monte Carlo coverage interval endpoints are themselves numerically stable to within \(\delta\) (JCGM 101:2008, 7.9). Therefore, the standard deviation of the full-sample endpoints, \(s_{\mathrm{MC}} = \mathrm{sd}(\mathrm{block\;endpoints})/\sqrt{h}\), is taken from the \(h\) blocks of propagate (sim.se: \(h = 10\) equal parts of the results, or the blocks of the adaptive procedure; the larger of the two values, for the lower and the upper endpoint, is reported as s.MC). If sim.se is not available, the stability cannot be assessed, and a note is printed. If \(2 s_{\mathrm{MC}} > \delta\), the result is reported as "INCONCLUSIVE", and nsim should be increased. If \(u(y) = 0\), which happens when all sensitivity coefficients vanish at the input means, the GUF is reported as not validated. If inputs are sampled with \(\nu \le 4\) degrees of freedom, the output has no finite fourth moment: u.MC and the Monte Carlo correlations then converge slowly, and a note recommends the coverage interval endpoints as the reliable Monte Carlo results.
Derivatives and contributions:
The symbolic gradient \(\nabla\), also often termed \(\mathbf{G}\), \(\mathbf{J}\) or \(\mathbf{C_x}\) (JCGM 102:2011).
The evaluated gradient, also known as "sensitivity coefficients", \(c_i = \partial f / \partial x_i\), evaluated at the input means.
The symbolic Hessian matrix \(\mathbf{H}\), also often termed \(\nabla^2\).
The evaluated Hessian matrix \(h_{ij} = \partial^2 f / \partial x_i \partial x_j\), evaluated at the input means.
Checks of the Monte Carlo input (extending JCGM 101:2008):
The input covariance matrix \(\mathbf{\Sigma}\) of the GUF, as provided from (i) the uncertainties in the 2nd row of the input data, (ii) estimation from the raw data, or (iii) an externally supplied one. Also often termed \(\mathbf{V}\) or \(\mathbf{U_x}\) (JCGM 102:2011).
The target covariance matrix \(\mathbf{\Sigma}^{\text{target}}\) that the Monte Carlo input should have under the chosen mc.margins model (covMat.target, see propagate), printed only if it differs from \(\mathbf{\Sigma}\). For "S1", its variances are larger than those of \(\Sigma\) by the factors \(\nu_i/(\nu_i-2)\), for "mvt" by \((n-1)/(n-N-2)\).
The covariance matrix estimated from the Monte Carlo input, \(\hat{\mathbf{\Sigma}}^\text{MC} = \mathrm{cov}(\mathbf{X}^\text{MC})\).
The relative Frobenius norm \(\varepsilon_F = \|\hat{\mathbf{\Sigma}}^\text{MC} - \mathbf{\Sigma}^{\text{target}}\|_F / \|\mathbf{\Sigma}^{\text{target}}\|_F\), together with its ratio to the Monte Carlo noise level (see below). "not applicable" is printed if the target covariance does not exist (\(\nu \le 2\)).
For correlated inputs only: the correlation matrix \(\mathbf{R}\) derived from \(\mathbf{\Sigma}\) (by cov2cor), used as the parameter of the Gaussian copula for "S1".
For correlated inputs only: the relative Frobenius norm \(\varepsilon_F = \|\hat{\mathbf{R}}^\text{MC} - \mathbf{R}\|_F / \|\mathbf{R}\|_F\) of the correlation matrix \(\hat{\mathbf{R}}^\text{MC} = \mathrm{cor}(\mathbf{X}^\text{MC})\) of the Monte Carlo input, together with its ratio to the Monte Carlo noise level. Being scale-free, it detects a wrong dependence structure also between inputs with small variances, which hardly affect the covariance norm.
A comparison of the input parameters (means, standard uncertainties \(u_i\), df.MC) with those of scaled and shifted t-distributions fitted to each margin of the Monte Carlo input (only if check = TRUE in propagate). The fitted scale is reported on the scale of \(u_i\).
Monte Carlo noise level: \(\varepsilon_F\) is never exactly zero, because \(\hat{\mathbf{\Sigma}}^\text{MC}\) and \(\hat{\mathbf{R}}^\text{MC}\) are estimates from nsim draws. The expected size of this sampling noise depends on the number of inputs, on nsim and on the tails of the margins, and is estimated empirically: the Monte Carlo input is split into two halves, the matrix is estimated from each half (\(\hat{\mathbf{\Sigma}}^{(1)}, \hat{\mathbf{\Sigma}}^{(2)}\)), and
$$\varepsilon_{\mathrm{noise}} = \frac{\|\hat{\mathbf{\Sigma}}^{(1)} - \hat{\mathbf{\Sigma}}^{(2)}\|_F}{2\,\|\mathbf{\Sigma}^{\text{target}}\|_F}$$
approximates the sampling noise of the full-sample estimate, since each half-sample estimate has twice the sampling variance of the full-sample one. The ratio \(\varepsilon_F / \varepsilon_{\mathrm{noise}}\) is printed, with "OK" if it is at most 3 and "CHECK" otherwise; the same is done for the correlation matrices. A ratio above 3 indicates a systematic deviation of the Monte Carlo input from the target. Being based on a single split, \(\varepsilon_{\mathrm{noise}}\) is itself a rough estimate, and occasional ratios slightly above 3 can occur by chance. For correlated inputs joined by the Gaussian copula ("S1") with small \(\nu_i\), a small systematic deviation is expected, because the Pearson correlation of t-margins is slightly smaller in absolute value than the copula parameter; it becomes visible as "CHECK" for large nsim. If any \(\nu_i \le 4\), the t-margins have infinite kurtosis, covariance estimates converge slowly and erratically, and a note is printed that the covariance check is unreliable.
Shape of the Monte Carlo output distribution:
Skewness and excess kurtosis of the model evaluations \(Y = f(\mathbf{X}^\text{MC})\), see skewness.
If normality = TRUE: Shapiro-Wilk test (2000 repeats, 5000 samples).
If normality = TRUE: Anderson-Darling test (1000 repeats, 10000 samples).
If normality = TRUE: Lilliefors test (1000 repeats, 10000 samples).
For each test, the median p-value over the repeats is reported (see Thode 2002 for the tests, and Royston 1992 for the Shapiro-Wilk test of large samples). Non-finite model evaluations are removed before subsampling. With large subsamples, these tests detect even small, practically irrelevant deviations from normality; the validation of the GUF (see above) is the more relevant criterion for deciding whether the GUF coverage interval can be used.
Evaluation of measurement data - Guide to the expression of uncertainty in measurement.
JCGM 100:2008. Joint Committee for Guides in Metrology (2008).
tools:::Rd_expr_doi("10.59161/JCGM100-2008E")
Evaluation of measurement data - Supplement 1 to the "Guide to the expression of uncertainty in measurement" - Propagation of distributions using a Monte Carlo method.
JCGM 101:2008. Joint Committee for Guides in Metrology (2008).
tools:::Rd_expr_doi("10.59161/JCGM101-2008")
Evaluation of measurement data - Supplement 2 to the "Guide to the expression of uncertainty in measurement" - Extension to any number of output quantities.
JCGM 102:2011. Joint Committee for Guides in Metrology (2011).
tools:::Rd_expr_doi("10.59161/JCGM102-2011")
Testing for Normality.
Thode Jr. HC.
Marcel Dekker (2002), New York.
Approximating the Shapiro-Wilk W-test for non-normality.
Royston P.
Stat Comp (1992), 2: 117-119.
For further references, see propagate.
propagate, plot.propagate, WelchSatter, report, sobol.
## GUM Supplement 1, 9.4, comparison loss with x1 = 0.010:
## asymmetric output distribution. The GUF is not validated
## (JCGM 101:2008, 8), and the shortest coverage interval
## differs from the probabilistically symmetric one.
## Use nsim = 1E6 in practice
EXPR1 <- expression(X1^2 + X2^2)
X1 <- c(0.010, 0.005)
X2 <- c(0, 0.005)
DF1 <- cbind(X1, X2)
RES1 <- propagate(expr = EXPR1, data = DF1, type = "stat",
nsim = 1E5, check = FALSE)
summary(RES1)
Run the code above in your browser using DataLab