From uncertainty budget to conformity decision: a complete workflow with propagate

Andrej-Nikolai Spiess

2026-11-10

1 Overview

This vignette follows one measurement through the complete evaluation with propagate, from the uncertainty budget to the decision whether the measured item conforms to its specification. Each step uses one function of the package:

Step Function Purpose Reference
1 makeDat() input quantities from an uncertainty budget JCGM 101:2008, 6.4
2 propagate() GUM uncertainty framework (GUF) and Monte Carlo method (MCM) JCGM 100:2008, JCGM 101:2008
3 summary() nonlinearity, uncertainty budgets, validation of the GUF JCGM 101:2008, 8
4 sobol() variance-based sensitivity of the Monte Carlo output Saltelli et al. (2010)
5 plot() distribution of the Monte Carlo output
6 fitDistr() a closed-form distribution for the output
7 conformity() conformance probability and decision rules JCGM 106:2012, ISO 14253-1
8 report() PDF report of the evaluation

The example is deliberately nonlinear: one input quantity has a sensitivity coefficient of zero, so that the first-order GUF misses its contribution completely, and the GUF and the Monte Carlo method lead to different conformity decisions. Section 10 then shows three further functions on their own examples: predictNLS(), propagateMV() and doGauge().

2 The measurement

The electrical power of a heating element with a rated power of 1000 W is determined from the voltage \(V\), the current \(I\) and the phase angle \(\varphi\) between them:

\[ P = V\,(I + \delta I)\,\cos\varphi . \]

The phase angle makes the model nonlinear. Its sensitivity coefficient at the estimate \(\varphi = 0\) is \(\partial P/\partial\varphi = -V(I + \delta I)\sin\varphi = 0\), so the first-order GUF assigns it no contribution at all. Yet \(\cos\varphi \le 1\): every phase angle different from zero reduces the power, so that the distribution of \(P\) is skewed to the left and its expectation lies below \(V I\).

The customer specifies the power as 1000 W \(\pm\) 0.6 %, i.e. the tolerance interval is [994, 1006] W.

3 Step 1: input quantities with makeDat()

The ten current readings give the Type A estimate and its standard uncertainty:

I.obs <- c(4.358, 4.349, 4.353, 4.361, 4.347, 4.356, 4.351, 4.359, 4.352, 4.355)
c(mean = mean(I.obs), u = sd(I.obs)/sqrt(length(I.obs)))
#>       mean          u 
#> 4.35410000 0.00142556

The uncertainty budget is written as a data frame with one row per input quantity, as in a spreadsheet. makeDat() understands the usual columns of an uncertainty budget: the estimate (mean), a standard uncertainty (u), or an expanded uncertainty with its coverage factor (U, k), or the half-width of a rectangular, triangular or arcsine distribution (a), the degrees of freedom (df) and the distribution (dist). The same table can be read from a .csv or .xlsx file with makeDat(file = ...).

BUD <- data.frame(name = c("V", "I", "dI", "phi"),
                  mean = c(230.1, mean(I.obs), 0, 0),
                  u    = c(NA, sd(I.obs)/sqrt(10), NA, NA),
                  U    = c(0.4, NA, NA, NA),
                  k    = c(2, NA, NA, NA),
                  a    = c(NA, NA, 0.01, 0.12),
                  df   = c(Inf, 9, Inf, Inf),
                  dist = c("normal", "t", "rect", "rect"))
BUD
#>   name     mean          u   U  k    a  df   dist
#> 1    V 230.1000         NA 0.4  2   NA Inf normal
#> 2    I   4.3541 0.00142556  NA NA   NA   9      t
#> 3   dI   0.0000         NA  NA NA 0.01 Inf   rect
#> 4  phi   0.0000         NA  NA NA 0.12 Inf   rect

With output = "sim", makeDat() draws nsim samples from the stated distribution of each input quantity, which are then used by the Monte Carlo method. The budget itself is kept as the attribute "budget", so that the GUF uses the stated estimates and standard uncertainties (not the sample moments):

