Standard Bayesian Workflow

The eight steps in full, with code: data preparation, model building, prior predictive check, fitting, diagnostics, posterior predictive check, analyzing results, and prediction.

Standard Bayesian Workflow

Follow this workflow for building and validating Bayesian models:

1. Data Preparation

import pymc as pm
import arviz as az
import numpy as np

# Load and prepare data
X = ...  # Predictors
y = ...  # Outcomes

# Standardize predictors for better sampling
X_mean = X.mean(axis=0)
X_std = X.std(axis=0)
X_scaled = (X - X_mean) / X_std

Key practices: - Standardize continuous predictors (improves sampling efficiency) - Center outcomes when possible - Handle missing data explicitly (treat as parameters) - Use named dimensions with coords for clarity

2. Model Building

coords = {
    'predictors': ['var1', 'var2', 'var3'],
    'obs_id': np.arange(len(y))
}

with pm.Model(coords=coords) as model:
    # Mutable data container so prediction data can be swapped later
    X_data = pm.Data('X_scaled', X_scaled, dims=('obs_id', 'predictors'))

    # Priors
    alpha = pm.Normal('alpha', mu=0, sigma=1)
    beta = pm.Normal('beta', mu=0, sigma=1, dims='predictors')
    sigma = pm.HalfNormal('sigma', sigma=1)

    # Linear predictor
    mu = alpha + pm.math.dot(X_data, beta)

    # Tie the observed variable's shape to X_data for out-of-sample prediction
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y, shape=X_data.shape[0], dims='obs_id')

Key practices: - Use weakly informative priors (not flat priors) - Use HalfNormal or Exponential for scale parameters - Use named dimensions (dims) instead of shape when possible - Use pm.Data() for values that will be updated for predictions

3. Prior Predictive Check

Always validate priors before fitting:

with model:
    prior_pred = pm.sample_prior_predictive(draws=1000, random_seed=42)

# Visualize
az.plot_ppc(prior_pred, group='prior')

Check: - Do prior predictions span reasonable values? - Are extreme values plausible given domain knowledge? - If priors generate implausible data, adjust and re-check

4. Fit Model

with model:
    # Optional: Quick exploration with ADVI
    # approx = pm.fit(n=20000)

    # Full MCMC inference
    idata = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        target_accept=0.9,
        random_seed=42,
        idata_kwargs={'log_likelihood': True}  # For model comparison
    )

Key parameters: - draws=2000: Number of samples per chain - tune=1000: Warmup samples (discarded) - chains=4: Run 4 chains for convergence checking - target_accept=0.9: Higher for difficult posteriors (0.95-0.99) - Include log_likelihood=True for model comparison - If using PyMC 6 sampler-specific kwargs, avoid deprecated nuts_sampler_kwargs; pass explicit NUTS kwargs through nuts={...} when needed

5. Check Diagnostics

Use the diagnostic script:

from scripts.model_diagnostics import check_diagnostics

results = check_diagnostics(idata, var_names=['alpha', 'beta', 'sigma'])

Check: - R-hat < 1.01: Chains have converged - ESS > 400: Sufficient effective samples - No divergences: NUTS sampled successfully - Trace plots: Chains should mix well (fuzzy caterpillar)

If issues arise: - Divergences → Increase target_accept=0.95, use non-centered parameterization - Low ESS → Sample more draws, reparameterize to reduce correlation - High R-hat → Run longer, check for multimodality

6. Posterior Predictive Check

Validate model fit:

with model:
    pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=42)

# Visualize
az.plot_ppc(idata)

Check: - Do posterior predictions capture observed data patterns? - Are systematic deviations evident (model misspecification)? - Consider alternative models if fit is poor

7. Analyze Results

# Summary statistics
print(az.summary(idata, var_names=['alpha', 'beta', 'sigma']))

# Posterior distributions
az.plot_posterior(idata, var_names=['alpha', 'beta', 'sigma'])

# Coefficient estimates
az.plot_forest(idata, var_names=['beta'], combined=True)

8. Make Predictions

X_new = ...  # New predictor values
X_new_scaled = (X_new - X_mean) / X_std

with model:
    pm.set_data({'X_scaled': X_new_scaled}, coords={'obs_id': np.arange(len(X_new_scaled))})
    post_pred = pm.sample_posterior_predictive(
        idata,
        var_names=['y_obs'],
        predictions=True,
        random_seed=42
    )

# Extract prediction intervals
y_pred_mean = post_pred.predictions['y_obs'].mean(dim=['chain', 'draw'])
y_pred_hdi = az.hdi(post_pred.predictions, var_names=['y_obs'])