R/propagateMV.R

In propagate: Propagation of Uncertainty

## propagateMV() -- propagate 1.3-0
## Uncertainty propagation for several output quantities (JCGM 102:2011), built on propagate():
##  - one propagate() run on a "carrier" model generates the joint Monte Carlo input sample of all
##    input quantities, with any input model of propagate (mc.margins, cov, type);
##  - all output models are evaluated on this same sample, giving the joint distribution of the outputs;
##  - one propagate() run per output gives its sensitivity coefficients; the GUF covariance matrix of
##    the outputs is U_y = C_x U_x C_x' (JCGM 102:2011, 6.2.1.3, Eq. 3).
## Coverage regions: hyper-ellipsoidal (k_p) and hyper-rectangular (k_q), JCGM 102:2011, 6.5 and 7.7.

propagateMV <- function(expr, data, type = c("stat", "raw", "sim"), cov = c("diag", "full"),
                        nsim = 1000000, alpha = 0.05, mc.margins = c("S1", "mvt"), ...)
{
  type <- match.arg(type)
  if (identical(mc.margins, "GUF")) stop("propagateMV: mc.margins = \"GUF\" was removed in 1.3-0; use \"S1\" or \"mvt\" (see ?propagate).")
  mc.margins <- match.arg(mc.margins)
  if (is.list(expr) && !is.expression(expr)) expr <- as.expression(lapply(expr, function(e) if (is.expression(e)) e[[1]] else e))
  if (!is.expression(expr) || length(expr) < 2) stop("propagateMV: 'expr' must be an expression with at least two elements (output quantities)!")
  m <- length(expr)
  NAMES <- names(expr); if (is.null(NAMES) || any(NAMES == "")) NAMES <- paste0("Y", seq_len(m))
  EXPR <- lapply(seq_len(m), function(l) as.expression(expr[[l]]))
  VARS <- unique(unlist(lapply(EXPR, all.vars)))
  if (!all(VARS %in% colnames(data))) stop("propagateMV: variable(s) ", paste(setdiff(VARS, colnames(data)), collapse = ", "), " not in 'data'!")
  p <- 1 - alpha
  quiet <- function(...) { out <- NULL; invisible(capture.output(out <- suppressMessages(propagate(...)))); out }
  
  ## 1) joint Monte Carlo input sample: carrier model = first output + 0 * (all other inputs)
  rest <- setdiff(VARS, all.vars(EXPR[[1]]))
  CARRIER <- if (length(rest)) as.expression(call("+", EXPR[[1]][[1]], str2lang(paste0("0 * (", paste(rest, collapse = " + "), ")")))) else EXPR[[1]]
  DAT <- data[, intersect(colnames(data), VARS), drop = FALSE]
  if (!is.null(attr(data, "budget"))) attr(DAT, "budget") <- attr(data, "budget")
  CAR <- quiet(CARRIER, DAT, type = type, cov = cov, nsim = nsim, alpha = alpha, mc.margins = mc.margins, check = FALSE, ...)
  Ux <- CAR$covMat
  
  ## 2) Monte Carlo: all outputs on the same input sample
  Y <- sapply(EXPR, function(e) eval(e[[1]], as.data.frame(CAR$datSIM)))
  Y <- matrix(Y, ncol = m, dimnames = list(NULL, NAMES))
  ok <- rowSums(!is.finite(Y)) == 0
  if (any(!ok)) { message("propagateMV: ", sum(!ok), " Monte Carlo trials with non-finite outputs removed."); Y <- Y[ok, , drop = FALSE] }
  yMC <- colMeans(Y); UMC <- cov(Y)
  Z <- sweep(Y, 2, yMC)
  iU <- tryCatch(solve(UMC), error = function(e) NULL)
  if (is.null(iU)) message("propagateMV: the covariance matrix of the outputs is singular (linearly dependent outputs); k_p is not available.")
  d2 <- if (!is.null(iU)) rowSums((Z %*% iU) * Z) else NA                 ## squared Mahalanobis distances
  kpMC <- if (!is.null(iU)) sqrt(sort(d2)[ceiling(p * nrow(Y))]) else NA  ## 7.7.2
  kqMC <- sort(apply(abs(sweep(Z, 2, sqrt(diag(UMC)), "/")), 1, max))[ceiling(p * nrow(Y))]   ## 7.7.3
  ciMC <- sapply(seq_len(m), function(l) .covIntJCGM(Y[, l], p))
  dimnames(ciMC) <- list(paste0(c(alpha/2, 1 - alpha/2) * 100, "%"), NAMES)
  
  ## 3) GUF: sensitivity matrix C_x from one propagate() run per output, U_y = C_x U_x C_x'
  RUNS <- lapply(EXPR, function(e) { v <- intersect(colnames(DAT), all.vars(e))
    Dv <- DAT[, v, drop = FALSE]; attr(Dv, "budget") <- attr(DAT, "budget")     ## GUF from the budget (makeDat)
    quiet(e, Dv, type = type, cov = if (is.matrix(cov)) cov[v, v, drop = FALSE] else cov,
          nsim = 1E4, alpha = alpha, mc.margins = mc.margins, check = FALSE, ...) })
  Cx <- t(sapply(RUNS, function(r) { g <- setNames(rep(0, ncol(Ux)), colnames(Ux)); g[colnames(r$covMat)] <- as.numeric(r$evalGrad); g }))
  Cx <- matrix(Cx, nrow = m, dimnames = list(NAMES, colnames(Ux)))
  yGUF <- setNames(sapply(RUNS, function(r) unname(r$prop["mean.1"])), NAMES)
  UGUF <- Cx %*% Ux %*% t(Cx)
  kpGUF <- sqrt(qchisq(p, m))                                             ## 6.5.2.2 (Gaussian)
  kqGUF <- qnorm((1 + p^(1/m))/2)                                         ## hyper-rectangle for independent Gaussian components
  
  OUT <- list(guf = list(y = yGUF, U = UGUF, u = sqrt(diag(UGUF)), r = cov2cor(UGUF), kp = kpGUF, kq = kqGUF, Cx = Cx),
              mc = list(y = yMC, U = UMC, u = sqrt(diag(UMC)), r = cov2cor(UMC), kp = kpMC, kq = kqMC, ci = ciMC),
              resSIM = Y, datSIM = CAR$datSIM, covMat = Ux, expr = expr, names = NAMES, alpha = alpha, nsim = nrow(Y),
              type = type, mc.margins = if (type == "sim") NA else mc.margins, df.MC = CAR$df.MC, runs = RUNS,
              input = list(data = DAT, type = type, cov = cov))
  class(OUT) <- "propagateMV"
  OUT
}

