Nothing
R/report.R
In propagate: Propagation of Uncertainty
Defines functions
report
Documented in
report
report <- function(object, file = tempfile(fileext = ".pdf"), title = "Uncertainty evaluation", measurand = NULL, unit = "",
compare = TRUE, nsim.compare = 1E5, fitDistr = FALSE, fit.nmax = 1E4,
gauge.mean = NULL, gauge.u = NULL, gauge.method = c("GUF", "MC"), gauge.dec = NULL, gauge.M = 1E6,
interval = c("symmetric", "shortest"), ndig = 2, sobol = TRUE, conformity = NULL, seed = NULL,
fontsize = 1)
{
## one object or a list of objects (one analysis each); per-analysis arguments are recycled
OBJ <- if (inherits(object, "propagate")) list(object) else object
if (!is.list(OBJ) || !length(OBJ) || !all(sapply(OBJ, inherits, "propagate")))
stop("report: 'object' must be a 'propagate' object or a list of them!")
nA <- length(OBJ)
per <- function(x, i) { if (is.null(x)) return(NULL); v <- rep_len(if (is.list(x)) x else as.list(x), nA)[[i]]
if (length(v) == 1 && !is.list(v) && is.na(v)) NULL else v }
## 'conformity': one specification list(lower, upper, ...) for all analyses, or a list of them (NULL/NA = none)
CONF <- if (is.null(conformity)) NULL else if (any(c("lower", "upper") %in% names(conformity))) list(conformity) else conformity
if (length(title) < nA && !is.null(names(OBJ)) && identical(title, "Uncertainty evaluation")) title <- names(OBJ)
interval <- match.arg(interval); gauge.method <- match.arg(gauge.method); fs <- fontsize
RED <- "#8b0000"; GRN <- "#006400"; VIO <- "#6a1fb0"; BLU <- "#1f4e79"; BAR <- "#4682b4"; ORA <- "#ff8c00"
f <- function(x, d) formatC(as.numeric(x), format = "f", digits = d)
g4 <- function(x, d = 4) formatC(signif(as.numeric(x), d), format = "g", digits = d)
num1 <- function(x) if (is.numeric(x) && length(x) && is.finite(x[1])) x[1] else NA_real_
PV <- tryCatch(as.character(utils::packageVersion("propagate")), error = function(e) "")
## ---------------------------------------------------------------- page layout helpers (shared)
if (!is.null(file)) { grDevices::cairo_pdf(file, width = 8.27, height = 11.69, onefile = TRUE); on.exit(grDevices::dev.off()) }
op <- par(mar = c(0, 0, 0, 0)); on.exit(par(op), add = TRUE, after = FALSE)
PAGE <- 0; Y <- 0; TITLE <- title[1]
newPage <- function() { plot.new(); plot.window(c(0, 1), c(0, 1)); PAGE <<- PAGE + 1
rect(0.04, 0.955, 0.96, 0.985, col = BLU, border = NA)
text(0.055, 0.97, TITLE, adj = c(0, 0.5), col = "white", font = 2, cex = 1.15 * fs)
text(0.945, 0.97, format(Sys.Date()), adj = c(1, 0.5), col = "white", cex = 0.75 * fs)
text(0.5, 0.02, paste0("propagate", if (nzchar(PV)) paste0(" ", PV), " - page ", PAGE), cex = 0.65 * fs, col = "grey40")
Y <<- 0.93 }
need <- function(h) if (Y - h < 0.035) newPage()
head <- function(t, x = 0.05) { need(0.05); text(x, Y, t, adj = c(0, 0.5), font = 2, col = BLU, cex = 0.95 * fs); Y <<- Y - 0.026 }
para <- function(t, cex = 0.62, col = "black", font = 1, family = "", x = 0.05, w = 0.9) {
nc <- floor(w/strwidth("m", cex = cex * fs, family = family) * if (nzchar(family)) 1 else 1.6)
for (l in unlist(lapply(t, strwrap, width = nc))) { need(0.015)
text(x, Y, l, adj = c(0, 0.5), cex = cex * fs, col = col, font = font, family = family); Y <<- Y - 0.0125 * fs }
Y <<- Y - 0.008 }
## table; 'bars' (and 'bars2'): shares drawn as horizontal bars in the last (empty) column, upper and lower half
tab <- function(D, w, x0 = 0.05, cex = 0.6, colr = NULL, fill = NULL, bars = NULL, bars2 = NULL, lh = 0.0145, align = NULL) {
lh <- lh * fs; xs <- x0 + c(0, cumsum(w)); nc <- ncol(D)
if (is.null(align)) align <- c("l", rep("r", nc - 1)) ## default: first column left-, others right-aligned
L <- align == "l"; xt <- ifelse(L, xs[-(nc + 1)] + 0.005, xs[-1] - 0.005)
row <- function(v, ...) { v <- as.character(v); if (any(L)) text(xt[L], Y, v[L], adj = c(0, 0.5), ...)
if (any(!L)) text(xt[!L], Y, v[!L], adj = c(1, 0.5), ...) }
hdr <- function() { rect(x0, Y - lh/2, max(xs), Y + lh/2, col = "#e2e2e2", border = NA)
row(names(D), font = 2, cex = cex * fs); Y <<- Y - lh }
bmax <- max(1, abs(as.numeric(c(bars, bars2))), na.rm = TRUE)
drawBar <- function(b, yc, hh, col) if (!is.null(b) && is.finite(b) && b != 0) {
wb <- (xs[nc + 1] - xs[nc] - 0.01) * abs(b)/bmax
rect(xs[nc] + 0.005, yc - hh, xs[nc] + 0.005 + wb, yc + hh, border = NA, col = if (b < 0) RED else col) }
need(2 * lh); hdr()
for (i in seq_len(nrow(D))) {
if (Y - lh < 0.035) { newPage(); hdr() }
if (!is.null(fill) && !is.na(fill[i])) rect(xs[nc], Y - lh/2, xs[nc + 1], Y + lh/2, col = fill[i], border = NA)
if (is.null(bars2)) drawBar(bars[i], Y, lh * 0.32, BAR) else {
drawBar(bars[i], Y + lh * 0.17, lh * 0.15, BAR); drawBar(bars2[i], Y - lh * 0.17, lh * 0.15, VIO) }
cl <- if (!is.null(colr) && !is.na(colr[i])) colr[i] else "black"
row(unlist(D[i, ]), cex = cex * fs, col = cl, font = if (cl == "black") 1 else 2)
segments(x0, Y - lh/2, max(xs), Y - lh/2, col = "#dddddd", lwd = 0.5); Y <<- Y - lh }
Y <<- Y - 0.008 }
## symbolic derivative as one-line string; shortened with "..." to width w
symStr <- function(e) if (is.language(e) || is.numeric(e)) paste(deparse(e, width.cutoff = 500L), collapse = "") else as.character(e)
symVec <- function(G) if (is.list(G)) vapply(G, symStr, "") else as.character(G)
fitTxt <- function(x, w, cex = 0.6) vapply(x, function(z) { if (strwidth(z, cex = cex * fs) <= w) return(z)
while (nchar(z) > 2 && strwidth(paste0(z, "..."), cex = cex * fs) > w) z <- substr(z, 1, nchar(z) - 1); paste0(z, "...") }, "", USE.NAMES = FALSE)
## a matrix with caption as table; cells are character
mtab <- function(M, cap, x0, cw, cex = 0.6, align = NULL) {
text(x0, Y, cap, adj = c(0, 0.5), cex = cex * fs, font = 3, col = "grey20"); Y <<- Y - 0.0135 * fs
D <- data.frame(rownames(M), matrix(M, nrow(M)), check.names = FALSE, stringsAsFactors = FALSE); names(D) <- c(" ", colnames(M))
tab(D, c(0.06, rep(cw, ncol(M))), x0 = x0, cex = cex, align = align) }
## several matrices side by side (all start at the same height)
mrow <- function(L, caps, x0s, cw, h, align = NULL) { need(h); y0 <- Y; ym <- Y
for (j in seq_along(L)) { Y <<- y0; mtab(L[[j]], caps[j], x0s[j], cw, align = align); ym <- min(ym, Y) }; Y <<- ym }
fmtM <- function(M) { M[M == 0] <- 0; F <- matrix(formatC(signif(M, 3), format = "g", digits = 3), nrow(M), dimnames = dimnames(M)); F }
## a sub-figure in user coordinates [x0, x1] x [y0, y1] of the page; back to the page with page()
ndc <- function(x, y) c(grconvertX(x, "user", "ndc"), grconvertY(y, "user", "ndc")) ## call while in page coordinates
subfig <- function(fig, mar) par(fig = fig, new = TRUE, mar = mar)
page <- function() { par(fig = c(0, 1, 0, 1), new = TRUE, mar = c(0, 0, 0, 0)); plot.new(); plot.window(c(0, 1), c(0, 1)) }
OUT <- vector("list", nA); names(OUT) <- names(OBJ)
for (ia in seq_len(nA)) {
O <- OBJ[[ia]]; TITLE <- per(title, ia); if (is.null(TITLE)) TITLE <- "Uncertainty evaluation"
MSi <- per(measurand, ia); UNi <- per(unit, ia); if (is.null(UNi)) UNi <- ""
GM <- per(gauge.mean, ia); GU <- per(gauge.u, ia); GD <- per(gauge.dec, ia)
SEED <- per(seed, ia); CF <- if (is.null(CONF)) NULL else per(CONF, ia)
u1 <- unname(O$prop["u.1"]); dec <- if (is.finite(u1) && u1 > 0) max(0, 3 - floor(log10(u1))) else 4
ciMC <- unname(if (interval == "shortest") O$ci.short else O$sim[5:6])
if (is.null(MSi)) { n <- names(O$expr); MSi <- if (length(n) && nzchar(n[1])) n[1] else "y" }
sym <- function(s) { e <- tryCatch(str2lang(gsub("_(\\w+)", "[\\1]", s)), error = function(e) NULL); if (is.null(e)) s else e }
MS <- sym(MSi)
UN <- { map <- c(Ohm = "Omega", degC = "degree*C", um = "mu*m", ug = "mu*g", uV = "mu*V", uA = "mu*A", "%" = "'%'")
if (!nzchar(UNi)) quote("") else if (UNi %in% names(map)) str2lang(map[[UNi]]) else
tryCatch({ e <- str2lang(UNi); if (is.name(e) || is.call(e)) e else UNi }, error = function(e) UNi) }
y <- O$resSIM[is.finite(O$resSIM)]
## ---------------------------------------------------------------- results of all methods
mcRow <- function(R, lab) { ci <- unname(if (interval == "shortest") R$ci.short else R$sim[5:6])
c(paste0("MC, ", lab, ", M = ", format(length(R$resSIM), big.mark = " ", scientific = FALSE)), f(R$sim["mean.MC"], dec),
if (any(R$df.MC <= 2)) "n/a" else f(R$sim["u.MC"], dec + 1), paste0("[", f(ci[1], dec), ", ", f(ci[2], dec), "]")) }
ciG <- .ciGUF(O$prop) ## v1.3-0: by name
RES <- rbind(c("GUF, 1st order", f(O$prop["mean.1"], dec), f(u1, dec + 1),
paste0("[", f(ciG[1], dec), ", ", f(ciG[2], dec), "] (k = ", f(num1(O$k), 3), ")")),
c("GUF, 2nd order", f(O$prop["mean.2"], dec), f(O$prop["u.2"], dec + 1), "-"),
if (is.finite(num1(O$prop["u.2m"]))) ## v1.3-0: second order with input moments
c("GUF, 2nd order, input moments", f(O$prop["mean.2"], dec), f(O$prop["u.2m"], dec + 1), "-"))
MCR <- list(); MCR[[if (is.na(O$mc.margins)) "samples" else O$mc.margins]] <- O
if (compare && O$type %in% c("stat", "raw") && !is.null(O$input$data))
for (m in setdiff(c("S1", "mvt"), O$mc.margins)) {
R <- NULL; invisible(capture.output(R <- tryCatch(suppressMessages(propagate(O$expr, O$input$data, type = O$type, cov = O$input$cov,
mc.margins = m, nsim = nsim.compare, alpha = O$alpha, ws.group = O$input$ws.group, check = FALSE)), error = function(e) NULL)))
if (!is.null(R)) MCR[[m]] <- R }
for (m in intersect(c("S1", "mvt", "samples"), names(MCR))) RES <- rbind(RES, mcRow(MCR[[m]], m))
RES <- as.data.frame(RES, stringsAsFactors = FALSE)
names(RES) <- c("Method", "Estimate", "u", paste0(signif(100 * (1 - O$alpha), 3), " % interval", if (interval == "shortest") " (shortest)" else ""))
resCol <- ifelse(grepl("^GUF, 1st|^MC, S1", RES$Method), RED, NA)
## ---------------------------------------------------------------- validation of the GUF by the MCM (JCGM 101:2008, 8)
VAL <- NULL; valCol <- NA
uG <- if (is.finite(num1(O$u.exp)) && is.finite(num1(O$k)) && O$k > 0) O$u.exp/O$k else NA
if (is.finite(uG) && uG > 0 && all(is.finite(ciG)) && all(is.finite(ciMC))) {
lv <- floor(log10(uG)) - (ndig - 1); delta <- 0.5 * 10^lv
dlo <- abs(ciG[1] - ciMC[1]); dhi <- abs(ciG[2] - ciMC[2])
SE <- O$sim.se
sE <- if (is.null(SE)) NA else if (interval == "shortest") max(SE[c("short.low", "short.high")]) else max(SE[c("sym.low", "sym.high")])
verdict <- if (!is.finite(sE)) "not assessable (no block results)" else if (2 * sE > delta) "INCONCLUSIVE: increase nsim" else
if (max(dlo, dhi) <= delta) "GUF VALIDATED" else "GUF NOT VALIDATED: use MC"
valCol <- if (grepl("^GUF VALIDATED", verdict)) GRN else if (grepl("NOT", verdict)) RED else ORA
cover <- mean(y >= ciG[1] & y <= ciG[2])
VAL <- data.frame(Item = c(paste0("numerical tolerance delta (", ndig, " significant digits of u)"),
"d_low = |low(GUF) - low(MC)|", "d_high = |high(GUF) - high(MC)|", "2 s of the MC endpoints",
"result", "MC coverage of the GUF interval"),
Value = c(g4(delta, 2), g4(dlo, 3), g4(dhi, 3), g4(2 * sE, 2), verdict,
paste0(f(100 * cover, 1), " % (nominal ", signif(100 * (1 - O$alpha), 3), " %)")),
stringsAsFactors = FALSE)
}
## ---------------------------------------------------------------- numerical accuracy of the MC results (JCGM 101:2008, 7.9)
ACC <- NULL
if (!is.null(O$sim.se)) {
SE <- O$sim.se; nmL <- if (interval == "shortest") c("short.low", "short.high") else c("sym.low", "sym.high")
v <- c(unname(O$sim["mean.MC"]), unname(O$sim["u.MC"]), ciMC); s <- unname(c(SE["mean.MC"], SE["u.MC"], SE[nmL]))
ok <- is.finite(v) & is.finite(s) & s > 0
dS <- ifelse(ok, floor(-log10(4 * s)), NA) ## 2 s <= 1/2 10^-d: d decimals are stable
sig <- ifelse(ok & v != 0, pmax(0, floor(log10(abs(v))) + dS + 1), NA)
rnd <- function(x, d) if (!is.finite(d)) "n/a" else if (d >= 0) formatC(round(x, d), format = "f", digits = d) else format(round(x, d))
ACC <- data.frame(Quantity = c("estimate", "u", "interval, low", "interval, high"),
"stable value" = mapply(rnd, v, dS), "2 s" = g4(2 * s, 2), digits = ifelse(is.na(sig), "-", sig),
check.names = FALSE, stringsAsFactors = FALSE)
if (any(O$df.MC <= 2)) ACC[2, 2:4] <- c("n/a", "-", "-")
}
## ---------------------------------------------------------------- input quantities: distribution and Type A/B
gr <- as.numeric(O$evalGrad); gr[gr == 0] <- 0; S <- as.matrix(O$covMat); nm <- colnames(S); N <- length(nm)
nu <- if (nrow(O$data) >= 3) as.numeric(O$data[3, nm]) else rep(Inf, N)
BG <- if (!is.null(O$input$budget)) O$input$budget else attr(O$input$data, "budget")
dname <- c(normal = "Gaussian", t = "t", rect = "rectangular", tri = "triangular", arcsine = "arcsine", ctrap = "curvil. trapezoid")
if (!is.null(BG) && all(nm %in% BG$name)) {
b <- BG[match(nm, BG$name), ]; dd <- as.character(b$dist); bdf <- as.numeric(b$df)
DIST <- ifelse(dd %in% names(dname), dname[dd], dd); DIST[is.na(DIST)] <- "sample"
TYP <- ifelse(dd == "t" & is.finite(bdf), "A", ifelse(is.na(dd), "-", "B"))
} else if (identical(O$type, "sim")) { DIST <- rep("sample", N); TYP <- rep("-", N)
} else {
mvt <- identical(O$mc.margins, "mvt"); A <- O$type == "raw" | (is.finite(nu) & nu < 1000)
DIST <- ifelse(A, if (mvt) "mvt" else "t", "Gaussian"); TYP <- ifelse(A, "A", "B")
}
## ---------------------------------------------------------------- uncertainty budget with first-order and Sobol' indices
Kc <- outer(gr, gr) * S; uy2 <- sum(Kc)
share <- diag(Kc)/uy2; cv <- uy2 - sum(diag(Kc))
indep <- abs(cv) <= 1E-12 * abs(uy2) && !identical(O$mc.margins, "mvt") &&
(!is.matrix(O$corMat) || all(abs(O$corMat[upper.tri(O$corMat)]) < 1E-12))
SOB <- if (isTRUE(sobol) && indep) tryCatch(suppressMessages(get("sobol", mode = "function")(O)), error = function(e) NULL) else NULL
Si <- if (!is.null(SOB)) pmax(0, as.numeric(SOB$S[nm])) else rep(NA, N)
BUD <- data.frame(Quantity = nm, Estimate = formatC(O$data[1, nm], format = "fg", digits = 8), "u(x_i)" = g4(sqrt(diag(S))),
nu_i = ifelse(!is.finite(nu) | nu >= 1000, "Inf", formatC(nu, format = "g")), "c_i" = g4(gr),
"|c_i| u(x_i)" = g4(abs(gr) * sqrt(diag(S))), Index = paste0(f(100 * share, 1), " %"),
"S_i (MC)" = if (is.null(SOB)) rep("-", N) else paste0(f(100 * Si, 1), " %"), " " = "",
check.names = FALSE, stringsAsFactors = FALSE)
BARS <- share; BARS2 <- if (is.null(SOB)) NULL else Si
if (abs(cv) > 1E-12 * abs(uy2)) { BUD <- rbind(BUD, c("covariance terms", rep("", 5), paste0(f(100 * cv/uy2, 1), " %"), "-", ""))
BARS <- c(BARS, cv/uy2); if (!is.null(BARS2)) BARS2 <- c(BARS2, NA) }
if (!is.null(SOB) && 1 - sum(Si) > 0.005) { BUD <- rbind(BUD, c("interactions (MC)", rep("", 5), "-", paste0(f(100 * (1 - sum(Si)), 1), " %"), ""))
BARS <- c(BARS, NA); BARS2 <- c(BARS2, 1 - sum(Si)) }
BUD <- rbind(BUD, c("u(y), first order", rep("", 4), g4(sqrt(uy2)), "100.0 %", if (is.null(SOB)) "-" else "100.0 %", ""))
BARS <- c(BARS, NA); if (!is.null(BARS2)) BARS2 <- c(BARS2, NA)
H <- O$evalHess
zeroC <- if (is.matrix(H) && all(dim(H) == length(gr))) nm[gr == 0 & rowSums(abs(H)) > 0] else character(0)
## ---------------------------------------------------------------- nonlinearity of the model and shape of the output
P2 <- O$prop; NL <- NULL
yy <- y - mean(y); m2 <- mean(yy^2); SH <- c(skewness = mean(yy^3)/m2^1.5, kurtosis = mean(yy^4)/m2^2 - 3)
heavy <- (!is.null(O$df.MC) && any(O$df.MC <= 4))
lev <- function(x, a, b) if (!is.finite(x)) "-" else if (abs(x) < a) "negligible" else if (abs(x) < b) "moderate" else "large"
uS <- if (is.finite(num1(P2["u.2m"]))) "u.2m" else "u.2" ## v1.3-0: with input moments if available
if (all(is.finite(P2[c("mean.1", "mean.2", "u.1", uS)])) && P2[uS] > 0) {
nlShare <- unname((P2[uS]^2 - P2["u.1"]^2)/P2[uS]^2)
nlShift <- if (P2["u.1"] > 0) unname((P2["mean.2"] - P2["mean.1"])/P2["u.1"]) else NA
qy <- quantile(y, c(0.025, 0.5, 0.975)) ## quantile skewness: robust, exists for any nu
sk <- unname((qy[3] + qy[1] - 2 * qy[2])/(qy[3] - qy[1]))
NL <- data.frame(Measure = c(paste0("share of the second-order terms in ", uS, "^2"), "shift of the estimate (mean.2 - mean.1)/u.1",
"quantile skewness of the MC output, (q97.5 + q2.5 - 2 q50)/(q97.5 - q2.5)",
"skewness of the MC output (moments; Gaussian: 0)", "excess kurtosis of the MC output (moments; Gaussian: 0)"),
Value = c(paste0(formatC(100 * nlShare, format = "fg", digits = 3), " %"), formatC(nlShift, format = "fg", digits = 3),
formatC(c(sk, SH), format = "fg", digits = 3)),
Assessment = c(lev(nlShare, 0.01, 0.1), lev(nlShift, 0.05, 0.2), lev(sk, 0.05, 0.2),
if (heavy) c("-", "-") else c(lev(SH[1], 0.1, 0.5), lev(SH[2], 0.2, 1))), stringsAsFactors = FALSE)
}
SHP <- data.frame(Measure = c("skewness", "excess kurtosis"), Value = formatC(SH, format = "fg", digits = 4), stringsAsFactors = FALSE)
## ---------------------------------------------------------------- comparison with gauge values
CHK <- NULL
if (!is.null(GM) || !is.null(GU)) {
decOf <- function(x) max(sapply(x, function(v) { s <- format(v, scientific = FALSE, digits = 15, drop0trailing = TRUE)
if (grepl(".", s, fixed = TRUE)) nchar(sub("^[^.]*\\.", "", s)) else 0 }))
dd <- rep_len(if (is.null(GD)) c(if (is.null(GM)) 0 else decOf(GM), if (is.null(GU)) 0 else decOf(GU)) else GD, 2)
isMC <- gauge.method == "MC"; se <- O$sim.se; fac <- sqrt(1 + O$nsim/gauge.M)
obs <- if (isMC) unname(O$sim[c("mean.MC", "u.MC")]) else c(unname(O$prop["mean.1"]), u1)
tol <- 0.5 * 10^-dd + if (isMC && !is.null(se)) 4 * fac * unname(se[c("mean.MC", "u.MC")]) else 0
pub <- c(if (is.null(GM)) NA else GM, if (is.null(GU)) NA else GU)
CHK <- data.frame(Item = paste0(if (isMC) "Monte Carlo" else "GUF, 1st order", c(": estimate", ": u")),
Gauge = c(f(pub[1], dd[1]), f(pub[2], dd[2])), propagate = c(f(obs[1], dd[1] + 1), f(obs[2], dd[2] + 1)),
Tolerance = g4(tol, 2), Result = ifelse(abs(obs - pub) <= tol, "PASS", "FAIL"), stringsAsFactors = FALSE)[!is.na(pub), ]
}
## ---------------------------------------------------------------- distribution fit, conformity assessment
FIT <- if (isTRUE(fitDistr)) tryCatch(suppressMessages(get("fitDistr", mode = "function")(y, nmax = fit.nmax, verbose = FALSE, plot = "none")),
error = function(e) NULL) else NULL
CA <- NULL
if (!is.null(CF)) {
CF <- modifyList(list(lower = -Inf, upper = Inf, method = "MC", rule = c("stringent", "simple"), p = 1 - O$alpha, interval = interval), CF)
CA <- tryCatch(get("conformity", mode = "function")(O, lower = CF$lower, upper = CF$upper, method = CF$method, rule = CF$rule,
p = CF$p, interval = CF$interval), error = function(e) { message("report: conformity: ", conditionMessage(e)); NULL })
}
## ---------------------------------------------------------------- reproducibility record (ISO/IEC 17025:2017, 7.5, 7.8)
CALL <- if (!is.null(O$call)) paste(deparse(O$call, width.cutoff = 500L), collapse = " ") else
paste0("propagate(expr = ", paste(deparse(O$expr, width.cutoff = 500L), collapse = ""), ", data = <", O$type, " data>, type = \"", O$type,
"\", cov = ", if (is.matrix(O$input$cov)) "<matrix>" else paste0("\"", O$input$cov, "\""), ", nsim = ", format(O$nsim, scientific = TRUE),
", alpha = ", O$alpha, if (!is.na(O$mc.margins)) paste0(", mc.margins = \"", O$mc.margins, "\""),
if (!is.null(O$input$exp.uc)) paste0(", exp.uc = \"", O$input$exp.uc, "\""), ") [reconstructed]")
if (nchar(CALL) > 400) CALL <- paste0(substr(CALL, 1, 400), " ...")
RK <- RNGkind()
REPRO <- c(Call = CALL,
"Monte Carlo" = paste0("M = ", format(length(O$resSIM), big.mark = " ", scientific = FALSE), " trials",
if (!is.null(O$adaptive)) paste0(" (adaptive, ", if (O$adaptive$converged) "converged" else "not converged", ")"),
"; alpha = ", O$alpha, "; input model: ", if (is.na(O$mc.margins)) "user-supplied samples" else O$mc.margins,
"; cov = ", if (is.matrix(O$input$cov)) "user matrix" else O$input$cov,
"; coverage interval: ", interval),
"RNG, software" = paste0("seed ", if (is.null(SEED)) "not recorded (use report(..., seed = ))" else paste(SEED, collapse = ", "),
"; RNG ", RK[1], "/", RK[2], "/", RK[3], "; propagate ", if (nzchar(PV)) PV else "(version unknown)",
"; R ", R.version$major, ".", R.version$minor, " (", R.version$platform, "); ", format(Sys.time(), "%Y-%m-%d %H:%M %Z")))
## ================================================================ page 1
newPage()
head("Measurement model")
text(0.08, Y - 0.006, bquote(italic(.(MS)) == .(O$expr[[1]])), adj = c(0, 0.5), cex = 0.95 * fs); Y <- Y - 0.042
head("Input data")
D <- O$input$data; y0 <- Y; corDone <- FALSE
IN <- data.frame(Quantity = nm, Estimate = BUD$Estimate[seq_len(N)], "u(x_i)" = BUD[["u(x_i)"]][seq_len(N)], nu_i = BUD$nu_i[seq_len(N)],
Distribution = DIST, Type = TYP, check.names = FALSE, stringsAsFactors = FALSE)
if (identical(O$type, "raw") && !is.null(D) && nrow(D) <= 15) {
RAW <- data.frame(k = seq_len(nrow(D)), apply(D[, nm, drop = FALSE], 2, function(v) formatC(v, format = "g", digits = 6)), check.names = FALSE)
tab(RAW, c(0.03, rep(min(0.09, 0.3/N), N)), lh = 0.0135); y1 <- Y; Y <- y0
names(IN)[5] <- "Distr."
tab(IN, c(0.07, 0.1, 0.09, 0.04, 0.07, 0.04), x0 = 0.5, align = c("l", "r", "r", "r", "l", "l"), lh = 0.0135)
if (N > 1 && N <= 6 && any(abs(cov2cor(S)[upper.tri(S)]) > 1E-12)) {
CR <- data.frame(r = nm, matrix(f(cov2cor(S), 3), N), check.names = FALSE); names(CR)[-1] <- nm
tab(CR, c(0.07, rep(min(0.07, 0.34/N), N)), x0 = 0.5, lh = 0.0135) }
Y <- min(Y, y1); corDone <- TRUE
} else {
tab(IN, c(0.14, 0.16, 0.13, 0.07, 0.14, 0.06), align = c("l", "r", "r", "r", "l", "l"), lh = if (N > 12) 0.0125 else 0.0135)
}
if (!corDone && N > 1 && N <= 6 && any(abs(cov2cor(S)[upper.tri(S)]) > 1E-12)) {
CR <- data.frame(r = nm, matrix(f(cov2cor(S), 3), N), check.names = FALSE); names(CR)[-1] <- nm
tab(CR, c(0.08, rep(0.08, N)), lh = 0.0135)
}
head("Result statement")
dU <- function(u) max(0, 1 - floor(log10(abs(u))))
fU <- function(u) if (is.finite(u) && u > 0) formatC(signif(u, 2), format = "f", digits = dU(u)) else "n/a"
fY <- function(y, u) if (is.finite(u) && u > 0) formatC(round(y, dU(u)), format = "f", digits = dU(u)) else format(signif(y, 6))
uM <- if (any(O$df.MC <= 2)) NA else unname(O$sim["u.MC"]); uR <- if (is.finite(uM)) uM else diff(ciMC)/4
PP <- paste0(signif(100 * (1 - O$alpha), 3), " %")
nuS <- if (is.na(num1(O$ws.df))) "-" else if (O$ws.df >= 1000) "infinity" else as.character(O$ws.df)
need(0.12); top <- Y + 0.012; bh <- 0.1 * fs
fit1 <- function(x, y, e, cex) text(x, y, e, cex = min(cex, 0.4/strwidth(e, cex = 1)))
text(0.17, top - 0.014, "GUM uncertainty framework", adj = c(0, 0.5), cex = 0.8 * fs, col = "grey20")
fit1(0.28, top - 0.045, bquote(italic(.(MS)) ~ "=" ~ .(fY(O$prop["mean.1"], u1)) ~ .(UN) * "," ~~ italic(u) ~ "=" ~ .(fU(u1)) ~ .(UN)), 1.05 * fs)
fit1(0.28, top - 0.078, bquote(.(PP) ~ "interval:" ~ group("[", .(fY(ciG[1], u1)) * "," ~ .(fY(ciG[2], u1)), "]") ~ .(UN) ~~
group("(", italic(k) ~ "=" ~ .(f(num1(O$k), 2)) * "," ~ nu[eff] ~ "=" ~ .(nuS), ")")), 0.8 * fs)
text(0.64, top - 0.014, "Monte Carlo method", adj = c(0, 0.5), cex = 0.8 * fs, col = "grey20")
fit1(0.72, top - 0.045, bquote(italic(.(MS)) ~ "=" ~ .(fY(O$sim["mean.MC"], uR)) ~ .(UN) * "," ~~ italic(u) ~ "=" ~ .(fU(uM)) ~ .(if (is.finite(uM)) UN else quote(""))), 1.05 * fs)
fit1(0.72, top - 0.078, bquote(.(PP) ~ .(interval) ~ "interval:" ~ group("[", .(fY(ciMC[1], uR)) * "," ~ .(fY(ciMC[2], uR)), "]") ~ .(UN)), 0.8 * fs)
Y <- top - bh - 0.016
head("Results of all methods")
tab(RES, c(0.32, 0.15, 0.13, 0.3), colr = resCol, lh = 0.0135)
## validation (left) and numerical accuracy (right), side by side
if (!is.null(VAL) || !is.null(ACC)) {
need(0.13); y0 <- Y; yL <- yR <- Y
if (!is.null(VAL)) { head("Validation of the GUF (JCGM 101:2008, 8)")
tab(VAL, c(0.29, 0.15), colr = c(NA, NA, NA, NA, valCol, NA), cex = 0.56, lh = 0.0135); yL <- Y }
Y <- y0
if (!is.null(ACC)) { head("Numerical accuracy of MC (7.9)", x = 0.52)
tab(ACC, c(0.12, 0.12, 0.08, 0.06), x0 = 0.52, cex = 0.56, lh = 0.0135); yR <- Y }
Y <- min(yL, yR)
para(paste0("GUF validated if d_low, d_high <= delta and the MC endpoints are stable (2 s <= delta). s: numerical standard uncertainty ",
"from 10 blocks of the MC trials; a digit is stable if 2 s <= half a unit of it; digits: stable significant digits."), cex = 0.54, col = "#3a3a3a")
}
if (!is.null(NL)) {
head("Nonlinearity of the model and shape of the output")
tab(NL, c(0.56, 0.12, 0.14), colr = c(NA, GRN, ORA, RED)[match(NL$Assessment, c("-", "negligible", "moderate", "large"))], lh = 0.0135)
para(c(paste0("All zero for a model linear at the estimates with Gaussian inputs (limits 1 %/10 %, 0.05/0.2, 0.05/0.2, 0.1/0.5, 0.2/1); ",
"large values: the first-order GUF may be inadequate, prefer the MC results (JCGM 101:2008, 5.8)."),
if (heavy) "Inputs with DOF <= 4: the moments of the output may not exist; skewness and kurtosis are not assessed.",
if (length(zeroC)) paste0("Note: ", paste(zeroC, collapse = ", "), " with c_i = 0 contribute(s) only through second-order terms ",
"(first-order index 0 %, see the Sobol' index on page 2).")), cex = 0.54, col = "#3a3a3a")
}
if (!is.null(CHK)) {
head(paste0("Comparison with the gauge values: ", sum(CHK$Result == "PASS"), " PASS, ", sum(CHK$Result == "FAIL"), " FAIL"))
tab(CHK, c(0.34, 0.14, 0.14, 0.12, 0.1), fill = ifelse(CHK$Result == "PASS", "#d9efd0", "#f4c7c3"))
para(paste0("Tolerance: rounding of the gauge values", if (gauge.method == "MC") paste0(" plus four Monte Carlo standard errors of both ",
"calculations (gauge with ", format(gauge.M, big.mark = " ", scientific = FALSE), " trials)") else "", "."), cex = 0.56)
}
## ================================================================ (page 2) -- blocks flow on and move to the next page if they do not fit
Y <- Y - 0.01; need(0.03 + (N + 3) * 0.0155 + 0.06)
head("Uncertainty budget: first-order law of propagation and Monte Carlo (Sobol') indices")
tab(BUD, c(0.13, 0.12, 0.08, 0.04, 0.09, 0.1, 0.07, 0.08, 0.19), bars = BARS, bars2 = BARS2,
colr = ifelse(BUD$Quantity == "covariance terms" & cv < 0, RED, NA), lh = if (N > 12) 0.0135 else 0.0155)
para(c(paste0("c_i = df/dx_i; Index = c_i^2 u^2(x_i)/u^2(y) (blue); covariance terms = sum_(i != j) c_i c_j u(x_i, x_j) (negative: dark red). ",
"S_i (MC) = first-order Sobol' index: share of Var(Y) due to x_i alone, with all nonlinear effects (violet); interactions = 1 - sum S_i."),
if (is.null(SOB)) paste0("Sobol' indices: not computed (", if (!indep) "they require independent input quantities)." else
if (!isTRUE(sobol)) "sobol = FALSE)." else "computation failed)."),
if (length(zeroC)) paste0("Inputs with c_i = 0 but non-zero second derivatives (", paste(zeroC, collapse = ", "),
") are not represented in the first-order index.")), cex = 0.54)
## histogram (left) with Q-Q plot and convergence trace (right)
need(0.225); yT <- Y - 0.005; yB <- max(0.07, Y - 0.225)
yM <- (yT + yB)/2
F1 <- ndc(c(0.05, 0.6), c(yB, yT)); F2 <- ndc(c(0.64, 0.95), c(yM + 0.005, yT)); F3 <- ndc(c(0.64, 0.95), c(yB, yM - 0.005))
subfig(F1, c(2.4, 1, 1.6, 1))
q <- quantile(y, c(5E-4, 1 - 5E-4)); h <- hist(y[y >= q[1] & y <= q[2]], breaks = 70, plot = FALSE)
band <- if (!is.null(GM) && !is.null(GU)) GM + c(-1, 1) * num1(O$k) * GU else NULL
plot(h, freq = FALSE, col = NA, border = NA, main = paste0("Monte Carlo output (", interval, " interval)"), xlab = "",
ylab = "", yaxt = "n", xlim = range(c(h$breaks, ciMC, ciG, band)), cex.main = 0.8 * fs, font.main = 3, cex.axis = 0.62 * fs)
if (!is.null(band)) rect(band[1], 0, band[2], max(h$density) * 1.05, col = "#dcebcf", border = NA)
plot(h, freq = FALSE, add = TRUE, col = "#c9b8de", border = "#b3a0cc")
abline(v = ciMC, col = VIO, lwd = 2); abline(v = ciG, col = GRN, lwd = 2, lty = 2)
if (!is.null(FIT)) { xs <- seq(min(h$breaks), max(h$breaks), length.out = 500)
lines(xs, do.call(FIT$dens[[1]], c(list(xs), as.list(FIT$bestpar))), col = RED, lwd = 2.5) }
legend("topright", bty = "n", cex = 0.58 * fs, legend = c("MC interval", "GUF interval", if (!is.null(band)) "gauge y +/- k u",
if (!is.null(FIT)) paste("fit:", FIT$stat$Distribution[1])), col = c(VIO, GRN, if (!is.null(band)) "#dcebcf", if (!is.null(FIT)) RED),
lty = c(1, 2, if (!is.null(band)) 1, if (!is.null(FIT)) 1), lwd = c(2, 2, if (!is.null(band)) 6, if (!is.null(FIT)) 2.5))
## Q-Q plot: quantiles of the MC output against those of the GUF distribution (t with nu_eff, or Gaussian)
subfig(F2, c(2, 2.4, 1.3, 0.4))
pp <- c(seq(0.001, 0.01, length.out = 10), seq(0.02, 0.98, by = 0.02), seq(0.99, 0.999, length.out = 10))
nuQ <- if (is.finite(num1(O$ws.df)) && O$ws.df < 1000) O$ws.df else Inf
qG <- unname(O$prop["mean.1"]) + u1 * (if (is.finite(nuQ)) qt(pp, nuQ) else qnorm(pp)); qM <- quantile(y, pp, names = FALSE)
plot(qG, qM, pch = 16, cex = 0.35, col = VIO, xlab = "", ylab = "", main = paste0("Q-Q: MC vs. GUF (", if (is.finite(nuQ)) paste0("t, nu = ", nuQ) else "Gaussian", ")"),
cex.main = 0.62 * fs, font.main = 3, axes = FALSE)
for (sd_ in 1:2) { at <- pretty(if (sd_ == 1) qG else qM, 4); axis(sd_, at = at, labels = formatC(at, format = "g", digits = 4), cex.axis = 0.5 * fs,
mgp = c(1.2, 0.3, 0), tcl = -0.2) }; box()
abline(0, 1, col = GRN, lty = 2)
a2 <- O$alpha/2; points(ciG, quantile(y, c(a2, 1 - a2)), pch = 4, col = RED, cex = 0.9, lwd = 1.5)
mtext("GUF quantile", 1, 1, cex = 0.45 * fs); mtext("MC quantile", 2, 1.4, cex = 0.45 * fs)
## convergence trace: running mean and u relative to the final values, in units of u
subfig(F3, c(2, 2.4, 1.3, 0.4))
yo <- O$resSIM[is.finite(O$resSIM)]; nT <- length(yo); ii <- unique(round(exp(seq(log(min(100, nT)), log(nT), length.out = 60))))
cs <- cumsum(yo); cs2 <- cumsum(yo^2); rm_ <- cs[ii]/ii; rs <- sqrt(pmax(0, (cs2[ii] - ii * rm_^2)/(ii - 1)))
uF <- sd(yo); mF <- mean(yo)
dm <- (rm_ - mF)/uF; du <- rs/uF - 1
plot(ii, dm, type = "l", log = "x", col = BLU, lwd = 1.5, ylim = range(c(dm, du, -0.05, 0.05), finite = TRUE), xlab = "", ylab = "",
main = "Convergence of mean and u", cex.main = 0.62 * fs, font.main = 3, cex.axis = 0.5 * fs, mgp = c(1.2, 0.3, 0), tcl = -0.2)
lines(ii, du, col = ORA, lwd = 1.5); abline(h = 0, col = "grey50", lty = 3)
mtext("number of trials M", 1, 1, cex = 0.45 * fs); mtext("relative to final, in u", 2, 1.4, cex = 0.45 * fs)
legend("topright", c("(mean_M - mean)/u", "u_M/u - 1"), col = c(BLU, ORA), lwd = 1.5, bty = "n", cex = 0.45 * fs)
page(); Y <- yB - 0.02
## conformity assessment (optional)
if (!is.null(CA)) {
head(paste0("Conformity assessment (JCGM 106:2012; decision rules ISO 14253-1), ", CA$method))
para(paste0("Tolerance interval [", format(CA$lower), ", ", format(CA$upper), "]; estimate y = ", g4(CA$estimate, 7), "; conformance probability p_c = ",
formatC(CA$pc, format = "f", digits = 4), " (1 - p_c = ", g4(CA$pnc, 3), "); expanded uncertainty U- = ", g4(CA$U[1], 3),
", U+ = ", g4(CA$U[2], 3), " (p = ", CA$p, ")."), cex = 0.6, font = 2)
RU <- CA$rules
if (!is.null(RU) && nrow(RU)) {
lim <- c(CA$lower, CA$upper, RU$RL, RU$RU, RU$AL, RU$AU, CA$estimate); lim <- lim[is.finite(lim)]
xr <- range(lim) + c(-0.1, 0.1) * diff(range(lim)); X <- function(v) 0.36 + 0.44 * (pmin(pmax(v, xr[1]), xr[2]) - xr[1])/diff(xr)
rh <- 0.017 * fs; need(rh * (nrow(RU) + 1.5))
for (i in seq_len(nrow(RU))) {
r <- RU[i, ]; yc <- Y; hb <- rh * 0.3
text(0.05, yc, fitTxt(r$description, 0.3, 0.52), adj = c(0, 0.5), cex = 0.52 * fs)
rect(X(xr[1]), yc - hb, X(xr[2]), yc + hb, col = "#C8102E", border = NA) ## reject
rect(X(r$RL), yc - hb, X(r$RU), yc + hb, col = "grey75", border = NA) ## uncertainty zone
if (r$AL < r$AU) rect(X(r$AL), yc - hb, X(r$AU), yc + hb, col = "#2E9D48", border = NA) ## accept
segments(X(c(CA$lower, CA$upper)), yc - hb * 1.4, X(c(CA$lower, CA$upper)), yc + hb * 1.4, lwd = 1.5)
points(X(CA$estimate), yc + hb * 1.9, pch = 25, bg = "black", cex = 0.6)
dcol <- switch(r$decision, accept = "#2E9D48", reject = "#C8102E", "grey35")
text(0.82, yc, toupper(r$decision), adj = c(0, 0.5), cex = 0.6 * fs, font = 2, col = dcol)
Y <- Y - rh }
tk <- pretty(xr); tk <- tk[tk >= xr[1] & tk <= xr[2]]
segments(X(xr[1]), Y + rh * 0.4, X(xr[2]), Y + rh * 0.4, col = "grey40"); segments(X(tk), Y + rh * 0.4, X(tk), Y + rh * 0.15, col = "grey40")
text(X(tk), Y - rh * 0.05, format(tk), cex = 0.5 * fs, col = "grey30")
text(X(c(CA$lower, CA$upper)), Y - rh * 0.55, c("LSL", "USL"), cex = 0.5 * fs, font = 2)
Y <- Y - rh * 1.2
para("Green: acceptance zone; grey: uncertainty zone (no decision); red: rejection; triangle: estimate; vertical lines: tolerance limits.", cex = 0.52, col = "#3a3a3a")
}
}
need(0.04 + (N + 2) * 0.0135); head("Sensitivity coefficients and Hessian matrix")
gS <- symVec(O$gradient); if (length(gS) != N) gS <- rep("", N)
lhD <- if (N > 12) 0.0125 else 0.0135
DER <- data.frame(Quantity = nm, "c_i = df/dx_i, symbolic" = fitTxt(gS, if (N <= 6) 0.33 else 0.5), "c_i, evaluated" = g4(gr),
check.names = FALSE, stringsAsFactors = FALSE)
if (is.matrix(H) && N <= 6) { ## evaluated Hessian as extra columns of the gradient table
cw <- min(0.08, 0.34/N); HE <- fmtM(matrix(H, N, dimnames = list(nm, nm)))
DH <- cbind(DER, as.data.frame(HE, stringsAsFactors = FALSE)); names(DH)[-(1:3)] <- paste0("H[ , ", nm, "]")
tab(DH, c(0.1, 0.35, 0.11, rep(cw, N)), align = c("l", "l", rep("r", N + 1)), lh = lhD, cex = 0.56)
hS <- symVec(O$hessian)
if (length(hS) == N^2 && any(H != 0) && N > 3)
para("The symbolic Hessian matrix (more than 3 input quantities) is printed by summary().", cex = 0.54)
if (length(hS) == N^2 && any(H != 0) && N <= 3) { ## symbolic Hessian (lower triangle; H is symmetric)
hS <- matrix(hS, N, dimnames = list(nm, nm)); hS[upper.tri(hS)] <- ""
cw2 <- min(0.28, 0.84/N); hS[] <- fitTxt(hS, cw2 - 0.01)
mrow(list(hS), "Hessian H = d^2 f/(dx_i dx_j), symbolic (lower triangle)", 0.05, cw2, (N + 3) * lhD, align = rep("l", N + 1))
}
if (all(H == 0)) para("H = 0: the model is linear in all input quantities at the estimates.", cex = 0.54)
} else {
tab(DER, c(0.14, 0.52, 0.16), align = c("l", "l", "r"), lh = lhD, cex = 0.56)
if (is.matrix(H)) para(paste0("Hessian matrix (", N, " x ", N, ") not shown; non-zero elements: ", sum(H != 0), "."), cex = 0.54)
}
## covariance check of the Monte Carlo input samples (as in summary.propagate)
need(0.04 + (min(N, 8) + 3) * 0.0135); head("Covariance check of the Monte Carlo input samples")
TGT <- if (!is.null(O$covMat.target)) O$covMat.target else S; XS <- O$datSIM
noiseRef <- function(FUN, target) { if (is.null(XS) || nrow(XS) < 20 || !all(is.finite(target))) return(NA)
n <- nrow(XS) %/% 2; M1 <- FUN(XS[1:n, , drop = FALSE], use = "pairwise.complete.obs")
M2 <- FUN(XS[(n + 1):(2 * n), , drop = FALSE], use = "pairwise.complete.obs"); norm(M1 - M2, "F")/(2 * norm(target, "F")) }
judge <- function(eps, noise) { if (!is.finite(eps)) return(list("not applicable (expected covariance does not exist)", "grey30"))
if (!is.finite(noise) || noise == 0) return(list(format(signif(eps, 3)), "black")); rr <- eps/noise
list(paste0(signif(eps, 3), " (", signif(rr, 3), " x MC noise level => ", if (rr > 3) "CHECK" else "OK", ")"), if (rr > 3) RED else GRN) }
CM <- list(S, TGT, if (is.matrix(O$covMat.MC)) O$covMat.MC else S * NA)
CM <- lapply(CM, function(M) { dimnames(M) <- list(nm, nm); fmtM(M) })
caps <- c("Input covariance (GUF)", "MC input covariance, expected", "MC design covariance, observed")
if (N <= 4) mrow(CM, caps, c(0.05, 0.36, 0.67), min(0.07, 0.23/N), (N + 3) * 0.0135)
else if (N <= 8) for (j in 1:3) mrow(CM[j], caps[j], 0.05, 0.09, (N + 3) * 0.0135)
JC <- judge(num1(O$frobCOV), noiseRef(cov, TGT))
para(paste("Relative Frobenius norm, expected vs. observed covariance:", JC[[1]]), cex = 0.58, col = JC[[2]], font = 2)
R0 <- O$corMat
if (is.matrix(R0) && ncol(R0) > 1 && any(abs(R0[upper.tri(R0)]) > 0)) {
JR <- judge(num1(O$frobCOR), noiseRef(cor, R0))
para(paste("Relative Frobenius norm, input vs. observed correlation:", JR[[1]]), cex = 0.58, col = JR[[2]], font = 2) }
if (heavy) para("Note: DOF <= 4 => the t-margins have infinite kurtosis; covariance estimates converge slowly and the check is unreliable.", cex = 0.54)
if (O$type %in% c("stat", "raw")) para(paste0("MC input models for Type A quantities: S1 = scaled and shifted t-distributions, scale u(x), nu = n - 1 ",
"(JCGM 101:2008, 6.4.9), correlated inputs joined by a Gaussian copula (Possolo 2010); mvt = multivariate t, nu = n - N (JCGM 102:2011, 5.3.2)."),
cex = 0.54, col = "#3a3a3a")
## reproducibility record
head("Reproducibility (ISO/IEC 17025:2017, 7.5, 7.8)")
for (k in names(REPRO)) {
need(0.015); text(0.05, Y, paste0(k, ":"), adj = c(0, 0.5), cex = 0.56 * fs, font = 2)
para(REPRO[[k]], cex = 0.54, family = if (k == "Call") "mono" else "", x = 0.16, w = 0.79); Y <- Y + 0.006 }
OUT[[ia]] <- list(results = RES, validation = VAL, accuracy = ACC, budget = BUD[, -ncol(BUD)], sobol = SOB, nonlinearity = NL,
derivatives = DER, shape = SHP, checks = CHK, conformity = CA, fit = FIT, reproducibility = REPRO)
}
invisible(c(list(file = file), if (nA == 1) OUT[[1]] else list(analyses = OUT)))
}
Try the propagate package in your browser
Any scripts or data that you put into this service are public.
propagate documentation built on Nov. 10, 2026, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions report
Documented in report
report <- function(object, file = tempfile(fileext = ".pdf"), title = "Uncertainty evaluation", measurand = NULL, unit = "",
compare = TRUE, nsim.compare = 1E5, fitDistr = FALSE, fit.nmax = 1E4,
gauge.mean = NULL, gauge.u = NULL, gauge.method = c("GUF", "MC"), gauge.dec = NULL, gauge.M = 1E6,
interval = c("symmetric", "shortest"), ndig = 2, sobol = TRUE, conformity = NULL, seed = NULL,
fontsize = 1)
{
## one object or a list of objects (one analysis each); per-analysis arguments are recycled
OBJ <- if (inherits(object, "propagate")) list(object) else object
if (!is.list(OBJ) || !length(OBJ) || !all(sapply(OBJ, inherits, "propagate")))
stop("report: 'object' must be a 'propagate' object or a list of them!")
nA <- length(OBJ)
per <- function(x, i) { if (is.null(x)) return(NULL); v <- rep_len(if (is.list(x)) x else as.list(x), nA)[[i]]
if (length(v) == 1 && !is.list(v) && is.na(v)) NULL else v }
## 'conformity': one specification list(lower, upper, ...) for all analyses, or a list of them (NULL/NA = none)
CONF <- if (is.null(conformity)) NULL else if (any(c("lower", "upper") %in% names(conformity))) list(conformity) else conformity
if (length(title) < nA && !is.null(names(OBJ)) && identical(title, "Uncertainty evaluation")) title <- names(OBJ)
interval <- match.arg(interval); gauge.method <- match.arg(gauge.method); fs <- fontsize
RED <- "#8b0000"; GRN <- "#006400"; VIO <- "#6a1fb0"; BLU <- "#1f4e79"; BAR <- "#4682b4"; ORA <- "#ff8c00"
f <- function(x, d) formatC(as.numeric(x), format = "f", digits = d)
g4 <- function(x, d = 4) formatC(signif(as.numeric(x), d), format = "g", digits = d)
num1 <- function(x) if (is.numeric(x) && length(x) && is.finite(x[1])) x[1] else NA_real_
PV <- tryCatch(as.character(utils::packageVersion("propagate")), error = function(e) "")
## ---------------------------------------------------------------- page layout helpers (shared)
if (!is.null(file)) { grDevices::cairo_pdf(file, width = 8.27, height = 11.69, onefile = TRUE); on.exit(grDevices::dev.off()) }
op <- par(mar = c(0, 0, 0, 0)); on.exit(par(op), add = TRUE, after = FALSE)
PAGE <- 0; Y <- 0; TITLE <- title[1]
newPage <- function() { plot.new(); plot.window(c(0, 1), c(0, 1)); PAGE <<- PAGE + 1
rect(0.04, 0.955, 0.96, 0.985, col = BLU, border = NA)
text(0.055, 0.97, TITLE, adj = c(0, 0.5), col = "white", font = 2, cex = 1.15 * fs)
text(0.945, 0.97, format(Sys.Date()), adj = c(1, 0.5), col = "white", cex = 0.75 * fs)
text(0.5, 0.02, paste0("propagate", if (nzchar(PV)) paste0(" ", PV), " - page ", PAGE), cex = 0.65 * fs, col = "grey40")
Y <<- 0.93 }
need <- function(h) if (Y - h < 0.035) newPage()
head <- function(t, x = 0.05) { need(0.05); text(x, Y, t, adj = c(0, 0.5), font = 2, col = BLU, cex = 0.95 * fs); Y <<- Y - 0.026 }
para <- function(t, cex = 0.62, col = "black", font = 1, family = "", x = 0.05, w = 0.9) {
nc <- floor(w/strwidth("m", cex = cex * fs, family = family) * if (nzchar(family)) 1 else 1.6)
for (l in unlist(lapply(t, strwrap, width = nc))) { need(0.015)
text(x, Y, l, adj = c(0, 0.5), cex = cex * fs, col = col, font = font, family = family); Y <<- Y - 0.0125 * fs }
Y <<- Y - 0.008 }
## table; 'bars' (and 'bars2'): shares drawn as horizontal bars in the last (empty) column, upper and lower half
tab <- function(D, w, x0 = 0.05, cex = 0.6, colr = NULL, fill = NULL, bars = NULL, bars2 = NULL, lh = 0.0145, align = NULL) {
lh <- lh * fs; xs <- x0 + c(0, cumsum(w)); nc <- ncol(D)
if (is.null(align)) align <- c("l", rep("r", nc - 1)) ## default: first column left-, others right-aligned
L <- align == "l"; xt <- ifelse(L, xs[-(nc + 1)] + 0.005, xs[-1] - 0.005)
row <- function(v, ...) { v <- as.character(v); if (any(L)) text(xt[L], Y, v[L], adj = c(0, 0.5), ...)
if (any(!L)) text(xt[!L], Y, v[!L], adj = c(1, 0.5), ...) }
hdr <- function() { rect(x0, Y - lh/2, max(xs), Y + lh/2, col = "#e2e2e2", border = NA)
row(names(D), font = 2, cex = cex * fs); Y <<- Y - lh }
bmax <- max(1, abs(as.numeric(c(bars, bars2))), na.rm = TRUE)
drawBar <- function(b, yc, hh, col) if (!is.null(b) && is.finite(b) && b != 0) {
wb <- (xs[nc + 1] - xs[nc] - 0.01) * abs(b)/bmax
rect(xs[nc] + 0.005, yc - hh, xs[nc] + 0.005 + wb, yc + hh, border = NA, col = if (b < 0) RED else col) }
need(2 * lh); hdr()
for (i in seq_len(nrow(D))) {
if (Y - lh < 0.035) { newPage(); hdr() }
if (!is.null(fill) && !is.na(fill[i])) rect(xs[nc], Y - lh/2, xs[nc + 1], Y + lh/2, col = fill[i], border = NA)
if (is.null(bars2)) drawBar(bars[i], Y, lh * 0.32, BAR) else {
drawBar(bars[i], Y + lh * 0.17, lh * 0.15, BAR); drawBar(bars2[i], Y - lh * 0.17, lh * 0.15, VIO) }
cl <- if (!is.null(colr) && !is.na(colr[i])) colr[i] else "black"
row(unlist(D[i, ]), cex = cex * fs, col = cl, font = if (cl == "black") 1 else 2)
segments(x0, Y - lh/2, max(xs), Y - lh/2, col = "#dddddd", lwd = 0.5); Y <<- Y - lh }
Y <<- Y - 0.008 }
## symbolic derivative as one-line string; shortened with "..." to width w
symStr <- function(e) if (is.language(e) || is.numeric(e)) paste(deparse(e, width.cutoff = 500L), collapse = "") else as.character(e)
symVec <- function(G) if (is.list(G)) vapply(G, symStr, "") else as.character(G)
fitTxt <- function(x, w, cex = 0.6) vapply(x, function(z) { if (strwidth(z, cex = cex * fs) <= w) return(z)
while (nchar(z) > 2 && strwidth(paste0(z, "..."), cex = cex * fs) > w) z <- substr(z, 1, nchar(z) - 1); paste0(z, "...") }, "", USE.NAMES = FALSE)
## a matrix with caption as table; cells are character
mtab <- function(M, cap, x0, cw, cex = 0.6, align = NULL) {
text(x0, Y, cap, adj = c(0, 0.5), cex = cex * fs, font = 3, col = "grey20"); Y <<- Y - 0.0135 * fs
D <- data.frame(rownames(M), matrix(M, nrow(M)), check.names = FALSE, stringsAsFactors = FALSE); names(D) <- c(" ", colnames(M))
tab(D, c(0.06, rep(cw, ncol(M))), x0 = x0, cex = cex, align = align) }
## several matrices side by side (all start at the same height)
mrow <- function(L, caps, x0s, cw, h, align = NULL) { need(h); y0 <- Y; ym <- Y
for (j in seq_along(L)) { Y <<- y0; mtab(L[[j]], caps[j], x0s[j], cw, align = align); ym <- min(ym, Y) }; Y <<- ym }
fmtM <- function(M) { M[M == 0] <- 0; F <- matrix(formatC(signif(M, 3), format = "g", digits = 3), nrow(M), dimnames = dimnames(M)); F }
## a sub-figure in user coordinates [x0, x1] x [y0, y1] of the page; back to the page with page()
ndc <- function(x, y) c(grconvertX(x, "user", "ndc"), grconvertY(y, "user", "ndc")) ## call while in page coordinates
subfig <- function(fig, mar) par(fig = fig, new = TRUE, mar = mar)
page <- function() { par(fig = c(0, 1, 0, 1), new = TRUE, mar = c(0, 0, 0, 0)); plot.new(); plot.window(c(0, 1), c(0, 1)) }
OUT <- vector("list", nA); names(OUT) <- names(OBJ)
for (ia in seq_len(nA)) {
O <- OBJ[[ia]]; TITLE <- per(title, ia); if (is.null(TITLE)) TITLE <- "Uncertainty evaluation"
MSi <- per(measurand, ia); UNi <- per(unit, ia); if (is.null(UNi)) UNi <- ""
GM <- per(gauge.mean, ia); GU <- per(gauge.u, ia); GD <- per(gauge.dec, ia)
SEED <- per(seed, ia); CF <- if (is.null(CONF)) NULL else per(CONF, ia)
u1 <- unname(O$prop["u.1"]); dec <- if (is.finite(u1) && u1 > 0) max(0, 3 - floor(log10(u1))) else 4
ciMC <- unname(if (interval == "shortest") O$ci.short else O$sim[5:6])
if (is.null(MSi)) { n <- names(O$expr); MSi <- if (length(n) && nzchar(n[1])) n[1] else "y" }
sym <- function(s) { e <- tryCatch(str2lang(gsub("_(\\w+)", "[\\1]", s)), error = function(e) NULL); if (is.null(e)) s else e }
MS <- sym(MSi)
UN <- { map <- c(Ohm = "Omega", degC = "degree*C", um = "mu*m", ug = "mu*g", uV = "mu*V", uA = "mu*A", "%" = "'%'")
if (!nzchar(UNi)) quote("") else if (UNi %in% names(map)) str2lang(map[[UNi]]) else
tryCatch({ e <- str2lang(UNi); if (is.name(e) || is.call(e)) e else UNi }, error = function(e) UNi) }
y <- O$resSIM[is.finite(O$resSIM)]
## ---------------------------------------------------------------- results of all methods
mcRow <- function(R, lab) { ci <- unname(if (interval == "shortest") R$ci.short else R$sim[5:6])
c(paste0("MC, ", lab, ", M = ", format(length(R$resSIM), big.mark = " ", scientific = FALSE)), f(R$sim["mean.MC"], dec),
if (any(R$df.MC <= 2)) "n/a" else f(R$sim["u.MC"], dec + 1), paste0("[", f(ci[1], dec), ", ", f(ci[2], dec), "]")) }
ciG <- .ciGUF(O$prop) ## v1.3-0: by name
RES <- rbind(c("GUF, 1st order", f(O$prop["mean.1"], dec), f(u1, dec + 1),
paste0("[", f(ciG[1], dec), ", ", f(ciG[2], dec), "] (k = ", f(num1(O$k), 3), ")")),
c("GUF, 2nd order", f(O$prop["mean.2"], dec), f(O$prop["u.2"], dec + 1), "-"),
if (is.finite(num1(O$prop["u.2m"]))) ## v1.3-0: second order with input moments
c("GUF, 2nd order, input moments", f(O$prop["mean.2"], dec), f(O$prop["u.2m"], dec + 1), "-"))
MCR <- list(); MCR[[if (is.na(O$mc.margins)) "samples" else O$mc.margins]] <- O
if (compare && O$type %in% c("stat", "raw") && !is.null(O$input$data))
for (m in setdiff(c("S1", "mvt"), O$mc.margins)) {
R <- NULL; invisible(capture.output(R <- tryCatch(suppressMessages(propagate(O$expr, O$input$data, type = O$type, cov = O$input$cov,
mc.margins = m, nsim = nsim.compare, alpha = O$alpha, ws.group = O$input$ws.group, check = FALSE)), error = function(e) NULL)))
if (!is.null(R)) MCR[[m]] <- R }
for (m in intersect(c("S1", "mvt", "samples"), names(MCR))) RES <- rbind(RES, mcRow(MCR[[m]], m))
RES <- as.data.frame(RES, stringsAsFactors = FALSE)
names(RES) <- c("Method", "Estimate", "u", paste0(signif(100 * (1 - O$alpha), 3), " % interval", if (interval == "shortest") " (shortest)" else ""))
resCol <- ifelse(grepl("^GUF, 1st|^MC, S1", RES$Method), RED, NA)
## ---------------------------------------------------------------- validation of the GUF by the MCM (JCGM 101:2008, 8)
VAL <- NULL; valCol <- NA
uG <- if (is.finite(num1(O$u.exp)) && is.finite(num1(O$k)) && O$k > 0) O$u.exp/O$k else NA
if (is.finite(uG) && uG > 0 && all(is.finite(ciG)) && all(is.finite(ciMC))) {
lv <- floor(log10(uG)) - (ndig - 1); delta <- 0.5 * 10^lv
dlo <- abs(ciG[1] - ciMC[1]); dhi <- abs(ciG[2] - ciMC[2])
SE <- O$sim.se
sE <- if (is.null(SE)) NA else if (interval == "shortest") max(SE[c("short.low", "short.high")]) else max(SE[c("sym.low", "sym.high")])
verdict <- if (!is.finite(sE)) "not assessable (no block results)" else if (2 * sE > delta) "INCONCLUSIVE: increase nsim" else
if (max(dlo, dhi) <= delta) "GUF VALIDATED" else "GUF NOT VALIDATED: use MC"
valCol <- if (grepl("^GUF VALIDATED", verdict)) GRN else if (grepl("NOT", verdict)) RED else ORA
cover <- mean(y >= ciG[1] & y <= ciG[2])
VAL <- data.frame(Item = c(paste0("numerical tolerance delta (", ndig, " significant digits of u)"),
"d_low = |low(GUF) - low(MC)|", "d_high = |high(GUF) - high(MC)|", "2 s of the MC endpoints",
"result", "MC coverage of the GUF interval"),
Value = c(g4(delta, 2), g4(dlo, 3), g4(dhi, 3), g4(2 * sE, 2), verdict,
paste0(f(100 * cover, 1), " % (nominal ", signif(100 * (1 - O$alpha), 3), " %)")),
stringsAsFactors = FALSE)
}
## ---------------------------------------------------------------- numerical accuracy of the MC results (JCGM 101:2008, 7.9)
ACC <- NULL
if (!is.null(O$sim.se)) {
SE <- O$sim.se; nmL <- if (interval == "shortest") c("short.low", "short.high") else c("sym.low", "sym.high")
v <- c(unname(O$sim["mean.MC"]), unname(O$sim["u.MC"]), ciMC); s <- unname(c(SE["mean.MC"], SE["u.MC"], SE[nmL]))
ok <- is.finite(v) & is.finite(s) & s > 0
dS <- ifelse(ok, floor(-log10(4 * s)), NA) ## 2 s <= 1/2 10^-d: d decimals are stable
sig <- ifelse(ok & v != 0, pmax(0, floor(log10(abs(v))) + dS + 1), NA)
rnd <- function(x, d) if (!is.finite(d)) "n/a" else if (d >= 0) formatC(round(x, d), format = "f", digits = d) else format(round(x, d))
ACC <- data.frame(Quantity = c("estimate", "u", "interval, low", "interval, high"),
"stable value" = mapply(rnd, v, dS), "2 s" = g4(2 * s, 2), digits = ifelse(is.na(sig), "-", sig),
check.names = FALSE, stringsAsFactors = FALSE)
if (any(O$df.MC <= 2)) ACC[2, 2:4] <- c("n/a", "-", "-")
}
## ---------------------------------------------------------------- input quantities: distribution and Type A/B
gr <- as.numeric(O$evalGrad); gr[gr == 0] <- 0; S <- as.matrix(O$covMat); nm <- colnames(S); N <- length(nm)
nu <- if (nrow(O$data) >= 3) as.numeric(O$data[3, nm]) else rep(Inf, N)
BG <- if (!is.null(O$input$budget)) O$input$budget else attr(O$input$data, "budget")
dname <- c(normal = "Gaussian", t = "t", rect = "rectangular", tri = "triangular", arcsine = "arcsine", ctrap = "curvil. trapezoid")
if (!is.null(BG) && all(nm %in% BG$name)) {
b <- BG[match(nm, BG$name), ]; dd <- as.character(b$dist); bdf <- as.numeric(b$df)
DIST <- ifelse(dd %in% names(dname), dname[dd], dd); DIST[is.na(DIST)] <- "sample"
TYP <- ifelse(dd == "t" & is.finite(bdf), "A", ifelse(is.na(dd), "-", "B"))
} else if (identical(O$type, "sim")) { DIST <- rep("sample", N); TYP <- rep("-", N)
} else {
mvt <- identical(O$mc.margins, "mvt"); A <- O$type == "raw" | (is.finite(nu) & nu < 1000)
DIST <- ifelse(A, if (mvt) "mvt" else "t", "Gaussian"); TYP <- ifelse(A, "A", "B")
}
## ---------------------------------------------------------------- uncertainty budget with first-order and Sobol' indices
Kc <- outer(gr, gr) * S; uy2 <- sum(Kc)
share <- diag(Kc)/uy2; cv <- uy2 - sum(diag(Kc))
indep <- abs(cv) <= 1E-12 * abs(uy2) && !identical(O$mc.margins, "mvt") &&
(!is.matrix(O$corMat) || all(abs(O$corMat[upper.tri(O$corMat)]) < 1E-12))
SOB <- if (isTRUE(sobol) && indep) tryCatch(suppressMessages(get("sobol", mode = "function")(O)), error = function(e) NULL) else NULL
Si <- if (!is.null(SOB)) pmax(0, as.numeric(SOB$S[nm])) else rep(NA, N)
BUD <- data.frame(Quantity = nm, Estimate = formatC(O$data[1, nm], format = "fg", digits = 8), "u(x_i)" = g4(sqrt(diag(S))),
nu_i = ifelse(!is.finite(nu) | nu >= 1000, "Inf", formatC(nu, format = "g")), "c_i" = g4(gr),
"|c_i| u(x_i)" = g4(abs(gr) * sqrt(diag(S))), Index = paste0(f(100 * share, 1), " %"),
"S_i (MC)" = if (is.null(SOB)) rep("-", N) else paste0(f(100 * Si, 1), " %"), " " = "",
check.names = FALSE, stringsAsFactors = FALSE)
BARS <- share; BARS2 <- if (is.null(SOB)) NULL else Si
if (abs(cv) > 1E-12 * abs(uy2)) { BUD <- rbind(BUD, c("covariance terms", rep("", 5), paste0(f(100 * cv/uy2, 1), " %"), "-", ""))
BARS <- c(BARS, cv/uy2); if (!is.null(BARS2)) BARS2 <- c(BARS2, NA) }
if (!is.null(SOB) && 1 - sum(Si) > 0.005) { BUD <- rbind(BUD, c("interactions (MC)", rep("", 5), "-", paste0(f(100 * (1 - sum(Si)), 1), " %"), ""))
BARS <- c(BARS, NA); BARS2 <- c(BARS2, 1 - sum(Si)) }
BUD <- rbind(BUD, c("u(y), first order", rep("", 4), g4(sqrt(uy2)), "100.0 %", if (is.null(SOB)) "-" else "100.0 %", ""))
BARS <- c(BARS, NA); if (!is.null(BARS2)) BARS2 <- c(BARS2, NA)
H <- O$evalHess
zeroC <- if (is.matrix(H) && all(dim(H) == length(gr))) nm[gr == 0 & rowSums(abs(H)) > 0] else character(0)
## ---------------------------------------------------------------- nonlinearity of the model and shape of the output
P2 <- O$prop; NL <- NULL
yy <- y - mean(y); m2 <- mean(yy^2); SH <- c(skewness = mean(yy^3)/m2^1.5, kurtosis = mean(yy^4)/m2^2 - 3)
heavy <- (!is.null(O$df.MC) && any(O$df.MC <= 4))
lev <- function(x, a, b) if (!is.finite(x)) "-" else if (abs(x) < a) "negligible" else if (abs(x) < b) "moderate" else "large"
uS <- if (is.finite(num1(P2["u.2m"]))) "u.2m" else "u.2" ## v1.3-0: with input moments if available
if (all(is.finite(P2[c("mean.1", "mean.2", "u.1", uS)])) && P2[uS] > 0) {
nlShare <- unname((P2[uS]^2 - P2["u.1"]^2)/P2[uS]^2)
nlShift <- if (P2["u.1"] > 0) unname((P2["mean.2"] - P2["mean.1"])/P2["u.1"]) else NA
qy <- quantile(y, c(0.025, 0.5, 0.975)) ## quantile skewness: robust, exists for any nu
sk <- unname((qy[3] + qy[1] - 2 * qy[2])/(qy[3] - qy[1]))
NL <- data.frame(Measure = c(paste0("share of the second-order terms in ", uS, "^2"), "shift of the estimate (mean.2 - mean.1)/u.1",
"quantile skewness of the MC output, (q97.5 + q2.5 - 2 q50)/(q97.5 - q2.5)",
"skewness of the MC output (moments; Gaussian: 0)", "excess kurtosis of the MC output (moments; Gaussian: 0)"),
Value = c(paste0(formatC(100 * nlShare, format = "fg", digits = 3), " %"), formatC(nlShift, format = "fg", digits = 3),
formatC(c(sk, SH), format = "fg", digits = 3)),
Assessment = c(lev(nlShare, 0.01, 0.1), lev(nlShift, 0.05, 0.2), lev(sk, 0.05, 0.2),
if (heavy) c("-", "-") else c(lev(SH[1], 0.1, 0.5), lev(SH[2], 0.2, 1))), stringsAsFactors = FALSE)
}
SHP <- data.frame(Measure = c("skewness", "excess kurtosis"), Value = formatC(SH, format = "fg", digits = 4), stringsAsFactors = FALSE)
## ---------------------------------------------------------------- comparison with gauge values
CHK <- NULL
if (!is.null(GM) || !is.null(GU)) {
decOf <- function(x) max(sapply(x, function(v) { s <- format(v, scientific = FALSE, digits = 15, drop0trailing = TRUE)
if (grepl(".", s, fixed = TRUE)) nchar(sub("^[^.]*\\.", "", s)) else 0 }))
dd <- rep_len(if (is.null(GD)) c(if (is.null(GM)) 0 else decOf(GM), if (is.null(GU)) 0 else decOf(GU)) else GD, 2)
isMC <- gauge.method == "MC"; se <- O$sim.se; fac <- sqrt(1 + O$nsim/gauge.M)
obs <- if (isMC) unname(O$sim[c("mean.MC", "u.MC")]) else c(unname(O$prop["mean.1"]), u1)
tol <- 0.5 * 10^-dd + if (isMC && !is.null(se)) 4 * fac * unname(se[c("mean.MC", "u.MC")]) else 0
pub <- c(if (is.null(GM)) NA else GM, if (is.null(GU)) NA else GU)
CHK <- data.frame(Item = paste0(if (isMC) "Monte Carlo" else "GUF, 1st order", c(": estimate", ": u")),
Gauge = c(f(pub[1], dd[1]), f(pub[2], dd[2])), propagate = c(f(obs[1], dd[1] + 1), f(obs[2], dd[2] + 1)),
Tolerance = g4(tol, 2), Result = ifelse(abs(obs - pub) <= tol, "PASS", "FAIL"), stringsAsFactors = FALSE)[!is.na(pub), ]
}
## ---------------------------------------------------------------- distribution fit, conformity assessment
FIT <- if (isTRUE(fitDistr)) tryCatch(suppressMessages(get("fitDistr", mode = "function")(y, nmax = fit.nmax, verbose = FALSE, plot = "none")),
error = function(e) NULL) else NULL
CA <- NULL
if (!is.null(CF)) {
CF <- modifyList(list(lower = -Inf, upper = Inf, method = "MC", rule = c("stringent", "simple"), p = 1 - O$alpha, interval = interval), CF)
CA <- tryCatch(get("conformity", mode = "function")(O, lower = CF$lower, upper = CF$upper, method = CF$method, rule = CF$rule,
p = CF$p, interval = CF$interval), error = function(e) { message("report: conformity: ", conditionMessage(e)); NULL })
}
## ---------------------------------------------------------------- reproducibility record (ISO/IEC 17025:2017, 7.5, 7.8)
CALL <- if (!is.null(O$call)) paste(deparse(O$call, width.cutoff = 500L), collapse = " ") else
paste0("propagate(expr = ", paste(deparse(O$expr, width.cutoff = 500L), collapse = ""), ", data = <", O$type, " data>, type = \"", O$type,
"\", cov = ", if (is.matrix(O$input$cov)) "<matrix>" else paste0("\"", O$input$cov, "\""), ", nsim = ", format(O$nsim, scientific = TRUE),
", alpha = ", O$alpha, if (!is.na(O$mc.margins)) paste0(", mc.margins = \"", O$mc.margins, "\""),
if (!is.null(O$input$exp.uc)) paste0(", exp.uc = \"", O$input$exp.uc, "\""), ") [reconstructed]")
if (nchar(CALL) > 400) CALL <- paste0(substr(CALL, 1, 400), " ...")
RK <- RNGkind()
REPRO <- c(Call = CALL,
"Monte Carlo" = paste0("M = ", format(length(O$resSIM), big.mark = " ", scientific = FALSE), " trials",
if (!is.null(O$adaptive)) paste0(" (adaptive, ", if (O$adaptive$converged) "converged" else "not converged", ")"),
"; alpha = ", O$alpha, "; input model: ", if (is.na(O$mc.margins)) "user-supplied samples" else O$mc.margins,
"; cov = ", if (is.matrix(O$input$cov)) "user matrix" else O$input$cov,
"; coverage interval: ", interval),
"RNG, software" = paste0("seed ", if (is.null(SEED)) "not recorded (use report(..., seed = ))" else paste(SEED, collapse = ", "),
"; RNG ", RK[1], "/", RK[2], "/", RK[3], "; propagate ", if (nzchar(PV)) PV else "(version unknown)",
"; R ", R.version$major, ".", R.version$minor, " (", R.version$platform, "); ", format(Sys.time(), "%Y-%m-%d %H:%M %Z")))
## ================================================================ page 1
newPage()
head("Measurement model")
text(0.08, Y - 0.006, bquote(italic(.(MS)) == .(O$expr[[1]])), adj = c(0, 0.5), cex = 0.95 * fs); Y <- Y - 0.042
head("Input data")
D <- O$input$data; y0 <- Y; corDone <- FALSE
IN <- data.frame(Quantity = nm, Estimate = BUD$Estimate[seq_len(N)], "u(x_i)" = BUD[["u(x_i)"]][seq_len(N)], nu_i = BUD$nu_i[seq_len(N)],
Distribution = DIST, Type = TYP, check.names = FALSE, stringsAsFactors = FALSE)
if (identical(O$type, "raw") && !is.null(D) && nrow(D) <= 15) {
RAW <- data.frame(k = seq_len(nrow(D)), apply(D[, nm, drop = FALSE], 2, function(v) formatC(v, format = "g", digits = 6)), check.names = FALSE)
tab(RAW, c(0.03, rep(min(0.09, 0.3/N), N)), lh = 0.0135); y1 <- Y; Y <- y0
names(IN)[5] <- "Distr."
tab(IN, c(0.07, 0.1, 0.09, 0.04, 0.07, 0.04), x0 = 0.5, align = c("l", "r", "r", "r", "l", "l"), lh = 0.0135)
if (N > 1 && N <= 6 && any(abs(cov2cor(S)[upper.tri(S)]) > 1E-12)) {
CR <- data.frame(r = nm, matrix(f(cov2cor(S), 3), N), check.names = FALSE); names(CR)[-1] <- nm
tab(CR, c(0.07, rep(min(0.07, 0.34/N), N)), x0 = 0.5, lh = 0.0135) }
Y <- min(Y, y1); corDone <- TRUE
} else {
tab(IN, c(0.14, 0.16, 0.13, 0.07, 0.14, 0.06), align = c("l", "r", "r", "r", "l", "l"), lh = if (N > 12) 0.0125 else 0.0135)
}
if (!corDone && N > 1 && N <= 6 && any(abs(cov2cor(S)[upper.tri(S)]) > 1E-12)) {
CR <- data.frame(r = nm, matrix(f(cov2cor(S), 3), N), check.names = FALSE); names(CR)[-1] <- nm
tab(CR, c(0.08, rep(0.08, N)), lh = 0.0135)
}
head("Result statement")
dU <- function(u) max(0, 1 - floor(log10(abs(u))))
fU <- function(u) if (is.finite(u) && u > 0) formatC(signif(u, 2), format = "f", digits = dU(u)) else "n/a"
fY <- function(y, u) if (is.finite(u) && u > 0) formatC(round(y, dU(u)), format = "f", digits = dU(u)) else format(signif(y, 6))
uM <- if (any(O$df.MC <= 2)) NA else unname(O$sim["u.MC"]); uR <- if (is.finite(uM)) uM else diff(ciMC)/4
PP <- paste0(signif(100 * (1 - O$alpha), 3), " %")
nuS <- if (is.na(num1(O$ws.df))) "-" else if (O$ws.df >= 1000) "infinity" else as.character(O$ws.df)
need(0.12); top <- Y + 0.012; bh <- 0.1 * fs
fit1 <- function(x, y, e, cex) text(x, y, e, cex = min(cex, 0.4/strwidth(e, cex = 1)))
text(0.17, top - 0.014, "GUM uncertainty framework", adj = c(0, 0.5), cex = 0.8 * fs, col = "grey20")
fit1(0.28, top - 0.045, bquote(italic(.(MS)) ~ "=" ~ .(fY(O$prop["mean.1"], u1)) ~ .(UN) * "," ~~ italic(u) ~ "=" ~ .(fU(u1)) ~ .(UN)), 1.05 * fs)
fit1(0.28, top - 0.078, bquote(.(PP) ~ "interval:" ~ group("[", .(fY(ciG[1], u1)) * "," ~ .(fY(ciG[2], u1)), "]") ~ .(UN) ~~
group("(", italic(k) ~ "=" ~ .(f(num1(O$k), 2)) * "," ~ nu[eff] ~ "=" ~ .(nuS), ")")), 0.8 * fs)
text(0.64, top - 0.014, "Monte Carlo method", adj = c(0, 0.5), cex = 0.8 * fs, col = "grey20")
fit1(0.72, top - 0.045, bquote(italic(.(MS)) ~ "=" ~ .(fY(O$sim["mean.MC"], uR)) ~ .(UN) * "," ~~ italic(u) ~ "=" ~ .(fU(uM)) ~ .(if (is.finite(uM)) UN else quote(""))), 1.05 * fs)
fit1(0.72, top - 0.078, bquote(.(PP) ~ .(interval) ~ "interval:" ~ group("[", .(fY(ciMC[1], uR)) * "," ~ .(fY(ciMC[2], uR)), "]") ~ .(UN)), 0.8 * fs)
Y <- top - bh - 0.016
head("Results of all methods")
tab(RES, c(0.32, 0.15, 0.13, 0.3), colr = resCol, lh = 0.0135)
## validation (left) and numerical accuracy (right), side by side
if (!is.null(VAL) || !is.null(ACC)) {
need(0.13); y0 <- Y; yL <- yR <- Y
if (!is.null(VAL)) { head("Validation of the GUF (JCGM 101:2008, 8)")
tab(VAL, c(0.29, 0.15), colr = c(NA, NA, NA, NA, valCol, NA), cex = 0.56, lh = 0.0135); yL <- Y }
Y <- y0
if (!is.null(ACC)) { head("Numerical accuracy of MC (7.9)", x = 0.52)
tab(ACC, c(0.12, 0.12, 0.08, 0.06), x0 = 0.52, cex = 0.56, lh = 0.0135); yR <- Y }
Y <- min(yL, yR)
para(paste0("GUF validated if d_low, d_high <= delta and the MC endpoints are stable (2 s <= delta). s: numerical standard uncertainty ",
"from 10 blocks of the MC trials; a digit is stable if 2 s <= half a unit of it; digits: stable significant digits."), cex = 0.54, col = "#3a3a3a")
}
if (!is.null(NL)) {
head("Nonlinearity of the model and shape of the output")
tab(NL, c(0.56, 0.12, 0.14), colr = c(NA, GRN, ORA, RED)[match(NL$Assessment, c("-", "negligible", "moderate", "large"))], lh = 0.0135)
para(c(paste0("All zero for a model linear at the estimates with Gaussian inputs (limits 1 %/10 %, 0.05/0.2, 0.05/0.2, 0.1/0.5, 0.2/1); ",
"large values: the first-order GUF may be inadequate, prefer the MC results (JCGM 101:2008, 5.8)."),
if (heavy) "Inputs with DOF <= 4: the moments of the output may not exist; skewness and kurtosis are not assessed.",
if (length(zeroC)) paste0("Note: ", paste(zeroC, collapse = ", "), " with c_i = 0 contribute(s) only through second-order terms ",
"(first-order index 0 %, see the Sobol' index on page 2).")), cex = 0.54, col = "#3a3a3a")
}
if (!is.null(CHK)) {
head(paste0("Comparison with the gauge values: ", sum(CHK$Result == "PASS"), " PASS, ", sum(CHK$Result == "FAIL"), " FAIL"))
tab(CHK, c(0.34, 0.14, 0.14, 0.12, 0.1), fill = ifelse(CHK$Result == "PASS", "#d9efd0", "#f4c7c3"))
para(paste0("Tolerance: rounding of the gauge values", if (gauge.method == "MC") paste0(" plus four Monte Carlo standard errors of both ",
"calculations (gauge with ", format(gauge.M, big.mark = " ", scientific = FALSE), " trials)") else "", "."), cex = 0.56)
}
## ================================================================ (page 2) -- blocks flow on and move to the next page if they do not fit
Y <- Y - 0.01; need(0.03 + (N + 3) * 0.0155 + 0.06)
head("Uncertainty budget: first-order law of propagation and Monte Carlo (Sobol') indices")
tab(BUD, c(0.13, 0.12, 0.08, 0.04, 0.09, 0.1, 0.07, 0.08, 0.19), bars = BARS, bars2 = BARS2,
colr = ifelse(BUD$Quantity == "covariance terms" & cv < 0, RED, NA), lh = if (N > 12) 0.0135 else 0.0155)
para(c(paste0("c_i = df/dx_i; Index = c_i^2 u^2(x_i)/u^2(y) (blue); covariance terms = sum_(i != j) c_i c_j u(x_i, x_j) (negative: dark red). ",
"S_i (MC) = first-order Sobol' index: share of Var(Y) due to x_i alone, with all nonlinear effects (violet); interactions = 1 - sum S_i."),
if (is.null(SOB)) paste0("Sobol' indices: not computed (", if (!indep) "they require independent input quantities)." else
if (!isTRUE(sobol)) "sobol = FALSE)." else "computation failed)."),
if (length(zeroC)) paste0("Inputs with c_i = 0 but non-zero second derivatives (", paste(zeroC, collapse = ", "),
") are not represented in the first-order index.")), cex = 0.54)
## histogram (left) with Q-Q plot and convergence trace (right)
need(0.225); yT <- Y - 0.005; yB <- max(0.07, Y - 0.225)
yM <- (yT + yB)/2
F1 <- ndc(c(0.05, 0.6), c(yB, yT)); F2 <- ndc(c(0.64, 0.95), c(yM + 0.005, yT)); F3 <- ndc(c(0.64, 0.95), c(yB, yM - 0.005))
subfig(F1, c(2.4, 1, 1.6, 1))
q <- quantile(y, c(5E-4, 1 - 5E-4)); h <- hist(y[y >= q[1] & y <= q[2]], breaks = 70, plot = FALSE)
band <- if (!is.null(GM) && !is.null(GU)) GM + c(-1, 1) * num1(O$k) * GU else NULL
plot(h, freq = FALSE, col = NA, border = NA, main = paste0("Monte Carlo output (", interval, " interval)"), xlab = "",
ylab = "", yaxt = "n", xlim = range(c(h$breaks, ciMC, ciG, band)), cex.main = 0.8 * fs, font.main = 3, cex.axis = 0.62 * fs)
if (!is.null(band)) rect(band[1], 0, band[2], max(h$density) * 1.05, col = "#dcebcf", border = NA)
plot(h, freq = FALSE, add = TRUE, col = "#c9b8de", border = "#b3a0cc")
abline(v = ciMC, col = VIO, lwd = 2); abline(v = ciG, col = GRN, lwd = 2, lty = 2)
if (!is.null(FIT)) { xs <- seq(min(h$breaks), max(h$breaks), length.out = 500)
lines(xs, do.call(FIT$dens[[1]], c(list(xs), as.list(FIT$bestpar))), col = RED, lwd = 2.5) }
legend("topright", bty = "n", cex = 0.58 * fs, legend = c("MC interval", "GUF interval", if (!is.null(band)) "gauge y +/- k u",
if (!is.null(FIT)) paste("fit:", FIT$stat$Distribution[1])), col = c(VIO, GRN, if (!is.null(band)) "#dcebcf", if (!is.null(FIT)) RED),
lty = c(1, 2, if (!is.null(band)) 1, if (!is.null(FIT)) 1), lwd = c(2, 2, if (!is.null(band)) 6, if (!is.null(FIT)) 2.5))
## Q-Q plot: quantiles of the MC output against those of the GUF distribution (t with nu_eff, or Gaussian)
subfig(F2, c(2, 2.4, 1.3, 0.4))
pp <- c(seq(0.001, 0.01, length.out = 10), seq(0.02, 0.98, by = 0.02), seq(0.99, 0.999, length.out = 10))
nuQ <- if (is.finite(num1(O$ws.df)) && O$ws.df < 1000) O$ws.df else Inf
qG <- unname(O$prop["mean.1"]) + u1 * (if (is.finite(nuQ)) qt(pp, nuQ) else qnorm(pp)); qM <- quantile(y, pp, names = FALSE)
plot(qG, qM, pch = 16, cex = 0.35, col = VIO, xlab = "", ylab = "", main = paste0("Q-Q: MC vs. GUF (", if (is.finite(nuQ)) paste0("t, nu = ", nuQ) else "Gaussian", ")"),
cex.main = 0.62 * fs, font.main = 3, axes = FALSE)
for (sd_ in 1:2) { at <- pretty(if (sd_ == 1) qG else qM, 4); axis(sd_, at = at, labels = formatC(at, format = "g", digits = 4), cex.axis = 0.5 * fs,
mgp = c(1.2, 0.3, 0), tcl = -0.2) }; box()
abline(0, 1, col = GRN, lty = 2)
a2 <- O$alpha/2; points(ciG, quantile(y, c(a2, 1 - a2)), pch = 4, col = RED, cex = 0.9, lwd = 1.5)
mtext("GUF quantile", 1, 1, cex = 0.45 * fs); mtext("MC quantile", 2, 1.4, cex = 0.45 * fs)
## convergence trace: running mean and u relative to the final values, in units of u
subfig(F3, c(2, 2.4, 1.3, 0.4))
yo <- O$resSIM[is.finite(O$resSIM)]; nT <- length(yo); ii <- unique(round(exp(seq(log(min(100, nT)), log(nT), length.out = 60))))
cs <- cumsum(yo); cs2 <- cumsum(yo^2); rm_ <- cs[ii]/ii; rs <- sqrt(pmax(0, (cs2[ii] - ii * rm_^2)/(ii - 1)))
uF <- sd(yo); mF <- mean(yo)
dm <- (rm_ - mF)/uF; du <- rs/uF - 1
plot(ii, dm, type = "l", log = "x", col = BLU, lwd = 1.5, ylim = range(c(dm, du, -0.05, 0.05), finite = TRUE), xlab = "", ylab = "",
main = "Convergence of mean and u", cex.main = 0.62 * fs, font.main = 3, cex.axis = 0.5 * fs, mgp = c(1.2, 0.3, 0), tcl = -0.2)
lines(ii, du, col = ORA, lwd = 1.5); abline(h = 0, col = "grey50", lty = 3)
mtext("number of trials M", 1, 1, cex = 0.45 * fs); mtext("relative to final, in u", 2, 1.4, cex = 0.45 * fs)
legend("topright", c("(mean_M - mean)/u", "u_M/u - 1"), col = c(BLU, ORA), lwd = 1.5, bty = "n", cex = 0.45 * fs)
page(); Y <- yB - 0.02
## conformity assessment (optional)
if (!is.null(CA)) {
head(paste0("Conformity assessment (JCGM 106:2012; decision rules ISO 14253-1), ", CA$method))
para(paste0("Tolerance interval [", format(CA$lower), ", ", format(CA$upper), "]; estimate y = ", g4(CA$estimate, 7), "; conformance probability p_c = ",
formatC(CA$pc, format = "f", digits = 4), " (1 - p_c = ", g4(CA$pnc, 3), "); expanded uncertainty U- = ", g4(CA$U[1], 3),
", U+ = ", g4(CA$U[2], 3), " (p = ", CA$p, ")."), cex = 0.6, font = 2)
RU <- CA$rules
if (!is.null(RU) && nrow(RU)) {
lim <- c(CA$lower, CA$upper, RU$RL, RU$RU, RU$AL, RU$AU, CA$estimate); lim <- lim[is.finite(lim)]
xr <- range(lim) + c(-0.1, 0.1) * diff(range(lim)); X <- function(v) 0.36 + 0.44 * (pmin(pmax(v, xr[1]), xr[2]) - xr[1])/diff(xr)
rh <- 0.017 * fs; need(rh * (nrow(RU) + 1.5))
for (i in seq_len(nrow(RU))) {
r <- RU[i, ]; yc <- Y; hb <- rh * 0.3
text(0.05, yc, fitTxt(r$description, 0.3, 0.52), adj = c(0, 0.5), cex = 0.52 * fs)
rect(X(xr[1]), yc - hb, X(xr[2]), yc + hb, col = "#C8102E", border = NA) ## reject
rect(X(r$RL), yc - hb, X(r$RU), yc + hb, col = "grey75", border = NA) ## uncertainty zone
if (r$AL < r$AU) rect(X(r$AL), yc - hb, X(r$AU), yc + hb, col = "#2E9D48", border = NA) ## accept
segments(X(c(CA$lower, CA$upper)), yc - hb * 1.4, X(c(CA$lower, CA$upper)), yc + hb * 1.4, lwd = 1.5)
points(X(CA$estimate), yc + hb * 1.9, pch = 25, bg = "black", cex = 0.6)
dcol <- switch(r$decision, accept = "#2E9D48", reject = "#C8102E", "grey35")
text(0.82, yc, toupper(r$decision), adj = c(0, 0.5), cex = 0.6 * fs, font = 2, col = dcol)
Y <- Y - rh }
tk <- pretty(xr); tk <- tk[tk >= xr[1] & tk <= xr[2]]
segments(X(xr[1]), Y + rh * 0.4, X(xr[2]), Y + rh * 0.4, col = "grey40"); segments(X(tk), Y + rh * 0.4, X(tk), Y + rh * 0.15, col = "grey40")
text(X(tk), Y - rh * 0.05, format(tk), cex = 0.5 * fs, col = "grey30")
text(X(c(CA$lower, CA$upper)), Y - rh * 0.55, c("LSL", "USL"), cex = 0.5 * fs, font = 2)
Y <- Y - rh * 1.2
para("Green: acceptance zone; grey: uncertainty zone (no decision); red: rejection; triangle: estimate; vertical lines: tolerance limits.", cex = 0.52, col = "#3a3a3a")
}
}
need(0.04 + (N + 2) * 0.0135); head("Sensitivity coefficients and Hessian matrix")
gS <- symVec(O$gradient); if (length(gS) != N) gS <- rep("", N)
lhD <- if (N > 12) 0.0125 else 0.0135
DER <- data.frame(Quantity = nm, "c_i = df/dx_i, symbolic" = fitTxt(gS, if (N <= 6) 0.33 else 0.5), "c_i, evaluated" = g4(gr),
check.names = FALSE, stringsAsFactors = FALSE)
if (is.matrix(H) && N <= 6) { ## evaluated Hessian as extra columns of the gradient table
cw <- min(0.08, 0.34/N); HE <- fmtM(matrix(H, N, dimnames = list(nm, nm)))
DH <- cbind(DER, as.data.frame(HE, stringsAsFactors = FALSE)); names(DH)[-(1:3)] <- paste0("H[ , ", nm, "]")
tab(DH, c(0.1, 0.35, 0.11, rep(cw, N)), align = c("l", "l", rep("r", N + 1)), lh = lhD, cex = 0.56)
hS <- symVec(O$hessian)
if (length(hS) == N^2 && any(H != 0) && N > 3)
para("The symbolic Hessian matrix (more than 3 input quantities) is printed by summary().", cex = 0.54)
if (length(hS) == N^2 && any(H != 0) && N <= 3) { ## symbolic Hessian (lower triangle; H is symmetric)
hS <- matrix(hS, N, dimnames = list(nm, nm)); hS[upper.tri(hS)] <- ""
cw2 <- min(0.28, 0.84/N); hS[] <- fitTxt(hS, cw2 - 0.01)
mrow(list(hS), "Hessian H = d^2 f/(dx_i dx_j), symbolic (lower triangle)", 0.05, cw2, (N + 3) * lhD, align = rep("l", N + 1))
}
if (all(H == 0)) para("H = 0: the model is linear in all input quantities at the estimates.", cex = 0.54)
} else {
tab(DER, c(0.14, 0.52, 0.16), align = c("l", "l", "r"), lh = lhD, cex = 0.56)
if (is.matrix(H)) para(paste0("Hessian matrix (", N, " x ", N, ") not shown; non-zero elements: ", sum(H != 0), "."), cex = 0.54)
}
## covariance check of the Monte Carlo input samples (as in summary.propagate)
need(0.04 + (min(N, 8) + 3) * 0.0135); head("Covariance check of the Monte Carlo input samples")
TGT <- if (!is.null(O$covMat.target)) O$covMat.target else S; XS <- O$datSIM
noiseRef <- function(FUN, target) { if (is.null(XS) || nrow(XS) < 20 || !all(is.finite(target))) return(NA)
n <- nrow(XS) %/% 2; M1 <- FUN(XS[1:n, , drop = FALSE], use = "pairwise.complete.obs")
M2 <- FUN(XS[(n + 1):(2 * n), , drop = FALSE], use = "pairwise.complete.obs"); norm(M1 - M2, "F")/(2 * norm(target, "F")) }
judge <- function(eps, noise) { if (!is.finite(eps)) return(list("not applicable (expected covariance does not exist)", "grey30"))
if (!is.finite(noise) || noise == 0) return(list(format(signif(eps, 3)), "black")); rr <- eps/noise
list(paste0(signif(eps, 3), " (", signif(rr, 3), " x MC noise level => ", if (rr > 3) "CHECK" else "OK", ")"), if (rr > 3) RED else GRN) }
CM <- list(S, TGT, if (is.matrix(O$covMat.MC)) O$covMat.MC else S * NA)
CM <- lapply(CM, function(M) { dimnames(M) <- list(nm, nm); fmtM(M) })
caps <- c("Input covariance (GUF)", "MC input covariance, expected", "MC design covariance, observed")
if (N <= 4) mrow(CM, caps, c(0.05, 0.36, 0.67), min(0.07, 0.23/N), (N + 3) * 0.0135)
else if (N <= 8) for (j in 1:3) mrow(CM[j], caps[j], 0.05, 0.09, (N + 3) * 0.0135)
JC <- judge(num1(O$frobCOV), noiseRef(cov, TGT))
para(paste("Relative Frobenius norm, expected vs. observed covariance:", JC[[1]]), cex = 0.58, col = JC[[2]], font = 2)
R0 <- O$corMat
if (is.matrix(R0) && ncol(R0) > 1 && any(abs(R0[upper.tri(R0)]) > 0)) {
JR <- judge(num1(O$frobCOR), noiseRef(cor, R0))
para(paste("Relative Frobenius norm, input vs. observed correlation:", JR[[1]]), cex = 0.58, col = JR[[2]], font = 2) }
if (heavy) para("Note: DOF <= 4 => the t-margins have infinite kurtosis; covariance estimates converge slowly and the check is unreliable.", cex = 0.54)
if (O$type %in% c("stat", "raw")) para(paste0("MC input models for Type A quantities: S1 = scaled and shifted t-distributions, scale u(x), nu = n - 1 ",
"(JCGM 101:2008, 6.4.9), correlated inputs joined by a Gaussian copula (Possolo 2010); mvt = multivariate t, nu = n - N (JCGM 102:2011, 5.3.2)."),
cex = 0.54, col = "#3a3a3a")
## reproducibility record
head("Reproducibility (ISO/IEC 17025:2017, 7.5, 7.8)")
for (k in names(REPRO)) {
need(0.015); text(0.05, Y, paste0(k, ":"), adj = c(0, 0.5), cex = 0.56 * fs, font = 2)
para(REPRO[[k]], cex = 0.54, family = if (k == "Call") "mono" else "", x = 0.16, w = 0.79); Y <- Y + 0.006 }
OUT[[ia]] <- list(results = RES, validation = VAL, accuracy = ACC, budget = BUD[, -ncol(BUD)], sobol = SOB, nonlinearity = NL,
derivatives = DER, shape = SHP, checks = CHK, conformity = CA, fit = FIT, reproducibility = REPRO)
}
invisible(c(list(file = file), if (nA == 1) OUT[[1]] else list(analyses = OUT)))
}
Try the propagate package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.