R/conformity.R

In propagate: Propagation of Uncertainty

conformity <- function(object, lower = -Inf, upper = Inf, method = c("MC", "GUF"), prob = NULL,
                       um = NULL, accept = NULL, risk = NULL, support = c(-Inf, Inf),
                       rule = NULL, p = 0.95, interval = c("symmetric", "shortest"))
{
  ## Conformity assessment with measurement uncertainty (JCGM 106:2012)
  ## 1) 'um' = NULL: specific conformance probability of one measured item, from a 'propagate' object
  ##    (Monte Carlo sample or GUF), a vector c(y, u[, nu]) or a sample; optional guarded acceptance
  ##    limits for a required conformance probability 'prob' (JCGM 106, 7 and 8).
  ## 2) 'um' given: global consumer's and producer's risks for a production process, whose distribution
  ##    is given by 'object' as a density function or as a sample, inspected by a measuring system with
  ##    standard uncertainty 'um' (normal); optional acceptance limits for a target consumer's risk
  ##    'risk' (JCGM 106, 9.5).
  ## 3) v1.3-0: 'rule': decision rules with the expanded uncertainty U (coverage probability 'p'),
  ##    ISO 14253-1, ASME B89.7.3.1, VDA 5: acceptance, rejection and uncertainty zones. For a single
  ##    item the decision; for a production process the accepted/rejected/undecided fractions and the
  ##    consumer's and producer's risks of each rule.
  method <- match.arg(method); interval <- match.arg(interval)
  if (lower >= upper) stop("conformity: 'lower' must be smaller than 'upper'.")
  if (!is.null(rule)) {
    rule <- .ruleNames(rule)
    if (p <= 0 || p >= 1) stop("conformity: 'p' must be in (0, 1).")
  }
  inC <- function(x) x >= lower & x <= upper

  ## ---------------------------------- 1) specific conformance probability
  if (is.null(um)) {
    if (inherits(object, "propagate")) {
      y <- unname(object$prop["mean.1"]); u <- unname(object$prop["u.1"])
      nu <- if (is.numeric(object$ws.df) && length(object$ws.df) && is.finite(object$ws.df)) object$ws.df else Inf
      S <- object$resSIM[is.finite(object$resSIM)]
    } else if (is.numeric(object) && length(object) %in% 2:3) {
      y <- object[1]; u <- object[2]; nu <- if (length(object) == 3) object[3] else Inf; S <- NULL
      method <- "GUF"
    } else if (is.numeric(object) && length(object) > 3) {
      S <- object[is.finite(object)]; y <- mean(S); u <- sd(S); nu <- Inf; method <- "MC"
    } else stop("conformity: 'object' must be a 'propagate' object, c(y, u), c(y, u, nu) or a sample.")
    if (method == "MC" && is.null(S)) stop("conformity: method = \"MC\" needs a Monte Carlo sample.")
    if (method == "MC") { y <- mean(S); u <- sd(S) }                  ## report the Monte Carlo estimate and u
    if (!is.finite(u) || u <= 0) stop("conformity: the standard uncertainty must be positive.")
    tdist <- is.finite(nu) && nu < 1000
    P <- function(q) if (tdist) pt(q, nu) else pnorm(q)               ## standardized GUF distribution
    Q <- function(p) if (tdist) qt(p, nu) else qnorm(p)
    pc <- if (method == "MC") mean(inC(S)) else P((upper - y)/u) - P((lower - y)/u)
   
    ## acceptance limits: a measured value at the limit has conformance probability 'prob' (one-sided,
    ## the distribution of the measurand shifted with the measured value)
    AL <- AU <- NA
    if (!is.null(prob)) {
      if (prob <= 0 || prob >= 1) stop("conformity: 'prob' must be in (0, 1).")
      D <- if (method == "MC") S - mean(S) else NULL
      q <- function(p) if (method == "MC") unname(quantile(D, p)) else u * Q(p)
      if (is.finite(upper)) AU <- upper - q(prob)
      if (is.finite(lower)) AL <- lower - q(1 - prob)
    }
    
    ## v1.3-0: coverage interval for 'p' (asymmetric for MC) and density of the measurand (for plot)
    CI <- if (method == "MC") .covIntJCGM(S, p, shortest = interval == "shortest") else y + c(-1, 1) * Q((1 + p)/2) * u
    dens <- if (method == "MC") density(S, n = 512) else {
      xs <- seq(y - 5 * u, y + 5 * u, length.out = 512)
      list(x = xs, y = (if (tdist) dt((xs - y)/u, nu) else dnorm((xs - y)/u))/u) }
    RULE <- NULL
    if (!is.null(rule)) {
      RULE <- .zoneLimits(lower, upper, y - CI[1], CI[2] - y, rule)
      RULE$decision <- with(RULE, ifelse(AL <= y & y <= AU, "accept", ifelse(y < RL | y > RU, "reject", "no decision")))
    }
    OUT <- list(type = "specific", method = method, estimate = y, u = u, nu = nu, lower = lower, upper = upper,
                pc = pc, pnc = 1 - pc, prob = prob, accept = c(AL = AL, AU = AU),
                z = c(lower = (y - lower)/u, upper = (upper - y)/u),
                p = p, U = c("U-" = y - CI[1], "U+" = CI[2] - y), interval = CI, rules = RULE, density = dens)
    class(OUT) <- "conformity"
    return(OUT)
  }

  ## ------------------ 2) global risks of a production process
  if (!is.numeric(um) || um <= 0) stop("conformity: 'um' must be a positive standard uncertainty.")
  if (!is.null(rule) && (!is.null(accept) || !is.null(risk)))
    stop("conformity: use either 'rule' or 'accept'/'risk' to define the acceptance limits.")
  if (inherits(object, "propagate")) object <- object$resSIM
  isFun <- is.function(object)
  if (!isFun && !(is.numeric(object) && length(object) > 3))
    stop("conformity: with 'um', 'object' must be a density function or a sample of the process.")
  if (!isFun) object <- object[is.finite(object)]
  pA <- function(x, A) pnorm((A[2] - x)/um) - pnorm((A[1] - x)/um)   ## P(measured value in A | eta)
  
  ## integral of f(eta) * g0(eta) over [a, b], within the support of the density
  int <- function(f, a, b) {
    a <- max(a, support[1]); b <- min(b, support[2]); if (a >= b) return(0)
    integrate(function(x) f(x) * object(x), a, b, subdivisions = 1000L, rel.tol = 1E-10)$value
  }
  E <- function(f, region) {                                       ## expectation over the process
    if (!isFun) return(mean(f(object) * switch(region, all = 1, C = inC(object), notC = !inC(object))))
    switch(region, all = int(f, -Inf, lower) + int(f, lower, upper) + int(f, upper, Inf),
           C = int(f, lower, upper), notC = int(f, -Inf, lower) + int(f, upper, Inf))
  }
  risks <- function(A) c(RC = E(function(x) pA(x, A), "notC"), RP = E(function(x) 1 - pA(x, A), "C"))
  
  ## acceptance limits: given, or solved for a target consumer's risk (symmetric guard band w on the
  ## finite tolerance limits, w > 0: guarded acceptance), or simple acceptance (A = T)
  w <- NA
  if (is.null(accept)) {
    if (!is.null(risk)) {
      Aw <- function(w) c(if (is.finite(lower)) lower + w else -Inf, if (is.finite(upper)) upper - w else Inf)
      f <- function(w) risks(Aw(w))["RC"] - risk
      w <- uniroot(f, c(-6 * um, if (all(is.finite(c(lower, upper)))) (upper - lower)/2 - 1E-9 else 6 * um),
                   tol = 1E-10 * um)$root
      accept <- Aw(w)
    } else accept <- c(lower, upper)
  }
  R <- risks(accept)
  pc0 <- E(function(x) 1, "C"); pAcc <- E(function(x) pA(x, accept), "all")
  
  ## process location and spread, capability indices (normal-theory definitions)
  I3 <- function(f) int(f, -Inf, lower) + int(f, lower, upper) + int(f, upper, Inf)   ## split at the limits
  m <- if (isFun) I3(function(x) x) else mean(object)
  s <- if (isFun) sqrt(I3(function(x) (x - m)^2)) else sd(object)
  Cp <- if (all(is.finite(c(lower, upper)))) (upper - lower)/(6 * s) else NA
  Cpk <- min(c((upper - m)/(3 * s), (m - lower)/(3 * s)), na.rm = TRUE)
  
  ## v1.3-0: decision rules with U = k um: accepted, rejected and undecided fractions, and the risks
  ## R_C = P(accepted, not conforming), R_P = P(rejected, conforming) of each rule
  RULE <- NULL; U <- NA
  if (!is.null(rule)) {
    U <- qnorm((1 + p)/2) * um
    RULE <- .zoneLimits(lower, upper, U, U, rule)
    stat <- t(sapply(seq_len(nrow(RULE)), function(i) {
      A <- c(RULE$AL[i], RULE$AU[i]); Rj <- c(RULE$RL[i], RULE$RU[i])
      acc <- if (RULE$AL[i] > RULE$AU[i]) 0 else E(function(x) pA(x, A), "all")
      rej <- E(function(x) 1 - pA(x, Rj), "all")
      und <- 1 - acc - rej; if (abs(und) < 1E-10) und <- 0           ## numerical integration noise
      c(accepted = acc, rejected = rej, undecided = max(0, und),
        RC = if (RULE$AL[i] > RULE$AU[i]) 0 else E(function(x) pA(x, A), "notC"),
        RP = E(function(x) 1 - pA(x, Rj), "C"))
    }))
    RULE <- cbind(RULE, as.data.frame(stat))
  }
  
  ## density of the process (for plot)
  dens <- if (isFun) {
    xr <- c(max(support[1], m - 5 * s), min(support[2], m + 5 * s))
    xs <- seq(xr[1], xr[2], length.out = 512); list(x = xs, y = object(xs))
  } else density(object, n = 512)
  OUT <- list(type = "global", method = if (isFun) "integration" else "MC", lower = lower, upper = upper, um = um,
              accept = c(AL = accept[1], AU = accept[2]), w = w, r = w/(2 * um),
              pc = pc0, pnc = 1 - pc0, RC = unname(R["RC"]), RP = unname(R["RP"]), accepted = pAcc,
              RC.accepted = unname(R["RC"])/pAcc, process = c(mean = m, sd = s, Cp = Cp, Cpk = Cpk, DPMO = 1E6 * (1 - pc0)),
              p = p, U = U, rules = RULE, density = dens)
  class(OUT) <- "conformity"
  OUT
}