set.seed(123)
EXPR <- expression(P = V * (I + dI) * cos(phi))
DAT  <- makeDat(EXPR, data = BUD, layout = "budget", output = "sim", nsim = M)
#> makeDat: 100000 x 4 data => use propagate(..., type = "sim").
head(DAT, 3)
#>             V        I           dI        phi
#> [1,] 229.9879 4.354375  0.007310058 0.06146650
#> [2,] 230.0540 4.354009 -0.008325165 0.08831624
#> [3,] 230.4117 4.352753 -0.006967282 0.11089883
attr(DAT, "budget")[, c("name", "mean", "u", "df", "dist")]
#>   name     mean           u  df   dist
#> 1    V 230.1000 0.200000000 Inf normal
#> 2    I   4.3541 0.001425560   9      t
#> 3   dI   0.0000 0.005773503 Inf   rect
#> 4  phi   0.0000 0.069282032 Inf   rect

The standard uncertainties were derived as \(u(V) = U/k = 0.2\) V, \(u(\delta I) = a/\sqrt{3} = 0.0058\) A and \(u(\varphi) = a/\sqrt{3} = 0.069\) rad.

Two other ways of creating the input data are described in ?makeDat and ?propagate: for repeated observations, the raw data (one column per quantity, type = "raw"); for correlated input quantities with arbitrary distributions, corDist() (Gaussian copula).

4 Step 2: propagation with propagate()

PROP <- propagate(EXPR, DAT, type = "sim", check = FALSE)
#> I see data with 1000 <= rows ... type = 'sim' seems correct!
#> propagate: the GUF uses the estimates, standard uncertainties and degrees of freedom of the uncertainty budget (attribute 'budget').

PROP$prop holds the results of the GUF: the estimate and standard uncertainty from the first-order (mean.1, u.1) and second-order (mean.2, u.2) Taylor expansion, the second-order uncertainty with the moments of the input distributions (u.2m, see below), and the coverage interval of the first-order GUF. PROP$sim holds the Monte Carlo results with the probabilistically symmetric 95 % coverage interval, PROP$ci.short the shortest coverage interval:

PROP$prop
#>      mean.1      mean.2         u.1         u.2        u.2m        2.5% 
#> 1001.878410  999.473902    1.621972    3.767508    2.693719  998.698688 
#>       97.5% 
#> 1005.058132
PROP$sim
#>     mean.MC        u.MC   median.MC      mad.MC        2.5%       97.5% 
#>  999.480401    2.696131  999.711982    2.821138  993.893076 1004.090294
PROP$ci.short
#>      2.5%     97.5% 
#>  994.1038 1004.2650

The first-order GUF gives \(P = 1001.9\) W with \(u = 1.6\) W and a symmetric interval. The Monte Carlo method gives an estimate 2.4 W lower, a standard uncertainty two-thirds larger, and an asymmetric interval that extends further to low powers. The second-order Taylor expansion (mean.2) already reproduces the shift of the estimate: for a rectangular \(\varphi\) on \([-a, a]\), \(E(\cos\varphi) \approx 1 - a^2/6\), i.e. the expected power is lower by about \(P\,a^2/6\):

c(GUF.2nd = unname(PROP$prop["mean.2"] - PROP$prop["mean.1"]),
  MCM     = unname(PROP$sim["mean.MC"] - PROP$prop["mean.1"]),
  approx  = unname(-PROP$prop["mean.1"] * 0.12^2/6))
#>   GUF.2nd       MCM    approx 
#> -2.404508 -2.398009 -2.404508

For final results, nsim = 1E6 or the adaptive Monte Carlo procedure of JCGM 101:2008, 7.9, are recommended, e.g. propagate(EXPR, DAT, type = "sim", adaptive = TRUE), which redraws the input samples from the budget block by block until the results are numerically stable.

5 Step 3: summary() – nonlinearity, budgets and validation

summary() collects all diagnostics of the evaluation:

