Merge branch 'master' into diff-diagnostics

Author: jgabry

R/psis.R

@@ -98,13 +98,14 @@ psis <- function(log_ratios, ...) UseMethod("psis")
 #' @template array
 #'
 psis.array <-
-  function(log_ratios, ...,
-           r_eff = 1,
-           cores = getOption("mc.cores", 1)) {
-  importance_sampling.array(log_ratios = log_ratios, ...,
-                            r_eff = r_eff,
-                            cores = cores,
-                            method = "psis")
+  function(log_ratios, ..., r_eff = 1, cores = getOption("mc.cores", 1)) {
+    importance_sampling.array(
+      log_ratios = log_ratios,
+      ...,
+      r_eff = r_eff,
+      cores = cores,
+      method = "psis"
+    )
   }
 
 
@@ -113,15 +114,14 @@ psis.array <-
 #' @template matrix
 #'
 psis.matrix <-
-  function(log_ratios,
-           ...,
-           r_eff = 1,
-           cores = getOption("mc.cores", 1)) {
-    importance_sampling.matrix(log_ratios,
-                               ...,
-                               r_eff = r_eff,
-                               cores = cores,
-                               method = "psis")
+  function(log_ratios, ..., r_eff = 1, cores = getOption("mc.cores", 1)) {
+    importance_sampling.matrix(
+      log_ratios,
+      ...,
+      r_eff = r_eff,
+      cores = cores,
+      method = "psis"
+    )
   }
 
 #' @export
@@ -130,9 +130,12 @@ psis.matrix <-
 #'
 psis.default <-
   function(log_ratios, ..., r_eff = 1) {
-    importance_sampling.default(log_ratios = log_ratios, ...,
-                                r_eff = r_eff,
-                                method = "psis")
+    importance_sampling.default(
+      log_ratios = log_ratios,
+      ...,
+      r_eff = r_eff,
+      method = "psis"
+    )
   }
 
 
@@ -149,25 +152,26 @@ is.psis <- function(x) {
 #' @noRd
 #' @seealso importance_sampling_object
 psis_object <-
-  function(unnormalized_log_weights,
-           pareto_k,
-           tail_len,
-           r_eff) {
-    importance_sampling_object(unnormalized_log_weights = unnormalized_log_weights,
-                               pareto_k = pareto_k,
-                               tail_len = tail_len,
-                               r_eff = r_eff,
-                               method = "psis")
+  function(unnormalized_log_weights, pareto_k, tail_len, r_eff) {
+    importance_sampling_object(
+      unnormalized_log_weights = unnormalized_log_weights,
+      pareto_k = pareto_k,
+      tail_len = tail_len,
+      r_eff = r_eff,
+      method = "psis"
+    )
   }
 
 
 #' @noRd
 #' @seealso do_importance_sampling
-do_psis <- function(log_ratios, r_eff, cores, method){
-  do_importance_sampling(log_ratios = log_ratios,
-                         r_eff = r_eff,
-                         cores = cores,
-                         method = "psis")
+do_psis <- function(log_ratios, r_eff, cores, method) {
+  do_importance_sampling(
+    log_ratios = log_ratios,
+    r_eff = r_eff,
+    cores = cores,
+    method = "psis"
+  )
 }
 
 #' Extract named components from each list in the list of lists obtained by
@@ -181,7 +185,9 @@ do_psis <- function(log_ratios, r_eff, cores, method){
 #' @return Numeric vector or matrix.
 #'
 psis_apply <- function(x, item, fun = c("[[", "attr"), fun_val = numeric(1)) {
-  if (!is.list(x)) stop("Internal error ('x' must be a list for psis_apply)")
+  if (!is.list(x)) {
+    stop("Internal error ('x' must be a list for psis_apply)")
+  }
   vapply(x, FUN = match.arg(fun), FUN.VALUE = fun_val, item)
 }
 
@@ -212,7 +218,7 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) {
     ord <- sort.int(lw_i, index.return = TRUE)
     tail_ids <- seq(S - tail_len_i + 1, S)
     lw_tail <- ord$x[tail_ids]
-    if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps/100) {
+    if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps / 100) {
       warning(
         "Can't fit generalized Pareto distribution ",
         "because all tail values are the same.",
@@ -252,11 +258,11 @@ psis_smooth_tail <- function(x, cutoff) {
   k <- fit$k
   sigma <- fit$sigma
   if (is.finite(k)) {
-      p <- (seq_len(len) - 0.5) / len
-      qq <- qgpd(p, k, sigma) + exp_cutoff
-      tail <- log(qq)
+    p <- (seq_len(len) - 0.5) / len
+    qq <- qgpd(p, k, sigma) + exp_cutoff
+    tail <- log(qq)
   } else {
-      tail <- x
+    tail <- x
   }
   list(tail = tail, k = k)
 }
@@ -322,7 +328,8 @@ throw_tail_length_warnings <- function(tail_lengths) {
     if (length(tail_lengths) == 1) {
       warning(
         "Not enough tail samples to fit the generalized Pareto distribution.",
-        call. = FALSE, immediate. = TRUE
+        call. = FALSE,
+        immediate. = TRUE
       )
     } else {
       bad <- which(tail_len_bad)
@@ -332,7 +339,11 @@ throw_tail_length_warnings <- function(tail_lengths) {
         "in some or all columns of matrix of log importance ratios. ",
         "Skipping the following columns: ",
         paste(if (Nbad <= 10) bad else bad[1:10], collapse = ", "),
-        if (Nbad > 10) paste0(", ... [", Nbad - 10, " more not printed].\n") else "\n",
+        if (Nbad > 10) {
+          paste0(", ... [", Nbad - 10, " more not printed].\n")
+        } else {
+          "\n"
+        },
         call. = FALSE,
         immediate. = TRUE
       )
@@ -352,17 +363,21 @@ throw_tail_length_warnings <- function(tail_lengths) {
 #' * If `r_eff` is `NA` then `rep(1, len)` is returned.
 #' * If `r_eff` is a scalar then `rep(r_eff, len)` is returned.
 #' * If `r_eff` is not a scalar but the length is not `len` then an error is thrown.
-#' * If `r_eff` has length `len` but has `NA`s then an error is thrown.
+#' * If `r_eff` has length `len` but has `NA`s then `NA`s are filled in with `1`s.
 #'
 prepare_psis_r_eff <- function(r_eff, len) {
   if (isTRUE(is.null(r_eff) || all(is.na(r_eff)))) {
     r_eff <- rep(1, len)
   } else if (length(r_eff) == 1) {
     r_eff <- rep(r_eff, len)
   } else if (length(r_eff) != len) {
-    stop("'r_eff' must have one value or one value per observation.", call. = FALSE)
+    stop(
+      "'r_eff' must have one value or one value per observation.",
+      call. = FALSE
+    )
   } else if (anyNA(r_eff)) {
-    stop("Can't mix NA and not NA values in 'r_eff'.", call. = FALSE)
+    message("Replacing NAs in `r_eff` with 1s")
+    r_eff[is.na(r_eff)] <- 1
   }
   r_eff
 }
@@ -390,4 +405,3 @@ throw_psis_r_eff_warning <- function() {
     call. = FALSE
   )
 }
-

README.md

@@ -13,21 +13,27 @@ __loo__ is an R package that allows users to compute efficient approximate
 leave-one-out cross-validation for fitted Bayesian models, as well as model
 weights that can be used to average predictive distributions. 
 The __loo__ package package implements the fast and stable computations for 
-approximate LOO-CV and WAIC from
+approximate LOO-CV
 
 * Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model 
 evaluation using leave-one-out cross-validation and WAIC. 
-_Statistics and Computing_. 27(5), 1413--1432. 
-doi:10.1007/s11222-016-9696-4. [Online](https://link.springer.com/article/10.1007/s11222-016-9696-4), 
-[arXiv preprint arXiv:1507.04544](https://arxiv.org/abs/1507.04544).
+_Statistics and Computing_. 27(5): 1413-1432. 
+[Journal](https://dx.doi.org/10.1007/s11222-016-9696-4), 
+[arXiv preprint](https://arxiv.org/abs/1507.04544)
+
+* Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2024).
+Pareto smoothed importance sampling. *Journal of Machine Learning Research*,
+25(72): 1-58. 
+[Journal](https://jmlr.org/papers/v25/19-556.html),
+[arXiv preprint](https://arxiv.org/abs/1507.02646)
 
 and computes model weights as described in
 
 * Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018). Using
-stacking to average Bayesian predictive distributions. In Bayesian
-Analysis, doi:10.1214/17-BA1091. 
-[Online](https://projecteuclid.org/euclid.ba/1516093227),
-[arXiv preprint arXiv:1704.02030](https://arxiv.org/abs/1704.02030).
+stacking to average Bayesian predictive distributions. *Bayesian Analysis* 
+13(3): 917-1007. 
+[Journal](https://dx.doi.org/10.1214/17-BA1091),
+[arXiv preprint](https://arxiv.org/abs/1704.02030)
 
 From existing posterior simulation draws, we compute approximate LOO-CV using
 Pareto smoothed importance sampling (PSIS), a new procedure for regularizing

inst/CITATION

@@ -26,6 +26,18 @@ bibentry(bibtype = "Article",
          header = "To cite the loo paper:"
 )
 
+bibentry(bibtype = "Article",
+         title = "Uncertainty in Bayesian leave-one-out cross-validation based model comparison",
+         author = c(person("Tuomas", "Sivula"),
+                    person("Måns", "Magnusson"),
+                    person("Asael Alonzo", "Matamoros"),
+    		    person("Aki", "Vehtari")),
+         journal = "Bayesian Analysis",
+         year = "2025",
+         note = "accepted for publication",
+         header = "To cite when using loo_compare():"
+)
+
 bibentry(bibtype = "Article",
          title = "Using stacking to average Bayesian predictive distributions",
          author = c(person("Yuling", "Yao"),
@@ -40,3 +52,71 @@ bibentry(bibtype = "Article",
          doi = "10.1214/17-BA1091",
          header = "To cite the stacking paper:"
 )
+
+bibentry(
+  title = "Pareto smoothed importance sampling",
+  bibtype = "Article",
+  author = c(
+    person("Aki", "Vehtari"),
+    person("Daniel", "Simpson"),
+    person("Andrew", "Gelman"),
+    person("Yuling", "Yao"),
+    person("Jonah", "Gabry")
+  ),
+  journal = "Journal of Machine Learning Research",
+  year = 2024,
+  volume = 25,
+  number = 72,
+  pages = "1-58",
+  header = "To cite Pareto-k diagnostics:"
+)
+
+bibentry(
+  bibtype = "Article",
+  author = c(
+    person(given = "Topi", family = "Paananen"),
+    person(given = "Juho", family = "Piironen"),
+    person(given = "Paul-Christian", family = "Buerkner"),
+    person(given = "Aki", family = "Vehtari")
+  ),
+  title = "Implicitly adaptive importance sampling",
+  journal = "Statistics and Computing",
+  volume = 31,
+  pages = "16",
+  year = 2021,
+  header = "To cite moment matching:"
+)
+
+bibentry(
+  bibtype = "InProceedings",
+  author = c(
+    person(given = "Måns", family = "Magnusson"),
+    person(given = "Michael Riis", family = "Andersen"),
+    person(given = "Johan", family = "Jonasson"),
+    person(given = "Aki", family = "Vehtari")
+  ),
+  title = "Leave-One-Out Cross-Validation for Large Data",
+  booktitle = "Thirty-sixth International Conference on Machine Learning",
+  publisher = "PMLR",
+  volume = "97",
+  pages = "4244-4253",
+  year = 2019,
+  header = "To cite subsampling loo:"
+)
+
+bibentry(
+  bibtype = "InProceedings",
+  author = c(
+    person(given = "Måns", family = "Magnusson"),
+    person(given = "Michael Riis", family = "Andersen"),
+    person(given = "Johan", family = "Jonasson"),
+    person(given = "Aki", family = "Vehtari")
+  ),
+  title = "Leave-One-Out Cross-Validation for Model Comparison in Large Data",
+  booktitle = "Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS)",
+  publisher = "PMLR",
+  volume = "108",
+  pages = "341-351",
+  year = 2019,
+  header = "To cite subsampling loo:"
+)

man-roxygen/loo-and-compare-references.R

@@ -10,10 +10,12 @@
 #' 25(72):1-58.
 #' [PDF](https://jmlr.org/papers/v25/19-556.html)
 #'
-#' Sivula, T, Magnusson, M., Matamoros A. A., and Vehtari, A. (2022).
+#' Sivula, T, Magnusson, M., Matamoros A. A., and Vehtari, A. (2025).
 #' Uncertainty in Bayesian leave-one-out cross-validation based model
-#' comparison. [preprint arXiv:2008.10296v3.](https://arxiv.org/abs/2008.10296v3).
+#' comparison. *Bayesian Analysis*, accepted for publication.
+#' [preprint arXiv:2008.10296v5.](https://arxiv.org/abs/2008.10296v5).
 #' 
-#' McLatchie, Y., and Vehtari, A. (2023).  Efficient estimation and
+#' McLatchie, Y., and Vehtari, A. (2024).  Efficient estimation and
 #' correction of selection-induced bias with order statistics.
-#' [preprint arXiv:2309.03742](https://arxiv.org/abs/2309.03742)
+#' *Statistics and Computing*. 34(132). 
+#' [doi:10.1007/s11222-024-10442-4](https://doi.org/10.1007/s11222-024-10442-4)

man-roxygen/loo-uncertainty-reference.R

@@ -1,4 +1,5 @@
-#' @references Sivula, T, Magnusson, M., Matamoros A. A., and Vehtari,
-#'   A. (2022).  Uncertainty in Bayesian leave-one-out
-#'   cross-validation based model comparison. [preprint
-#'   arXiv:2008.10296v3.](https://arxiv.org/abs/2008.10296v3).
+#' @references
+#' Sivula, T, Magnusson, M., Matamoros A. A., and Vehtari, A. (2025).
+#' Uncertainty in Bayesian leave-one-out cross-validation based model
+#' comparison. *Bayesian Analysis*, accepted for publication.
+#' [preprint arXiv:2008.10296v5.](https://arxiv.org/abs/2008.10296v5).