pymc

name	pymc
description	Use when "PyMC", "Bayesian", "MCMC", "probabilistic programming", or asking about "Bayesian regression", "hierarchical model", "NUTS sampler", "posterior distribution", "prior predictive", "credible intervals", "uncertainty quantification"
version	1.0.0

PyMC Bayesian Modeling

Probabilistic programming with MCMC - build and fit Bayesian models.

When to Use

• Bayesian regression and classification
• Hierarchical/multilevel models
• Uncertainty quantification
• Prior/posterior predictive checks
• Model comparison (LOO, WAIC)

Quick Start

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

# Build model
with pm.Model() as model:
    # Priors
    alpha = pm.Normal('alpha', mu=0, sigma=1)
    beta = pm.Normal('beta', mu=0, sigma=1)
    sigma = pm.HalfNormal('sigma', sigma=1)

    # Likelihood
    mu = alpha + beta * X
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

    # Sample
    idata = pm.sample(2000, tune=1000, chains=4)

# Analyze
az.summary(idata)
az.plot_posterior(idata)

Model Building

coords = {'predictors': ['var1', 'var2', 'var3']}

with pm.Model(coords=coords) as model:
    # Priors with dimensions
    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 model
    mu = alpha + pm.math.dot(X, beta)

    # Likelihood
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

Common Models

Linear Regression

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

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

Logistic Regression

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

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

Hierarchical Model

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

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

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

Sampling

with model:
    # MCMC (default NUTS)
    idata = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        target_accept=0.9,
        random_seed=42
    )

    # Variational inference (faster, approximate)
    approx = pm.fit(n=20000, method='advi')

Diagnostics

import arviz as az

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

# R-hat should be < 1.01
# ESS should be > 400

# Trace plots
az.plot_trace(idata)

# Check divergences
print(f"Divergences: {idata.sample_stats.diverging.sum().values}")

Predictive Checks

with model:
    # Prior predictive (before fitting)
    prior_pred = pm.sample_prior_predictive(1000)
    az.plot_ppc(prior_pred, group='prior')

    # Posterior predictive (after fitting)
    pm.sample_posterior_predictive(idata, extend_inferencedata=True)
    az.plot_ppc(idata)

Model Comparison

# Fit models with log_likelihood
with model1:
    idata1 = pm.sample(idata_kwargs={'log_likelihood': True})

with model2:
    idata2 = pm.sample(idata_kwargs={'log_likelihood': True})

# Compare using LOO
comparison = az.compare({'model1': idata1, 'model2': idata2})
print(comparison)

Predictions

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

# Extract intervals
y_pred_mean = post_pred.posterior_predictive['y'].mean(dim=['chain', 'draw'])
y_pred_hdi = az.hdi(post_pred.posterior_predictive)

Best Practices

1. Standardize predictors for better sampling
2. Use weakly informative priors (not flat)
3. Non-centered parameterization for hierarchical models
4. Check diagnostics before interpretation (R-hat, ESS)
5. Prior predictive checks before fitting
6. Posterior predictive checks after fitting

Troubleshooting

• Divergences → Increase target_accept=0.95, use non-centered
• Low ESS → More draws, reparameterize
• High R-hat → Run longer chains

Resources

• Docs: https://www.pymc.io/projects/docs/
• Examples: https://www.pymc.io/projects/examples/