summary(PROP)
#> 
#> Monte Carlo input model:
#> user-supplied samples (type = "sim")
#> 
#> Results from uncertainty propagation:
#>   mean.1   mean.2      u.1      u.2     u.2m     2.5%    97.5% 
#> 1002.000  999.500    1.622    3.768    2.694  998.700 1005.000 
#> 
#> Results from Monte Carlo simulation (probabilistically symmetric coverage interval, JCGM 101:2008, 7.7.1):
#>   mean.MC      u.MC median.MC    mad.MC      2.5%     97.5% 
#>   999.500     2.696   999.700     2.821   993.900  1004.000 
#> 
#> Shortest coverage interval from Monte Carlo simulation (JCGM 101:2008, 7.7.2):
#>   2.5%  97.5% 
#>  994.1 1004.0 
#> 
#> Numerical standard uncertainties of the Monte Carlo results (block standard deviations / sqrt(h)):
#>    mean.MC       u.MC    sym.low   sym.high  short.low short.high 
#>    0.00753    0.00644    0.01630    0.01230    0.05390    0.04010 
#> 
#> Standard uncertainty: Taylor expansion (1st order, 2nd order Gaussian, 2nd order input moments) vs. Monte Carlo:
#>                        u.1     u.2   u.2m   u.MC
#> u                    1.622   3.768  2.694  2.696
#> difference to u.MC -39.8 % +39.7 % -0.1 % +0.0 %
#> u.2m: second order with the skewness and excess kurtosis of each input (object$moments), exact for quadratic
#> models with independent inputs; u.2m - u.MC = higher-order terms of the model. mean.2 does not depend on these moments.
#> 
#> Welch-Satterthwaite degrees of freedom: 5380 
#> 
#> Coverage factor (k): 1.96 
#> 
#> Expanded uncertainty: 3.17972 
#> 
#> Nonlinearity of the model (second-order Taylor terms):
#> Share of the second-order terms in u.2m^2: 63.7 %
#> Shift of the mean (mean.2 - mean.1)/u.1:    -1.48
#> => large values indicate that the first-order GUF may be inadequate (see the validation below)
#> 
#> Uncertainty budget (first order):
#> Quantity    Value  Std. unc.  Distribution  DOF  Sens. coeff.  Contribution    Index                           
#> ---------------------------------------------------------------------------------------------------------------
#> V         230.100        0.2  sample        Inf         4.354        0.8708   28.8 %  ███████                  
#> I         4.35410   0.001426  sample          9         230.1         0.328    4.1 %  █                        
#> dI        0.00000   0.005774  sample        Inf         230.1         1.328   67.1 %  █████████████████        
#> phi        0.0000    0.06928  sample        Inf             0             0    0.0 %                           
#> ---------------------------------------------------------------------------------------------------------------
#> Contribution = |c_i| u(x_i); Index = c_i^2 u^2(x_i)/u^2(y)
#> Note: phi with c_i = 0 contribute(s) only through second-order terms (see the second-order budget below; for MC-based shares use sobol()).
#> 
#> Uncertainty budget (second order):
#> Quantity                  1st order    2nd order     Variance      Index   Index 1st                           
#> ---------------------------------------------------------------------------------------------------------------
#> V                         7.583e-01            0    7.583e-01     10.5 %      28.8 %  ███                      
#> I                         1.076e-01            0    1.076e-01      1.5 %       4.1 %                           
#> dI                        1.765e+00            0    1.765e+00     24.3 %      67.1 %  ██████                   
#> phi                               0    4.625e+00    4.625e+00     63.7 %       0.0 %  ████████████████         
#> other interactions (2)            0    1.415e-06    1.415e-06      0.0 %                                       
#> ---------------------------------------------------------------------------------------------------------------
#> 1st order = c_i^2 u^2(x_i); 2nd order = 1/4 H_ii^2 (g2_i + 2) u^4(x_i) + c_i H_ii g1_i u^3(x_i), interactions H_ij^2 u^2(x_i) u^2(x_j);
#> g1, g2 = skewness and excess kurtosis of x_i (object$moments; 0 for Gaussian inputs, then 2nd order = 1/2 H_ii^2 u^4(x_i));
#> Index = Variance/u.2m^2 (u.2m^2 = 7.256e+00); Index 1st = first-order index of the table above
#> other interactions = 2 interaction term(s) each < 0.05 % of u.2m^2
#> Note: exact for quadratic models with independent inputs; correlated pairs use Gaussian joint moments. For other models, use the Monte Carlo results.
#> 
#> Validation of GUM uncertainty framework by MCM (JCGM 101:2008, 8), symmetric MC interval:
#>  d.low d.high  delta   s.MC 
#> 4.8100 0.9680 0.0500 0.0163 
#> => GUF NOT validated for 2 significant digits of u(y); use the MCM results
#> 
#> Symbolic gradient matrix: (I + dI) * cos(phi) V * cos(phi) V * cos(phi) -(V * (I + dI) * sin(phi))
#> Evaluated gradient matrix (Sensitivities): 4.3541 230.1 230.1 0 
#> 
#> Symbolic Hessian matrix:
#> [1] "0"                      "cos(phi)"               "cos(phi)"              
#> [4] "-((I + dI) * sin(phi))"
#> [1] "cos(phi)"        "0"               "0"               "-(V * sin(phi))"
#> [1] "cos(phi)"        "0"               "0"               "-(V * sin(phi))"
#> [1] "-((I + dI) * sin(phi))"     "-(V * sin(phi))"           
#> [3] "-(V * sin(phi))"            "-(V * (I + dI) * cos(phi))"
#> 
#> Evaluated Hessian matrix:
#> [1] 0 1 1 0
#> [1] 1 0 0 0
#> [1] 1 0 0 0
#> [1]     0.000     0.000     0.000 -1001.878
#> 
#> Input covariance (GUF):
#>        V        I       dI    phi
#> V   0.04 0.00e+00 0.00e+00 0.0000
#> I   0.00 2.03e-06 0.00e+00 0.0000
#> dI  0.00 0.00e+00 3.33e-05 0.0000
#> phi 0.00 0.00e+00 0.00e+00 0.0048
#> MC design covariance, observed:
#>             V         I       dI       phi
#> V    4.00e-02  2.20e-07 5.82e-06 -5.38e-05
#> I    2.20e-07  2.63e-06 1.75e-08 -5.14e-07
#> dI   5.82e-06  1.75e-08 3.33e-05  6.38e-07
#> phi -5.38e-05 -5.14e-07 6.38e-07  4.80e-03
#> Relative Frobenius norm, expected vs. observed:
#> 0.00197 (1.54 x MC noise level => OK) 
#> 
#> Skewness / Excess Kurtosis of MC evaluations:
#> -0.312876 / -0.4211226

