---

name: "pymc-bayesian-modeling" description: "Bayesian modeling with PyMC 5: priors, likelihood, NUTS/ADVI sampling, diagnostics (R-hat, ESS), LOO/WAIC comparison, prediction. Hierarchical, logistic, GP variants; predictive checks." license: "Apache-2.0"

PyMC Bayesian Modeling

Overview

PyMC is a Python library for Bayesian statistical modeling and probabilistic programming. It provides an expressive syntax for defining probabilistic models and efficient inference via MCMC (NUTS) and variational methods (ADVI). This skill covers the full Bayesian modeling cycle from model specification through diagnostics, comparison, and prediction.

When to Use

• Estimating parameters with full uncertainty quantification (credible intervals, not just point estimates)
• Fitting hierarchical/multilevel models to grouped or nested data
• Performing prior and posterior predictive checks to validate model assumptions
• Comparing candidate models using information criteria (LOO-CV, WAIC)
• Building regression models (linear, logistic, Poisson) in a Bayesian framework
• Handling missing data or measurement error as latent parameters
• Modeling time series with autoregressive or random walk priors
• Generating posterior predictions for new observations with uncertainty bounds
• Use Stan/PyStan instead for compiled, more scalable Bayesian inference on large models; use statsmodels for frequentist statistical tests

Prerequisites

• Python packages: pymc >= 5.0, arviz, numpy, matplotlib
• Data: NumPy arrays or pandas DataFrames with numeric columns
• Environment: CPU sufficient for most models; GPU via JAX backend for large models

pip install pymc arviz numpy matplotlib
# Optional: JAX backend for GPU acceleration
pip install pymc[jax]

Quick Start

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

# Simulate data
np.random.seed(42)
X = np.random.randn(100)
y = 2.5 + 1.3 * X + np.random.randn(100) * 0.5

# Build and fit model
with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=0, sigma=5)
    beta = pm.Normal("beta", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=1)
    mu = alpha + beta * X
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)
    idata = pm.sample(1000, tune=1000, chains=4, random_seed=42)

print(az.summary(idata, var_names=["alpha", "beta", "sigma"]))
# Expected: alpha ~ 2.5, beta ~ 1.3, sigma ~ 0.5

Workflow

Step 1: Prepare Data

Standardize continuous predictors for better sampling efficiency. Use named coordinates for readable models and ArviZ integration.

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

# Load data
X = np.random.randn(200, 3)  # 200 obs, 3 predictors
y = X @ np.array([1.0, -0.5, 0.3]) + np.random.randn(200) * 0.8

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

# Define coordinates for named dimensions
coords = {
    "predictors": ["var1", "var2", "var3"],
    "obs_id": np.arange(len(y)),
}
print(f"Data shape: X={X_scaled.shape}, y={y.shape}")

Step 2: Define Model and Set Priors

Specify the model structure inside a pm.Model() context. Use weakly informative priors, dims for named dimensions, and HalfNormal or Exponential for scale parameters.

with pm.Model(coords=coords) as model:
    # Priors — weakly informative, not flat
    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")

# Inspect model variables
print(model.basic_RVs)  # Lists: [alpha, beta, sigma, y_obs]

Step 3: Prior Predictive Check

Validate that priors produce plausible data ranges before fitting. Adjust priors if simulated data is unreasonable.

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

# Check prior-implied data range
prior_y = prior_pred.prior_predictive["y_obs"].values.flatten()
print(f"Prior predictive range: [{prior_y.min():.1f}, {prior_y.max():.1f}]")
print(f"Observed data range:    [{y.min():.1f}, {y.max():.1f}]")

az.plot_ppc(prior_pred, group="prior", num_pp_samples=100)

Step 4: Sample Posterior (MCMC)

Run NUTS sampling with multiple chains. Include log_likelihood=True if you plan model comparison later.

with model:
    idata = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        target_accept=0.9,
        random_seed=42,
        idata_kwargs={"log_likelihood": True},
    )

print(f"Posterior shape: {idata.posterior['beta'].shape}")
# Expected: (4 chains, 2000 draws, 3 predictors)

Step 5: Diagnose Sampling

Check convergence before interpreting results. All three diagnostics (R-hat, ESS, divergences) must pass.

# Summary with convergence diagnostics
summary = az.summary(idata, var_names=["alpha", "beta", "sigma"])
print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat", "ess_bulk"]])

# R-hat convergence check
bad_rhat = summary[summary["r_hat"] > 1.01]
if len(bad_rhat) > 0:
    print(f"WARNING: {len(bad_rhat)} parameters with R-hat > 1.01")
    print(bad_rhat[["r_hat"]])

# Effective sample size check
low_ess = summary[summary["ess_bulk"] < 400]
if len(low_ess) > 0:
    print(f"WARNING: {len(low_ess)} parameters with ESS < 400")

# Divergence check
n_div = idata.sample_stats.diverging.sum().item()
total = len(idata.posterior.draw) * len(idata.posterior.chain)
print(f"Divergences: {n_div}/{total} ({n_div / total * 100:.2f}%)")

# Visual diagnostics — trace plots and rank plots
az.plot_trace(idata, var_names=["alpha", "beta", "sigma"])
az.plot_rank(idata, var_names=["alpha", "beta", "sigma"])

Step 6: Posterior Predictive Check

Validate model fit by comparing simulated data from the posterior to observed data.

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

az.plot_ppc(idata, num_pp_samples=100)
# Blue = observed data, grey = posterior simulations
# Systematic deviations indicate model misspecification

Step 7: Compare Models