## v1.3-0: decision rules (acceptance criterion . rejection criterion), ISO 14253-1, ASME B89.7.3.1, VDA 5
.RULES <- c(stringent = "stringent.stringent", simple = "simple.simple", stringent.relaxed = "stringent.relaxed",
            relaxed.stringent = "relaxed.stringent", stringent.simple = "stringent.simple",
            simple.stringent = "simple.stringent")
.RULELABEL <- c(stringent = "Stringent acceptance - stringent rejection (ISO 14253-1)",
                simple = "Simple acceptance - simple rejection",
                stringent.relaxed = "Stringent acceptance - relaxed rejection",
                relaxed.stringent = "Relaxed acceptance - stringent rejection",
                stringent.simple = "Stringent acceptance - simple rejection",
                simple.stringent = "Simple acceptance - stringent rejection")
.ruleNames <- function(rule) {
  if (identical(rule, "all")) return(names(.RULES))
  bad <- setdiff(rule, names(.RULES))
  if (length(bad)) stop("conformity: unknown rule(s) '", paste(bad, collapse = "', '"), "'; use ",
                        paste(names(.RULES), collapse = ", "), " or \"all\".")
  rule
}

## acceptance interval [AL, AU] and rejection limits RL, RU (reject if y < RL or y > RU) of each rule,
## with U- = y - y_low and U+ = y_high - y of the coverage interval
.zoneLimits <- function(lower, upper, Um, Up, rule) {
  Z <- t(sapply(rule, function(r) {
    a <- strsplit(.RULES[[r]], ".", fixed = TRUE)[[1]]
    acc <- switch(a[1], stringent = c(lower + Um, upper - Up), simple = c(lower, upper), relaxed = c(lower - Up, upper + Um))
    rej <- switch(a[2], stringent = c(lower - Up, upper + Um), simple = c(lower, upper), relaxed = c(lower + Um, upper - Up))
    c(acc, rej)
  }))
  data.frame(rule = rule, description = unname(.RULELABEL[rule]), AL = Z[, 1], AU = Z[, 2], RL = Z[, 3], RU = Z[, 4],
             stringsAsFactors = FALSE, row.names = NULL)
}

