learning-bayesian-statistics/baygent-skills/bayesian-workflow

Bayesian Workflow

Workflow overview

Every Bayesian analysis follows this sequence. Do not skip steps -- especially model criticism.

1. Formulate — Define the generative story. What underlying process, that we're precisely trying to model, created the data?
2. Specify priors — See references/priors.md
3. Implement in PyMC — Write the model. Prefer PyMC 5+ syntax. Use the latest version possible.
4. Run prior predictive checks — pm.sample_prior_predictive(). Verify priors produce plausible data ranges before fitting
5. Inference — pm.sample(nuts_sampler="nutpie"). Always use nutpie for speed (the nutpie python package provides cutting-edge sampling). Don't hardcode the number of chains — let the sampler pick the best default for the platform.
6. Diagnose convergence — Use arviz_stats.diagnose(idata) as the first check (requires arviz-stats >= 1.0.0). It covers R-hat, ESS, divergences, tree depth, and E-BFMI in one call. See references/diagnostics.md
7. Criticize the model — See references/model-criticism.md
8. Check prior sensitivity — Run psense_summary(idata) to verify conclusions are robust to prior choices. Visualize with plot_psense_dist(idata) from arviz_plots. Requires log_likelihood and log_prior in the InferenceData — compute them after sampling if needed. See references/sensitivity.md
9. Compare models (if applicable) — See references/model-comparison.md
10. Report results — Generate <slug>/report.md using the canonical template in references/reporting.md. Run scripts/check_diagnostics.py to turn raw diagnostics into qualitative ratings + an ordered next-steps list, and use that output to fill the Assessment lines and Suggested Next Steps section. When the user mentions a non-technical audience or is new to Bayesian stats, additionally adapt the prose to plain language and include a glossary — but keep the canonical report structure as the audit trail.

Installation

Prefer conda-forge / mamba-forge to install PyMC and its dependencies — pip can cause issues with compiled backends (nutpie, JAX). Example:

mamba install -c conda-forge pymc nutpie arviz arviz-stats preliz

See Stack compatibility below for the PyMC 5.x vs 6.x / ArviZ 0.23 vs 1.x notes.

Stack compatibility (PyMC 5.x and 6.x)