Use LOO-CV or WAIC to compare candidate models. Lower information criterion is better.

# Fit multiple models with log_likelihood=True, then compare
# Example: compare linear vs a second model
idatas = {"linear": idata}  # add more fitted models here

comparison = az.compare(idatas, ic="loo")
print(comparison[["rank", "elpd_loo", "p_loo", "d_loo", "weight"]])

# Check LOO reliability via Pareto-k diagnostics
loo_result = az.loo(idata, pointwise=True)
high_k = (loo_result.pareto_k > 0.7).sum().item()
print(f"Observations with Pareto-k > 0.7: {high_k}")
# Interpretation: Dloo < 2 = similar models; Dloo > 10 = strong evidence

az.plot_compare(comparison)

Step 8: Generate Predictions

Produce posterior predictions for new data with full uncertainty propagation.

X_new = np.array([[0.5, -1.0, 0.2]])
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
    )

y_pred = post_pred.posterior_predictive["y_obs"]
print(f"Predicted mean: {y_pred.mean().item():.3f}")
print(f"94% HDI: {az.hdi(y_pred, hdi_prob=0.94).values}")

Key Parameters

| Parameter | Default | Range / Options | Effect | |-----------|---------|-----------------|--------| | draws | 1000 | 500-10000 | Number of posterior samples per chain | | tune | 1000 | 500-5000 | Warmup iterations (discarded); increase for complex posteriors | | chains | 4 | 2-8 | Number of independent chains; minimum 4 for reliable R-hat | | cores | all CPUs | 1-N | Parallel chains; set equal to chains for full parallelism | | target_accept | 0.8 | 0.8-0.99 | NUTS acceptance rate; increase to reduce divergences | | init | "auto" | "adapt_diag", "jitter+adapt_diag", "advi" | Initialization strategy for sampler | | random_seed | None | any int | Seed for reproducibility | | idata_kwargs | {} | {"log_likelihood": True} | Store log-likelihood for LOO/WAIC model comparison | | method (pm.fit) | "advi" | "advi", "fullrank_advi", "svgd" | Variational inference algorithm | | n (pm.fit) | 10000 | 5000-100000 | VI optimization iterations | | samples (prior pred) | 500 | 100-5000 | Prior predictive samples for validation |

Key Concepts

Prior/Distribution Selection Guide

| Distribution | Use When | Key Parameters | |-------------|----------|----------------| | Normal(mu, sigma) | Unbounded real-valued parameter (standardized data) | mu: center, sigma: spread | | HalfNormal(sigma) | Scale/standard deviation parameter (positive) | sigma: spread of positive half | | Exponential(lam) | Scale parameter, alternative to HalfNormal | lam: rate (1/mean) | | StudentT(nu, mu, sigma) | Robust alternative to Normal (outlier-resistant) | nu: degrees of freedom (<10 = heavier tails) | | Beta(alpha, beta) | Probability or proportion in [0,1] | alpha=beta=2: weakly informative | | Gamma(alpha, beta) | Positive parameter (rate, concentration) | alpha: shape, beta: rate | | LogNormal(mu, sigma) | Positive parameter with multiplicative effects | mu, sigma: of underlying Normal | | LKJCorr(n, eta) | Correlation matrix prior | eta=1: uniform; eta>1: prefer identity | | Dirichlet(a) | Probability vector (sums to 1) | a: concentration; uniform if all equal | | Bernoulli(p / logit_p) | Binary outcome likelihood | Use logit_p for numerical stability | | Poisson(mu) | Count data (equidispersed) | mu: rate; use NegBinomial if overdispersed | | NegativeBinomial(mu, alpha) | Overdispersed count data | alpha: dispersion (smaller = more overdispersion) |

Diagnostic Thresholds

| Metric | Threshold | Interpretation | Action if Failed | |--------|-----------|---------------|------------------| | R-hat | < 1.01 | Chains converged | Run longer chains; check multimodality | | ESS bulk | > 400 | Sufficient independent samples | Increase draws; reparameterize | | ESS tail | > 400 | Reliable tail estimates | Increase draws | | Divergences | 0 | NUTS explored successfully | Increase target_accept; non-centered param. | | Pareto-k (LOO) | < 0.7 | LOO estimate reliable | Use WAIC or k-fold CV | | Max tree depth | < 10 | No trajectory truncation | Reparameterize or increase max_treedepth |

Model Variants Overview

| Problem Type | Recipe | Likelihood | Key Feature | |-------------|--------|------------|-------------| | Grouped/nested data | Hierarchical Model | Normal (varies) | Non-centered parameterization, partial pooling | | Binary outcome | Logistic Regression | Bernoulli | logit_p link function | | Nonlinear/spatial | Gaussian Process | Normal | Kernel-based covariance, flexible shape | | Count data | (use Poisson in Workflow) | Poisson / NegBinomial | Log link; NegBinomial for overdispersion | | Time series | (see references) | AR / GaussianRandomWalk | Autoregressive coefficients | | Mixture/clustering | (see references) | Mixture / NormalMixture | Component weights via Dirichlet |

Common Recipes

Recipe: Hierarchical Model

When to use: data has natural grouping (patients within hospitals, students within schools). Non-centered parameterization avoids divergences from funnel geometry.

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

n_groups, n_per_group = 5, 30
group_idx = np.repeat(np.arange(n_groups), n_per_group)
group_names = [f"group_{i}" for i in range(n_groups)]

# Simulated grouped data
true_alphas = np.random.normal(3.0, 1.5, n_groups)
y_obs = np.random.normal(true_alphas[group_idx], 0.5)

with pm.Model(co