print.conformity <- function(x, digits = 4, ...)
{
  f <- function(v) format(signif(v, digits))
  g <- function(v) format(signif(v, 10))                            ## limits and estimates: full precision
  lim <- paste0("[", g(x$lower), ", ", g(x$upper), "]")
  if (x$type == "specific") {
    cat("Conformity assessment of a measured item (JCGM 106:2012, 7), ", x$method, "\n", sep = "")
    cat("Estimate y = ", g(x$estimate), ", u = ", f(x$u), if (x$method == "GUF") paste0(", nu = ", f(x$nu)),
        ", tolerance interval ", lim, "\n", sep = "")
    cat("Conformance probability pc = ", f(x$pc), ", non-conformance 1 - pc = ", f(x$pnc), "\n", sep = "")
    if (!is.null(x$prob)) cat("Acceptance limits for pc >= ", x$prob, " (guarded acceptance, JCGM 106, 8.3): AL = ",
                             g(x$accept["AL"]), ", AU = ", g(x$accept["AU"]), "\n", sep = "")
    if (!is.null(x$rules)) {
      cat("Decision rules with the expanded uncertainty U- = ", f(x$U[1]), ", U+ = ", f(x$U[2]), " (p = ", x$p, "):\n", sep = "")
      .printRules(x$rules, x$U, digits, cols = c("decision"))
    }
  } else {
    cat("Global risks of a production process (JCGM 106:2012, 9.5), ", x$method, "\n", sep = "")
    cat("Tolerance interval ", lim, ", measurement uncertainty um = ", f(x$um), ", acceptance interval [",
        g(x$accept[1]), ", ", g(x$accept[2]), "]", if (is.finite(x$w)) paste0(" (guard band w = ", f(x$w), ", r = ", f(x$r), ")"), "\n", sep = "")
    cat("Conforming fraction of the process pc = ", f(x$pc), " (DPMO = ", f(x$process["DPMO"]), ")\n", sep = "")
    cat("Consumer's risk RC = ", f(x$RC), ", producer's risk RP = ", f(x$RP), ", accepted fraction = ", f(x$accepted), "\n", sep = "")
    cat("Non-conforming fraction among accepted items = ", f(x$RC.accepted), "\n", sep = "")
    cat("Process: mean = ", g(x$process["mean"]), ", sd = ", f(x$process["sd"]), ", Cp = ", f(x$process["Cp"]),
        ", Cpk = ", f(x$process["Cpk"]), "\n", sep = "")
    if (!is.null(x$rules)) {
      cat("Decision rules with the expanded uncertainty U = ", f(x$U), " (k um, p = ", x$p, "):\n", sep = "")
      .printRules(x$rules, x$U, digits, cols = c("accepted", "rejected", "undecided", "RC", "RP"))
    }
  }
  invisible(x)
}