This skill teaches the latest PyMC 6 / ArviZ 1.x idioms and stays runnable on PyMC 5.x during the transition (regulated/corporate environments can't always upgrade freely). The scripts are verified on both stacks.

Most code is identical across versions. Where an API genuinely diverges, prefer the form that runs on both:

arviz_plots (the ArviZ 1.x plotting package, imported as azp) installs on both stacks, so leading with azp.* plots is the most portable choice. Bare az.summary(idata) works on both, but there is no single interval kwarg that does — ci_prob/ci_kind is ArviZ 1.x and hdi_prob is ArviZ 0.23 (the bare default also differs: 89% ETI vs 94% HDI). arviz_stats.diagnose, az.loo, az.plot_forest/energy/pair/khat, and idata.to_netcdf() are unchanged on both; InferenceData and DataTree are both xarray-backed, so downstream .sel() / az.* calls are identical.

PyMC model template

import pymc as pm
import arviz as az
import numpy as np
RANDOM_SEED = sum(map(ord, "churn-logistic-v1"))
rng = np.random.default_rng(RANDOM_SEED)
# always use dimensions and coordinates in PyMC models
with pm.Model(coords=coords) as model:
    # use Data containers when working on a PyMC model
    data = pm.Data("data", df["y"].to_numpy(), dims="obs")
    # --- Priors ---
    # Always document WHY each prior was chosen
    mu = pm.Normal("mu", mu=0, sigma=10)  # Weakly informative: allows wide range
    # --- Data model ---
    pm.Normal("obs", mu=mu, sigma=1, observed=data, dims="obs")
    # --- Prior predictive check ---
    prior_pred = pm.sample_prior_predictive(random_seed=rng)
    # --- Inference ---
    idata = pm.sample(nuts_sampler="nutpie", random_seed=rng)
    idata.extend(prior_pred)
    # --- Posterior predictive check ---
    idata.extend(pm.sample_posterior_predictive(idata, random_seed=rng))
    # --- Compute log-likelihood and log-prior for sensitivity checks & LOO ---
    pm.compute_log_likelihood(idata, model=model)
    pm.compute_log_prior(idata, model=model)
    # --- Save immediately after sampling ---
    # Late crashes can destroy valid results. Save to disk before any post-processing.
    idata.to_netcdf("model_output.nc")

Critical rules

• Always run prior predictive checks before sampling. If prior predictions span implausible ranges, fix priors first. If you have issues or doubts for some parameters, use the PreliZ package to elicit priors from the user.
• Always check convergence before interpreting results. R-hat > 1.01 or ESS < 100 * nbr_chains means the results are unreliable.
• Always run posterior predictive checks. A model that fits well numerically but cannot reproduce the data is useless.
• Always run calibration checks (PIT / coverage). Use ArviZ's plot_ppc_pit for this — it handles all data types (continuous, binary, count) correctly. See references/model-criticism.md.
• Document every prior choice with a brief justification in a code comment.
• Never report point estimates alone. Always include credible intervals — a 94% HDI is a fine default, but no interval width is magic (see references/reporting.md).
• Use `arviz_stats.diagnose(idata)` as the first diagnostic on every model (arviz-stats >= 1.0.0). It checks R-hat, ESS, divergences, tree depth saturation, and E-BFMI in one call. Follow up with az.plot_trace(idata, var_names=[...]) for visual inspection, or az.plot_rank(idata, var_names=[...]) for rank-based convergence views (both run on both stacks). Pass var_names to focus on the parameters — ArviZ 1.x errors if the auto-selected set (e.g. a vector Deterministic over an obs dim) exceeds its subplot cap. The older az.plot_trace(idata, kind="rank_vlines") is ArviZ-0.23-only.
• Don't hardcode number of chains. Let PyMC / nutpie choose the optimal default for the user's platform. Just call pm.sample() without specifying chains.
• Use reproducible, descriptive seeds. Never use magic numbers like 42. Instead, derive a seed from the analysis name: RANDOM_SEED = sum(map(ord, "my-analysis-name")). Pass it to pm.sample(random_seed=rng), pm.sample_prior_predictive(random_seed=rng), and numpy via rng = np.random.default_rng(RANDOM_SEED).
• Save InferenceData immediately after sampling with idata.to_netcdf("model_output.nc"). Late crashes or kernel restarts can destroy valid MCMC results — save before any post-processing.
• Use ArviZ for all plots and calibration. Don't write custom plotting code when ArviZ already handles it — including for binary data, count data, and calibration. ArviZ developers have thought through edge cases so you don't have to.
• Prefer xarray over numpy for InferenceData operations. InferenceData and DataTree objects are backed by xarray — use xarray's labeled indexing (.sel(), .mean(dim=...), etc.) instead of converting to numpy arrays. This preserves dimension labels, avoids shape bugs, and makes code more readable. Fall back to numpy only when xarray can't do what you need.
• Always generate `<slug>/report.md` after a full analysis run. Store all artifacts (inference_data.nc, trace.png, forest.png, posterior_predictive.png, pit_ecdf.png, summary.csv, diagnostics.json, calibration.json) in a slug-named results folder, and produce report.md from the canonical template in references/reporting.md. Code without an interpreted, fixed-shape report is incomplete.
• Use `scripts/check_diagnostics.py` to interpret diagnostics, not hand-rolled prose. Pipe the JSON outputs of diagnose_model.py (and optionally calibration_check.py and psense_summary) into check_diagnostics.py to get per-section qualitative ratings and an ordered, actionable Suggested Next Steps list. Use those outputs verbatim in the report's Assessment lines; expand only with problem-specific context.
• Always use the posterior mean (not median) for predictive probabilities. The proper Bayesian predictive distribution averages over the posterior: P(Y=k|x) = (1/S) Σ P(Y=k|x,θₛ). This is the mean, not the median. The median does not correspond to the posterior predictive distribution, can violate probability coherence (probabilities may not sum to 1), and biases calibration due to Jensen's inequality. In code: use np.mean(probs, axis=sample_axis), never np.median(...).
• Use `pm.set_data()` + `pm.sample_posterior_predictive()` for out-of-sample predictions. Don't manually extract posterior samples and recompute predictions — let PyMC propagate uncertainty properly. Define predictors as pm.Data(...) during model building, then swap in new data:

# After fitting the model:
with model:
    pm.set_data({"X": X_new, "group_idx": group_idx_new})
    oos_preds = pm.sample_posterior_predictive(idata, predictions=True, random_seed=rng)

• Check model identifiability before interpreting components. If two model components always appear together in the likelihood (e.g., a league intercept and a home advantage term when every observation is from home perspective), their individual posteriors reflect prior assumptions, not data signal — only their sum is identified. Use az.plot_pair() to check for strong posterior correlations between components. If correlation is near ±1, the components are not separately identifiable — either merge them or restructure the data.

Common model families