fix merge conflicts for rvar_weights

.github/workflows/covr.yml

@@ -35,7 +35,7 @@ jobs:
         shell: Rscript {0}
 
       - name: Cache R packages
-        uses: actions/cache@v2
+        uses: actions/cache@v4
         with:
           path: ${{ env.R_LIBS_USER }}
           key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}

.github/workflows/rcmdcheck.yml

@@ -58,7 +58,7 @@ jobs:
 
       - name: Cache R packages
         if: runner.os != 'Windows'
-        uses: actions/cache@v2
+        uses: actions/cache@v4
         with:
           path: ${{ env.R_LIBS_USER }}
           key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
@@ -98,7 +98,7 @@ jobs:
 
       - name: Upload check results
         if: failure()
-        uses: actions/upload-artifact@main
+        uses: actions/upload-artifact@v4
         with:
           name: ${{ runner.os }}-r${{ matrix.config.r }}-results
           path: check

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: posterior
 Title: Tools for Working with Posterior Distributions
-Version: 1.5.0.9000
-Date: 2023-10-31
+Version: 1.6.1
+Date: 2025-02-27
 Authors@R: c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com", role = c("aut", "cre")),
              person("Jonah", "Gabry", email = "jsg2201@columbia.edu", role = c("aut")),
              person("Matthew", "Kay", email = "mjskay@northwestern.edu", role = c("aut")),
@@ -11,7 +11,8 @@ Authors@R: c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.c
              person("Ben", "Lambert", role = c("ctb")),
              person("Ozan", "Adıgüzel", role = c("ctb")), 
              person("Jacob", "Socolar", role = c("ctb")),
-             person("Noa", "Kallioinen", role = c("ctb")))
+             person("Noa", "Kallioinen", role = c("ctb")),
+             person("Teemu", "Säilynoja", role = c("ctb")))
 Description: Provides useful tools for both users and developers of packages 
   for fitting Bayesian models or working with output from Bayesian models. 
   The primary goals of the package are to: 
@@ -56,5 +57,5 @@ LazyData: false
 URL: https://mc-stan.org/posterior/, https://discourse.mc-stan.org/
 BugReports: https://github.com/stan-dev/posterior/issues
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.0
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr

NAMESPACE

@@ -245,6 +245,9 @@ S3method(pareto_min_ss,rvar)
 S3method(pareto_smooth,default)
 S3method(pareto_smooth,rvar)
 S3method(pillar_shaft,rvar)
+S3method(pit,default)
+S3method(pit,draws_matrix)
+S3method(pit,rvar)
 S3method(print,draws_array)
 S3method(print,draws_df)
 S3method(print,draws_list)
@@ -480,6 +483,11 @@ export(pareto_khat)
 export(pareto_khat_threshold)
 export(pareto_min_ss)
 export(pareto_smooth)
+export(pit)
+export(ps_convergence_rate)
+export(ps_khat_threshold)
+export(ps_min_ss)
+export(ps_tail_length)
 export(quantile2)
 export(r_scale)
 export(rdo)

NEWS.md

@@ -1,11 +1,33 @@
-# posterior (development version)
+# posterior 1.6.1
+
+### Bug Fixes
+
+* Fix a test issue that led to an R CMD check failure on R devel.
 
 ### Enhancements
 
-* Add `pareto_smooth` option to `weight_draws`, to Pareto smooth
-  weights before adding to a draws object.
+* Convert lists of matrices to `draws_array` objects.
+* Improve the documentation in various places.
+* Implement `pit()` for draws and rvar objects. LOO-PIT can be computed using `weights`.
 * Add support for applying weights to individual `rvar` objects.
 * Add `log_weights()` function for easy access to raw internal weights.
