R/corDist.R

In propagate: Propagation of Uncertainty

corDist <- function(mean, u, dist = "normal", cor = NULL, df = Inf, d = NA, nsim = 1E6,
                    cor.type = c("pearson", "copula"))
{
  ## Correlated Monte Carlo samples of input quantities with stated distributions, joined by a
  ## Gaussian copula (Sklar's theorem; Possolo 2010), for propagate(..., type = "sim", cov = "full").
  ## Margins (expectation 'mean', standard deviation 'u'): "normal", "rect", "tri", "arcsine",
  ## "ctrap" (inexactness 'd', JCGM 101:2008, 6.4.3), "t" (scale u, 'df' degrees of freedom, 6.4.9);
  ## or any sampler function(n) (e.g. the rDistr samplers or rgamma) or a numeric sample, whose sorted
  ## values are assigned in the rank order of the copula (rank reordering; mean and u from the sample).
  ## cor.type = "pearson": the copula correlation is adjusted so that the Pearson correlation of the
  ## samples equals 'cor' (NORTA, Cario & Nelson 1997).
  cor.type <- match.arg(cor.type)
  nm <- names(mean); p <- length(mean)
  if (is.null(nm) || any(nm == "")) stop("corDist: 'mean' must be a named vector (names = variables of the expression).")
  rec <- function(v) { v <- if (is.list(v)) v else as.list(v); rep_len(v, p) }
  u <- rep_len(u, p); dist <- rec(dist); df <- rep_len(df, p); d <- rep_len(d, p)
  if (is.null(cor)) cor <- diag(p)
  if (!isTRUE(all.equal(dim(cor), c(p, p))) || any(abs(diag(cor) - 1) > 1E-12)) stop("corDist: 'cor' must be a p x p correlation matrix.")
  ## samplers and samples: sorted values (empirical quantile function); mean and u from the sample
  SRT <- lapply(seq_len(p), function(i) {
    s <- dist[[i]]
    if (is.function(s)) { v <- s(nsim); if (length(v) != nsim) stop("corDist: sampler for '", nm[i], "' must return n values.") }
    else if (is.numeric(s) && length(s) > 1) v <- if (length(s) == nsim) s else sample(s, nsim, replace = TRUE)
    else return(NULL)
    if (any(!is.finite(v))) stop("corDist: non-finite values in the sample of '", nm[i], "'.")
    sort(v) })
  isS <- !sapply(SRT, is.null)
  mean[isS] <- vapply(SRT[isS], base::mean, numeric(1)); u[isS] <- vapply(SRT[isS], sd, numeric(1)); df[isS] <- Inf
  ## quantile functions of the margins
  QF <- lapply(seq_len(p), function(i) {
    s <- dist[[i]]; x <- mean[i]; ui <- u[i]
    if (isS[i]) { S <- SRT[[i]]; return(function(v) S[pmin(nsim, pmax(1, ceiling(v * nsim)))]) }
    if (!is.character(s) || length(s) != 1) stop("corDist: 'dist' of '", nm[i], "' must be a name, a sampler function(n) or a sample.")
    switch(s,
      normal  = function(v) x + ui * qnorm(v),
      t       = function(v) x + ui * qt(v, df[i]),
      rect    = function(v) x + sqrt(3) * ui * (2 * v - 1),
      tri     = function(v) x + sqrt(6) * ui * ifelse(v < 0.5, sqrt(2 * v) - 1, 1 - sqrt(2 * (1 - v))),
      arcsine = function(v) x - sqrt(2) * ui * cos(pi * v),
      ctrap   = { di <- d[i]; a2 <- 3 * (ui^2 - di^2/9)
                  if (!is.finite(di) || a2 <= di^2) stop("corDist: 'ctrap' needs d with u^2 > 4 d^2/9.")
                  a <- sqrt(a2); zg <- seq(0, a + di, length.out = 2001)          ## P(|X - x| <= z), inverted
                  Pz <- sapply(zg, function(z) integrate(function(s) pmin(1, z/s), a - di, a + di)$value/(2 * di))
                  function(v) x + sign(v - 0.5) * approx(Pz, zg, xout = abs(2 * v - 1), ties = "ordered", rule = 2)$y },
      stop("corDist: unknown distribution '", s, "'."))
  })
  ## copula correlation: Pearson matching by root finding with 2D Gauss-Hermite quadrature
  RZ <- cor
  if (cor.type == "pearson" && p > 1 && any(cor[upper.tri(cor)] != 0)) {
    n <- 60; J <- matrix(0, n, n); J[cbind(1:(n - 1), 2:n)] <- J[cbind(2:n, 1:(n - 1))] <- sqrt(1:(n - 1))
    E <- eigen(J, symmetric = TRUE); z <- E$values; w <- E$vectors[1, ]^2
    Phi <- function(q) pmin(pmax(pnorm(q), 1E-15), 1 - 1E-15)
    G <- lapply(QF, function(q) { g <- q(Phi(z)); m <- sum(w * g); (g - m)/sqrt(sum(w * (g - m)^2)) })
    ok <- !(sapply(dist, identical, "t") & df <= 4)                ## no finite 4th moment: copula value kept
    for (i in 1:(p - 1)) for (j in (i + 1):p) {
      if (cor[i, j] == 0 || !ok[i] || !ok[j]) next
      mj <- sum(w * QF[[j]](Phi(z))); sj <- sqrt(sum(w * (QF[[j]](Phi(z)) - mj)^2))
      rho <- function(rz) sum(outer(w, w) * G[[i]] * (QF[[j]](Phi(outer(rz * z, sqrt(1 - rz^2) * z, "+"))) - mj)/sj)
      lim <- c(rho(-0.9999), rho(0.9999))
      if (cor[i, j] < lim[1] || cor[i, j] > lim[2])
        stop("corDist: correlation ", cor[i, j], " of '", nm[i], "' and '", nm[j], "' is not attainable with these margins (",
             round(lim[1], 3), " to ", round(lim[2], 3), ").")
      RZ[i, j] <- RZ[j, i] <- uniroot(function(rz) rho(rz) - cor[i, j], c(-0.9999, 0.9999), tol = 1E-8)$root
    }
  }
  L <- try(chol(RZ), silent = TRUE)
  if (inherits(L, "try-error")) stop("corDist: the (copula) correlation matrix is not positive definite.")
  Z <- matrix(rnorm(nsim * p), ncol = p) %*% L
  X <- sapply(seq_len(p), function(i)
    if (isS[i]) SRT[[i]][rank(Z[, i], ties.method = "first")] else QF[[i]](pnorm(Z[, i])))   ## rank reordering for samples
  X <- matrix(X, ncol = p, dimnames = list(NULL, nm))
  ## the uncertainty budget makes propagate() use the stated estimates, u and DOF for the GUF
  ## (with cov = "full": sample correlation and budget u)
  attr(X, "budget") <- data.frame(name = nm, mean = unname(mean), u = u, df = df,
                                  dist = ifelse(isS, "sample", sapply(dist, function(s) if (is.character(s)) s else "sample")), stringsAsFactors = FALSE)
  attr(X, "corMat.copula") <- RZ
  X
}