## table of decision rules (internal): rule as acceptance/rejection criterion
.printRules <- function(Z, U, digits, cols) {
  Ur <- min(U[U > 0]); d <- if (is.finite(Ur)) max(0, 1 - floor(log10(Ur))) + 1 else 4
  fz <- function(v) ifelse(is.finite(v), formatC(round(v, d), format = "f", digits = d), ifelse(v > 0, "Inf", "-Inf"))
  lab <- gsub(".", "/", .RULES[Z$rule], fixed = TRUE); lab[Z$rule == "stringent"] <- "stringent/stringent (ISO 14253-1)"
  T <- data.frame("acceptance/rejection" = lab, accept = paste0("[", fz(Z$AL), ", ", fz(Z$AU), "]"),
                  "reject below | above" = paste0(fz(Z$RL), " | ", fz(Z$RU)), stringsAsFactors = FALSE, check.names = FALSE)
  for (cc in cols) T[[cc]] <- if (is.numeric(Z[[cc]])) trimws(formatC(ifelse(abs(Z[[cc]]) < 1E-10, 0, Z[[cc]]), digits = digits, format = "g")) else toupper(Z[[cc]])
  W <- pmax(nchar(names(T)), apply(T, 2, function(z) max(nchar(z))))
  line <- function(v) paste(mapply(function(s, w) formatC(s, width = -w), v, W), collapse = "  ")
  cat(line(names(T)), "\n", sep = "")
  for (i in seq_len(nrow(T))) cat(line(unlist(T[i, ])), "\n", sep = "")
}