+* Add support for weights in `ess` and `mcse` functions.
+
+# posterior 1.6.0
+
+### Enhancements
+
+* Add `exclude` option to `subset_draws()`, which can be used to exclude
+  the matched selection.
+* Add `are_log_weights` option to `pareto_smooth()`, which is necessary
+  for correct Pareto smoothing computation if the input vector
+  consists of log weights.
+* Add `pareto_smooth` option to `weight_draws()`, to Pareto smooth
+  weights before adding to a draws object.
+* Add individual Pareto diagnostic functions (`pareto_khat()`,
+  `pareto_khat_threshold()`, `pareto_min_ss()`, `pareto_convergence_rate()`)
+* `thin_draws()` now automatically thins draws based on ESS by default,
+  and non-integer thinning is possible.
 * Matrix multiplication of `rvar`s can now be done with the base matrix 
   multiplication operator (`%*%`) instead of `%**%` in R >= 4.3.
 * `variables()`, `variables<-()`, `set_variables()`, and `nvariables()` now

R/as_draws.R

@@ -76,10 +76,11 @@ closest_draws_format <- function(x) {
     out <- "rvars"
   } else if (is_draws_list_like(x)) {
     out <- "list"
-  }
-  else {
-    stop_no_call("Don't know how to transform an object of class ",
-          "'", class(x)[1L], "' to any supported draws format.")
+  } else {
+    stop_no_call(
+      "Don't know how to transform an object of class '",
+      class(x)[1L], "' to any supported draws format."
+    )
   }
   paste0("draws_", out)
 }

R/as_draws_array.R

