domains / science / bayesian-modeling / pymc

pymc

Skill card

combined

The K-Dense-AI copy is a highly specific, version-locked (PyMC 6.3.2) procedure that includes critical modern API details (ArviZ 1.3 changes, `nuts_sampler` vs `nuts_sampler_kwargs`, `PlotCollection` handling) and a rigorous scientific workflow (prior sensitivity, parameter recovery, estimand definition). The other two copies (Vibe-Skills and OpenClaw) are nearly identical to each other, targeting a generic 'PyMC 5.x+' API with standard code snippets and a simpler 8-step workflow. They lack the version-specific constraints and the deeper diagnostic requirements of the K-Dense copy. However, the Vibe/OpenClaw copies provide explicit code templates for common patterns (Linear/Logistic Regression) and a 'Routing Boundary' section that the K-Dense copy omits. Combining them allows for a skill that is both version-aware and practically templated.

Repositories

Other copies

Combined skill

---
name: pymc
description: Builds and checks Bayesian models with PyMC, including hierarchical models, NUTS MCMC, variational inference, mutable-data predictions, posterior predictive checks, diagnostics, and PSIS-LOO model comparison. Use for probabilistic modeling and uncertainty inference in PyMC.
allowed-tools: Read Write Edit Bash
compatibility: Requires Python 3.12+ with PyMC 6.3.2, PyTensor 3.3.2 and ArviZ 1.3-compatible dependencies; NumPy, pandas, Matplotlib, h5netcdf and h5py for bundled helpers/artifacts. Network access for installation only. Optional nutpie, NumPyro and BlackJAX samplers need separate dependencies.
license: Apache License, Version 2.0
metadata:
  version: "2.1"
  skill-author: K-Dense Inc.
  last-reviewed: "2026-10-01"
---

# PyMC Bayesian Modeling

## Version and scope

Targets **PyMC 6.3.2**, **PyTensor 3.3.2**, and **ArviZ 1.3.0**. These are local Python APIs, with no service endpoint, credentials, or remote inference required. The native tests use the PyMC NUTS sampler on small synthetic models. Optional nutpie/JAX samplers, SMC, and production-scale convergence are not validated by those tests. Record the complete resolved environment, backend, seeds, preprocessing, and data provenance for a scientific run.

Local macOS validation used `PYTENSOR_FLAGS="cxx="` after the C linker failed with `library 'd64' not found`; this disables PyTensor C compilation for that process. The Python library workflows passed under that setting; it does not establish that the C backend works on this host.

ArviZ 1 splits into base/stats/plots packages, re-exported by `import arviz as az`. PyMC returns an xarray `DataTree` despite the legacy `return_inferencedata` option name. `az.summary(..., round_to="none")` returns numerical columns; the default formats for display. Use `az.plot_ppc_dist`, `az.plot_dist`, and `az.plot_trace_dist`; save their returned `PlotCollection`, not a separate pyplot figure. `az.hdi(..., prob=.95)` uses `prob`, with bounds along `ci_bound`.

## Routing Boundary

Use this skill only for PyMC, probabilistic programming, Bayesian hierarchical models, NUTS/MCMC sampling, posterior predictive checks, and PyMC model diagnostics. Do not use it for generic regression, scikit-learn modeling, causal analysis, pymoo optimization, or materials-science prompts.

## Workflow

1. Define the estimand and predictive unit: another measurement, a new patient or group, or a future time block. Preserve units, missingness assumptions and grouping. Fit any scaling only on training data; reject constant columns.
2. Build a generative model with scientifically calibrated priors and named dimensions. Check likelihood support, finite initial log probability, design rank, confounding and symmetries before sampling.
3. Run a **prior predictive check** and inspect domain-relevant summaries before fitting. A broad prior is not automatically noninformative or plausible.
4. Fit several independent chains with explicit seeds. For reproducible backend selection use `nuts_sampler="pymc"`; the automatic choice can prefer installed nutpie. `nuts={...}` replaces deprecated `nuts_sampler_kwargs`.
5. Inspect rank R-hat, bulk/tail ESS, estimand-specific MCSE, divergences, energy BFMI and sampler-specific tree depth. Fix geometry and initialization before simply increasing draws. A clean screen does **not** prove convergence.
6. Generate posterior predictive replicates and inspect residual structure, tails, zeros, group effects and time dependence relevant to the question. In-sample agreement is not external validation or parameter identifiability.
7. Compare prior sensitivity and parameter recovery on synthetic data. Weakly informed scales, correlated parameters, separated logistic data, or mixture label switching can survive good MCMC diagnostics. A proper prior may make a posterior finite without the likelihood identifying its parameters.
8. Predict using the correct training transform and uncertainty levels. Keep `predictions` separate from in-sample `posterior_predictive`. Restore model data before recomputing training likelihoods.

See [standard_workflow.md](references/standard_workflow.md) for a complete bounded API smoke example, [workflows.md](references/workflows.md) for missing values, scaling, predictive scoring and serialization, and [model_patterns.md](references/model_patterns.md) for model structures.

## Runnable templates

From the skill directory, run the synthetic linear or hierarchical template:

```bash
uv run --isolated --with "pymc==6.3.2" --with "arviz==1.3.0" --with "pytensor==3.3.2" --with matplotlib --with pandas --with h5netcdf --with h5py python assets/linear_regression_template.py --draws 80 --tune 80 --chains 2 --output-dir linear_smoke
uv run --isolated --with "pymc==6.3.2" --with "arviz==1.3.0" --with "pytensor==3.3.2" --with matplotlib --with pandas --with h5netcdf --with h5py python assets/hierarchical_model_template.py --draws 80 --tune 80 --chains 2 --output-dir hierarchical_smoke
```

These short commands intentionally test execution and artifacts; they are not adequate evidence of converged scientific estimates. For a real analysis, inspect priors before fitting and budget chains/draws around diagnostics and required MCSE. The templates default to four chains and longer runs, but those counts are not acceptance criteria by themselves.

The hierarchical template predicts an unseen group by drawing group effects from every hyperposterior draw, sharing each group effect across that group's new rows, then adding observation noise. Using only population means omits between-group uncertainty. Centered, non-centered, and partial parameterizations can each be appropriate; non-centering is not universally superior.

## Diagnostic and comparison helpers

Import the modules from this skill's `scripts/` directory (add that directory to `sys.path` when working elsewhere). Example fragments below assume fitted data:

```python
from scripts.model_diagnostics import check_diagnostics, create_diagnostic_report
result = check_diagnostics(idata, var_names=["alpha", "beta", "sigma"])
create_diagnostic_report(idata, var_names=["alpha", "beta", "sigma"], output_dir="diagnostics")
```

`check_diagnostics` flags nonfinite statistics and provides a structured report. Use `create_diagnostic_report` to save plots and summaries.

## Common Model Patterns

### Linear Regression

For continuous outcomes with linear relationships:

```python
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:

```python
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)
    p = pm.math.sigmoid(logit_p)
    y = pm.Bernoulli('y', p=p, observed=y_obs)
```

**Use template:** `assets/logistic_regression_template.py`

### Hierarchical Model

For grouped data:

```python
with pm.Model() as hier_model:
    # Hyperpriors
    mu_alpha = pm.Normal('mu_alpha', mu=0, sigma=10)
    sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=1)

    # Group effects
    alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_alpha, shape=n_groups)

    # Likelihood
    mu = alpha[group_idx] + X @ beta
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)
```

**Use template:** `assets/hierarchical_model_template.py`