The output shows, from top to bottom:

6 Step 4: sensitivity analysis with sobol()

The second-order budget uses the moments of the input distributions: the variance of \(\varphi^2\) is \(4a^4/45\) for a rectangular distribution, but \(2u^4 = 2a^4/9\) for a Gaussian one with the same standard uncertainty, a ratio of 0.4 (the Gaussian assumption of u.2 would give the phase angle 81 %). The variance-based (Sobol’) sensitivity indices of the Monte Carlo output are independent of the Taylor expansion:

S <- sobol(PROP)
#> sobol: 'sim' data => using rows 1..50000 as A and 50001..1e+05 as B (rows must be independent draws).
round(rbind(first = S$S, total = S$ST), 3)
#>           V     I    dI   phi
#> first 0.107 0.024 0.247 0.632
#> total 0.104 0.019 0.242 0.635

The phase angle accounts for about 63 % of the variance of \(P\), the clamp meter correction for 25 %. This agrees with the exact second-order contribution of the rectangular \(\varphi\), \(P^2\,\mathrm{Var}(\varphi^2/2) = P^2 a^4/45\):

P0 <- unname(PROP$prop["mean.1"]); a <- 0.12
Vphi <- P0^2 * a^4/45                         # exact for a rectangular phase angle
c(share.phi = Vphi/(unname(PROP$prop["u.1"])^2 + Vphi), MC.variance = unname(PROP$sim["u.MC"])^2,
  first.order.plus.phi = unname(PROP$prop["u.1"])^2 + Vphi)
#>            share.phi          MC.variance first.order.plus.phi 
#>            0.6374381            7.2691205            7.2561202

The three budgets thus tell a consistent story: the first-order budget ignores the phase angle (0 %), the second-order budget with the input moments detects it (64 %), and the Monte Carlo indices confirm this for the stated rectangular distribution (63 %).

7 Step 5: plots

plot() shows the distribution of the Monte Carlo output together with the coverage intervals:

plot(PROP)

8 Step 6: a distribution for the output with fitDistr()

To pass the result on to a further evaluation, e.g. as an input quantity of another model, a closed-form distribution is often convenient. fitDistr() fits all distributions of the package to the Monte Carlo output and ranks them by the Bayesian information criterion (BIC). nmax limits the number of values used for the fit:

FIT <- fitDistr(PROP, nmax = 5000, verbose = FALSE)
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in base::log((w + d)/(w - d)): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced
#> Warning in dt(z, df = nu, log = TRUE): NaNs produced

