Practical Python Pseudocode Patterns for Bayesian Workflows

Why “pseudocode patterns” matter in Bayesian Python work

In real projects, Bayesian analysis is less about a single model and more about a repeatable workflow: define data contracts, build a model, fit it, validate it, generate predictions, and turn those predictions into decisions. “Pseudocode patterns” are reusable code shapes that keep this workflow consistent across problems and teams. They help you avoid one-off notebooks that cannot be reproduced, and they make it easier to swap modeling choices (different likelihoods, priors, or link functions) without rewriting everything. This chapter focuses on practical Python-style pseudocode you can adapt to PyMC, Stan (via CmdStanPy), NumPyro, or other probabilistic programming tools, without re-teaching Bayesian fundamentals already covered earlier.

Pattern 1: A project skeleton that separates concerns

A common failure mode is mixing data cleaning, modeling, plotting, and decision logic in one place. A simple skeleton prevents that by making each step explicit and testable. The key idea is to treat Bayesian work like software: inputs and outputs are defined, and each stage can be rerun.

Step-by-step: minimal folder and module layout

• data/: raw extracts and immutable snapshots
• src/: reusable code (data prep, model builders, evaluation)
• reports/: generated figures/tables
• notebooks/: exploration only; final runs happen via scripts
• configs/: YAML/JSON configs for model and sampler settings

# src/pipeline.py (pseudocode) def run_pipeline(config_path: str):     cfg = load_config(config_path)     raw = load_raw_data(cfg.data)     df = clean_and_validate(raw, cfg.schema)     features = build_features(df, cfg.features)     model = build_model(features, cfg.model)     fit = fit_model(model, cfg.sampler)     checks = run_model_checks(model, fit, features, cfg.checks)     preds = generate_predictions(model, fit, features, cfg.prediction)     decision = compute_decision_metrics(preds, cfg.decision)     save_artifacts(fit, checks, preds, decision, cfg.output)     return decision

This pattern forces you to define a “contract” between steps: what columns exist, what shapes arrays have, and what artifacts are saved. It also makes it easy to rerun the same analysis with a different config, which is essential when stakeholders ask “what if we change X?”

Pattern 2: Data contracts and validation before modeling

Bayesian models are sensitive to data encoding mistakes: wrong units, swapped labels, missingness treated as zeros, or leakage from future information. A data contract pattern validates assumptions before sampling begins. The goal is to fail fast with clear errors rather than silently fitting a model to corrupted inputs.

Step-by-step: schema checks and shape checks

# src/validation.py (pseudocode) def clean_and_validate(raw, schema):     df = standardize_column_names(raw)     assert set(schema.required_columns).issubset(df.columns)     for col, rule in schema.rules.items():         if rule.type == "nonnegative":             assert (df[col] >= 0).all()         if rule.type == "in_range":             lo, hi = rule.lo, rule.hi             assert ((df[col] >= lo) & (df[col] <= hi)).all()         if rule.type == "categorical":             assert df[col].isin(rule.allowed).all()     df = handle_missingness(df, schema.missingness)     return df

For modeling code, also validate array shapes explicitly. For example, if you create an index array for groups, ensure it is integer-coded from 0 to G-1, and that the number of observations matches the length of the outcome vector. These checks prevent subtle bugs that can look like “bad convergence” but are actually data issues.

Continue in our app.

You can listen to the audiobook with the screen off, receive a free certificate for this course, and also have access to 5,000 other free online courses.

Or continue reading below...

Download the app

Pattern 3: A “model builder” function with a stable interface

Instead of writing a model directly inside a notebook cell, define a function that takes a feature dictionary and returns a model object. The stable interface is the pattern: the rest of the pipeline does not care whether the model is logistic regression, a time series, or a hierarchical model. It only needs the same return types: a model and a mapping of named variables you will later extract.

Step-by-step: model builder signature

# src/models.py (pseudocode) def build_model(data: dict, cfg: dict):     # data contains arrays like: y, X, group_idx, offsets, etc.     # cfg contains modeling choices: likelihood, link, priors, etc.     if cfg.family == "gaussian":         return build_gaussian_regression(data, cfg)     if cfg.family == "poisson":         return build_poisson_regression(data, cfg)     if cfg.family == "logistic":         return build_logistic_regression(data, cfg)     raise ValueError("Unknown family")

