pymc-bayesian-modeling

name	pymc-bayesian-modeling
description	Bayesian modeling with PyMC. Build hierarchical models, MCMC (NUTS), variational inference, LOO/WAIC comparison, posterior checks, for probabilistic programming and inference.
license	Apache License, Version 2.0
metadata	{"skill-author":"K-Dense Inc."}
verified	false
lastVerifiedAt	"2026-02-19T05:29:09.098Z"
source	builtin
trust_score	100
provenance_sha	0747c4743ebe53b5

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:
    
    alpha = pm.Normal(, mu=, sigma=)
    beta = pm.Normal(, mu=, sigma=, dims=)
    sigma = pm.HalfNormal(, sigma=)

    
    mu = alpha + pm.math.dot(X_scaled, beta)

    
    y_obs = pm.Normal(, mu=mu, sigma=sigma, observed=y, dims=)