head(FIT$stat[, c("Distribution", "k", "BIC", "dBIC", "weight")], 5)
#>     Distribution k      BIC       dBIC       weight
#> 1 Normal mixture 5 23964.39  0.0000000 5.184729e-01
#> 2     Johnson SB 4 23964.53  0.1484051 4.813936e-01
#> 3        4P Beta 4 23980.92 16.5369789 1.329740e-04
#> 4     3P Weibull 3 23991.83 27.4464837 5.685882e-07
#> 5  Skewed-normal 3 24011.45 47.0592735 3.132814e-11

The skewed, bounded output is best described by a Johnson SB or a two-component normal mixture; the normal distribution, implicitly assumed by the GUF, ranks far behind. The parameters of the best fit are in FIT$bestpar; plot = "qq" draws a QQ plot instead of the histogram.

9 Step 7: conformity assessment with conformity()

The tolerance interval of the customer is [994, 1006] W. conformity() computes the conformance probability \(p_c\), the probability that the true power lies within the tolerance interval (JCGM 106:2012, 7), and, with rule, the decisions of the decision rules of ISO 14253-1 and JCGM 106 (rule = "all" gives all six combinations of stringent, simple and relaxed acceptance and rejection):

CMC <- conformity(PROP, lower = 994, upper = 1006, rule = "all")
CMC
#> Conformity assessment of a measured item (JCGM 106:2012, 7), MC
#> Estimate y = 999.4804, u = 2.696, tolerance interval [994, 1006]
#> Conformance probability pc = 0.9718, non-conformance 1 - pc = 0.02816
#> Decision rules with the expanded uncertainty U- = 5.587, U+ = 4.61 (p = 0.95):
#> acceptance/rejection               accept             reject below | above  decision   
#> stringent/stringent (ISO 14253-1)  [999.59, 1001.39]  989.39 | 1011.59      NO DECISION
#> simple/simple                      [994.00, 1006.00]  994.00 | 1006.00      ACCEPT     
#> stringent/relaxed                  [999.59, 1001.39]  999.59 | 1001.39      REJECT     
#> relaxed/stringent                  [989.39, 1011.59]  989.39 | 1011.59      ACCEPT     
#> stringent/simple                   [999.59, 1001.39]  994.00 | 1006.00      NO DECISION
#> simple/stringent                   [994.00, 1006.00]  989.39 | 1011.59      ACCEPT

From the Monte Carlo results, \(p_c = 0.97\). The stringent acceptance zone of ISO 14253-1 is the tolerance interval reduced by the expanded uncertainty on each side; because the Monte Carlo interval is asymmetric (\(U_- = 5.6\) W, \(U_+ = 4.6\) W), the zone is [999.6, 1001.4] W. The estimate of 999.5 W lies just outside of it, and as it also lies inside the stringent rejection limits, ISO 14253-1 gives no decision: the heater can neither be proven to conform nor to not conform.

plot(CMC)

From the first-order GUF, the same assessment reaches a different conclusion:

CGUF <- conformity(PROP, lower = 994, upper = 1006, method = "GUF", rule = c("stringent", "simple"))
CGUF
#> Conformity assessment of a measured item (JCGM 106:2012, 7), GUF
#> Estimate y = 1001.878, u = 1.622, nu = 5380, tolerance interval [994, 1006]
#> Conformance probability pc = 0.9945, non-conformance 1 - pc = 0.005526
#> Decision rules with the expanded uncertainty U- = 3.179, U+ = 3.179 (p = 0.95):
#> acceptance/rejection               accept             reject below | above  decision
#> stringent/stringent (ISO 14253-1)  [997.18, 1002.82]  990.82 | 1009.18      ACCEPT  
#> simple/simple                      [994.00, 1006.00]  994.00 | 1006.00      ACCEPT

The GUF places the estimate at 1001.9 W with a smaller, symmetric uncertainty, so that \(p_c = 0.99\) and the heater is accepted under the stringent rule. Since the GUF was not validated (Section 5), this decision is not reliable: the nonlinearity of the model changes the outcome of the conformity assessment.

10 Step 8: the report with report()

