proof of concept for posterior_pit

[avehtari]

Summary

Fixes #1855

• Extend posterior_predict() beyond random draws (PIT / CDF / density / quantile) for calibration and predictive checks.
• More detailed description of this PR and its status can be found in notes/pr-summary.md

Current API

• new argument output with values: "random", "probability", "pit", "density", "quantile" (default is "random")
• additional support of q / p / lower.tail / log.p / log.

# example usage
pits <- posterior_predict(fit, output = "pit")
ppc_pit_ecdf(pit = colMeans(pits), method = "correlated")

What’s implemented now

• Shared helpers; broad family support (see for overview: proof of concept for posterior_pit #1857 (comment) and proof of concept for posterior_pit #1857 (comment))
  • R/posterior_predict.R: output API, continuous/discrete helpers, per-family predict methods
  • R/distributions.R: New/extended d/p/q/r helpers (hurdle, ZI, ordinal/categorical, numerical quantiles, etc.)
• Tests in tests/testthat
• Vignettes (vignettes/brms_posterior_predict.Rmdand vignettes/brms_distribution_output_cheatsheet.Rmd; see for graphical overview here: proof of concept for posterior_pit #1857 (comment))

Known limitations

• Some families still "random" only (cox none; custom only if implemented)
• Discrete truncated "random": rejection sampling

TODOs

• validity checks
• test code in different case-studies
• remove "note" folder before merging
• update vignettes and decide what we want to do with them
• update News.md

Original PR text (click to expand)

This is a proof of concept branch and far from ready to be merged, but easier to show a possible way to solve issue #1855

To avoid duplication of code in brms, the idea is to use posterior_predict(), but add a new argument type with default value "r" referring to random draws. New type="p" would instead return CDFs / PITs.

PITs = posterior_predict(fit, type = "p")
PIT_post = colMeans(PITs)
PIT_loo = E_loo(PITs, loo(fit)$psis_object)$value

These PIT values are more useful than the ones based on ranks.

The proof of concept works with families "gaussian" and "student_t" without truncation. If the general idea is acceptable then support for the rest of the families and truncation can be added.

Adding here that natural other types are "d" for posterior predictive densities

ppd = posterior_predict(fit, type = "d")

and
"q" for quantiles that are more accurate than current rank based quantiles (used e.g. by predictive_interval())

ppqs = posterior_predict(fit, type = "q", probs = c(0.05, 0.95))
ppq = apply(ppqs, 1, mean)

[paul-buerkner]

Thanks! I like the idea. Please find some comments to code structure and naming in my review.

[paul-buerkner]

Yes this is exactly what I meant! Florence Bockting ***@***.***> schrieb am Di. 24. Feb. 2026 um 20:06:

…

***@***.**** commented on this pull request. ------------------------------ In R/posterior_predict.R <#1857 (comment)>: > mu <- get_dpar(prep, "mu", i = i) sigma <- get_dpar(prep, "sigma", i = i) sigma <- add_sigma_se(sigma, prep, i = i) - rcontinuous( - n = prep$ndraws, dist = "norm", - mean = mu, sd = sigma, - lb = prep$data$lb[i], ub = prep$data$ub[i], - ntrys = ntrys + switch(type, How about something like the following? We have a helper function that prepares the switching and can be called within the corresponding distributions. In the following a (non-polished) example for a helper function and its function call in the two example distributions (gaussian and student). (Note I changed the argument type to output in the following code... just to get an idea how it would look like) .predict_continuous_helper <- function(output, prep, i, dist, ntrys, ...) { lb <- prep$data$lb[i] ub <- prep$data$ub[i] if (output == "probability") { q <- prep$data$Y[i] return(pcontinuous( q = q, dist = dist, lb = lb, ub = ub, ntrys = ntrys, ndraws = prep$ndraws, ... )) } else if (output == "random") { return(rcontinuous( n = prep$ndraws, dist = dist, lb = lb, ub = ub, ntrys = ntrys, ... )) } } 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) .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, 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) .predict_continuous_helper( output = output, prep = prep, i = i, ntrys = ntrys, dist = "student_t", df = nu, mu = mu, sigma = sigma ) } — Reply to this email directly, view it on GitHub <#1857 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCW2AAA5QELHJGMCXTWVBL4NSOLVAVCNFSM6AAAAACV2RC4P6VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZTQNBZHE4TMMJZGE> . You are receiving this because you commented.Message ID: ***@***.***>

[florence-bockting]

I have now updated the proof of concept for 4 functions:

• posterior_predict_gaussian()
• posterior_predict_student()
• posterior_predict_binomial()
• posterior_predict_poisson()

They include two additional new arguments:

• output (which can take currently the values random or probability) and
• randomized (which can be NULL or a boolean).

The current implementation follows the idea that we can call:

• rpred <- posterior_predict(fit_norm)
  • computes random draws from the normal predictive
  • output = "random" and randomized = NULL by default
• ppred <- posterior_predict(fit_norm, output = "probability")
  • computes cdf of the normal predictive
  • randomized = NULL by default
• rpred <- posterior_predict(fit_pois)
  • computes random draws from the poisson predictive
  • output = "random" and randomized = NULL by default
• ppred <- posterior_predict(fit_pois, output = "probability")
  • computes randomized PITs from the poisson predictive
  • randomized = TRUE by default
• ppred <- posterior_predict(fit_pois, output = "probability", randomized = FALSE)
  • computes the non-randomized PITs from the poisson predictive, which are currently defined as the mid points (F(x) +F(x-1)) * 0.5

[avehtari]

computes the non-randomized PITs from the poisson predictive, which are currently defined as the mid points (F(x) +F(x-1)) * 0.5

I think these should be just F(x) to follow what ppois etc. do. For the same reason, I think that it might be more logical to have non-randomized as the default

[avehtari]

Or we could have separately output options probability (not randomized F(x)) and pit randomized for discrete

[florence-bockting]

Or we could have separately output options probability (not randomized F(x)) and pit randomized for discrete

Okay, one option could be to have for outcome the values random, probability, density, quantile, and pit. Whereby pit is only a valid option for discrete cases. Cons are that there is a conceptual overlap between pit and probability and pit is only a valid choice for a subclass of inputs. Pro is that we have all alternatives in one function and no additional argument (like randomized).

An alternative is to outsource the pit case entirely and to have an additional function such as posterior_pit(fit_pois) that applies only to discrete cases and additionally, we have posterior_predict(fit_pois, outcome = "probability") supporting the default implementations like ppois(). In my opinion, this option is from a design point relatively clean, as we don't have the problem that there are "cases" which work/don't work. However, we would have to introduce an additional function.

[paul-buerkner]

I would argue putting everything in one function. The code overhead is otherwise way too high. I like the probability + randomized argument approach. The only question is what the default for randomized is? Aki can you elaborate which use cases you had for much choice of randomized? Florence Bockting ***@***.***> schrieb am Mo. 16. März 2026 um 07:50:

…

*florence-bockting* left a comment (paul-buerkner/brms#1857) <#1857 (comment)> Or we could have separately output options probability (not randomized F(x)) and pit randomized for discrete Okay, one option could be to have for outcome the values random, probability, density, quantile, and pit. Whereby pit is only a valid option for discrete cases. Cons are that there is a conceptual overlap between pit and probability and pit is only a valid choice for a subclass of inputs. Pro is that we have all alternatives in one function and no additional argument (like randomized). An alternative is to outsource the pit case entirely and to have an additional function such as posterior_pit(fit_pois) that applies only to discrete cases and additionally, we have posterior_predict(fit_pois, outcome = "probability") supporting the default implementations like ppois(). In my opinion, this option is from a design point relatively clean, as we don't have the problem that there are "cases" which work/don't work. However, we would have to introduce an additional function. — Reply to this email directly, view it on GitHub <#1857 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCW2ACTEEEYGBYRD56TRK34Q6P4VAVCNFSM6AAAAACV2RC4P6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DANRVGQ3DMNZVHA> . You are receiving this because you commented.Message ID: ***@***.***>

[avehtari]

Okay, one option could be to have for outcome the values random, probability, density, quantile, and pit. Whereby pit is only a valid option for discrete cases

In this case pit for continuous would be valid, but the result would be the same as with probability

• PITs (with randomization for discrete) are needed fot PIT and LOO-PIT uniformity checks (pp_check(..., type = "pit_ecdf|loo_pit_ecdf"))
• Probabilities (non-randomized) are needed for calibration plots (reliabilitydiag type plots), for which we have PR in bayesplot PPC Calibration plots stan-dev/bayesplot#352
• Probabilities are useful in general for making predictions and modeling visualization
• If the outcome options are random, probability, density, quantile, and pit, we don't need separate argument for randomization, and pit computes such PIT values which are useful for uniformity tests both for continuous and discrete case

[paul-buerkner]

I am fine with option random, probability, density, quantile, and pit (or whatever subset we go with first), as long all all is realized in the same function (posterior_predict). Otherwise code overhead is too high. If we want, we can still have lightweight wrappers around posterior_predict with output argument fixed.

[florence-bockting]

I am fine with option random, probability, density, quantile, and pit (or whatever subset we go with first), as long all all is realized in the same function (posterior_predict).

Alright. We definitely stick with one function - the idea of an extra function is out.
I will change the implementation to having posterior_predict with output = (random, probability, density, quantile, pit). Whereby pit is the randomized pit for discrete distributions and the "normal" cdf for continuous distributions.

[florence-bockting]

Overview of current output support in R/posterior_predict.R (click to expand)

Legend: Y - yes,
column* - newly implemented

method	random	probability*	density*	pit*	quantile*	notes
gaussian	Y	Y	Y	Y	Y
student	Y	Y	Y	Y	Y
lognormal	Y	Y	Y	Y	Y
shifted_lognormal	Y	Y	Y	Y	Y
skew_normal	Y	Y	Y	Y	Y
gaussian_mv	Y					no support yet
student_mv	Y					no support yet
gaussian_time	Y					no support yet
student_time	Y					no support yet
gaussian_lagsar	Y					no support yet
student_lagsar	Y					no support yet
gaussian_errorsar	Y					no support yet
student_errorsar	Y					no support yet
gaussian_fcor	Y					no support yet
student_fcor	Y					no support yet
binomial	Y	Y	Y	Y	Y
beta_binomial	Y	Y	Y	Y	Y
bernoulli	Y	Y	Y	Y	Y
poisson	Y	Y	Y	Y	Y
negbinomial	Y	Y	Y	Y	Y
negbinomial2	Y	Y	Y	Y	Y
geometric	Y	Y	Y	Y	Y
discrete_weibull	Y	Y	Y	Y	Y
com_poisson	Y	Y	Y	Y	Y
exponential	Y	Y	Y	Y	Y
gamma	Y	Y	Y	Y	Y
weibull	Y	Y	Y	Y	Y
frechet	Y	Y	Y	Y	Y
gen_extreme_value	Y	Y	Y	Y	Y
inverse.gaussian	Y	Y	Y	Y	Y
exgaussian	Y	Y	Y	Y	Y
wiener	Y
beta	Y	Y	Y	Y	Y
xbeta	Y	Y	Y	Y	Y
von_mises	Y	Y	Y	Y	Y
asym_laplace	Y	Y	Y	Y	Y
zero_inflated_asym_laplace	Y	Y	Y	Y	Y
cox						explicitly unsupported
hurdle_poisson	Y	Y	Y	Y	Y	no truncation for "random" implemented
hurdle_negbinomial	Y	Y	Y	Y	Y	no truncation for "random" implemented
hurdle_gamma	Y	Y	Y	Y	Y	no truncation for "random" implemented
hurdle_lognormal	Y	Y	Y	Y	Y	no truncation for "random" implemented
hurdle_cumulative	Y	Y	Y	Y	Y
zero_inflated_beta	Y	Y	Y	Y	Y
zero_one_inflated_beta	Y	Y	Y	Y	Y
zero_inflated_poisson	Y	Y	Y	Y	Y
zero_inflated_negbinomial	Y	Y	Y	Y	Y
zero_inflated_binomial	Y	Y	Y	Y	Y
zero_inflated_beta_binomial	Y	Y	Y	Y	Y
categorical	Y	Y	Y	Y	Y
multinomial	Y					no support yet
dirichlet_multinomial	Y					no support yet
dirichlet	Y					no support yet
dirichlet2	Y					no support yet
logistic_normal	Y					no support yet
cumulative	Y	Y	Y	Y	Y
sratio	Y	Y	Y	Y	Y
cratio	Y	Y	Y	Y	Y
acat	Y	Y	Y	Y	Y
ordinal	Y	Y	Y	Y	Y
custom						delegated to custom family method
mixture	depends	depends	depends	depends	depends	delegated to component families

"no support yet" usually because of different output structure.

[florence-bockting]

Overview of current distribution helpers (R/distributions.R; click to expand)

Legend: Y - yes, Y* - yes, newly implemented

Family	ddist	pdist	qdist	rdist	notes
acat	Y
asym_laplace	Y	Y	Y	Y
beta_binomial	Y	Y	Y*	Y
categorical	Y	Y	Y*	Y*
com_poisson	Y	Y	Y*	Y	refactored rdist to reuse qdist code
cox	Y	Y
cratio	Y
cumulative	Y
dirichlet	Y			Y
dirichletmultinomial	Y
discrete_weibull	Y	Y	Y	Y
exgaussian	Y	Y	Y*	Y	check robustness of qdist implementation
frechet	Y	Y	Y	Y
gen_extreme_value	Y	Y	Y	Y
hurdle_gamma	Y	Y	Y*
hurdle_lognormal	Y	Y	Y*
hurdle_negbinomial	Y	Y	Y*
hurdle_poisson	Y	Y	Y*
inv_gaussian	Y	Y	Y*	Y	check robustness of qdist implementation
invgamma		Y
logistic_normal	Y			Y
multi_normal	Y			Y
multi_student_t	Y			Y
multinomial	Y
ordinal	Y*	Y	Y*	Y*
shifted	Y	Y	Y	Y
shifted_lnorm	Y	Y	Y	Y
skew_normal	Y	Y	Y	Y
sratio	Y
student_t	Y	Y	Y	Y
von_mises	Y	Y	Y*	Y
wiener	Y			Y
xbeta	Y	Y	Y	Y
zero_inflated_asym_laplace	Y	Y	Y*
zero_inflated_beta	Y	Y	Y*
zero_inflated_beta_binomial	Y	Y	Y*
zero_inflated_binomial	Y	Y	Y*
zero_inflated_negbinomial	Y	Y	Y*
zero_inflated_poisson	Y	Y	Y*

[florence-bockting]

Graphics from the knitted vignettes\brms_distribution_output_cheatsheet.Rmd:

output = "random" (click to expand)

output = "probability" (click to expand)

output = "pit" (click to expand)

output = "density" (click to expand)

output = "quantile" (click to expand)

[florence-bockting]

PR is currently failing because of some add_subsampling_vars_to_pointwise(pointwise = x$pointwise, idxs, elpd_loo_approx), which seems to be unrelated to my changes. Any idea @paul-buerkner where this comes from?

[paul-buerkner]

Issue in loo. Should be fixed soon Florence Bockting ***@***.***> schrieb am Mo. 20. Juli 2026 um 15:54:

…

*florence-bockting* left a comment (paul-buerkner/brms#1857) <#1857 (comment)> PR is currently failing because of some add_subsampling_vars_to_pointwise(pointwise = x$pointwise, idxs, elpd_loo_approx), which seems to be unrelated to my changes. Any idea @paul-buerkner <https://github.com/paul-buerkner> where this comes from? — Reply to this email directly, view it on GitHub <#1857?email_source=notifications&email_token=ADCW2AAXHE7A5XYHNYPF3PL5FYQCHA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBSGMYDGNBUGE32M4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-5023034417>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCW2ADXZCCN3TLEGR5BKTT5FYQCHAVCNFSNUABEKJSXA33TNF2G64TZHMZTONRQGY4DSNR3JFZXG5LFHMZTSNRYGY2DQMBTHGQXMAQ> . Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS <https://github.com/notifications/mobile/ios/ADCW2AG6OD22DZTZ4EE5KHD5FYQCHA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBSGMYDGNBUGE32M4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJKTGN5XXIZLSL5UW64Y> and Android <https://github.com/notifications/mobile/android/ADCW2ADJSPMTT4G56QHIAFT5FYQCHA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMBSGMYDGNBUGE32M4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLTGN5XXIZLSL5QW4ZDSN5UWI>. Download it today! You are receiving this because you were mentioned.Message ID: ***@***.***>

[florence-bockting]

@paul-buerkner When you have time, you can probably take a look at the PR.
I implemented not with the support of posterior_predict for most families.
See PR description for reference to relevant files.