--- title: "From uncertainty budget to conformity decision: a complete workflow with propagate" author: "Andrej-Nikolai Spiess" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{From uncertainty budget to conformity decision: a complete workflow with propagate} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4.5, out.width = "90%") library(propagate) M <- 1E5 # Monte Carlo trials; use 1E6 (or adaptive = TRUE) for final results ``` # 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 . \] * \(V\): reading of a calibrated voltmeter, 230.1 V; calibration certificate: \(U = 0.4\) V with \(k = 2\) (Gaussian). * \(I\): mean of ten repeated readings of a clamp meter (Type A evaluation, scaled and shifted \(t\)-distribution with \(\nu = 9\), JCGM 101:2008, 6.4.9). * \(\delta I\): correction for the clamp meter, zero; manufacturer's accuracy statement \(\pm 0.01\) A (rectangular). * \(\varphi\): the load is nominally resistive, so \(\varphi = 0\); the phase angle is only known to lie within \(\pm 0.12\) rad (rectangular). 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: ```{r typeA} 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))) ``` 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 = ...)`. ```{r budget} 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 ``` 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): ```{r makeDat} set.seed(123) EXPR <- expression(P = V * (I + dI) * cos(phi)) DAT <- makeDat(EXPR, data = BUD, layout = "budget", output = "sim", nsim = M) head(DAT, 3) attr(DAT, "budget")[, c("name", "mean", "u", "df", "dist")] ``` 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()` ```{r propagate} PROP <- propagate(EXPR, DAT, type = "sim", check = FALSE) ``` `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: ```{r results} PROP$prop PROP$sim PROP$ci.short ``` 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\): ```{r shift} 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)) ``` 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: ```{r summary} summary(PROP) ``` The output shows, from top to bottom: * **Taylor expansion vs. Monte Carlo.** The standard uncertainty grows from \(u_1 = 1.62\) W (first order) to \(u_2 = 3.77\) W (second order with Gaussian inputs), and the second order with the actual input moments gives \(u_{2m} = 2.69\) W, equal to the Monte Carlo value 2.70 W. \(u_2\) overshoots because a Gaussian phase angle with the same standard uncertainty has a larger fourth moment than the rectangular one; with the skewness and excess kurtosis of each input (here \(-1.2\) for the rectangular \(\varphi\) and \(\delta I\)), the second-order terms are exact for this quadratic dependence on \(\varphi\). * **Nonlinearity of the model.** The second-order terms make up 64 % of \(u_{2m}^2\), and the second-order estimate is shifted by 1.5 first-order standard uncertainties: the model is strongly nonlinear at the estimates. * **First-order uncertainty budget.** The index \(c_i^2 u^2(x_i)/u^2(y)\) attributes the uncertainty to \(\delta I\), \(V\) and \(I\); the phase angle has index 0 % because \(c_\varphi = 0\). A note points out that it contributes through second-order terms. * **Second-order uncertainty budget.** It decomposes \(u_{2m}^2(y)\) into the first- and second-order terms of each input quantity, with the input moments, and the phase angle now dominates (64 %). * **Validation of the GUF by the MCM** (JCGM 101:2008, 8). The endpoints of the GUF interval differ from those of the Monte Carlo interval by up to 4.8 W, much more than the numerical tolerance \(\delta = 0.05\) W: the GUF is **not validated**, and the Monte Carlo results are to be used. * **Derivatives, covariance checks and the shape of the output**: the input samples have the stated covariance matrix (relative Frobenius norm at the level of the Monte Carlo noise), and the output distribution is skewed to the left. # 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: ```{r sobol} S <- sobol(PROP) round(rbind(first = S$S, total = S$ST), 3) ``` 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\): ```{r sobolcheck} 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) ``` 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: ```{r plot} 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: ```{r fitDistr} FIT <- fitDistr(PROP, nmax = 5000, verbose = FALSE) head(FIT$stat[, c("Distribution", "k", "BIC", "dBIC", "weight")], 5) ``` 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): ```{r conformityMC} CMC <- conformity(PROP, lower = 994, upper = 1006, rule = "all") CMC ``` 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. ```{r conformityPlot, fig.height = 3.5} plot(CMC) ``` From the first-order GUF, the same assessment reaches a different conclusion: ```{r conformityGUF} CGUF <- conformity(PROP, lower = 994, upper = 1006, method = "GUF", rule = c("stringent", "simple")) CGUF ``` 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`. ```{r report} 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 ``` 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`: ```{r 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\): ```{r predictNLS1, results = "hide", message = FALSE} 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) ``` ```{r predictNLS2} round(rbind("u(conc) = 0" = PR1$summary[1, ], "u(conc) = 0.2" = PR2$summary[1, ]), 4) ``` 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: ```{r predictNLS3, results = "hide", message = FALSE} 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 ``` ```{r predictNLS4} 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\)): ```{r MV} 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 ``` 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: ```{r MVplot, fig.height = 6.5} 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: ```{r doGauge, eval = FALSE} 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: * `summary()` flags the nonlinearity (64 % second-order share), shows how \(u_1\), \(u_2\) and \(u_{2m}\) approach the Monte Carlo value, and reports the failed validation of the GUF; * the second-order budget and `sobol()` reveal the input quantity that the first-order budget misses entirely; * the Monte Carlo output is skewed, and `fitDistr()` provides a closed-form description of it; * `conformity()` with the ISO 14253-1 decision rule gives "no decision" from the Monte Carlo results, but "accept" from the (invalid) GUF; * `report()` documents the whole evaluation on two pages. # 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.