PyMC Bayesian Modeling

Overview

PyMC is a Python library for Bayesian modeling and probabilistic programming. Build, fit, validate, and compare Bayesian models using PyMC's modern API (version 5.x+), including hierarchical models, MCMC sampling (NUTS), variational inference, and model comparison (LOO, WAIC).

When to Use This Skill

This skill should be used when:

• Building Bayesian models (linear/logistic regression, hierarchical models, time series, etc.)

• Performing MCMC sampling or variational inference

• Conducting prior/posterior predictive checks

• Diagnosing sampling issues (divergences, convergence, ESS)

• Comparing multiple models using information criteria (LOO, WAIC)

• Implementing uncertainty quantification through Bayesian methods

• Working with hierarchical/multilevel data structures

• Handling missing data or measurement error in a principled way

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: # 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_scaled, beta)

# Likelihood
y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y, 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(samples=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

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}) post_pred = pm.sample_posterior_predictive( idata.posterior, var_names=['y_obs'], random_seed=42 )

Extract prediction intervals

y_pred_mean = post_pred.posterior_predictive['y_obs'].mean(dim=['chain', 'draw']) y_pred_hdi = az.hdi(post_pred.posterior_predictive, var_names=['y_obs'])

Common Model Patterns

Linear Regression

For continuous outcomes with linear relationships:

with pm.Model() as linear_model: alpha = pm.Normal('alpha', mu=0, sigma=10) beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors) sigma = pm.HalfNormal('sigma', sigma=1)

mu = alpha + pm.math.dot(X, beta)
y = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)

Use template: assets/linear_regression_template.py

Logistic Regression

For binary outcomes:

with pm.Model() as logistic_model: alpha = pm.Normal('alpha', mu=0, sigma=10) beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)

logit_p = alpha + pm.math.dot(X, beta)
y = pm.Bernoulli('y', logit_p=logit_p, observed=y_obs)

Hierarchical Models

For grouped data (use non-centered parameterization):

with pm.Model(coords={'groups': group_names}) as hierarchical_model: # Hyperpriors mu_alpha = pm.Normal('mu_alpha', mu=0, sigma=10) sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=1)

# Group-level (non-centered)
alpha_offset = pm.Normal('alpha_offset', mu=0, sigma=1, dims='groups')
alpha = pm.Deterministic('alpha', mu_alpha + sigma_alpha * alpha_offset, dims='groups')

# Observation-level
mu = alpha[group_idx]
sigma = pm.HalfNormal('sigma', sigma=1)
y = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)

Use template: assets/hierarchical_model_template.py

Critical: Always use non-centered parameterization for hierarchical models to avoid divergences.

Poisson Regression

For count data:

with pm.Model() as poisson_model: alpha = pm.Normal('alpha', mu=0, sigma=10) beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)

log_lambda = alpha + pm.math.dot(X, beta)
y = pm.Poisson('y', mu=pm.math.exp(log_lambda), observed=y_obs)

For overdispersed counts, use NegativeBinomial instead.

Time Series

For autoregressive processes:

with pm.Model() as ar_model: sigma = pm.HalfNormal('sigma', sigma=1) rho = pm.Normal('rho', mu=0, sigma=0.5, shape=ar_order) init_dist = pm.Normal.dist(mu=0, sigma=sigma)

y = pm.AR('y', rho=rho, sigma=sigma, init_dist=init_dist, observed=y_obs)

Model Comparison

Comparing Models

Use LOO or WAIC for model comparison:

from scripts.model_comparison import compare_models, check_loo_reliability

Fit models with log_likelihood

models = { 'Model1': idata1, 'Model2': idata2, 'Model3': idata3 }

Compare using LOO

comparison = compare_models(models, ic='loo')

Check reliability

check_loo_reliability(models)

Interpretation:

• Δloo < 2: Models are similar, choose simpler model

• 2 < Δloo < 4: Weak evidence for better model

• 4 < Δloo < 10: Moderate evidence

• Δloo > 10: Strong evidence for better model

Check Pareto-k values:

• k < 0.7: LOO reliable

• k > 0.7: Consider WAIC or k-fold CV

Model Averaging

When models are similar, average predictions:

from scripts.model_comparison import model_averaging

averaged_pred, weights = model_averaging(models, var_name='y_obs')

Distribution Selection Guide

For Priors

Scale parameters (σ, τ):

• pm.HalfNormal('sigma', sigma=1)

• Default choice

• pm.Exponential('sigma', lam=1)

• Alternative

• pm.Gamma('sigma', alpha=2, beta=1)

• More informative

Unbounded parameters:

• pm.Normal('theta', mu=0, sigma=1)

• For standardized data

• pm.StudentT('theta', nu=3, mu=0, sigma=1)

• Robust to outliers

Positive parameters:

• pm.LogNormal('theta', mu=0, sigma=1)

• pm.Gamma('theta', alpha=2, beta=1)

Probabilities:

• pm.Beta('p', alpha=2, beta=2)

• Weakly informative

• pm.Uniform('p', lower=0, upper=1)

• Non-informative (use sparingly)

Correlation matrices:

• pm.LKJCorr('corr', n=n_vars, eta=2)
• eta=1 uniform, eta>1 prefers identity

For Likelihoods

Continuous outcomes:

• pm.Normal('y', mu=mu, sigma=sigma)

• Default for continuous data

• pm.StudentT('y', nu=nu, mu=mu, sigma=sigma)

• Robust to outliers

Count data:

• pm.Poisson('y', mu=lambda)

• Equidispersed counts

