Vectorized Bayesian bootstrap

R/loo_model_weights.R

@@ -329,22 +329,16 @@ pseudobma_weights <-
       return(wts)
     }
 
-    temp <- matrix(NA, BB_n, K)
     BB_weighting <- dirichlet_rng(BB_n, rep(alpha, N))
-    for (bb in 1:BB_n) {
-      z_bb <- BB_weighting[bb, ] %*% lpd_point * N
-      uwts <- exp(z_bb - max(z_bb))
-      temp[bb, ] <- uwts / sum(uwts)
-    }
-    wts <- structure(
-      colMeans(temp),
+    z <- BB_weighting %*% lpd_point * N
+    uwts <- exp(z - matrixStats::rowMaxs(z))
+    structure(
+      colMeans(uwts / rowSums(uwts)),
       names = paste0("model", 1:K),
       class = "pseudobma_bb_weights"
     )
-    return(wts)
   }
 
-
 #' Generate dirichlet simulations, rewritten version
 #' @importFrom stats rgamma
 #' @noRd

tests/testthat/_snaps/model_weighting.md

@@ -42,3 +42,7 @@
       model2 1.000 
       model3 0.000 
 
+# Bayesian bootstrap gives expected result
+
+    c(0.188359351057998, 0.309962881543203, 0.501677767398798)
+

tests/testthat/test_model_weighting.R

@@ -111,6 +111,43 @@ test_that("loo_model_weights (stacking and pseudo-BMA) gives expected result", {
   expect_identical(w3, w3_b)
 })
 
+test_that("Bayesian bootstrap gives expected result", {
+  lpd_point <- matrix(
+    c(
+      -0.2,
+      -0.8,
+      -1.1,
+      -1.4,
+      -0.3,
+      -0.5,
+      -0.6,
+      -0.7,
+      -0.1,
+      -1.0,
+      -1.2,
+      -0.4,
+      -0.9,
+      -0.5,
+      -0.8
+    ),
+    ncol = 3,
+    byrow = TRUE
+  )
+  BB_n <- 25
+  alpha <- 0.7
+
+  set.seed(0)
+  weights <- pseudobma_weights(lpd_point, BB = TRUE, BB_n = BB_n, alpha = alpha)
+
+  expect_s3_class(weights, "pseudobma_bb_weights")
+  expect_named(weights, paste0("model", 1:3))
+  expect_equal(sum(weights), 1)
+  expect_snapshot_value(
+    unname(as.numeric(weights)),
+    style = "deparse"
+  )
+})
+
 test_that("stacking_weights and pseudobma_weights throw correct errors", {
   xx <- cbind(rnorm(10))
   expect_error(stacking_weights(xx), "two models are required")