print.propagateMV <- function(x, digits = 7, ...)
{
  for (meth in c("guf", "mc")) {
    R <- x[[meth]]
    cat(if (meth == "guf") "GUM uncertainty framework (JCGM 102:2011, 6):\n" else
        paste0("Monte Carlo method (JCGM 102:2011, 7), M = ", format(x$nsim, big.mark = " "), ", input model ", x$mc.margins, ":\n"))
    print(signif(rbind(estimate = R$y, u = R$u), digits))
    cat("Correlation matrix:\n"); print(round(R$r, 4))
    cat(sprintf("Coverage factors for %g %% coverage regions: k_p (ellipsoid) = %.3f, k_q (rectangle) = %.3f\n\n", 100 * (1 - x$alpha), R$kp, R$kq))
  }
  invisible(x)
}

summary.propagateMV <- function(object, ...) print.propagateMV(object, ...)

## pairwise plots: Monte Carlo sample with the GUF (dashed) and Monte Carlo (solid) coverage ellipses
plot.propagateMV <- function(x, pairs = NULL, npoints = 3000, ...)
{
  m <- length(x$names)
  if (is.null(pairs)) pairs <- utils::combn(m, 2, simplify = FALSE)
  op <- par(mfrow = c(1, length(pairs)), mar = c(4, 4, 2, 1)); on.exit(par(op))
  ell <- function(y, U, k) { e <- eigen(U, symmetric = TRUE); a <- seq(0, 2 * pi, length.out = 200)
    t(y + k * e$vectors %*% diag(sqrt(pmax(e$values, 0))) %*% rbind(cos(a), sin(a))) }
  i <- unique(round(seq(1, nrow(x$resSIM), length.out = min(npoints, nrow(x$resSIM)))))
  for (pr in pairs) {
    G <- ell(x$guf$y[pr], x$guf$U[pr, pr], sqrt(qchisq(1 - x$alpha, 2)))
    Mm <- ell(x$mc$y[pr], x$mc$U[pr, pr], sqrt(qchisq(1 - x$alpha, 2)))
    plot(x$resSIM[i, pr], pch = 16, cex = 0.3, col = "grey60", xlab = x$names[pr[1]], ylab = x$names[pr[2]],
         xlim = range(c(G[, 1], Mm[, 1], x$resSIM[i, pr[1]])), ylim = range(c(G[, 2], Mm[, 2], x$resSIM[i, pr[2]])), ...)
    lines(G, lty = 2, lwd = 2, col = "#7570b3"); lines(Mm, lwd = 2, col = "#1b7837")
  }
  invisible(NULL)
}