@@ -131,7 +131,11 @@ as_draws_array.mcmc.list <- function(x, ...) {
 
 # try to convert any R object into a 'draws_array' object
 .as_draws_array <- function(x) {
-  x <- as.array(x)
+  if (is_matrix_list_like(x)) {
+    x <- as_array_matrix_list(x)
+  } else {
+    x <- as.array(x)
+  }
   new_dimnames <- list(iteration = NULL, chain = NULL, variable = NULL)
   if (!is.null(dimnames(x)[[3]])) {
     new_dimnames[[3]] <- dimnames(x)[[3]]
@@ -178,7 +182,14 @@ is_draws_array <- function(x) {
 
 # is an object looking like a 'draws_array' object?
 is_draws_array_like <- function(x) {
-  is.array(x) && length(dim(x)) == 3L
+  is.array(x) && length(dim(x)) == 3L ||
+    is_matrix_list_like(x)
+}
+
+# is an object likely a list of matrices?
+# such an object can be easily converted to a draws_array
+is_matrix_list_like <- function(x) {
+  is.list(x) && length(dim(x[[1]])) == 2L
 }
 
 #' Extract parts of a `draws_array` object
@@ -217,15 +228,8 @@ variance.draws_array <- function(x, ...) {
 # convert a list of matrices to an array
 as_array_matrix_list <- function(x) {
   stopifnot(is.list(x))
-  if (length(x) == 1) {
-    tmp <- dimnames(x[[1]])
-    x <- x[[1]]
-    dim(x) <- c(dim(x), 1)
-    dimnames(x) <- tmp
-  } else {
-    x <- abind::abind(x, along = 3L)
-  }
-  x <- aperm(x, c(1, 3, 2))
+  x <- abind::abind(x, along = 3L)
+  aperm(x, c(1, 3, 2))
 }
 
 # create an empty draws_array object
@@ -246,4 +250,3 @@ empty_draws_array <- function(variables = character(0), nchains = 0,
   class(out) <- class_draws_array()
   out
 }
-

R/convergence.R

@@ -16,11 +16,14 @@
 #' | [ess_basic()] | Basic version of effective sample size |
 #' | [ess_bulk()] | Bulk effective sample size |
 #' | [ess_tail()] | Tail effective sample size |
+#' | [ess_mean()] | Effective sample sizes for the mean |
+#' | [ess_median()] | Effective sample sizes for the median |
 #' | [ess_quantile()] | Effective sample sizes for quantiles |
-#' | [ess_sd()] | Effective sample sizes for standard deviations |
+#' | [ess_sd()] | Effective sample sizes for the standard deviation |
 #' | [mcse_mean()] | Monte Carlo standard error for the mean |
+#' | [mcse_median()] | Monte Carlo standard error for the median |
 #' | [mcse_quantile()] | Monte Carlo standard error for quantiles |
-#' | [mcse_sd()] | Monte Carlo standard error for standard deviations |
+#' | [mcse_sd()] | Monte Carlo standard error for the standard deviation |
 #' | [pareto_khat()] | Pareto khat diagnostic for tail(s) |
 #' | [pareto_diags()] | Additional diagnostics related to Pareto khat |
 #' | [rhat_basic()] | Basic version of Rhat |
@@ -88,6 +91,7 @@ rhat_basic.rvar <- function(x, split = TRUE, ...) {
 #' @family diagnostics
 #' @template args-conv
 #' @template args-conv-split
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv
 #' @template ref-gelman-bda-2013
@@ -178,6 +182,7 @@ rhat.rvar <- function(x, ...) {
 #'
 #' @family diagnostics
 #' @template args-conv
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv
 #' @template ref-vehtari-rhat-2021
@@ -231,6 +236,7 @@ ess_bulk.rvar <- function(x, ...) {
 #'
 #' @family diagnostics
 #' @template args-conv
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv
 #' @template ref-vehtari-rhat-2021
@@ -270,6 +276,7 @@ ess_tail.rvar <- function(x, ...) {
 #' @family diagnostics
 #' @template args-conv
 #' @template args-conv-quantile
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv-quantile
 #' @template ref-vehtari-rhat-2021
@@ -360,6 +367,7 @@ ess_median <- function(x, ...) {
 #' estimate of a single variable.
 #'
 #' @template args-conv
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv
 #' @template ref-gelman-bda-2013
@@ -374,7 +382,7 @@ ess_median <- function(x, ...) {
 #' @export
 ess_mean <- function(x, ...) UseMethod("ess_mean")
 
-#' @rdname ess_quantile
+#' @rdname ess_mean
 #' @export
 ess_mean.default <- function(x, weights = NULL, ...) {
 
@@ -410,6 +418,7 @@ ess_mean.rvar <- function(x, ...) {
 #'
 #' @family diagnostics
 #' @template args-conv
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv
 #' @template ref-vehtari-rhat-2021
@@ -428,6 +437,9 @@ ess_sd <- function(x, ...) UseMethod("ess_sd")
 #' @export
 ess_sd.default <- function(x, weights = NULL, ...) {
   if (is.null(weights)) {
+    # var/sd are not a simple expectation of g(X), e.g. variance
+    # has (X-E[X])^2. The following ESS is based on a relevant quantity
+    # in the computation and is empirically a good choice.
     .ess(.split_chains(abs(x - mean(x))))
   } else {
 
@@ -449,7 +461,7 @@ ess_sd.rvar <- function(x, ...) {
   summarise_rvar_by_element_with_chains(x, ess_sd, weights = weights, ...)
 }
 
-# TODO: ess_weights
+# TODO: ess_weights function
 
 #' Monte Carlo standard error for quantiles
 #'
@@ -460,6 +472,7 @@ ess_sd.rvar <- function(x, ...) {
 #' @family diagnostics
 #' @template args-conv
 #' @template args-conv-quantile
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv-quantile
 #' @template ref-vehtari-rhat-2021
@@ -555,6 +568,7 @@ mcse_median <- function(x, ...) {
 #'
 #' @family diagnostics
 #' @template args-conv
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv
 #' @template ref-gelman-bda-2013
@@ -583,7 +597,7 @@ mcse_mean.default <- function(x, weights = NULL, ...) {
     x <- as.matrix(x)
 
     r_eff <- .ess(.split_chains(x)) / (nrow(x) * ncol(x))
-    .mcse_weighted(x, weights, r_eff, ...)
+    .mcse_weighted(x, weights = weights, r_eff = r_eff, ...)
   }
 }
 
@@ -602,6 +616,7 @@ mcse_mean.rvar <- function(x, ...) {
 #'
 #' @family diagnostics
 #' @template args-conv
+#' @template args-conv-weights
 #' @template args-methods-dots
 #' @template return-conv
 #' @template ref-vehtari-rhat-2021
@@ -632,7 +647,6 @@ mcse_sd.default <- function(x, weights = NULL, ...) {
     # which doesn't assume normality of sims.
     Evar <- mean(sims_c^2)
     varvar <- (mean(sims_c^4) - Evar^2) / ess # (Equation 6.20)
-
     # 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
@@ -656,7 +670,7 @@ mcse_sd.default <- function(x, weights = NULL, ...) {
     second_moment_weighted <- weighted.mean(x_centered^2, w = weights)
     fourth_moment_weighted <- weighted.mean(x_centered^4, w = weights)
 
-    r_eff <- .ess(x_centered^2) / (nrow(x) * ncol(x))
+    r_eff <- .ess(.split_chains(x_centered^2)) / (nrow(x) * ncol(x))
     weighted_ess <- .ess_weighted(x_centered^2, weights = weights, r_eff = r_eff)
 
     # Kenney and Keeping (1951, eq 6.20)
@@ -681,6 +695,7 @@ mcse_sd.rvar <- function(x, ...) {
 #' with other summary functions in the \pkg{posterior} package.
 #'
 #' @template args-conv
+#' @template args-conv-weights
 #' @template args-conv-quantile
 #' @param na.rm (logical) Should `NA` and `NaN` values be removed from `x` prior
 #'   to computing quantiles? The default is `FALSE`.
@@ -697,19 +712,19 @@ mcse_sd.rvar <- function(x, ...) {
 #' quantile2(mu)
 #'
 #' @export
-quantile2 <- function(x, probs = c(0.05, 0.95), na.rm = FALSE, ...) {
+quantile2 <- function(x, ...) {
   UseMethod("quantile2")
 }
 
 #' @rdname quantile2
 #' @export
 quantile2.default <- function(
-  x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, ...
+  x, probs = c(0.05, 0.95), na.rm = FALSE, names = TRUE, weights = NULL, ...
 ) {
   names <- as_one_logical(names)
   na.rm <- as_one_logical(na.rm)
 
-  out <- weighted_quantile(x, probs = probs, na.rm = na.rm, ...)
+  out <- weighted_quantile(x, probs = probs, na.rm = na.rm, weights = weights, ...)
 
   if (names) {
     names(out) <- paste0("q", probs * 100)
@@ -726,8 +741,8 @@ quantile2.rvar <- function(
 ) {
   weights <- weights(x)
   summarise_rvar_by_element(x, function(draws) {
-    quantile2(
-      draws, probs = probs, weights = weights, na.rm = na.rm, names = names, ...
+    quantile2.default(
+      draws, probs = probs, na.rm = na.rm, names = names, weights = weights, ...
     )
   })
 }
@@ -777,6 +792,7 @@ autocovariance <- function(x) {
 autocorrelation <- function(x) {
   ac <- autocovariance(x)
   ac <- ac / ac[1]
+  ac
 }
 
 #' Rank normalization
@@ -807,8 +823,10 @@ z_scale <- function(x, c = 3/8) {
 #'
 #' Compute rank uniformization for a numeric array. First replace each value by
 #' its rank. Average rank for ties are used to conserve the number of unique
-#' values of discrete quantities. Second, uniformize ranks to the scale
-#' `[1/(2S), 1-1/(2S)]`, where `S` is the number of values.
+#' values of discrete quantities. Second, uniformize ranks using formula
+#' `(r - c) / (S - 2 * c + 1)`, where `r` is a rank, `S` is the number of
+#' values, and `c` is a fractional offset which defaults to c = 3/8 as
+#' recommend by Blom (1958).
 #'
 #' @template args-scale
 #' @template args-frac-offset

R/draws-index.R

@@ -68,6 +68,10 @@ iteration_ids.draws_df <- function(x) {
 
 #' @export
 iteration_ids.draws_list <- function(x) {
+  if (!length(x) || !length(x[[1]])) {
+    # no chains or no variables within chains
+    return(integer(0))
+  }
   seq_along(x[[1]][[1]])
 }