pymc-bayesian-modeling

•

24.4k23993 months ago

Bayesian modeling with PyMC. Build hierarchical models, MCMC (NUTS), variational inference, LOO/WAIC comparison, posterior checks, for probabilistic programming and inference.

Quick Install

npxskills add davila7/claude-code-templates--skill pymc-bayesian-modeling

Instructions

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

Loading…

Related Skills

skill-creator

Create or update AgentSkills. Use when designing, structuring, or packaging skills with scripts, references, and assets.

367,360

skill-creator

Create or update AgentSkills. Use when designing, structuring, or packaging skills with scripts, references, and assets.

367,309

1password

Set up and use 1Password CLI (op). Use when installing the CLI, enabling desktop app integration, signing in (single or multi-account), or reading/injecting/running secrets via op.

365,915

himalaya

CLI to manage emails via IMAP/SMTP. Use `himalaya` to list, read, write, reply, forward, search, and organize emails from the terminal. Supports multiple accounts and message composition with MML (MIME Meta Language).

365,885

Tags & Topics

anthropicanthropic-claudeclaudeclaude-code

PyMC Bayesian Modeling

Overview

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

pymc-bayesian-modeling

PyMC Bayesian Modeling

Overview

When to Use This Skill

Standard Bayesian Workflow

1. Data Preparation

2. Model Building

3. Prior Predictive Check

4. Fit Model

5. Check Diagnostics

6. Posterior Predictive Check

7. Analyze Results

8. Make Predictions

Common Model Patterns

Linear Regression

Logistic Regression

Hierarchical Models

Poisson Regression

Time Series

Model Comparison

Comparing Models

Model Averaging

Distribution Selection Guide

For Priors

For Likelihoods

Sampling and Inference

MCMC with NUTS

Variational Inference

Diagnostic Scripts

Comprehensive Diagnostics

Quick Diagnostic Check

Common Issues and Solutions

Divergences

Low Effective Sample Size

High R-hat

Slow Sampling

Best Practices

Model Building

Sampling

Validation

Workflow

Resources

References (references/)

Scripts (scripts/)

Templates (assets/)

Quick Reference

Model Building

Sampling

Diagnostics

Model Comparison

Predictions

Additional Notes

Related Skills

skill-creator

skill-creator

1password

himalaya

PyMC Bayesian Modeling

Overview

When to Use This Skill

Standard Bayesian Workflow

1. Data Preparation

2. Model Building

3. Prior Predictive Check

4. Fit Model

5. Check Diagnostics

6. Posterior Predictive Check

7. Analyze Results

8. Make Predictions

Common Model Patterns

Linear Regression

Logistic Regression

Hierarchical Models

Poisson Regression

Time Series

Model Comparison

Comparing Models

Model Averaging

Distribution Selection Guide

For Priors

References (`references/`)

Scripts (`scripts/`)

Templates (`assets/`)

References (`references/`)

Scripts (`scripts/`)

Templates (`assets/`)