Skip to main content Skills Marketplace Discover and explore AI skills built by the community.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
Show prompt details A direct command skips the review prompt. Inspect the source before running it.
npx skills add https://github.com/synthetic-sciences/openscience --skill pymc-bayesian-modelingThe command stays on one line. Scroll horizontally to inspect it before copying.
Prefer a local copy? Download the files currently available to SkillsMP.
Download Zip Downloading... More from this repository
Related occupations SOC
Based on SOC occupation classification
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. category coding license Apache License, Version 2.0 metadata {"skill-author":"Synthetic Sciences"} version 1.0.0 author Synthetic Sciences tags ["Bayesian","MCMC","Probabilistic Programming","Statistics"] dependencies ["pymc>=5.10.0","arviz>=0.17.0","numpy>=1.25.0"]
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
X = ...
y = ...
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'
0
1
'beta'
0
1
'predictors'
'sigma'
1
'y_obs'
'obs_id'
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 )
az.plot_ppc(prior_pred, group='prior' )
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:
idata = pm.sample(
draws=2000 ,
tune=1000 ,
chains=4 ,
target_accept=0.9 ,
random_seed=42 ,
idata_kwargs={'log_likelihood' : True }
)
draws=2000: Number of samples per chaintune=1000: Warmup samples (discarded)chains=4: Run 4 chains for convergence checkingtarget_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' ])
R-hat < 1.01 : Chains have convergedESS > 400 : Sufficient effective samplesNo divergences : NUTS sampled successfullyTrace plots : Chains should mix well (fuzzy caterpillar)
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 with model:
pm.sample_posterior_predictive(idata, extend_inferencedata=True , random_seed=42 )
az.plot_ppc(idata)
Do posterior predictions capture observed data patterns?
Are systematic deviations evident (model misspecification)?
Consider alternative models if fit is poor
7. Analyze Results
print (az.summary(idata, var_names=['alpha' , 'beta' , 'sigma' ]))
az.plot_posterior(idata, var_names=['alpha' , 'beta' , 'sigma' ])
az.plot_forest(idata, var_names=['beta' ], combined=True )
8. Make Predictions X_new = ...
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_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 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:
mu_alpha = pm.Normal('mu_alpha' , mu=0 , sigma=10 )
sigma_alpha = pm.HalfNormal('sigma_alpha' , sigma=1 )
alpha_offset = pm.Normal('alpha_offset' , mu=0 , sigma=1 , dims='groups' )
alpha = pm.Deterministic('alpha' , mu_alpha + sigma_alpha * alpha_offset, dims='groups' )
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 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
models = {
'Model1' : idata1,
'Model2' : idata2,
'Model3' : idata3
}
comparison = compare_models(models, ic='loo' )
check_loo_reliability(models)
Δloo < 2 : Models are similar, choose simpler model2 < Δloo < 4 : Weak evidence for better model4 < Δloo < 10 : Moderate evidenceΔloo > 10 : Strong evidence for better model
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
pm.HalfNormal('sigma', sigma=1) - Default choicepm.Exponential('sigma', lam=1) - Alternativepm.Gamma('sigma', alpha=2, beta=1) - More informative
pm.Normal('theta', mu=0, sigma=1) - For standardized datapm.StudentT('theta', nu=3, mu=0, sigma=1) - Robust to outliers
pm.LogNormal('theta', mu=0, sigma=1)pm.Gamma('theta', alpha=2, beta=1)
pm.Beta('p', alpha=2, beta=2) - Weakly informativepm.Uniform('p', lower=0, upper=1) - Non-informative (use sparingly)
pm.LKJCorr('corr', n=n_vars, eta=2) - eta=1 uniform, eta>1 prefers identity
For Likelihoods
pm.Normal('y', mu=mu, sigma=sigma) - Default for continuous datapm.StudentT('y', nu=nu, mu=mu, sigma=sigma) - Robust to outliers
pm.Poisson('y', mu=lambda) - Equidispersed countspm.NegativeBinomial('y', mu=mu, alpha=alpha) - Overdispersed countspm.ZeroInflatedPoisson('y', psi=psi, mu=mu) - Excess zeros
pm.Bernoulli('y', p=p) or pm.Bernoulli('y', logit_p=logit_p)
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
)
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' )
start = approx.sample(return_inferencedata=False )[0 ]
idata = pm.sample(start=start)
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/'
)
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
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
Sample more draws: draws=5000
Reparameterize to reduce posterior correlation
Use QR decomposition for regression with correlated predictors
High R-hat
Run longer chains: tune=2000, draws=5000
Check for multimodality
Improve initialization with ADVI
Slow Sampling
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 samplingUse weakly informative priors (not flat)Use named dimensions (dims) for clarityNon-centered parameterization for hierarchical modelsCheck prior predictive before fitting
Sampling
Run multiple chains (at least 4) for convergenceUse target_accept=0.9 as baseline (higher if needed)Include log_likelihood=True for model comparisonSet random seed for reproducibility
Validation
Check diagnostics before interpretation (R-hat, ESS, divergences)Posterior predictive check for model validationCompare multiple models when appropriateReport 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
References (references/)
distributions.md
sampling_inference.md
workflows.md
Scripts (scripts/)
model_diagnostics.pycheck_diagnostics() for quick checks, create_diagnostic_report() for comprehensive analysis with plots.
model_comparison.pycompare_models(), check_loo_reliability(), model_averaging().
Templates (assets/)
linear_regression_template.py
hierarchical_model_template.py
Quick Reference
Model Building with pm.Model(coords={'var' : names}) as model:
param = pm.Normal('param' , mu=0 , sigma=1 , dims='var' )
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
