Skip to main content Skills Marketplace Découvrez et explorez les compétences IA créées par la communauté.
Installer avec Codex ou Claude Copiez ce prompt, collez-le dans Codex, Claude ou un autre assistant, puis laissez-le vérifier la page du skill et l'installer pour vous.
Afficher les détails du prompt Une commande directe contourne le prompt de vérification. Examinez la source avant de l'exécuter.
npx skills add https://github.com/synthetic-sciences/openscience --skill pymc-bayesian-modelingLa commande reste sur une seule ligne. Faites défiler horizontalement pour la vérifier avant de la copier.
Vous préférez une copie locale ? Téléchargez les fichiers actuellement disponibles dans SkillsMP.
Télécharger Zip Téléchargement... Métiers associés SOC
Basé sur la classification professionnelle SOC
Explorateur de fichiers
8 fichiers 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