This pattern makes experimentation safe: you can compare models by swapping cfg.family while keeping the rest of the workflow identical.

Pattern 4: Centralized sampler configuration and reproducibility

Sampling settings are part of your analysis. If they live in scattered notebook cells, you cannot reproduce results. Centralize them in config and log them with the fitted model artifacts. Also, set random seeds in a single place and record library versions.

Step-by-step: sampler config object

# configs/model_run.yaml (conceptual) sampler:   seed: 202501   chains: 4   draws: 2000   tune: 2000   target_accept: 0.9   max_treedepth: 12 output:   run_id: "2026-01-06_pricing_model_v3"   artifact_dir: "reports/runs/"

# src/fitting.py (pseudocode) def fit_model(model, sampler_cfg):     set_global_seed(sampler_cfg.seed)     fit = sample(         model=model,         chains=sampler_cfg.chains,         draws=sampler_cfg.draws,         tune=sampler_cfg.tune,         target_accept=sampler_cfg.target_accept,         max_treedepth=sampler_cfg.max_treedepth,     )     log_run_metadata(sampler_cfg, library_versions(), system_info())     return fit

Even if you later change the sampler (NUTS vs. SVI vs. importance sampling), the interface remains stable: fit_model returns a fit object with posterior draws and diagnostics.

Pattern 5: Posterior extraction as a first-class step

Downstream tasks—plots, predictions, decisions—should not directly poke around in the fit object in ad hoc ways. Create a single extraction function that returns a tidy structure (arrays or a dataframe) with consistent naming. This reduces mistakes like mixing transformed and untransformed parameters or using the wrong dimension ordering.

Step-by-step: extract draws into a standard format

# src/posterior.py (pseudocode) def extract_posterior(fit, var_names: list[str]) -> dict:     draws = {}     for name in var_names:         arr = fit.get_draws(name)  # shape: (n_draws, ...)         draws[name] = arr     draws["n_draws"] = infer_n_draws(draws)     return draws

Once you standardize extraction, you can write generic utilities like “compute summary table,” “compute posterior predictive,” or “compute decision metric,” all without caring which modeling library produced the draws.

Pattern 6: Posterior predictive as a pure function

Predictions should be generated by a function that takes (a) posterior draws, (b) new data, and (c) a clear definition of what is being predicted (mean response, full outcome distribution, counterfactual difference, etc.). Treat prediction as a pure function: no hidden global state, no reliance on notebook variables.

Step-by-step: prediction API for new data

# src/predict.py (pseudocode) def posterior_predict(draws: dict, new_data: dict, cfg: dict) -> dict:     # Example for regression-like models     beta = draws["beta"]          # (S, K)     intercept = draws["alpha"]    # (S,)     X_new = new_data["X"]         # (N, K)     mu = intercept[:, None] + beta @ X_new.T  # (S, N)     if cfg.link == "logit":         p = sigmoid(mu)             # (S, N)         if cfg.return == "mean":             return {"p": p}         if cfg.return == "y":             y_rep = bernoulli_rng(p) # (S, N)             return {"p": p, "y_rep": y_rep}     if cfg.link == "identity":         if cfg.return == "y":             sigma = draws["sigma"]   # (S,)             y_rep = normal_rng(mu, sigma[:, None])             return {"mu": mu, "y_rep": y_rep}     raise ValueError("Unsupported link/return")

This pattern makes it straightforward to produce both point summaries (like posterior mean predictions) and full predictive distributions (simulated outcomes). It also encourages you to keep the transformation logic (e.g., logit) in one place.

Pattern 7: Decision computation as a separate layer

In applied work, the model is not the deliverable; the decision is. Implement decision metrics in a separate module that consumes predictions (or parameter draws) and returns business-facing quantities: expected utility, expected loss, probability a metric exceeds a threshold, or the distribution of incremental impact. Keeping this separate prevents the model from being “baked into” the decision logic and makes it easier to reuse the same decision layer across different models.

Step-by-step: expected value and threshold probability

