Improve floating point accuracy

R/convergence.R

@@ -495,11 +495,12 @@ mcse_sd.default <- function(x, ...) {
   # has (X-E[X])^2. The following ESS is based on a relevant quantity
   # in the computation and is empirically a good choice.
   sims_c <- x - mean(x)
-  ess <- ess_mean((sims_c)^2)
+  sims_sq <- sims_c^2
+  ess <- ess_mean(sims_sq)
   # Variance of variance estimate by Kenney and Keeping (1951, p. 141),
   # which doesn't assume normality of sims.
-  Evar <- mean(sims_c^2)
-  varvar <- (mean(sims_c^4) - Evar^2) / ess
+  Evar <- mean(sims_sq)
+  varvar <- mean((sims_sq - Evar)^2) / ess
   # The first order Taylor series approximation of variance of sd.
   # Kenney and Keeping (1951, p. 141) write "...since fluctuations of
   # any moment are of order N^{-1/2}, squares and higher powers of

R/discrete-summaries.R

@@ -148,7 +148,7 @@ dissent.default <- function(x) {
     out <- 0
   } else {
     x_i <- tab$x
-    out <- -sum(p * log2(1 - abs(x_i - mean(x)) / d))
+    out <- -sum(p * log2_one_minus(abs(x_i - mean(x)) / d))
   }
   out
 }

R/gpd.R

@@ -3,7 +3,9 @@
 #' Computes the quantile function for a generalized Pareto distribution
 #' with location `mu`, scale `sigma`, and shape `k`.
 #' 
-#' @param p Numeric vector of probabilities.
+#' @param p Numeric vector of probabilities, in `[0, 1]`, or of log
+#'   probabilities, in `[-Inf, 0]`, if `log.p` is `TRUE`. Values outside the
+#'   valid range return `NaN` with a warning.
 #' @param mu Location parameter.
 #' @param sigma Scale parameter (must be positive).
 #' @param k Shape parameter.
@@ -19,16 +21,23 @@ qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
   if (is.na(sigma) || sigma <= 0) {
     return(rep(NaN, length(p)))
   }
-  if (log.p) {
-    p <- exp(p)
+  # probabilities outside the valid range give NaN with a warning, as in the
+  # base R quantile functions. which() leaves NA and NaN inputs to propagate,
+  # and replacing them here avoids a second warning from log() or log1p().
+  invalid <- which(if (log.p) p > 0 else p < 0 | p > 1)
+  if (length(invalid)) {
+    warning_no_call("NaNs produced")
+    p[invalid] <- NaN
   }
-  if (!lower.tail) {
-    p <- 1 - p
+  if (log.p) {
+    log_survival <- if (lower.tail) log1m_exp(p) else p
+  } else {
+    log_survival <- if (lower.tail) log1p(-p) else log(p)
   }
   if (k == 0) {
-    q <-  mu - sigma * log1p(-p)
+    q <- mu - sigma * log_survival
   } else {
-    q <- mu + sigma * expm1(-k * log1p(-p)) / k
+    q <- mu + sigma * expm1(-k * log_survival) / k
   }
   q
 }
@@ -55,20 +64,17 @@ pgeneralized_pareto <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE,
     return(rep(NaN, length(q)))
   }
   z <- (q - mu) / sigma