report() writes the evaluation to a PDF file with two pages: the measurement model, the input data, the result statement of the GUF and the MCM, the results of all methods, the sensitivity coefficients and the Hessian matrix, the nonlinearity of the model, the uncertainty budget, a histogram of the Monte Carlo output, the covariance check and the shape of the output distribution. Published or certified values can be compared with gauge.mean and gauge.u.

f <- file.path(tempdir(), "heater-power.pdf")
REP <- report(PROP, file = f, title = "Power of a 1000 W heating element",
              measurand = "P", unit = "W")
REP$nonlinearity
#>                                                                     Measure
#> 1                                 share of the second-order terms in u.2m^2
#> 2                               shift of the estimate (mean.2 - mean.1)/u.1
#> 3 quantile skewness of the MC output, (q97.5 + q2.5 - 2 q50)/(q97.5 - q2.5)
#> 4                          skewness of the MC output (moments; Gaussian: 0)
#> 5                   excess kurtosis of the MC output (moments; Gaussian: 0)
#>    Value Assessment
#> 1 63.7 %      large
#> 2  -1.48      large
#> 3 -0.141   moderate
#> 4 -0.313   moderate
#> 5 -0.421   moderate

A list of propagate objects gives one report with two pages per evaluation, e.g. for all examples of a study: report(list("Heater A" = PROP1, "Heater B" = PROP2), file = ...).

11 Further functions

11.1 predictNLS(): uncertainty of predictions from a nonlinear regression

predictNLS() propagates the uncertainty of the fitted parameters of a nonlinear model (nls), and optionally of the predictor values, to the predicted response. The example is the logistic model of the DNase data of ?nls:

DNase1 <- subset(DNase, Run == 1)
fm <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = DNase1,
          start = list(Asym = 3, xmid = 0, scal = 1))

A prediction at a concentration of 2, first with an exactly known concentration, then with \(u(\mathrm{conc}) = 0.2\):

set.seed(123)
PR1 <- predictNLS(fm, newdata = data.frame(conc = 2), nsim = M)
PR2 <- predictNLS(fm, newdata = data.frame(conc = 2), newerror = data.frame(conc = 0.2), nsim = M)
round(rbind("u(conc) = 0" = PR1$summary[1, ], "u(conc) = 0.2" = PR2$summary[1, ]), 4)
#>               mean.1 mean.2    u.1    u.2   2.5%  97.5% mean.MC   u.MC
#> u(conc) = 0    0.748 0.7482 0.0080 0.0081 0.7307 0.7654  0.7481 0.0088
#> u(conc) = 0.2  0.748 0.7466 0.0496 0.0496 0.6509 0.8452  0.7466 0.0499
#>               median.MC mad.MC   2.5%  97.5%
#> u(conc) = 0      0.7481 0.0083 0.7306 0.7654
#> u(conc) = 0.2    0.7482 0.0495 0.6447 0.8403

The first-order mean equals predict(fm); the uncertainty of the concentration increases the standard uncertainty of the prediction about sixfold. Over a range of concentrations, the confidence band of the fitted curve (interval = "confidence") and the prediction band for a new observation (interval = "prediction") are obtained in the same way:

CONC <- exp(seq(log(0.05), log(12.5), length.out = 25))
CB <- predictNLS(fm, newdata = data.frame(conc = CONC), nsim = M)$summary
PB <- predictNLS(fm, newdata = data.frame(conc = CONC), interval = "prediction", nsim = M)$summary
plot(density ~ conc, data = DNase1, log = "x", pch = 16, xlab = "concentration", ylab = "optical density")
lines(CONC, CB[, "mean.1"], lwd = 2)
matlines(CONC, CB[, 11:12], lty = 2, col = 2)      # Monte Carlo confidence band
matlines(CONC, PB[, 11:12], lty = 3, col = 4)      # Monte Carlo prediction band
legend("topleft", c("fit", "95 % confidence band", "95 % prediction band"), lty = 1:3, col = c(1, 2, 4), bty = "n")

11.2 propagateMV(): several output quantities

propagateMV() evaluates several measurement models on the same input quantities and gives the joint distribution of the output quantities (JCGM 102:2011): their estimates, standard uncertainties and correlations, and the coverage factors of hyper-ellipsoidal (\(k_p\)) and hyper-rectangular (\(k_q\)) 95 % coverage regions. The example is JCGM 102:2011, 9.4: resistance \(R\), reactance \(X\) and impedance \(Z\) from six sets of simultaneous observations of voltage, current and phase angle, with the multivariate \(t\)-distribution of the inputs (JCGM 102:2011, 5.3.2, \(\nu = n - N = 3\)):