# src/decision.py (pseudocode) def prob_exceeds(preds: dict, threshold: float, key: str) -> float:     x = preds[key]  # e.g., (S, N) or (S,)     return (x > threshold).mean()  # Monte Carlo estimate  def expected_profit(preds: dict, unit_margin: float, traffic: float) -> dict:     # Example: profit = traffic * conversion_rate * unit_margin     p = preds["p"]  # (S, N) or (S,)     profit = traffic * p * unit_margin     return {"profit": profit, "profit_mean": profit.mean()} 

Notice the deliberate simplicity: decision functions are just Monte Carlo computations over draws. This keeps your workflow transparent and auditable.

Pattern 8: Scenario and counterfactual analysis via “design matrices”

A powerful Bayesian workflow pattern is to represent “what-if” questions as alternative input datasets. You create multiple versions of new_data (scenarios) and run the same posterior_predict function on each. This is cleaner than writing special-case code for each scenario and reduces the risk of mixing assumptions.

Step-by-step: build scenarios and compare

# src/scenarios.py (pseudocode) def make_scenarios(base_df):     scenarios = {}     scenarios["current"] = build_features(base_df)     df_discount = base_df.copy()     df_discount["price"] *= 0.9     scenarios["10pct_discount"] = build_features(df_discount)     df_ads = base_df.copy()     df_ads["ad_spend"] += 1000     scenarios["more_ads"] = build_features(df_ads)     return scenarios  def compare_scenarios(draws, scenarios, pred_cfg):     out = {}     for name, data in scenarios.items():         out[name] = posterior_predict(draws, data, pred_cfg)     # Example: incremental effect distribution     inc = out["10pct_discount"]["p"] - out["current"]["p"]     return {"preds": out, "incremental": inc}

This pattern generalizes to policy evaluation, pricing, staffing, and forecasting: define scenarios as data transformations, then compare predictive distributions.

Pattern 9: Vectorization and shape discipline for speed and clarity

Bayesian workflows often involve large arrays of draws. A common performance trap is looping over draws in Python. Instead, vectorize operations so they run in NumPy/JAX. The pattern is to standardize shapes and document them. For example, use S for number of posterior draws, N for number of observations, K for number of features, and keep arrays consistently shaped as (S, N) for predicted quantities.

Step-by-step: a shape convention checklist

• Parameters: (S,) for scalars, (S, K) for coefficients
• Linear predictor: (S, N)
• Predicted mean/probability: (S, N)
• Simulated outcomes: (S, N)

# src/shapes.py (pseudocode) def assert_shapes(draws, data):     S = draws["n_draws"]     assert draws["alpha"].shape == (S,)     assert draws["beta"].shape[0] == S     assert data["X"].ndim == 2

Shape discipline is not cosmetic. It prevents silent broadcasting errors that can produce plausible-looking but wrong results.

Pattern 10: A unified “summary table” generator for stakeholders

Stakeholders rarely want raw draws. They want a small table of key quantities with uncertainty. Create a summary function that takes an array of draws and returns a consistent set of statistics (mean, median, intervals, and tail probabilities). This is a pattern because you will reuse it everywhere: parameter summaries, predicted KPI summaries, incremental impact summaries, and risk summaries.

Step-by-step: generic Monte Carlo summarizer

# src/summarize.py (pseudocode) def summarize_draws(x, probs=(0.05, 0.5, 0.95)) -> dict:     # x can be (S,) or (S, N); summarize along S     q = quantile(x, probs, axis=0)     return {         "mean": x.mean(axis=0),         "sd": x.std(axis=0),         "q_low": q[0],         "median": q[1],         "q_high": q[2],     }  def summarize_increment(inc_draws, threshold=0.0) -> dict:     s = summarize_draws(inc_draws)     s["prob_gt_0"] = (inc_draws > threshold).mean(axis=0)     return s

Because this function is generic, you can apply it to any derived quantity. That encourages a workflow where you compute the quantity you care about (profit, uplift, risk) and then summarize it, rather than summarizing parameters and hoping they translate into decisions.

Pattern 11: Caching expensive computations and artifact management

Sampling and posterior predictive simulation can be expensive. A practical pattern is to cache intermediate artifacts keyed by a run ID and a hash of the config and data snapshot. This allows you to regenerate reports without resampling and to compare runs reliably.

