Examples
Worked examples
Each example is a narrated, executed notebook: the question, the data, the model written in mathematics and in PyMC, every check, and what the result means. Read them in the browser, then run and modify them locally. The foundations come first; after that, each series stands on its own.
E01–E06
Foundations
The core of applied Bayesian work: writing a model down, checking what the priors claim, reading a posterior, pooling information across groups, comparing models and handling censored data. Every challenge assumes these.
E01The Bayesian workflow in PyMC 6Howell !Kung census
coords/dims, pm.Data, prior & posterior predictive checks, DataTree, diagnostics, out-of-sample predictionE02Hierarchical modelsMinnesota radon surveypooling, shrinkage, divergences, centred vs non-centred, group-level predictorsE03GLMs and model comparisonBangladesh arsenic wellslogistic regression, interpretable priors, calibration checks, PSIS-LOOE04Changepoints and discrete latentsUK coal-mining disasterswhy NUTS cannot do discrete, compound steps, marginalisationE05Gaussian processes with HSGPMotorcycle crash testGP priors, exact vs Hilbert-space GPs, heteroskedastic noiseE06Survival analysis and censoringMastectomy survivalpm.Censored, Weibull AFT, survival curves, prior sensitivityD01–D03
Data work
Most models that refuse to sample are data problems in disguise. These cover the craft around the model: preparing arrays, treating missingness and measurement error as part of the model, and carrying an analysis through to a report.
D01From a raw table to model-ready arraysPalmer penguins field records, IBM Telcoauditing and the unit of analysis, index vs dummy vs sum-to-zero coding as prior choices, centring vs standardising (measured), back-transforming with uncertainty, scaler leakage, the retransformation trap, collinearity before the sampler finds it, …D02Imperfect data as part of the modelPalmer penguins, US divorce rates with standard errorsmissing outcomes and covariates (complete-case vs mean-fill vs joint model, under a controlled deletion), marginalising a missing discrete covariate, a missing-not-at-random experiment and what identifies a selection model, measurement error and shrinkage, …D03Bayesian data analysis, end to endsix datasets for EDA, then penguinsexploratory plots that choose the likelihood (mean-variance slopes), priors from domain knowledge with PreliZ, why "weak" priors are not weak, estimands on the outcome scale, ROPE, prior sensitivity, a multiverse table, …
E07–E10, E46–E51, E78–E80
Mechanistic and scientific models
When theory supplies the regression function. Earthquakes as a point process in space and time, the expansion of the universe from supernovae, differential equations with JAX gradients, a black-hole ringdown in real LIGO strain, chaotic dynamics (fitting, forecasting and controlling systems that amplify every error, from blowflies to Lorenz's toy atmospheres), jump-diffusion SDEs for market crashes, and reaction-diffusion PDEs for cell invasion and signaling gradients, designing heaters, insulation and sensors for a heat-equation plate, and finding planets in Kepler light curves and proving an interstellar comet is unbound; and cosmology from whole images: every pixel of a Planck CMB patch, and a Hubble Einstein ring modelled pixel by pixel.
E07Geo-temporal: induced earthquakes in OklahomaUSGS catalogue, every M3+ quake 2005-19log-Gaussian Cox process, 2-D HSGP, space-time interaction so a hot-spot can move, taming priors behind an exp link, maps, hold-out forecastE08Physics-based models: weighing the universeUnion2.1, 580 type Ia supernovaethe theory is the regression function, a differentiable integral inside the model, a perfect H0-M degeneracy, P(expansion is accelerating)E09Your own JAX code, with gradients, inside a modelHudson's Bay lynx-hare peltshand-written Op +
jax.vjp, pytensor.wrap_jax, whole-model JAX (nutpie/NumPyro), a JAX ODE solver, multimodal likelihoods and start pointsE10Signal in the noise: GW150914Real LIGO Hanford + Livingston strain (GWOSC)PSD estimation and whitening, exact time-domain likelihood, black-hole ringdown -> remnant mass and spin, Bayes factor signal vs noise, an honest chirp-mass failureE46Taming chaos: inference, prediction and control of chaotic dynamicsSimulated chaotic Ricker population (Wood 2010's benchmark); Nicholson's blowflies, 180 bidaily counts (gamair)bifurcation diagram and Lyapunov exponent, why trajectory matching fails (a comb of ~1000 likelihood spikes, gradients growing like e^(lambda t), frozen NUTS chains), multiple shooting, …E47Weather jiu-jitsu on the butterfly: Bayesian assimilation and regime control of Lorenz-63Simulated Lorenz-63 with process noise and multiplicative observation error (the setting of Liu, Huang & Lall 2026, Chaos, Solitons & Fractals)two-regime attractor and residence times, instantaneous vs finite-time Lyapunov exponents (and what the LLE does not predict), multiple-shooting centred latent states with RK4 in PyTensor recovering sigma/rho/beta, …E48Regimes and extremes in a seasonal chaotic atmosphere: a Bayesian NHMM for Lorenz-84 and a test of NHMM-triggered controlSimulated seasonally forced Lorenz-84 (Lorenz 1990), daily noisy observations (the setting of Liu, Huang & Lall 2026, Phys. Rev. E)nonautonomous chaos, choosing regime variables (eddy amplitude, not phase), nonhomogeneous HMM with seasonal softmax transitions and AR(1) emissions via scan, HMM multimodality and multi-start MAP initialisation, ordered transform vs a -inf Potential, …E49Jumps: Bayesian jump-diffusion SDEs for the S&P 500S&P 500 daily returns 2008-19compound-Poisson jump-diffusion simulated exactly and by Euler, the Poisson-mixture transition density with the jump count summed out, Merton model (jumps absorbing volatility clustering), non-centred stochastic volatility with and without jumps (SVJ), …E50Reaction-diffusion in biology: cell invasion in a scratch assay, and reading a signaling gradientPC-3 scratch assay, 38 columns x 5 times (Jin et al. 2016, counts via Simpson, Murphy & Maclaren 2024); simulated Nodal/Lefty-scale morphogen images and FRAPFisher-KPP by the method of lines (finite volumes, no-flux, explicit Euler in scan with a stability check), Poisson counts with an estimated carrying capacity, a start-up delay tanh(beta t) and the fourfold bias in D without it, …E51Designing a thermal gradient plate: heaters, insulation and sensors for the 2-D heat equationSimulated 40 x 20 cm plate, calibration experiment (3 test heaters, 12 thermocouples) and operating disturbancesfinite-volume 2-D steady heat equation as a linear system, influence maps and superposition, Bayesian calibration with pt.linalg.solve inside PyMC (a weakly identified edge loss whose interval misses the truth), …E78Finding exoplanets and an interstellar cometKepler-10 light curves (NASA/MAST, 15 quarters) + NASA Exoplanet Archive; MPC astrometry of 2I/Borisov + JPL Horizons/SBDBBLS transit search in Numba, Gumbel false-alarm posterior, limb-darkened transit model (tested quadrature, tabulated in PyMC), marginalised local baselines, b-density degeneracy and reparameterisation, odd/even + occultation as model expansion, …E79CMB field-level inference: the sky behind the LMCPlanck 2018 SMICA map patches (IRSA mirror of the Planck Legacy Archive), half-mission noise, common mask, best-fit TT spectrumGaussian random fields on a flat-sky patch, naive masked-power bias checked on simulations, joint sky + bandpower posterior in PyMC (Hartley basis), centred / non-centred / partial non-centring at low and high signal-to-noise, why a mask hole defeats HMC, …E80A gravitational lens, pixel by pixelHST/ACS image of the Einstein ring SDSS J1627-0053 (MAST), SLACS published modelspixel-level forward model (SIE + shear, Sérsic lens light, FFT PSF from field stars), drizzle-correlated noise whitened in Fourier space, linear amplitudes integrated out, staged vs random starts, wrap_jax + nutpie low-rank, …E52–E61
Cell biology and biophysics
Mechanistic models of living cells, fitted to what biologists measure: signals that travel through the membrane and the cytosol, the fluctuation test that showed mutations arise before selection, why a neuron fires (conductance models fitted to real patch-clamp recordings), transcriptional bursting read from allele-resolved single-cell RNA counts, and how bacteria control their size, measured through noisy segmentation, and single molecules tracked in living nuclei, and super-resolution (STORM) microscopy: localising, detecting and counting single blinking molecules below the diffraction limit, and kinesin's steps counted from MINFLUX traces, the forces cells exert on their substrate (traction force microscopy), and decisions from dynamics: Turing patterns and kinetic proofreading in T cells. See also E50 (reaction-diffusion in tissues).
E52Signals in two dimensions and one: bulk-surface reaction-diffusion in a cell, the cytosolic shortcut, and polaritySimulated 20 um round cell; photoactivation of a membrane-cytosol shuttling protein imaged at 6 times (64 membrane arcs, 961 cytosolic volumes)2-D cytosolic PDE coupled to a 1-D membrane PDE by a Robin (binding/unbinding) boundary condition, finite volumes on a disk, exact angular Fourier decomposition (33 problems of size 17), …E53The Luria-Delbrück fluctuation test: jackpots and Bayesian mutation ratesLuria & Delbrück 1943 fluctuation tables (E. coli / phage T1, 281 cultures; transcribed, cross-checked with the flan package) + simulated culturescompound-Poisson likelihood as a triangular solve, partial plating, censored/binned counts, estimator simulation study (mean vs p0 vs median vs MLE), PSIS-LOO induced vs spontaneous mutation, hierarchical rates, tail posterior predictive checks, …E54Why a neuron fires: Bayesian inference of a conductance-based model from patch-clamp recordingsAllen Cell Types Database, one mouse Sst interneuron (21 current steps)RC fit with two exponentials + stationary AR(1), Pospischil HH model, Rush-Larsen in Numba and scan, trace-matching pathology, pm.Simulator + SMC feature likelihood, bimodal posterior, prior-whitened sloppiness spectrum, …E55Transcriptional bursting: the telegraph model and what mRNA snapshots identifyLarsson et al. 2019 allele-resolved scRNA-seq UMIs (24 genes, 165 G1 fibroblasts) + simulated cellsGillespie SSA, chemical master equation / finite state projection, tridiagonal recurrence, pt.linalg.expm time course, NB-limit map, capture as binomial thinning, Poisson/NB/telegraph LOO, prior sensitivity, experiment design, hierarchical LKJ, wrap_jaxE56Sizer, adder or timer? Bacterial size control and regression dilutionTanouchi et al. 2017 mother machine, E. coli at 25/27/37 °C, 19k cyclesnoisy linear map, regression dilution, errors-in-variables as a closed-form 16-dim MvNormal over 8-generation windows, measurement noise identified from lineage autocorrelation, add-noise stress test, over-correction check, PSIS-LOO sizer/adder/timer, PPC, …E57Single-particle tracking: diffusion states, localisation error, anomalous diffusionRARA-HaloTag and HaloTag-NLS tracks in live U2OS nuclei (saSPT example data, Heckert et al. 2022; MIT)exact MA(1) track likelihood via a sine eigenbasis, MSD-fit bias, D-sigma trade-off, Spot-On-style state mixture with state-dependent track loss, LOO, hierarchy across cells, GPB1 switching HMM in scan, two-filter smoothing, fractional Brownian motion, …E58Seeing below the diffraction limit: Bayesian localisation, detection and counting in STORMReal raw dSTORM EMCCD frames of Unc-13 at a Drosophila NMJ (Dannhäuser et al. 2022, Zenodo CC BY) + simulated nuclear-pore ringsphoton-transfer camera calibration, erf-pixelated PSF, EMCCD excess noise, NUTS vs MLE vs least squares, Cramér-Rao/Mortensen bounds, batched Fisher scoring + Laplace evidence, SMC evidence check, emitter counting as model selection, …E59Counting steps you cannot see: kinesin stepping and dwell-time kineticsMINFLUX kinesin-1 step tables + raw traces (Wolff, Scheiderer et al. 2023, Zenodo CC BY 4.0); simulated traces for step detectionhypoexponential dwells, randomness parameter, penalised step finder vs 2-phase lattice HMM (forward algorithm in scan), missed-step bias, ordered rates and label switching, LOO over hidden steps, dead-time/hidden-step equivalence, …E60How hard does a cell pull? Traction force microscopy as a Bayesian inverse problemSimulated cell with 16 known adhesions + real PIV displacement field of a 7-cell colony (pyTFM example data, Bauer et al. 2021)FTTC Boussinesq Green's function, longitudinal/transverse Fourier diagonalisation, Tikhonov as a Gaussian prior, L-curve vs discrepancy vs marginal likelihood, hyperparameters in PyMC with the field integrated out, exact FFT posterior draws, …E61Decisions from dynamics: Turing patterns and kinetic proofreadingSimulated Schnakenberg patterns; Nodal/Lefty diffusion coefficients (Müller et al. 2012, transcribed); real 1G4 T-cell CD69 dose-response + SPR KDs (Pettmann et al. 2021, eLife CC BY)linear stability, FFT semi-implicit PDE solver, calibrated feature likelihood by importance sampling, exact Fourier-mode likelihood (Woodbury), P(Turing conditions), discrete N marginalised with logsumexp, reparameterising to remove multimodality, …
E62
When the model cannot answer
Most failed analyses do not crash: they answer a question nobody asked, report numbers the prior chose, or never finish. Three real problems fail in these ways - an MMM that cannot split credit between channels, a Gaussian process too large to fit, a trial whose hazard ratio is not the clinic's question - and each is followed by the pivot: a question the data can answer, a cheaper model checked against the exact one, and an answer audit.
E63–E66
Causal designs, doses, cascades and experiments
Workhorse models from four fields, each carried through to a decision. Difference-in-differences, regression discontinuity and instrumental variables on classic natural experiments, with each identifying assumption written as a prior; population pharmacokinetics from theophylline and warfarin to a patient's dose after two monitoring samples; Hawkes and ETAS aftershock forecasts for the 2019 Ridgecrest sequence; and thousands of real headline A/B tests: winner's curse, peeking, and bandits.
E63Bayesian quasi-experiments: difference-in-differences, regression discontinuity, instrumental variablesCard & Krueger NJ/PA minimum wage; Mississippi banks 1929-34; Carpenter & Dobkin drinking-age mortality; Card (1995) college proximity; US Senate marginsfirst-difference DiD, Student-t robustness, placebo and alternative control groups, hierarchical regional trends (few clusters), parallel-trends violation as a prior with a break-down value, event study with counterfactual forecast, …E64Population pharmacokinetics: hierarchical compartment models and dosing decisionsTheophylline (12 subjects, Boeckmann/Sheiner/Beal via R datasets) and warfarin PK/PD (32 subjects, O'Reilly et al. via nlmixr2data)one-compartment oral model in closed form, hierarchical log-normal CL/V/ka with LKJ correlation, flip-flop mirror-image non-identifiability and multi-start chains, ordering constraint vs physiological prior, combined error compared by LOO, …E65Self-exciting point processes: Hawkes and ETAS aftershock forecastingUSGS ComCat catalogue, 2019 Ridgecrest M6.4/M7.1 sequence (6,909 M2+ events, 2016-2019)conditional intensity and branching ratio, point-process likelihood with closed-form compensator, Omori-Utsu ETAS with Gutenberg-Richter b-value, O(n^2) vs O(n) recursion vs power law as a mixture of exponentials (timed), …E66Bayesian online experiments: A/B tests at scale, peeking, winner's curse, banditsUpworthy Research Archive (4,873 headline A/B tests, 2013-2015; 1,000-test subsample for the model)Beta-Binomial P(best) and expected loss, natural A/A check, hierarchical Student-t effects over 1,000 tests, centred vs non-centred, tail PPC, split-half replication of the winner's curse, population prior for a new test, …
E67–E73
Rankings, choices, raters, compositions, topics and decisions
Outcomes that are not one number on a line. Skills estimated from who beat whom in tennis and from whole Formula 1 finishing orders; what people will pay for, from a panel of stated choices, with a market simulator; the true label when pathologists or anaesthetists disagree and nobody knows the answer; data that are parts of a whole; and quick decisions, where a drift diffusion model explains both the choice and how long it took, and says how cautious a person should be; and topic models for anything that comes as a bag of features: the stories behind Upworthy's headlines, the palettes of a century of posters, the intents behind web visits and the missions behind shopping baskets - and the newer topic models built for LDA's problems, tested on the same data.
E67Ranking from head-to-head results: Bradley-Terry, dynamic skill ratings and Plackett-LuceATP tour matches 2021-2026 (Jeff Sackmann, via archival mirror); Formula 1 race and sprint results 2021-2025 (Jolpica-F1)Bradley-Terry, sum-to-zero vs anchoring identification, surface partial pooling, sufficient-statistic vs informative PPC, Gaussian random-walk skills vs Elo, held-out log/Brier scores with paired SEs, P(best) and rank posteriors, knockout simulation, …E68Discrete choice: multinomial logit, mixed logit, willingness to pay, market simulationTrain's electricity-supplier stated-preference survey (mlogit Electricity: 361 people x 12 choice tasks)random-utility logit on (task, alt, attr) dims, choice-scale prior predictive, MLE cross-check, IIA, person-level PPC, mixed logit with LKJ + non-centred person effects, lognormal price coefficient, WTP posteriors (median vs mean, individual), …E69Truth without a gold standard: latent class models for noisy raters and diagnostic tests7 pathologists x 118 carcinoma slides (poLCA); Dawid & Skene's 5 anaesthetists x 45 patients (rater); simulated Hui-Walter and known-truth checksmajority vote as a model, two-class latent class model with logsumexp marginalisation, parameter vs degrees-of-freedom counting, label switching and a Youden (better-than-chance) fix, score/pairwise-agreement PPC, …E70Compositional data: parts of a whole (Dirichlet, logistic-normal, Dirichlet-multinomial)Aitchison's Arctic lake sediments (sand/silt/clay vs depth); 2024 UK general election, 542 English constituenciesclosure and spurious correlation, ALR/CLR/ILR log-ratios, Dirichlet regression (mean-precision), logistic-normal with LKJ and covariate-dependent scale, LOO across scales via the log-ratio Jacobian, multinomial vs Dirichlet-multinomial overdispersion, …E71Drift diffusion models: how people trade speed for accuracy, and what they should doWagenmakers et al. (2008) lexical decision, speed vs accuracy instructions (rtdists speed_acc: 17 people, 31,351 trials; a quarter fitted, the rest held out)Euler simulation with a first-passage continuity correction, speed-accuracy closed forms and SPRT optimality, Wiener first-passage density in PyTensor (Navarro-Fuss small/large-time series with a NaN-safe switch), …E72Topic models beyond text: LDA for headlines, poster colours, web sessions and shopping basketsUpworthy Research Archive (4,873 headline tests; tests as documents); 750 public-domain Library of Congress posters (1890s magazine, WWI, WPA) as 64x48 thumbnails; msnbc.com anonymous web data (UCI, 60,000 sessions); UCI Online Retail (3,661 customers x 2,087 products)LDA with token topics summed out in PyMC, known-truth check, label switching and Hungarian alignment, why NUTS fails on a sparse vocabulary, collapsed Gibbs in Numba, topic stability across chains, held-out likelihood vs purpose-driven choice of K, …E73After LDA: newer topic models, problem by problemUpworthy headline tests (as E72); UCI Online Retail customers x productscommon yardsticks (NPMI coherence, diversity, top-word Jaccard stability across restarts, downstream CV R2), Dirichlet-multinomial mixture (GSDMM) and biterm topic model for short texts, anchor-word spectral recovery and anchor-started Gibbs, …
E74–E77
Colour, light and images
When the estimate is a colour, a spectrum, a region or a word. How the dyes of Japanese woodblock prints fade, and what a 1766 Harunobu print looked like when new; why two paints can match in the shop and differ at home (metamers), inverted from colour to spectrum with a prior learned from 1,269 Munsell chips; segmenting photographs with an unknown number of regions and checking the result, and its confidence, against the people who outlined them; and how 110 languages divide colour space, speaker by speaker - each with uncertainty displays that still work when colour cannot be used to show the uncertainty.
E74Faded colours: running a fading model backwards, and showing uncertainty when colour is the answerRijksmuseum/RCE light-ageing study of Edo-period woodblock-print colourants (Baines et al. 2026, Zenodo, CC BY: microfading and xenon-arc reflectance spectra of 19 colourants and mixtures, Blue Wool references); CIE 1964 observer and D65; two impressions of Harunobu's Descending Geese of the Koto Bridges (c. 1766, Art Institute of Chicago, CC0)spectra to CIELAB and CIEDE2000 in NumPy and PyTensor, Kubelka-Munk K/S additivity, rank-one check of "one chromophore fading", hierarchical stretched-exponential fading rendered through the observer and fitted to colour, …E75Image segmentation with an unknown number of regions, checked against peopleBerkeley Segmentation Dataset BSDS500 (Martin et al. 2001; Arbelaez et al. 2011): ten photographs at half resolution (two training, eight test) with the segmentations of 5-6 people eachSLIC superpixels in NumPy/Numba, a truncated stick-breaking DP mixture in PyMC and why NUTS fails on it (label switching, chains in different modes, a binding truncation), label-invariant summaries (co-assignment, boundary probability, posterior of K), …E76Metamers: one colour, many spectraReflectance spectra of the 1,269 Munsell matte chips (Univ. of Kuopio / Eastern Finland spectral database, Zenodo mirror); CIE 1964 observer with D65, A (from Planck's law), FL11 and LED-B2the observer as a 3 x 41 map and its null space, metameric blacks, metamer-mismatch bodies by linear programming, hat-function illuminant integration, three logit-reflectance priors (smoothness GP, learned covariance, …E77How languages carve up colour: naming models, speaker variation and universalsWorld Color Survey (Kay, Berlin, Maffi, Merrifield & Cook 2009; via the jvosten/wcs mirror): 110 unwritten languages, 2,616 speakers each naming 330 Munsell chips, plus best-example (focus) choicessoftmax naming model over quadratic CIELAB logits (Gaussian categories in natural form, zero-sum weights), grid vs hue-circle vs CIELAB features scored on held-out speakers, pooled vs hierarchical speaker effects with a per-speaker predictive check, …
E81–E83
Directions, experiments and quantiles
Three model families that need their own tricks. Angles wrap around, so means, priors and links that work on a line break on a circle (wind at Santa Barbara, ocelots and their prey); Bayesian optimisation chooses the next reaction on a fully measured chemistry screen and is scored against 50 real chemists; and quantile regression for birth weights shows why the asymmetric Laplace posterior is too sure of itself, and how to fix it.
E81Circular data: directions, times of day and seasonsSanta Barbara airport hourly wind 2021-23 (IEM ASOS); Barro Colorado camera traps (Rowcliffe et al. 2014, CC BY); Oklahoma quakes as a nullmean resultant vector, von Mises / wrapped / projected normal, tan-half vs canonical (natural-parameter) link, wrapped-prior and step-mode failures with zero divergences, harmonic circular predictors, a regime mixture, circular PPC (rose bands, PIT), …E82Bayesian optimisation: choosing the next experimentShields et al. 2021 direct-arylation screen: all 1,728 reactions (12 ligands x 4 bases x 4 solvents x 3 conc x 3 temp), plus 50 chemists' campaigns and EDBO runs (MIT)GP over categorical + continuous factors (per-factor correlation kernels, main effects + interactions), MAP refits of a compiled PyMC model vs NUTS, EI/PI/UCB/Thompson (Kronecker + Matheron draws), paired benchmark vs random/OFAT/chemists, encoding failure, …E83Bayesian quantile regression, and when it lies20,000 US singleton births 2023 (NCHS natality public-use file) + a known-truth simulationcheck loss and a Frisch-Newton LP solver, the quantile process, the asymmetric Laplace working likelihood (mode right, spread wrong, units trap), coverage simulation, Yang-Wang-He sandwich, Gibbs posterior with a bootstrap-calibrated learning rate, …
E11–E17
At the research frontier
Neural networks inside differential equations, inference beyond plain NUTS, state-space models, likelihood-free inference, BART, the display of uncertainty and privacy-preserving meta-analysis. Expect honest negative results next to the successes.
E11Scientific ML: a neural network inside a differential equationHudson's Bay lynx-hare peltsuniversal / neural ODEs in JAX, priors over functions, why weight-space r_hat is meaningless, did the network rediscover mass-action?, physics vs flexibility on a held-out forecastE12The inference frontier: when plain NUTS is not enoughS&P 500 latent volatility (from C07)a reference posterior, low-rank and normalizing-flow adapted NUTS, Laplace / ADVI / Pathfinder, Pareto k-hat as a trust diagnostic when you have no referenceE13State-space models: let the Kalman filter integrate out the statesUK road casualties 1969-84 and the 1983 seat-belt lawthe Kalman recursion in 10 lines,
pymc_extras.statespace, 3 sampled parameters instead of 400, decomposition, one-step-ahead checks, a counterfactual with a negative control, missing data and forecasting for free - and an honest speed comparisonE14Simulation-based inference: Bayes without a likelihoodS&P 500 returns; 114 years of Canadian lynxrejection ABC calibrated against a known posterior, tolerance and summary statistics as the two approximations, pm.Simulator + SMC, the g-and-k distribution, Wood's synthetic likelihood, a population model that cannot produce the cycleE15Bayesian additive regression treesCapital Bikeshare hourly rentals (subsampled); LaLonde job trainingthe BART prior over functions, pmb.BART inside an ordinary PyMC model, PGBART + NUTS and diagnosing in function space, partial dependence / ICE / variable importance, out-of-sample prediction, an honest comparison with a GLM, BART for causal g-computationE16Showing uncertainty: one posterior, twenty displaysMotorcycle crash test; Minnesota radon with county boundaries; Premier League 2024/25epistemic vs predictive bands, fan charts, why a band is not a set of curves, spaghetti, animated hypothetical outcome plots, quantile dotplots, exceedance curves, caterpillars and ridgelines, …E17Meta-analysis: a database of studies that never shares a patientmetadat benchmarks: 13 BCG trials, 16 magnesium trials (ISIS-4), 56 studies in 11 districtssufficient statistics and why aggregate tables are lossless, two-stage vs one-stage models, forest plots with prediction intervals, tau / I² posteriors, meta-regression, leave-one-out and cumulative meta-analysis, …E18–E23
Correlation structures and graphs
Dependence that lives on a map, a network or a tree: areal models, learning graphs and causal structure from data, network outcomes, many correlated outcomes at once, and regression along a phylogeny.
E18Areal models on a graph: ICAR, proper CAR, BYM2Scottish lip cancer 1975-80, 56 district polygons (GeoDa) + WinBUGS neighbour listcontiguity graph from polygons in NumPy, islands and disconnected components, graph Laplacian and the BYM2 scaling factor,
pm.ICAR / pm.CAR / BYM2, Moran's I residual check, exceedance-probability maps, PSIS-LOO vs refit K-fold, spatial confoundingE19Learning a graph: Bayesian Gaussian graphical modelsSachs et al. (2005) T-cell flow cytometry, 11 phosphoproteins + 20-edge consensus networkmarginal vs partial correlation, LKJCholeskyCov precision and partial-correlation deterministics, why a horseshoe on the precision matrix leaves the positive-definite cone, Bayesian neighbourhood selection (node-wise horseshoe, AND/OR), …E20Learning a directed graph: Bayesian causal discoverySachs et al. (2005), 7466 cells in 9 intervention conditionslinear-Gaussian SEMs, Markov equivalence seen through LOO and pm.do, NOTEARS acyclicity as a pm.Potential (and why it is multimodal and scale-dependent), fixed-order horseshoe DAGs, exact BGe + dynamic-programming averaging over all orders, …E21Network data as the outcome: latent-space and stochastic block modelsLazega law-firm coworker network (71 lawyers, 2485 dyads)dyad logistic vs sociality random effects, Hoff latent space + Procrustes, mixed-membership SBM with marginalised roles (degree-corrected or not), co-clustering matrix, network PPCs (degree, transitivity, shared partners, geodesics), …E22Many correlated outcomes: full covariance, factor models, copulasKen French 12 and 30 US industry portfolios, daily 2015-2024EWMA volatility standardisation, MvNormal + LKJCholeskyCov, low-rank + diagonal factor models, sign-flip / rotation non-identifiability and what r_hat does (not) see, anchored loadings vs Procrustes alignment, choosing K by held-out density, …E23Correlation from a tree: phylogenetic regressionPrimates301 brain / body / group size + 10kTrees phylogenyreading an R phylo object without R, covariance from shared branch length, Brownian motion / Pagel's lambda / OU-as-GP regression, the OU ridge and a reparameterisation, conditional LOO from the precision matrix (and two wrong shortcuts), whitening, …E24–E29
More model families
Hidden Markov models, joint longitudinal and survival models, Dirichlet-process mixtures, sparse regression, imperfect detection in ecology and Bayesian neural networks.
E24Hidden Markov modelsUSGS M7+ world earthquake counts 1900-2024; elk GPS tracks (Morales et al. 2004)forward algorithm in
pytensor.scan vs pymc_extras DiscreteMarkovChain + marginalize, label switching and ordering, forward-backward and Viterbi per draw, dwell times, pseudo-residuals, integrated LOO and leave-future-out, padded multi-track HMM, …E25Joint longitudinal-survival modelsMayo Clinic PBC trial: serial bilirubin + survivalmixed model + Weibull hazard sharing LKJ random effects, cumulative hazard by Gauss-Legendre quadrature, baseline / LOCF / two-stage biases, informative dropout, centred vs non-centred vs hierarchically centred, …E26Bayesian nonparametrics: Dirichlet process mixturesEnzyme activity of 245 people (Richardson & Green 1997); PBC bilirubin slopesstick-breaking and the implied prior on the number of clusters, pm.StickBreakingWeights + pm.NormalMixture, truncation check, label-free diagnostics, co-clustering and Binder / VI point clustering, sensitivity to alpha, …E27Sparse regression: horseshoe, regularised horseshoe, spike-and-slabProstate cancer (ESL split); LARS diabetes, 64 predictors on 120 patientsshrinkage factor kappa, ridge vs Bayesian lasso vs horseshoe, tau0 from a prior guess and the simulated m_eff prior, horseshoe divergences and their fixes, spike-and-slab via pymc_extras.marginalize vs exact enumeration, correlated predictors, …E28Imperfect detection: occupancy, N-mixture, capture-recaptureSwiss breeding-bird survey (crossbill, mallard); snowshoe hare capture historieslatent presence / abundance summed out in pm.CustomDist (checked against pymc_extras.marginalize), finite-sample occupancy, naive vs corrected trends, visits x prior confounding, N-mixture truncation check, Poisson vs NegBin, …E29Bayesian neural networks: what the posterior over weights buys youMotorcycle crash test; concrete compressive strength (Yeh 1998, UCI)an MLP in PyMC with dims, priors over functions and the 1/sqrt(width) GP limit, NUTS vs ADVI / Pathfinder / MAP ensemble / last-layer Laplace on held-out density, function-space r_hat, per-input prior scales (ARD), BNN vs HSGP vs GLM, PIT and coverage, …E30–E31
Random partitions and random features
Exchangeable partitions and latent binary features: Pólya urns, Pitman–Yor processes and unseen species, and the Indian buffet process.
E30Exchangeable partitions: Pólya urns, the Chinese restaurant process, Pitman-Yor and unseen speciesMoby-Dick word tokens (Project Gutenberg); Barro Colorado Island 50-ha tree censusurn simulation and de Finetti limits, the EPPF checked by brute force, NUTS on (alpha, d) with no latent assignments, a stable log rising factorial, frequency-of-frequencies and Heaps-curve PPCs, exact held-out predictive from EPPF ratios, …E31The Indian buffet process: latent binary featuresUCI optdigits handwritten digits (300 training / 200 test digits by other writers)the IBP as a restaurant, stick-breaking and finite beta-Bernoulli, collapsed Gibbs in NumPy checked by brute-force enumeration, noise level vs alpha as the driver of K+, stuck chains and split-merge, label-free summaries, …
E32–E35, E45, E84–E85
Media measurement and market structure
One open multi-brand e-commerce dataset and five questions a marketing team actually asks: a custom likelihood for spend that chases demand, a hierarchical marketing-mix model with a budget decision, a synthetic control, causal mediation, and which shape of response curve the data support (a Weibull transform written in JAX). Plus market structure from household panels: which margarines compete, read off a Bayesian choice map, and what that means for a merger screen; and the same question from the 6-month purchase-incidence tables and partial demographics most studies actually get.
E32Building a custom likelihood: media spend that chases demandConjura MMM data: a UK skincare brand, 884 days of new customers, discounts and Google + Meta spendthe data-generating story before the likelihood, identification strategies compared, a Gaussian copula with a negative-binomial margin derived for a discrete outcome (Park-Gupta endogeneity correction), stable log-space implementation (incomplete-beta tails, …E33Hierarchical multi-market MMM and a budget decision under uncertaintyConjura MMM data: Scandinavian apparel brand, 3 markets x Google/Meta, 81 weeksGeometric adstock + Hill saturation, priors on cost per incremental customer (Meridian ROI-style), partial pooling across markets (centred vs non-centred), holiday/trend/seasonal baseline, prior predictive in euros, LOO vs no/full pooling, …E34Bayesian synthetic control of a media natural experimentConjura MMM data: a brand's US Meta switch-off (Oct 2023) + 21 donor brandschange-point scan for interventions, simplex (Dirichlet) SC, augmented/structural SC (horseshoe + random-walk level integrated out as a GP, CausalImpact-style), own control series, placebos in time/space/outcome, donor-pool sensitivity, cost per lost customer, …E35Halo effects as Bayesian causal mediation: does Meta work through branded search?Conjura MMM data: one UK clothing brand, weekly new customers, branded-search clicks, Meta/Google spend, UK + USnatural direct/indirect effects, joint mediator + outcome negative binomial model with adstock and saturation, counterfactual propagation from posterior draws (plug-in and Monte Carlo), product-of-coefficients check, "search as control" vs "search left out", …E45One family for every response curve: a Weibull transform in JAX for media modelsConjura MMM data: a UK apparel brand, 180 weeks of new customers and Google + Meta spendthe saturation family as a prior, a normalised Weibull response $(1-e^{-cx^k})/(1-e^{-c})$ spanning power, concave and S-curves with a shape map and closed-form inflection, a float32-safe custom JAX function (
expm1 + series switch, …E84Market structure analysis: who competes with whom?Allenby & Rossi (1991) margarine scanner panel (bayesm): 516 households, 4,470 purchases of 10 products with all shelf pricesswitching matrix and switching lift, IIA (logit diversion proportional to share), nested-logit tests of form vs brand vs tier hierarchies scored on held-out last purchases, factor-structured hierarchical logit as a choice map, …E85Market structure from purchase incidence: combining panel, activity and demographic data84.51 Complete Journey (CC0): household x 14 cereals 0/1 for Jan-Jun 2017 (2,391 households), shopping trips, demographics for a third, Jul-Dec as a test waveco-purchase lift confounded by volume, Bernoulli-logit factor model ladder (activity, category propensity, centred-loading choice map), a naive map whose first axis is volume, …E36–E41
State of the art, explained for everyone
Six model families used in practice today, each with a plain-language opening and a closing section of uncertainty displays for readers with no statistics: MRP, extreme-event attribution, epidemic nowcasting, radiocarbon chronologies, election forecasting and growth charts.
E36Multilevel regression and poststratification (MRP): what every state thinks, from one national survey2018 CCES (5,000-person sample + 55,000 held out; abortion-coverage question) + ACS poststratification table + 2016 vote share (Lopez-Martin, Phillips & Gelman case study)why raw state means and raking fail in small states, binomial-cell multilevel logit with
ZeroSumNormal effects and interactions, random-walk prior on ordered groups, state-level predictors, per-draw poststratification (states, age within state), …E37Extreme-event attribution: how much did warming load the dice for the 2021 Pacific Northwest heatwave?NOAA GHCN-Daily annual maximum temperature at Sea-Tac and Portland; NASA GISTEMP global temperatureblock maxima and the GEV (pymc_extras GenExtreme checked against scipy), a geophysical shape prior, a reparameterisation that removes support-wall divergences, non-stationary GEV on smoothed global temperature, …E38Real-time epidemic tracking: is it growing right now? Nowcasting + a renewal-equation RGerman COVID-19 hospitalisation reporting triangle (RKI via the Hospitalization Nowcast Hub), winter 2021/22the reporting triangle and real data-as-of-then snapshots, a delay-hazard nowcast (reporting weekday + weekly drift, negative binomial cells), the renewal equation as a triangular solve, a spurious "wave is ending" from truncated data, …E39How old is it? Radiocarbon calibration and Bayesian chronologiesIntCal20 calibration curve; all 49 measurements from the 1988 Shroud of Turin dating (3 labs, shroud + 3 known-age controls, hand-transcribed from Damon et al. 1989); 31 dates from the Sluggan Bog peat core (Bchron)calibration through a wiggly, uncertain curve, why NUTS on a calendar date fails (r_hat 2.85) and summing the date out on a 1-year grid instead, lab offsets and lab-by-sample heterogeneity anchored by controls, …E40Who will win? Dynamic Bayesian poll aggregation and election forecasting2016 US presidential polls (HuffPost Pollster via the Economist's us-potus-model), 1976-2016 state results, fundamentals, ACS state demographicsLinzer reverse random walks from a fundamentals prior, correlated state innovations, house/mode/population effects, poll-level + shared polling error with a rotation for the unidentified truth/bias split, electoral-college simulation and tipping points, …E41Is my child normal? Distributional regression (GAMLSS/LMS) and growth centile chartsFourth Dutch Growth Study: BMI of 7,294 boys aged 0-21 (gamlss.data dbbmi)constant-spread vs smooth-spread vs skewed vs heavy-tailed models, Bayesian P-splines (centred RW2, low-rank mass matrix) for every distribution parameter, …E42–E44
Posteriors as labelled arrays
A posterior is a labelled array with named chain, draw and model dimensions. These notebooks do all their data and posterior work by name, and use it for derived quantities, time aggregation, binning and unusual plots.
E42xarray foundations for posteriors: dims, coords, broadcasting and timeCapital Bikeshare hourly rentals 2011-12 (17k hours)table ->
date x hour Dataset with non-dimension coords, stack/unstack round trip to model rows, DataTree navigation and thinning, az.extract, broadcasting counterfactual grids by name, apply_ufunc RNG, "sum the draws, …E43Transforming posteriors with coords: dates, exceedances, ranks and anomaliesNOAA GHCN-Daily TMAX, Sea-Tac + Portland 1948-2025long -> wide xarray, align join pitfalls, concat over a named station dim, a harmonic basis as a DataArray shared by model and posterior, date posteriors with idxmax + interp (summer onset and length), …E44Plotting posteriors creatively from labelled arraysPalmer penguins (raw field records)multivariate normal with one LKJ covariance per species, apply_ufunc for KDEs / Cholesky / rankdata / solve, xr.dot ellipses and correlated draws, pairwise contrast matrices by renaming a dim, stack/sortby, FacetGrids straight from the posterior, …