/pymc

name: pymc
description: Bayesian modeling with PyMC. Build hierarchical models, MCMC (NUTS), variational inference, LOO/WAIC comparison, posterior checks, for probabilistic programming and inference.
allowed-tools: Read Write Edit Bash
compatibility: Requires Python 3.12+ and PyMC 6.0.1-compatible dependencies. Install reproducible environments with `uv pip install "pymc[nutpie]==6.0.1"`; optional NumPyro or BlackJAX samplers require separately pinned JAX-compatible dependencies.
license: Apache License, Version 2.0
metadata:
  version: "1.4"
  skill-author: K-Dense Inc.

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 6.x+), including hierarchical models, MCMC sampling (NUTS), variational inference, posterior predictive checks, and model comparison (LOO, WAIC).

Current Version and Setup

PyMC 6.0.1 is the current stable release as of June 2026. It requires Python 3.12+, uses PyTensor 3 as the computational graph backend, and defaults to compiled backends such as Numba. For reproducible local environments, pin the version:

uv pip install "pymc[nutpie]==6.0.1"

The `nutpie` extra enables the faster Rust/Numba NUTS implementation. If using NumPyro or BlackJAX, install those optional sampler dependencies in the same environment and pin them in the project lockfile.

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

Never sample first and check later. The eight-step workflow — documented with code in [references/standard_workflow.md](references/standard_workflow.md) — is:

1. **Data preparation** — including standardizing predictors so priors are interpretable. 2. **Model building** — priors and likelihood in a `pm.Model` context. 3. **Prior predictive check** — confirm the priors imply plausible data *before* fitting. 4. **Fit model** — `pm.sample()` with an explicit seed. 5. **Check diagnostics** — R-hat, ESS, divergences. Divergences invalidate the fit; fix the model or reparameterize rather than raising `target_accept` and hoping. 6. **Posterior predictive check** — does the fitted model reproduce the observed data? 7. **Analyze results** — summaries and intervals from the posterior. 8. **Make predictions** — on new data via `pm.set_data` and posterior predictive sampling.

Reusable model structures and model comparison are in [references/model_patterns.md](references/model_patterns.md).

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.LKJCholeskyCov('chol', n=n_vars, eta=2, sd_dist=pm.HalfNormal.dist(1))` - Preferred covariance prior
• `pm.LKJCorr('corr', n=n_vars, eta=2)` - Correlation-only prior; 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
• `pm.HurdleNegativeBinomial('y', psi=psi, mu=mu, alpha=alpha)` - Excess zeros plus overdispersion

**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
    initvals = approx.sample(return_inferencedata=False)[0]
    idata = pm.sample(initvals=initvals)

**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
• Local ESS plots
• Summary statistics CSV

Quick Diagnostic Check

from scripts.model_diagnostics imp