Step-by-step: artifact store with run IDs

# src/artifacts.py (pseudocode) def artifact_paths(base_dir, run_id):     return {         "fit": f"{base_dir}/{run_id}/fit.nc",         "preds": f"{base_dir}/{run_id}/preds.parquet",         "checks": f"{base_dir}/{run_id}/checks.json",         "decision": f"{base_dir}/{run_id}/decision.json",         "config": f"{base_dir}/{run_id}/config.yaml",     }  def save_artifacts(fit, checks, preds, decision, output_cfg):     paths = artifact_paths(output_cfg.artifact_dir, output_cfg.run_id)     save_config(paths["config"])     save_fit(fit, paths["fit"])     save_json(checks, paths["checks"])     save_table(preds, paths["preds"])     save_json(decision, paths["decision"]) 

With this pattern, you can always answer: “Which model run produced this chart?” and “Can we reproduce it?” without relying on memory.

Pattern 12: A lightweight “model registry” for comparing candidates

In practice you will fit multiple candidate models: different feature sets, different likelihoods, different pooling structures, or different priors. A model registry pattern stores each run’s config, diagnostics, and key predictive metrics in a single table so you can compare candidates systematically.

Step-by-step: register and compare runs

# src/registry.py (pseudocode) def register_run(run_id, cfg, diagnostics, metrics, registry_path):     row = {         "run_id": run_id,         "timestamp": now(),         "family": cfg.model.family,         "features": cfg.features.name,         "sampler": cfg.sampler,         "diagnostics": diagnostics,         "metrics": metrics,     }     append_jsonl(registry_path, row)  def select_best(registry_rows, key="expected_loss"):     # Example: pick smallest expected loss among runs that pass diagnostics     ok = [r for r in registry_rows if r["diagnostics"]["pass"]]     return sorted(ok, key=lambda r: r["metrics"][key])[0]

This pattern encourages disciplined iteration: you do not just “try another model,” you create a comparable run with recorded settings and outcomes.

Pattern 13: Testing Bayesian code with deterministic checks

Bayesian inference is stochastic, but much of your pipeline can be tested deterministically: data validation, feature building, shape checks, and even some probabilistic functions by fixing seeds and using small synthetic datasets. A testing pattern reduces the risk of shipping a pipeline that silently breaks when data changes.

Step-by-step: unit tests for feature and prediction functions

# tests/test_features.py (pseudocode) def test_build_features_shapes():     df = make_tiny_dataframe()     data = build_features(df, cfg={"use_price": True})     assert data["X"].shape[0] == len(df)     assert "y" in data  # tests/test_predict.py (pseudocode) def test_posterior_predict_shapes():     draws = {"alpha": zeros(100), "beta": zeros((100, 3)), "n_draws": 100}     new_data = {"X": zeros((10, 3))}     preds = posterior_predict(draws, new_data, cfg={"link": "identity", "return": "mean"})     assert preds["mu"].shape == (100, 10)

These tests do not prove the model is “correct,” but they do ensure your workflow is stable and that refactors do not break core assumptions.

Pattern 14: A single “analysis report” function that consumes artifacts

To avoid re-running expensive steps just to make plots, generate reports from saved artifacts. The pattern is: the pipeline produces artifacts; the report reads artifacts and renders tables/figures. This separation also makes it easier to schedule nightly runs and publish updated dashboards.

Step-by-step: report generation from stored outputs

# src/report.py (pseudocode) def build_report(artifact_dir, run_id):     paths = artifact_paths(artifact_dir, run_id)     fit = load_fit(paths["fit"])     preds = load_table(paths["preds"])     decision = load_json(paths["decision"])     # Create stakeholder-ready outputs     fig1 = plot_key_predictions(preds)     fig2 = plot_risk_distribution(decision)     table = make_summary_table(preds, decision)     save_fig(fig1, f"reports/{run_id}_preds.png")     save_fig(fig2, f"reports/{run_id}_risk.png")     save_table(table, f"reports/{run_id}_summary.csv")

Once you adopt this pattern, your notebooks can focus on exploration, while the production path is scriptable, repeatable, and auditable.