plot.conformity <- function(x, density = TRUE, col = c(accept = "#2E9D48", reject = "#C8102E", zone = "grey70"), ...)
{
  ## top: density of the measurand (single item) or of the process (global), shaded inside/outside the
  ## tolerance interval; bottom: acceptance, rejection and uncertainty zones of the decision rules
  Z <- x$rules; n <- if (is.null(Z)) 0 else nrow(Z); L <- x$lower; U <- x$upper
  fin <- c(L, U)[is.finite(c(L, U))]
  D <- x$density
  if (x$type == "specific") { Um <- x$U[1]; Up <- x$U[2]; y <- x$estimate } else { Um <- Up <- x$U; y <- NA }
  lims <- c(fin, if (n) c(Z$AL, Z$AU, Z$RL, Z$RU), if (x$type == "specific") x$interval)
  lims <- lims[is.finite(lims)]
  span <- diff(range(c(lims, D$x[D$y > 1E-3 * max(D$y)])))
  xl <- range(c(lims, D$x[D$y > 1E-3 * max(D$y)])) + c(-0.08, 0.08) * span
  op <- par(no.readonly = TRUE); on.exit(par(op))
  density <- density || n == 0
  if (density && n) layout(matrix(1:2), heights = c(1.1, 0.35 + 0.55 * n))
  if (density) {
    par(mar = c(if (n) 1 else 4, 1, 3, 1))
    keep <- D$x >= xl[1] & D$x <= xl[2]
    plot(D$x[keep], D$y[keep], type = "n", xlim = xl, ylim = c(0, 1.08 * max(D$y)), axes = FALSE, xlab = "", ylab = "")
    if (!n) axis(1, cex.axis = 0.8)
    inside <- D$x >= L & D$x <= U
    for (part in list(list(keep & inside, col["accept"]), list(keep & D$x < L, col["reject"]), list(keep & D$x > U, col["reject"]))) {
      s <- part[[1]]
      if (sum(s) > 1) polygon(c(D$x[s][1], D$x[s], rev(D$x[s])[1]), c(0, D$y[s], 0), col = adjustcolor(part[[2]], 0.35), border = NA)
    }
    lines(D$x[keep], D$y[keep], lwd = 2)
    abline(v = fin, lty = 2); abline(h = 0, col = "grey50")
    if (x$type == "specific") {
      title(main = sprintf("Measurand: P(LSL <= Y <= USL) = %.4g  (%s)", x$pc, x$method), cex.main = 1)
    } else {
      title(main = sprintf("Process: conforming fraction %.4g, um = %.3g", x$pc, x$um), cex.main = 1)
      if (any(is.finite(x$accept))) abline(v = x$accept[is.finite(x$accept)], col = "grey40", lty = 3)
    }
  }
  if (!n) return(invisible(x))
  par(mar = c(3, 1, if (density) 1 else 3, 1))
  plot(NA, xlim = xl, ylim = c(0.4, n + 0.6), axes = FALSE, xlab = "", ylab = "")
  axis(1, cex.axis = 0.8)
  h <- 0.18
  for (i in seq_len(n)) {
    yy <- n + 1 - i
    brk <- sort(unique(pmin(pmax(c(xl, Z$AL[i], Z$AU[i], Z$RL[i], Z$RU[i]), xl[1]), xl[2])))
    for (j in seq_len(length(brk) - 1)) {
      m <- mean(brk[j:(j + 1)])
      cz <- if (Z$AL[i] <= m && m <= Z$AU[i]) col["accept"] else if (m < Z$RL[i] || m > Z$RU[i]) col["reject"] else col["zone"]
      rect(brk[j], yy - h, brk[j + 1], yy + h, col = cz, border = NA)
    }
    segments(fin, yy - h - 0.05, fin, yy + h + 0.05, lwd = 2)
    text(xl[1], yy + h + 0.22, Z$description[i], adj = c(0, 0), cex = 0.75)
    if (x$type == "specific") {
      arrows(y - Um, yy + h + 0.12, y + Up, yy + h + 0.12, angle = 90, code = 3, length = 0.04, lwd = 1.3)
      points(y, yy + h + 0.12, pch = 25, bg = "black", cex = 0.9)
      dc <- c(accept = col[["accept"]], reject = col[["reject"]], "no decision" = "grey30")[Z$decision[i]]
      text(xl[2], yy + h + 0.22, toupper(Z$decision[i]), adj = c(1, 0), cex = 0.8, font = 2, col = dc)
    } else {
      text(xl[2], yy + h + 0.22, sprintf("RC = %.3g, RP = %.3g, undecided = %.3g", Z$RC[i], Z$RP[i], Z$undecided[i]),
           adj = c(1, 0), cex = 0.75)
    }
  }
  mtext(c(if (is.finite(L)) "LSL", if (is.finite(U)) "USL"), side = 1, at = fin, line = 1.8, cex = 0.8, font = 2)
  invisible(x)
}