-  if (abs(k) < 1e-15) {
-    # for very small values of indistinguishable in floating point accuracy from the case k=0
-    p <- -expm1(-z)
+  if (k == 0) {
+    log_survival <- -z
   } else {
     # pmax handles values outside the support
-    p <- -expm1(log1p(pmax(k * z, -1)) / -k)
+    log_survival <- log1p(pmax(k * z, -1)) / -k
   }
-  # force to [0, 1] for values outside the support
-  p <- pmin(pmax(p, 0), 1)
-  if (!lower.tail) {
-    p <- 1 - p
-  }
-  if (log.p) {
-    p <- log(p)
+  log_survival <- pmin(log_survival, 0)
+  if (lower.tail) {
+    p <- if (log.p) log1m_exp(log_survival) else -expm1(log_survival)
+  } else {
+    p <- if (log.p) log_survival else exp(log_survival)
   }
   p
 }

R/misc-math.R

@@ -0,0 +1,41 @@
+# numerically stable version of log(sum(exp(x)))
+log_sum_exp <- function(x) {
+  x <- as.numeric(x)
+  if (length(x) == 0) {
+    return(-Inf)
+  }
+  max <- max(x)
+  if (max == -Inf) {
+    res <- -Inf
+  } else if (max == Inf) {
+    res <- Inf
+  } else {
+    sum <- sum(exp(x - max))
+    res <- max + log(sum)
+  }
+  res
+}
+
+# numerically stable version of exp(x) - exp(y)
+exp_x_minus_exp_y <- function(x, y) {
+  out <- -exp(x) * expm1(y - x)
+  # equal infinite inputs give 0 * NaN above; which() drops the NA comparisons
+  # that NA or NaN inputs would produce
+  out[which(x == y)] <- 0
+  out
+}
+
+# numerically stable version of log2(1 - x)
+log2_one_minus <- function(x) {
+  log1p(-x) / log(2)
+}
+
+# numerically stable version of log(1 - exp(x)) for x <= 0
+log1m_exp <- function(x) {
+  log(-expm1(x))
+}
+
+# numerically stable version of log(exp(x) / sum(exp(x)))
+log_normalize <- function(x) {
+  x - log_sum_exp(x)
+}

R/misc.R

@@ -219,20 +219,6 @@ escape_all <- function(x) {
   gsub(specials, "\\\\\\1", x)
 }
 
-# numerically stable version of log(sum(exp(x)))
-log_sum_exp <- function(x) {
-  max <- max(as.numeric(x), warnings = FALSE)
-  if (max == -Inf) {
-    res <- 0
-  } else if (max == Inf) {
-    res <- Inf
-  } else {
-    sum <- sum(exp(x - max))
-    res <- max + log(sum)
-  }
-  res
-}
-
 # simple version of destructuring assignment
 `%<-%` <- function(vars, values, envir = parent.frame()) {
   vars <- as.character(substitute(vars)[-1])

R/pareto_smooth.R

@@ -525,10 +525,12 @@ ps_tail <- function(x,
   max_tail <- max(draws_tail)
 
   if (are_log_weights) {
-    draws_tail <- exp(draws_tail)
+    tail_excesses <- exp_x_minus_exp_y(draws_tail, cutoff)
     cutoff <- exp(cutoff)
+  } else {
+    tail_excesses <- draws_tail - cutoff
   }
-  fit <- gpdfit(draws_tail - cutoff, sort_x = FALSE, ...)
+  fit <- gpdfit(tail_excesses, sort_x = FALSE, ...)
   k <- fit$k
   sigma <- fit$sigma
   if (is.finite(k) && smooth_draws) {
@@ -631,22 +633,31 @@ ps_khat_threshold <- function(ndraws, ...) {
 #' @return convergence rate
 #' @export
 ps_convergence_rate <- function(k, ndraws, ...) {
-  # allow array of k's
-  rate <- numeric(length(k))
-  # k<0 bounded distribution
-  rate[k < 0] <- 1
-  # k>0 non-finite mean
-  rate[k > 1] <- 0
-  # limit value at k=1/2
-  rate[k == 0.5] <- 1 - 1 / log(ndraws)
-  # smooth approximation for the rest (see Appendix B of PSIS paper)
-  ki <- (k > 0 & k < 1 & k != 0.5)
+  # allow array of k's; which() keeps NA and NaN shapes out of every branch,
+  # so they keep the NA initialization instead of silently reading as 0
+  rate <- rep(NA_real_, length(k))
+  # k<0 bounded distribution, k=0 exponential tail; both converge at rate 1
+  rate[which(k <= 0)] <- 1
+  # k>=1 non-finite mean
+  rate[which(k >= 1)] <- 0
+  # smooth approximation for the rest
+  # this is a numerically stable calculation of the equation on page 31, in Appendix B of PSIS paper
+  ki <- which(k > 0 & k < 1)
   kk <- k[ki]
-  rate[ki] <- pmax(
-    0,
-    (2 * (kk - 1) * ndraws^(2 * kk + 1) + (1 - 2 * kk) * ndraws^(2 * kk) + ndraws^2) /
-      ((ndraws - 1) * (ndraws - ndraws^(2 * kk)))
-  )
+  aa <- 2 * kk - 1
+  near_half <- abs(aa) <= abs(1 - aa)
+  rate_ki <- numeric(length(kk))
+  rate_ki[aa == 0] <- ndraws / (ndraws - 1) - 1 / log(ndraws)
+  use_half <- near_half & aa != 0
+  rate_ki[use_half] <-
+    ndraws / (ndraws - 1) -
+    aa[use_half] / -expm1(-aa[use_half] * log(ndraws))
+  use_one <- !near_half
+  bb <- 1 - aa[use_one]
+  rate_ki[use_one] <-
+    (bb * (ndraws - 1) - expm1(bb * log(ndraws))) /
+    ((ndraws - 1) * -expm1((bb - 1) * log(ndraws)))
+  rate[ki] <- pmax(0, rate_ki)
   rate
 }
 

R/pit.R

@@ -79,12 +79,12 @@ pit.default <- function(x, y, weights = NULL, log = FALSE, ...) {
 pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE, ...) {
   y <- validate_y(y, x)
   if (!is.null(weights)) {
-    weights <- sapply(seq_len(nvariables(x)), function(var_idx) {
-      validate_weights(weights[, var_idx], x[, var_idx], log)
-    })
-    weights <- normalize_log_weights(weights)
+    weights <- validate_pit_weights(weights, x, log)
   }
   pit <- vapply(seq_len(ncol(x)), function(j) {
+    if (!is.null(weights) && anyNA(weights[, j])) {
+      return(NA_real_)
+    }
     sel_min <- x[, j] < y[j]
     if (!any(sel_min)) {
       pit <- 0
@@ -112,7 +112,7 @@ pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE, ...) {
     pit
   }, FUN.VALUE = 1.0)
 
-  if (any(pit > 1 + 1e-10)) {
+  if (any(pit > 1 + 1e-10, na.rm = TRUE)) {
     warning_no_call(
       paste(
         "Some PIT values larger than 1. ",
@@ -241,10 +241,7 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
 
   # validate and normalize weights to log scale (same as pit.draws_matrix)
   if (!is.null(weights)) {
-    weights <- sapply(seq_len(nvariables(x)), function(var_idx) {
-      validate_weights(weights[, var_idx], x[, var_idx], log)
-    })
-    weights <- normalize_log_weights(weights)
+    weights <- validate_pit_weights(weights, x, log)
   }
 
   ndraws <- ndraws(x)
@@ -270,6 +267,9 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
   }
 
   pit_values <- vapply(seq_len(ncol(x)), function(j) {
+    if (!is.null(weights) && anyNA(weights[, j])) {
+      return(NA_real_)
+    }
     draws <- x[, j]
 
     # --- raw PIT (same logic as pit.draws_matrix) ---
@@ -312,9 +312,13 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
     }
 
     # --- right tail ---
+    # A tail carrying no probability mass cannot inform the fit: the raw PIT is
+    # already the whole answer, and gpdfit() would be handed all-zero weights.
+    right_tail_empty <- !is.null(log_wt_sorted) && tail_proportion == 0
+
     right_replaced <- FALSE
     right_tail <- sorted[tail_ids]
-    if (!is_constant(right_tail)) {
+    if (!right_tail_empty && !is_constant(right_tail)) {
       right_cutoff <- sorted[min(tail_ids) - 1]
       if (right_cutoff == right_tail[1]) {
         right_cutoff <- right_cutoff - .Machine$double.eps
@@ -344,8 +348,11 @@ pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
       left_sorted <- left_ord$x
       log_wt_left_sorted <- if (!is.null(weights)) weights[left_ord$ix, j] else NULL
 
+      left_tail_empty <- !is.null(log_wt_left_sorted) &&
+        log_sum_exp(log_wt_left_sorted[tail_ids]) == -Inf
+
       left_tail <- left_sorted[tail_ids]
-      if (!is_constant(left_tail)) {
+      if (!left_tail_empty && !is_constant(left_tail)) {
         left_cutoff <- left_sorted[min(tail_ids) - 1]
         if (left_cutoff == left_tail[1]) {
           left_cutoff <- left_cutoff - .Machine$double.eps
@@ -423,5 +430,44 @@ validate_y <- function(y, x = NULL) {
 }
 
 normalize_log_weights <- function(log_weights) {
-  apply(log_weights, 2, function(col) col - log_sum_exp(col))
+  apply(log_weights, 2, log_normalize)
+}
+
+#' Validate and normalize a per-variable weights matrix
+#'
+#' Columns whose weights carry no probability mass are reported and returned as
+#' `NA`, so that one degenerate variable does not take down the whole call. This
+#' matches [pareto_khat()], which warns and returns `NA` rather than erroring.
+#' Weights that are malformed for any other reason still raise an error.
+#'
+#' @noRd
+#' @param weights A matrix of weights, one column per variable in `x`.
+#' @param x A `draws_matrix`.
+#' @param log (logical) Are the weights already on the log scale?
+#' @return A matrix of normalized log weights, `NA` for degenerate columns.
+#'
+validate_pit_weights <- function(weights, x, log) {
+  vars <- variables(x)
+  cols <- lapply(seq_along(vars), function(var_idx) {
+    tryCatch(
+      validate_weights(weights[, var_idx], x[, var_idx], log),
+      posterior_degenerate_weights_error = function(cnd) NULL
+    )
+  })
+  degenerate <- vapply(cols, is.null, logical(1))
+  if (any(degenerate)) {
+    warning_no_call(
+      "All draws have zero weight for ",
+      paste0("'", vars[degenerate], "'", collapse = ", "),
+      ". Returning NA for ",
+      if (sum(degenerate) > 1) "those variables." else "that variable."
+    )
+    cols[degenerate] <- list(rep(NA_real_, ndraws(x)))
+  }
+  out <- do.call(cbind, cols)
+  # normalize only the columns that carry mass; log_sum_exp() has no NA handling
+  if (any(!degenerate)) {
+    out[, !degenerate] <- normalize_log_weights(out[, !degenerate, drop = FALSE])
+  }
+  out
 }

R/weight_draws.R

@@ -6,6 +6,13 @@
 #' `.log_weight`). See [weights.draws()] for details how to extract weights from
 #' `draws` objects.
 #'
+#' Because the stored log-weights are unnormalized and are normalized only when
+#' they are extracted, subsetting a weighted `draws` object conditions on the
+#' retained draws: the weights of the draws that remain are renormalized to sum
+#' to one. Any operation that drops draws has this effect, including
+#' [subset_draws()], [thin_draws()] and `[` indexing. Weights are therefore
+#' comparable only within one subset, not across subsets of the same object.
+#'
 #' @template args-methods-x
 #' @param weights (numeric vector) A vector of weights of length `ndraws(x)`.
 #'   Weights will be internally stored on the log scale (in a variable called
@@ -169,7 +176,12 @@ weights.draws <- function(object, log = FALSE, normalize = TRUE, ...) {
   }
   out <- extract_variable(object, ".log_weight")
   if (normalize) {
-    out <- out - log_sum_exp(out)
+    # weight_draws() rejects this at construction, but dropping draws from an
+    # object that was valid can still remove every weighted draw
+    if (!any(is.finite(out))) {
+      stop_no_call("All draws have zero weight.")
+    }
+    out <- log_normalize(out)
   }
   if (!log) {
     out <- exp(out)
@@ -193,6 +205,18 @@ validate_weights <- function(weights, draws, log = FALSE) {
     }
     weights <- log(weights)
   }
+  if (!any(is.finite(weights))) {
+    # classed so that callers working variable by variable, such as pit(), can
+    # handle one degenerate column without aborting the whole call
+    stop(errorCondition(
+      paste0(
+        "All weights are zero, so no draw carries any probability mass. ",
+        "If the weights underflowed to zero, pass them on the log scale ",
+        "with `log = TRUE`."
+      ),
+      class = "posterior_degenerate_weights_error"
+    ))
+  }
   weights
 }