From 1b046ebb0536c8b1a0ce23ffbfb993ad508222b7 Mon Sep 17 00:00:00 2001
From: florence-bockting <florence.bockting@tu-dortmund.de>
Date: Fri, 27 Feb 2026 08:02:38 +0200
Subject: [PATCH 1/5] switch 'type' to 'output'; add wrapper function

---
 R/posterior_predict.R | 79 ++++++++++++++++++++++++-------------------
 1 file changed, 45 insertions(+), 34 deletions(-)

diff --git a/R/posterior_predict.R b/R/posterior_predict.R
index cd948f94c..483734201 100644
--- a/R/posterior_predict.R
+++ b/R/posterior_predict.R
@@ -82,7 +82,7 @@
 posterior_predict.brmsfit <- function(
   object, newdata = NULL, re_formula = NULL, re.form = NULL,
   transform = NULL, resp = NULL, negative_rt = FALSE,
-  ndraws = NULL, draw_ids = NULL, sort = FALSE, ntrys = 5, type = "r",
+  ndraws = NULL, draw_ids = NULL, sort = FALSE, ntrys = 5, output = "random",
   cores = NULL, ...
 ) {
   cl <- match.call()
@@ -93,11 +93,11 @@ posterior_predict.brmsfit <- function(
   object <- restructure(object)
   prep <- prepare_predictions(
     object, newdata = newdata, re_formula = re_formula, resp = resp,
-    ndraws = ndraws, draw_ids = draw_ids, check_response = FALSE, type = type, ...
+    ndraws = ndraws, draw_ids = draw_ids, check_response = FALSE, ...
   )
   posterior_predict(
     prep, transform = transform, sort = sort, ntrys = ntrys,
-    negative_rt = negative_rt, cores = cores, summary = FALSE, type = type
+    negative_rt = negative_rt, cores = cores, summary = FALSE, output = output
   )
 }
 
@@ -119,7 +119,7 @@ posterior_predict.mvbrmsprep <- function(object, ...) {
 posterior_predict.brmsprep <- function(object, transform = NULL, sort = FALSE,
                                        summary = FALSE, robust = FALSE,
                                        probs = c(0.025, 0.975),
-                                       cores = NULL, type = "r", ...) {
+                                       cores = NULL, ...) {
   summary <- as_one_logical(summary)
   cores <- validate_cores_post_processing(cores)
   if (is.customfamily(object$family)) {
@@ -136,7 +136,7 @@ posterior_predict.brmsprep <- function(object, transform = NULL, sort = FALSE,
   pp_fun <- paste0("posterior_predict_", object$family$fun)
   pp_fun <- get(pp_fun, asNamespace("brms"))
   N <- choose_N(object)
-  out <- plapply(seq_len(N), pp_fun, .cores = cores, prep = object, type = type, ...)
+  out <- plapply(seq_len(N), pp_fun, .cores = cores, prep = object, output = output, ...)
   if (grepl("_mv$", object$family$fun)) {
     out <- do_call(abind, c(out, along = 3))
     out <- aperm(out, perm = c(1, 3, 2))
@@ -301,52 +301,63 @@ validate_pp_method <- function(method) {
   method
 }
 
+# Helper function to predict continuous distributions
+# @param output "probability" or "random"
+# @param prep A named list returned by prepare_predictions containing
+#   all required data and posterior draws
+# @param i index of the observation for which to compute pp values
+# @param dist name of the distribution
+# @param ntrys number of trys in rejection sampling for truncated models
+# @param ... additional arguments passed to the distribution functions
+# @return a vector of draws from the distribution
+.predict_continuous_helper <- function(output, prep, i, dist, ntrys, ...) {
+  lb <- prep$data$lb[i]
+  ub <- prep$data$ub[i]
+  
+  switch(output,
+    "probability" = {
+      q <- prep$data$Y[i]
+      pcontinuous(
+        q = q, dist = dist, lb = lb, ub = ub, ntrys = ntrys,
+        ndraws = prep$ndraws, ...
+      )
+    },
+    "random" = {
+      rcontinuous(
+        n = prep$ndraws, dist = dist, lb = lb, ub = ub, ntrys = ntrys, ...
+      )
+    }
+  )
+}
+
 # ------------------- family specific posterior_predict methods ---------------------
 # All posterior_predict_<family> functions have the same arguments structure
-# @param i index of the observatio for which to compute pp values
+# @param i index of the observation for which to compute pp values
 # @param prep A named list returned by prepare_predictions containing
 #   all required data and posterior draws
 # @param ... ignored arguments
 # @param A vector of length prep$ndraws containing draws
 #   from the posterior predictive distribution
-posterior_predict_gaussian <- function(i, prep, ntrys = 5, type = "r", ...) {
+posterior_predict_gaussian <- function(i, prep, ntrys = 5, output = "random", ...) {
   mu <- get_dpar(prep, "mu", i = i)
   sigma <- get_dpar(prep, "sigma", i = i)
   sigma <- add_sigma_se(sigma, prep, i = i)
-  switch(type,
-         r = rcontinuous(
-           n = prep$ndraws, dist = "norm",
-           mean = mu, sd = sigma,
-           lb = prep$data$lb[i], ub = prep$data$ub[i],
-           ntrys = ntrys
-         ),
-         p = pcontinuous(
-           n = prep$ndraws, dist = "norm",
-           q = prep$data$Y[i], mean = mu, sd = sigma,
-           lb = prep$data$lb[i], ub = prep$data$ub[i],
-           ntrys = ntrys
-         )
+
+  .predict_continuous_helper(
+    output = output, prep = prep, i = i, ntrys = ntrys,
+    dist = "norm", mean = mu, sd = sigma
   )
 }
 
-posterior_predict_student <- function(i, prep, ntrys = 5, type = "r", ...) {
+posterior_predict_student <- function(i, prep, ntrys = 5, output = "random", ...) {
   nu <- get_dpar(prep, "nu", i = i)
   mu <- get_dpar(prep, "mu", i = i)
   sigma <- get_dpar(prep, "sigma", i = i)
   sigma <- add_sigma_se(sigma, prep, i = i)
-  switch(type,
-         r =   rcontinuous(
-           n = prep$ndraws, dist = "student_t",
-           df = nu, mu = mu, sigma = sigma,
-           lb = prep$data$lb[i], ub = prep$data$ub[i],
-           ntrys = ntrys
-         ),
-         p = pcontinuous(
-           n = prep$ndraws, dist = "student_t",
-           q = prep$data$Y[i], df = nu, mu = mu, sigma = sigma,
-           lb = prep$data$lb[i], ub = prep$data$ub[i],
-           ntrys = ntrys
-         )
+  
+  .predict_continuous_helper(
+    output = output, prep = prep, i = i, ntrys = ntrys,
+    dist = "student_t", df = nu, mu = mu, sigma = sigma
   )
 }
 

From 51026133530a015be8f23420855a4e042707de1b Mon Sep 17 00:00:00 2001
From: florence-bockting <florence.bockting@tu-dortmund.de>
Date: Fri, 27 Feb 2026 09:46:09 +0200
Subject: [PATCH 2/5] add test for posterior_predict with output arg

---
 tests/testthat/tests.posterior_predict.R | 26 ++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/tests/testthat/tests.posterior_predict.R b/tests/testthat/tests.posterior_predict.R
index 533797b15..3e3c9c9d3 100644
--- a/tests/testthat/tests.posterior_predict.R
+++ b/tests/testthat/tests.posterior_predict.R
@@ -439,3 +439,29 @@ test_that("posterior_predict_custom runs without errors", {
   }
   expect_equal(length(brms:::posterior_predict_custom(sample(1:nobs, 1), prep)), ns)
 })
+
+test_that("posterior_predict for location shift models runs with 'output' argument without error", {
+  ns <- 30
+  nobs <- 10
+  prep <- structure(list(ndraws = ns), class = "brmsprep")
+  prep$dpars <- list(
+    mu = matrix(rnorm(ns * nobs), ncol = nobs),
+    sigma = rchisq(ns, 3), nu = rgamma(ns, 4)
+  )
+  prep$data <- list(Y  = rpred)
+  i <- sample(nobs, 1)
+
+  pnorm(prep$data$Y[i])
+
+  # probability
+  rpred <- brms:::posterior_predict_gaussian(i, prep = prep, output = "random")
+  expect_equal(length(rpred), ns)
+
+  qpred <- brms:::posterior_predict_gaussian(i, prep = prep, output = "probability")
+  
+  expect_equal(length(qpred), ns)
+
+  pred <- brms:::posterior_predict_student(i, prep = prep)
+  expect_equal(length(pred), ns)
+})
+

From bbfd28cf73be6f2d9d74f2351efaf83e2040d19b Mon Sep 17 00:00:00 2001
From: florence-bockting <florence.bockting@tu-dortmund.de>
Date: Fri, 27 Feb 2026 13:02:19 +0200
Subject: [PATCH 3/5] fix 'probability' method for posterior_predict

---
 R/posterior_predict.R | 68 ++++++++++++++++++++++---------------------
 1 file changed, 35 insertions(+), 33 deletions(-)

diff --git a/R/posterior_predict.R b/R/posterior_predict.R
index 483734201..a7c9c522a 100644
--- a/R/posterior_predict.R
+++ b/R/posterior_predict.R
@@ -301,35 +301,6 @@ validate_pp_method <- function(method) {
   method
 }
 
-# Helper function to predict continuous distributions
-# @param output "probability" or "random"
-# @param prep A named list returned by prepare_predictions containing
-#   all required data and posterior draws
-# @param i index of the observation for which to compute pp values
-# @param dist name of the distribution
-# @param ntrys number of trys in rejection sampling for truncated models
-# @param ... additional arguments passed to the distribution functions
-# @return a vector of draws from the distribution
-.predict_continuous_helper <- function(output, prep, i, dist, ntrys, ...) {
-  lb <- prep$data$lb[i]
-  ub <- prep$data$ub[i]
-  
-  switch(output,
-    "probability" = {
-      q <- prep$data$Y[i]
-      pcontinuous(
-        q = q, dist = dist, lb = lb, ub = ub, ntrys = ntrys,
-        ndraws = prep$ndraws, ...
-      )
-    },
-    "random" = {
-      rcontinuous(
-        n = prep$ndraws, dist = dist, lb = lb, ub = ub, ntrys = ntrys, ...
-      )
-    }
-  )
-}
-
 # ------------------- family specific posterior_predict methods ---------------------
 # All posterior_predict_<family> functions have the same arguments structure
 # @param i index of the observation for which to compute pp values
@@ -345,7 +316,7 @@ posterior_predict_gaussian <- function(i, prep, ntrys = 5, output = "random", ..
 
   .predict_continuous_helper(
     output = output, prep = prep, i = i, ntrys = ntrys,
-    dist = "norm", mean = mu, sd = sigma
+    dist = "norm", mean = mu, sd = sigma, ...
   )
 }
 
@@ -357,7 +328,7 @@ posterior_predict_student <- function(i, prep, ntrys = 5, output = "random", ...
   
   .predict_continuous_helper(
     output = output, prep = prep, i = i, ntrys = ntrys,
-    dist = "student_t", df = nu, mu = mu, sigma = sigma
+    dist = "student_t", df = nu, mu = mu, sigma = sigma, ...
   )
 }
 
@@ -1029,12 +1000,12 @@ rcontinuous <- function(n, dist, ..., lb = NULL, ub = NULL, ntrys = 5) {
   out
 }
 
-pcontinuous <- function(n, dist, ..., lb = NULL, ub = NULL, ntrys = 5) {
+pcontinuous <- function(q, dist, ..., lb = NULL, ub = NULL, ntrys = 5) {
   args <- list(...)
   if (is.null(lb) && is.null(ub)) {
     # sample as usual
     pdist <- paste0("p", dist)
-    out <- do_call(pdist, c(list(n), args))
+    out <- do_call(pdist, c(list(q), args))
   } else {
     error("not implemented yet")
     # sample from truncated distribution
@@ -1140,3 +1111,34 @@ check_discrete_trunc_bounds <- function(x, lb = NULL, ub = NULL, thres = 0.01) {
   }
   round(x)
 }
+
+# predict random numbers or probability values from continuous distributions
+# @param output "probability" or "random"
+# @param prep A named list returned by prepare_predictions containing
+#   all required data and posterior draws
+# @param i index of the observation for which to compute pp values
+# @param dist name of the distribution
+# @param ntrys number of trys in rejection sampling for truncated models
+# @param q optional custom quantile value; if NULL, the default is prep$data$Y[i] 
+# @param ... additional arguments passed to the distribution functions
+# @return a vector of draws
+.predict_continuous_helper <- function(output, prep, i, dist, ntrys, q = NULL, ...) {
+  lb <- prep$data$lb[i]
+  ub <- prep$data$ub[i]
+  
+  switch(output,
+    "probability" = {
+      if (is.null(q)) {
+        q <- prep$data$Y[i]
+      }
+      pcontinuous(
+        q = prep$data$Y[i], dist = dist, lb = lb, ub = ub, ntrys = ntrys, ...
+      )
+    },
+    "random" = {
+      rcontinuous(
+        n = prep$ndraws, dist = dist, lb = lb, ub = ub, ntrys = ntrys, ...
+      )
+    }
+  )
+}

From 780ad34a8fdaf07db11fe7152a1f3f753e73542b Mon Sep 17 00:00:00 2001
From: florence-bockting <florence.bockting@tu-dortmund.de>
Date: Fri, 27 Feb 2026 13:02:40 +0200
Subject: [PATCH 4/5] add test for posterior_predict_gaussian

---
 tests/testthat/tests.posterior_predict.R | 33 +++++++++++-------------
 1 file changed, 15 insertions(+), 18 deletions(-)

diff --git a/tests/testthat/tests.posterior_predict.R b/tests/testthat/tests.posterior_predict.R
index 3e3c9c9d3..9e2798647 100644
--- a/tests/testthat/tests.posterior_predict.R
+++ b/tests/testthat/tests.posterior_predict.R
@@ -440,28 +440,25 @@ test_that("posterior_predict_custom runs without errors", {
   expect_equal(length(brms:::posterior_predict_custom(sample(1:nobs, 1), prep)), ns)
 })
 
-test_that("posterior_predict for location shift models runs with 'output' argument without error", {
-  ns <- 30
-  nobs <- 10
-  prep <- structure(list(ndraws = ns), class = "brmsprep")
-  prep$dpars <- list(
-    mu = matrix(rnorm(ns * nobs), ncol = nobs),
-    sigma = rchisq(ns, 3), nu = rgamma(ns, 4)
-  )
-  prep$data <- list(Y  = rpred)
-  i <- sample(nobs, 1)
-
-  pnorm(prep$data$Y[i])
+test_that("posterior_predict_gaussian runs with various 'output' values without error", {
+  fit <- rename_pars(brms:::brmsfit_example3)
+  prep <- brms::prepare_predictions(fit)
+  model_fit <- fit$fit@sim
+  S <- model_fit$chains * (model_fit$iter - model_fit$warmup)
+  i <- 1
 
   # probability
   rpred <- brms:::posterior_predict_gaussian(i, prep = prep, output = "random")
-  expect_equal(length(rpred), ns)
+  expect_equal(length(rpred), S)
 
-  qpred <- brms:::posterior_predict_gaussian(i, prep = prep, output = "probability")
+  # compute PIT values (q = prep$data$Y[i])
+  PITs <- brms:::posterior_predict_gaussian(i, prep = prep, output = "probability")
+  expect_equal(length(PITs), S)
+  expect_true(all(PITs >= 0 & PITs <= 1))
   
-  expect_equal(length(qpred), ns)
-
-  pred <- brms:::posterior_predict_student(i, prep = prep)
-  expect_equal(length(pred), ns)
+  # compute cdf based on custom 'q'
+  qpred <- brms:::posterior_predict_gaussian(i, q = 15, prep = prep, output = "probability")
+  expect_equal(length(qpred), S)
+  expect_true(all(qpred >= 0 & qpred <= 1))
 })
 

From 3106b93188979dcafe78267d3ad210c4ad84db91 Mon Sep 17 00:00:00 2001
From: florence-bockting <florence.bockting@tu-dortmund.de>
Date: Fri, 27 Feb 2026 14:39:57 +0200
Subject: [PATCH 5/5] fix setting of q

---
 R/posterior_predict.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/posterior_predict.R b/R/posterior_predict.R
index a7c9c522a..e7de15c0f 100644
--- a/R/posterior_predict.R
+++ b/R/posterior_predict.R
@@ -1132,7 +1132,7 @@ check_discrete_trunc_bounds <- function(x, lb = NULL, ub = NULL, thres = 0.01) {
         q <- prep$data$Y[i]
       }
       pcontinuous(
-        q = prep$data$Y[i], dist = dist, lb = lb, ub = ub, ntrys = ntrys, ...
+        q = q, dist = dist, lb = lb, ub = ub, ntrys = ntrys, ...
       )
     },
     "random" = {