H6 <- cbind(V   = c(5.007, 4.994, 5.005, 4.990, 4.999, 4.999),
            I   = c(19.663, 19.639, 19.640, 19.685, 19.678, 19.661)/1000,
            phi = c(1.0456, 1.0438, 1.0468, 1.0428, 1.0433, 1.0445))
EXMV <- expression(R = (V/I) * cos(phi), X = (V/I) * sin(phi), Z = V/I)
set.seed(123)
MV <- propagateMV(EXMV, H6, type = "raw", cov = "full", mc.margins = "mvt", nsim = M)
MV
#> GUM uncertainty framework (JCGM 102:2011, 6):
#>                     R           X           Z
#> estimate 127.73070000 219.8474000 254.2597000
#> u          0.05804901   0.2413427   0.1929676
#> Correlation matrix:
#>         R       X       Z
#> R  1.0000 -0.5883 -0.4851
#> X -0.5883  1.0000  0.9925
#> Z -0.4851  0.9925  1.0000
#> Coverage factors for 95 % coverage regions: k_p (ellipsoid) = 2.795, k_q (rectangle) = 2.388
#> 
#> Monte Carlo method (JCGM 102:2011, 7), M = 100 000, input model mvt:
#>                    R           X           Z
#> estimate 127.7298000 219.8502000 254.2619000
#> u          0.1316133   0.5353225   0.4269184
#> Correlation matrix:
#>         R       X       Z
#> R  1.0000 -0.5958 -0.4922
#> X -0.5958  1.0000  0.9924
#> Z -0.4922  0.9924  1.0000
#> Coverage factors for 95 % coverage regions: k_p (ellipsoid) = 2.579, k_q (rectangle) = 2.238

The GUF and the MCM agree on the estimates and the correlations, but the Monte Carlo standard uncertainties are more than twice as large: the multivariate \(t\)-distribution with three degrees of freedom has the covariance matrix \(\frac{n-1}{n-N-2}\,\mathbf{S}/n = 5\,\mathbf{S}/n\), as in JCGM 102:2011, Table 11. plot(MV) shows the pairwise scatter plots of the output quantities with the coverage ellipses:

plot(MV)

11.3 doGauge(): validation against published results

doGauge() runs the demos of the package, which reproduce the worked examples of JCGM 100:2008 (GUM, Annex H), JCGM 101:2008, JCGM 102:2011, the EMUE compendium, NIST TN 1900 and the EURACHEM/CITAC guide QUAM, and compares each result with the published value. Every comparison is classified as PASS (agreement within the rounding of the published value, plus four Monte Carlo standard errors for Monte Carlo results), NOTE (a documented discrepancy, e.g. a misprint in the source) or DIFF. The output is a table with the counts per demo and one HTML or PDF file per demo with the full comparison, which can serve as a software validation record, e.g. for ISO/IEC 17025:2017, 7.11.2:

G <- doGauge(c("JCGM100", "JCGM101", "JCGM102"), output = "html", outdir = tempdir())
G$summary                     # PASS / NOTE / DIFF per demo
browseURL(G$files[1])         # the full comparison of the first demo

The demos use \(M = 10^6\) Monte Carlo trials and take a few minutes; they are therefore not run in this vignette.

12 Summary

For the heater, the complete workflow shows why a nonlinear model needs more than the first-order GUF:

References

Evaluation of measurement data – Guide to the expression of uncertainty in measurement. JCGM 100:2008.

Evaluation of measurement data – Supplement 1 to the GUM – Propagation of distributions using a Monte Carlo method. JCGM 101:2008.

Evaluation of measurement data – Supplement 2 to the GUM – Extension to any number of output quantities. JCGM 102:2011.

Evaluation of measurement data – The role of measurement uncertainty in conformity assessment. JCGM 106:2012.

ISO 14253-1:2017. Geometrical product specifications (GPS) – Inspection by measurement of workpieces and measuring equipment – Part 1: Decision rules for verifying conformity or nonconformity with specifications.

Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S (2010). Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer Physics Communications, 181, 259-270.