Bayesian Workflow: From Priors to Posterior Predictive Checks

The difficult part of Bayesian analysis is not Bayes' theorem. It is building a model that deserves to be updated. A practical Bayesian workflow is iterative:

This article focuses on that workflow rather than repeating the basic prior-likelihood-posterior definitions.

Start from a generative model

Suppose observations satisfy

with prior

The Bayesian model defines a joint distribution

That joint model should be interpretable as a data-generating story. For example, a count model might be

Each component should have a reason to exist.

Prior predictive checks

Before seeing the observed outcome, simulate

then

These draws come from the prior predictive distribution

The question is:

Does the prior imply data that could plausibly exist?

If a prior for a logistic-regression coefficient routinely produces probabilities numerically indistinguishable from 0 or 1 for ordinary predictor values, the prior may be too wide. A prior predictive check detects that before the data push the posterior into a strange geometry.

Weakly informative does not mean vague

A prior can be weakly informative relative to the domain while still excluding absurd parameter regions. For standardized predictors, a normal prior

may already imply a wide range of odds ratios. For a scale parameter,

a half-normal or half-t prior may be preferable to a prior that places substantial mass on impossibly large scales. The appropriate scale depends on the scientific units. “Flat” is not a scientific unit.

Fit the model, then diagnose the computation

Posterior simulation is an algorithm. It can fail. For multiple MCMC chains, diagnostics include:

• rank-normalized split-$\hat R$;
• bulk effective sample size;
• tail effective sample size;
• divergent transitions;
• maximum tree depth for HMC/NUTS;
• energy diagnostics;
• trace plots.

A posterior summary is not trustworthy merely because the software completed.

Effective sample size

Autocorrelated draws contain less information than independent draws. For a scalar quantity,

If $N=10{,}000$ but strong autocorrelation reduces

to 200, posterior Monte Carlo error behaves more like a few hundred independent draws than ten thousand. Report enough computation to make Monte Carlo error small relative to posterior uncertainty.

Divergences are not cosmetic warnings

Hamiltonian Monte Carlo can report divergent transitions when numerical trajectories fail to explore the posterior geometry reliably. Common causes include:

• funnel-shaped hierarchical posteriors;
• weakly identified scale parameters;
• extreme posterior curvature;
• poor parameterization.

Increasing the target acceptance rate can help in some cases. Reparameterizing the model is often more important. A non-centered parameterization can transform a hierarchical model from difficult geometry into a tractable one.

Posterior predictive checks

After fitting, simulate replicated data

Compare the replicated data with the observed data on features that matter. Possible checks include:

• mean;
• variance;
• zero frequency;
• tail quantiles;
• maximum value;
• autocorrelation;
• group-level spread;
• calibration;
• spatial pattern.

A model can fit average behavior while failing badly in the tails. One posterior predictive plot cannot validate every aspect of the model.

A p-value analogue is not the goal

Posterior predictive checking should not be reduced to another thresholded hypothesis test. The purpose is model criticism. If simulated data fail systematically to reproduce a relevant feature, ask which assumption is responsible. A useful check points toward model revision.

Hierarchical shrinkage

Suppose group-specific parameters satisfy

Groups with little data are partially pooled toward the population mean. The amount of pooling is learned through

This is one of the central strengths of Bayesian multilevel modeling. It avoids the extremes of:

• fitting every group independently;
• forcing every group to share one identical parameter.

The uncertainty in group-level effects is propagated naturally.

Small datasets do not make Bayesian analysis automatically reliable

Priors can stabilize weak-data problems. That does not create information from nowhere. When the likelihood is weak, the posterior can be dominated by the prior. That may be appropriate if the prior is genuinely informed. It should be visible. Prior sensitivity is therefore particularly important in small samples.

Sensitivity analysis

Refit the model under several defensible prior specifications. Compare:

• posterior means;
• interval widths;
• tail probabilities;
• decisions;
• predictive performance.

If a conclusion changes materially under modest prior changes, report that fragility. The right conclusion may be

the data do not strongly identify this quantity.

That is more informative than hiding the prior dependence.

Model comparison

Bayesian model comparison can mean several different things.

Predictive comparison

Approximate out-of-sample predictive performance with methods such as:

• leave-one-out cross-validation;
• K-fold cross-validation;
• expected log predictive density.

This asks which model predicts unseen data better under a specified scoring rule.

Bayes factors

A Bayes factor compares marginal likelihoods:

Marginal likelihoods average likelihood over the prior. They can be highly sensitive to prior mass in parameter regions that fit poorly. Bayes factors are therefore not “objective Bayesian p-values.” The prior is part of the comparison.

LOO is not magic either

Leave-one-out predictive comparison estimates performance when each observation is held out in turn. For hierarchical, grouped, or time-series data, the held-out unit must match the deployment target. Leaving out one row when deployment involves a new patient or future time period can still leak structure. The validation unit matters in Bayesian models just as it does elsewhere.

Sequential updating requires a stable model

Bayes' rule makes sequential updating algebraically coherent:

But if the data-generating process changes over time, simply carrying yesterday's posterior forward can accumulate model misspecification. Dynamic models need explicit state evolution. For example,

Sequential Bayes is not a substitute for modeling drift.

A/B testing and optional stopping

Bayesian monitoring does not inherit the same Type I error interpretation as a fixed-horizon frequentist test. But that does not mean stopping rules become irrelevant to decisions. If the organization repeatedly checks posterior probabilities and deploys when a threshold is crossed, the long-run operating characteristics of that policy still matter.

Simulation should evaluate:

• false deployment rate;
• expected sample size;
• regret;
• probability of choosing the worse variant.

A posterior probability is an inferential quantity. A decision rule is an operating policy.

Bayesian networks are not automatically causal

A Bayesian network factorizes a joint distribution according to a directed acyclic graph. For variables

the factorization is

That graph can represent conditional independence without having a causal interpretation. Causal claims require additional assumptions about interventions and graph semantics. A directed graph is not causal simply because it is Bayesian.

Reproducible workflow example

A simple Beta-Binomial model can illustrate prior predictive and posterior predictive checking.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
from __future__ import annotations

import numpy as np
from scipy import stats

rng = np.random.default_rng(2026)

alpha: float = 2.0
beta: float = 8.0
n: int = 100

prior_theta = rng.beta(
    alpha,
    beta,
    size=10_000,
)

prior_counts = rng.binomial(
    n=n,
    p=prior_theta,
)

print(
    np.quantile(
        prior_counts,
        [0.05, 0.5, 0.95],
    )
)

observed_successes: int = 24

post_alpha: float = (
    alpha + observed_successes
)

post_beta: float = (
    beta + n - observed_successes
)

posterior_theta = rng.beta(
    post_alpha,
    post_beta,
    size=10_000,
)

posterior_rep = rng.binomial(
    n=n,
    p=posterior_theta,
)

print(
    np.quantile(
        posterior_rep,
        [0.05, 0.5, 0.95],
    )
)

The prior predictive distribution asks whether the prior implied plausible counts. The posterior predictive distribution asks whether the fitted model reproduces plausible future counts given the observed data.

Conclusion

A Bayesian analysis is not complete when the posterior sampler returns draws. The workflow is

Bayesian inference is powerful because every uncertainty statement belongs to one coherent probability model. That is also the responsibility: if the model is wrong, the posterior can be precisely wrong.

References

• Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). CRC Press.
• Gelman, A., Vehtari, A., Simpson, D., et al. (2020). Bayesian workflow. arXiv:2011.01808.
• Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27, 1413–1432.
• Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.