• pm.NegativeBinomial('y', mu=mu, alpha=alpha)

• Overdispersed counts

• pm.ZeroInflatedPoisson('y', psi=psi, mu=mu)

• Excess zeros

Binary outcomes:

• pm.Bernoulli('y', p=p) or pm.Bernoulli('y', logit_p=logit_p)

Categorical outcomes:

• pm.Categorical('y', p=probs)

See: references/distributions.md for comprehensive distribution reference

Sampling and Inference

MCMC with NUTS

Default and recommended for most models:

idata = pm.sample( draws=2000, tune=1000, chains=4, target_accept=0.9, random_seed=42 )

Adjust when needed:

• Divergences → target_accept=0.95 or higher

• Slow sampling → Use ADVI for initialization

• Discrete parameters → Use pm.Metropolis() for discrete vars

Variational Inference

Fast approximation for exploration or initialization:

with model: approx = pm.fit(n=20000, method='advi')

# Use for initialization
start = approx.sample(return_inferencedata=False)[0]
idata = pm.sample(start=start)

Trade-offs:

• Much faster than MCMC

• Approximate (may underestimate uncertainty)

• Good for large models or quick exploration

See: references/sampling_inference.md for detailed sampling guide

Diagnostic Scripts

Comprehensive Diagnostics

from scripts.model_diagnostics import create_diagnostic_report

create_diagnostic_report( idata, var_names=['alpha', 'beta', 'sigma'], output_dir='diagnostics/' )

Creates:

• Trace plots

• Rank plots (mixing check)

• Autocorrelation plots

• Energy plots

• ESS evolution

• Summary statistics CSV

Quick Diagnostic Check

from scripts.model_diagnostics import check_diagnostics

results = check_diagnostics(idata)

Checks R-hat, ESS, divergences, and tree depth.

Common Issues and Solutions

Divergences

Symptom: idata.sample_stats.diverging.sum() > 0

Solutions:

• Increase target_accept=0.95 or 0.99

• Use non-centered parameterization (hierarchical models)

• Add stronger priors to constrain parameters

• Check for model misspecification

Low Effective Sample Size

Symptom: ESS < 400

Solutions:

• Sample more draws: draws=5000

• Reparameterize to reduce posterior correlation

• Use QR decomposition for regression with correlated predictors

High R-hat

Symptom: R-hat > 1.01

Solutions:

• Run longer chains: tune=2000, draws=5000

• Check for multimodality

• Improve initialization with ADVI

Slow Sampling

Solutions:

• Use ADVI initialization

• Reduce model complexity

• Increase parallelization: cores=8, chains=8

• Use variational inference if appropriate

Best Practices

Model Building

• Always standardize predictors for better sampling

• Use weakly informative priors (not flat)

• Use named dimensions (dims ) for clarity

• Non-centered parameterization for hierarchical models

• Check prior predictive before fitting

Sampling

• Run multiple chains (at least 4) for convergence

• Use target_accept=0.9 as baseline (higher if needed)

• Include log_likelihood=True for model comparison

• Set random seed for reproducibility

Validation

• Check diagnostics before interpretation (R-hat, ESS, divergences)

• Posterior predictive check for model validation

• Compare multiple models when appropriate

• Report uncertainty (HDI intervals, not just point estimates)

Workflow

• Start simple, add complexity gradually

• Prior predictive check → Fit → Diagnostics → Posterior predictive check

• Iterate on model specification based on checks

• Document assumptions and prior choices

Resources

This skill includes:

References (references/ )

distributions.md : Comprehensive catalog of PyMC distributions organized by category (continuous, discrete, multivariate, mixture, time series). Use when selecting priors or likelihoods.

sampling_inference.md : Detailed guide to sampling algorithms (NUTS, Metropolis, SMC), variational inference (ADVI, SVGD), and handling sampling issues. Use when encountering convergence problems or choosing inference methods.

workflows.md : Complete workflow examples and code patterns for common model types, data preparation, prior selection, and model validation. Use as a cookbook for standard Bayesian analyses.

Scripts (scripts/ )

model_diagnostics.py : Automated diagnostic checking and report generation. Functions: check_diagnostics() for quick checks, create_diagnostic_report() for comprehensive analysis with plots.

model_comparison.py : Model comparison utilities using LOO/WAIC. Functions: compare_models() , check_loo_reliability() , model_averaging() .

Templates (assets/ )

linear_regression_template.py : Complete template for Bayesian linear regression with full workflow (data prep, prior checks, fitting, diagnostics, predictions).

hierarchical_model_template.py : Complete template for hierarchical/multilevel models with non-centered parameterization and group-level analysis.

Quick Reference

Model Building

with pm.Model(coords={'var': names}) as model: # Priors param = pm.Normal('param', mu=0, sigma=1, dims='var') # Likelihood y = pm.Normal('y', mu=..., sigma=..., observed=data)

Sampling

idata = pm.sample(draws=2000, tune=1000, chains=4, target_accept=0.9)

Diagnostics

from scripts.model_diagnostics import check_diagnostics check_diagnostics(idata)

Model Comparison

from scripts.model_comparison import compare_models compare_models({'m1': idata1, 'm2': idata2}, ic='loo')

Predictions

with model: pm.set_data({'X': X_new}) pred = pm.sample_posterior_predictive(idata.posterior)

Additional Notes

• PyMC integrates with ArviZ for visualization and diagnostics

• Use pm.model_to_graphviz(model) to visualize model structure

• Save results with idata.to_netcdf('results.nc')

• Load with az.from_netcdf('results.nc')

• For very large models, consider minibatch ADVI or data subsampling