32.1 The Bayesian Workflow

1. Bayesian Models as Scientific Hypotheses

In 2015, physicists at LIGO faced a question no amount of raw data could answer on its own: given a faint ripple in a detector's signal, what was the probability that two black holes had just merged a billion light-years away, versus the probability that the instrument had glitched? The framework that let them assign a credible answer, and quantify exactly how uncertain that answer remained, is the same one you will build in this section.

When a clinical trial estimates a drug's effect without quantifying uncertainty, regulators cannot distinguish a genuine benefit from statistical noise, and patients bear the consequences of that ambiguity. The Bayesian workflow exists to make uncertainty explicit and actionable, turning vague confidence into calibrated probability statements that support real decisions.

Every Bayesian model consists of three components. The prior \(p(\theta)\) encodes what you know (or assume) about the parameters \(\theta\) before seeing data. The likelihood \(p(D \mid \theta)\) describes the probability of observing data \(D\) given parameter values. And the posterior \(p(\theta \mid D)\) is what you learn from combining the two through Bayes' theorem:

The posterior is a complete probability distribution over parameter values, conditioned on observed data. It tells you not just the best estimate of each parameter but the full range of plausible values and their relative credibilities. Point estimates alone can mislead when data are noisy or scarce; the posterior width directly quantifies how much the data constrain your hypothesis. The mechanism is multiplicative: at each candidate parameter value, you multiply the prior probability by the likelihood of the observed data under that value. The normalized product gives the posterior weight. Use posterior inference when you need calibrated uncertainty bands, want to incorporate prior knowledge from earlier experiments, or plan to compare competing models on the same data. Prefer maximum likelihood when the dataset is large enough that the prior's contribution is negligible and you need only a point prediction.

The Posterior in One Equation

The denominator \(p(D) = \int p(D \mid \theta) \, p(\theta) \, d\theta\) is the marginal likelihood (or evidence), a normalizing constant that ensures the posterior integrates to one. Computing this integral is the central computational challenge of Bayesian inference, and it is what motivates the Markov chain Monte Carlo (MCMC) and variational methods we cover in Section 32.2. In short: the posterior is just a weighted compromise between what you believed before (the prior) and what the data tell you (the likelihood), and every computational method in this chapter exists to evaluate that compromise efficiently.

This framework turns competing scientific hypotheses into competing models. If a chemical reaction rate varies with temperature, one hypothesis posits an Arrhenius relationship (an exponential dependence on inverse temperature, derived from activation energy theory) while another posits a power law; each defines its own likelihood function. You can fit both to the same data, compare them quantitatively (Section 32.3), and determine which the data support more strongly.

2. Choosing Priors

Prior selection is the most contested part of Bayesian analysis, and for good reason: the prior directly affects the posterior, especially when data are scarce. Three broad strategies cover most scientific applications.

Uninformative (Flat) Priors

A flat prior \(p(\theta) \propto 1\) says "all parameter values are equally plausible." This sounds objective, but it is not. On the real line, a flat prior is improper (it does not integrate to a finite value), which can cause computational problems. On a bounded interval, a uniform prior is proper but still embeds assumptions: a uniform prior on a rate constant \(k \in [0, 100]\) implies that \(k = 99\) is just as plausible as \(k = 0.01\), which no chemist would endorse. Flat priors are useful as baselines but rarely represent genuine scientific indifference.

Weakly Informative Priors

The recommendation from Gelman et al. (2020) is to use weakly informative priors that rule out clearly implausible values without strongly favoring any particular region. For a parameter you expect to be "around 1, give or take an order of magnitude," a \(\text{Normal}(0, 10)\) or \(\text{HalfNormal}(10)\) prior on the appropriate scale (log, logit (the log-odds transform, used for parameters bounded between 0 and 1), etc.) is weakly informative. The idea is to let the data dominate in regions where you have observations, while preventing the sampler from exploring physically meaningless territory.

Domain-Informed Priors

For discovery applications, the most valuable priors are those that encode genuine domain knowledge. If previous experiments measured a binding affinity of \(K_d = 50 \pm 15\) nM, using \(\text{Normal}(50, 15)\) as a prior for \(K_d\) in a new experiment on a related compound is not subjective bias; it is efficient use of accumulated scientific knowledge. The posterior from one experiment becomes the prior for the next, making Bayesian updating a natural model for the scientific literature.

Let us build a concrete example. A materials scientist is studying how the yield strength \(\sigma_y\) of an alloy depends on grain size \(d\). The Hall-Petch relationship predicts \(\sigma_y = \sigma_0 + k_{\text{HP}} \, d^{-1/2}\). We encode this as a Bayesian model in PyMC, a Python library for specifying and fitting probabilistic models using automatic posterior sampling:

import pymc as pm
import numpy as np
import arviz as az

# Simulated grain-size data (micrometers) and yield strength (MPa)
rng = np.random.default_rng(42)
d_obs = rng.uniform(5, 100, size=30)           # grain sizes
sigma_0_true, k_hp_true = 50.0, 500.0          # true parameters
sigma_y_obs = sigma_0_true + k_hp_true * d_obs**(-0.5) + rng.normal(0, 8, size=30)

with pm.Model() as hall_petch_model:
    # Weakly informative priors grounded in physical knowledge
    sigma_0 = pm.Normal("sigma_0", mu=0, sigma=100)     # friction stress (MPa)
    k_hp = pm.HalfNormal("k_hp", sigma=1000)            # Hall-Petch slope (positive)
    noise = pm.HalfNormal("noise", sigma=20)             # measurement noise

    # Deterministic Hall-Petch relationship
    mu = sigma_0 + k_hp * d_obs**(-0.5)

    # Likelihood: observed yield strengths
    y_obs = pm.Normal("y_obs", mu=mu, sigma=noise, observed=sigma_y_obs)

print(hall_petch_model)

3. Likelihoods for Common Scientific Data

With priors encoding what we believe before seeing data, the next question is how to formalize the connection between those beliefs and the measurements we actually collect.

The likelihood function connects your model's parameters to the data you actually observe. Choosing the right likelihood is as important as choosing the right prior, and the choice depends on the data-generating process. Table 32.1 summarizes the most common likelihood functions for scientific data types.

A frequent mistake is defaulting to a Normal likelihood for all continuous data. If your measurements are strictly positive (concentrations, intensities, durations), a Lognormal or Gamma likelihood respects that constraint and often fits better. If your data have occasional large outliers (common in spectroscopy or environmental monitoring), a Student-t likelihood with low degrees of freedom \(\nu\) provides robustness that a Normal likelihood cannot. The choice of likelihood is itself a scientific hypothesis about how noise enters your measurement process.

4. Prior Predictive Checks

Before fitting the model to real data, we should check whether the prior produces reasonable data. A prior predictive check samples parameters from the prior, generates synthetic data, and asks: "Does this look anything like data I might actually observe?" If your prior predictive distribution includes yield strengths of \(-10{,}000\) MPa or \(+10^{12}\) MPa, the prior is too vague. If it concentrates all mass within \(\pm 1\) MPa of a particular value, it is too informative and will overwhelm the data.

with hall_petch_model:
    prior_pred = pm.sample_prior_predictive(samples=500, random_seed=42)

# Inspect the prior predictive distribution
prior_y = az.extract(prior_pred.prior_predictive)["y_obs"].values
print(f"Prior predictive range: [{prior_y.min():.0f}, {prior_y.max():.0f}] MPa")
print(f"Prior predictive mean:  {prior_y.mean():.0f} MPa")
print(f"Prior predictive std:   {prior_y.std():.0f} MPa")

# Visualize with ArviZ
az.plot_ppc(prior_pred, group="prior", num_pp_samples=100)

The prior predictive check is the Bayesian analogue of a unit test. It catches modeling errors before you spend computational resources on inference. Gelman et al. (2020) recommend it as the first step in every Bayesian analysis, and modern tools make it nearly free: sample_prior_predictive requires no MCMC and runs in seconds. (A full MCMC posterior takes minutes to hours; a prior predictive check takes under one second on the same model, yet it catches the majority of specification errors before you spend any of that compute.)

5. The Iterative Model-Building Cycle

The Bayesian workflow is not "specify model, press run, report results." It is an iterative cycle of model building, checking, and revision. Gelman et al. (2020) describe this cycle explicitly. Figure 32.1 illustrates the six steps and the feedback loop that drives model improvement. Figure 32.1.1 illustrates the iterative Bayesian model-building cycle.

1. Specify the model. Choose priors and likelihood based on domain knowledge and the data-generating process.
2. Prior predictive check. Simulate data from the prior; verify that it covers the plausible range without extending into nonsense.
3. Fit the model. Run MCMC or variational inference (VI) to approximate the posterior (covered in Section 32.2).
4. Diagnose the fit. Check convergence (R-hat, a statistic that compares within-chain and between-chain variance to detect non-convergence, effective sample size (ESS), divergences, where a divergence is a numerical failure in the sampler that signals the posterior geometry is too difficult for the current configuration) to ensure the posterior samples are trustworthy.
5. Posterior predictive check. Simulate data from the fitted model; compare to the observed data. Systematic mismatches reveal model misspecification.
6. Revise and iterate. If the model fails checks, diagnose why (wrong likelihood? missing covariate? prior too tight?) and return to step 1.

Let us complete the cycle for our Hall-Petch model by running the fit and performing posterior predictive checks:

with hall_petch_model:
    # Step 3: Fit the model with NUTS (default in PyMC)
    trace = pm.sample(2000, tune=1000, chains=4, random_seed=42)

    # Step 4: Diagnose convergence
    summary = az.summary(trace, var_names=["sigma_0", "k_hp", "noise"])
    print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat", "ess_bulk"]])

    # Step 5: Posterior predictive check
    post_pred = pm.sample_posterior_predictive(trace, random_seed=42)

# Visualize the posterior predictive check
az.plot_ppc(post_pred, num_pp_samples=100)

The posterior predictive check is widely regarded as one of the most powerful diagnostics in the Bayesian toolkit. If the model is well-specified, data simulated from the posterior should be statistically indistinguishable from the observed data. Systematic deviations (for example, the model consistently underestimates yield strength at very small grain sizes) point to specific model failures and suggest specific fixes (in this case, adding an inverse Hall-Petch correction term for nanoscale grains).

6. From Single Models to Model Comparison

The workflow described above applies to a single model. But scientific discovery often involves comparing competing hypotheses, not just fitting one. Does the data favor an Arrhenius model or a power-law model? Is a linear dose-response adequate, or does a sigmoidal model fit better? The Bayesian framework handles model comparison naturally through the marginal likelihood, leave-one-out cross-validation (LOO-CV), and the widely applicable information criterion (WAIC), which we develop in Section 32.3.

Model comparison in the Bayesian framework is not about finding the model that "fits best" in the sense of minimizing residuals. A model with enough parameters can fit any dataset perfectly, but that does not make it scientifically useful. Bayesian model comparison penalizes complexity through the marginal likelihood. A simpler model that explains the data almost as well as a complex one scores a higher marginal likelihood, because its prior concentrates more mass on the parameter values that actually generated the data. This built-in Occam's razor (the principle that simpler explanations should be preferred when they account for the data equally well) is one of the strongest arguments for the Bayesian approach in scientific discovery.

7. When Not to Use Bayesian Methods

Model comparison gives the Bayesian framework its scientific teeth, but that power comes with real costs in computation and modeling effort that make it the wrong choice in some situations.

Bayesian methods are not always the right tool. They are computationally expensive compared to frequentist alternatives (though this gap is shrinking). They require explicit model specification, which can be difficult for complex systems where the data-generating process is poorly understood. And they can give misleading results when the model is misspecified and the analyst does not perform posterior predictive checks.

Concrete cases where frequentist or machine learning approaches may be preferable:

• Pure prediction with no need for uncertainty. If you only need a point prediction and do not care about confidence intervals, gradient-boosted trees or neural networks are faster and often more accurate. See Chapter 26.
• Very large datasets with simple models. When \(n > 10^6\) and the model is a standard regression, the posterior is essentially determined by the data, and the prior is irrelevant. Maximum likelihood gives the same answer much faster.
• Exploratory analysis before model specification. The exploratory discovery techniques of Chapter 25 are better for initial data understanding. Bayesian modeling comes after you have a hypothesis to test.

The Bayesian workflow shines when uncertainty quantification matters, when you have genuine prior information, when you need to compare competing hypotheses, and when the dataset is small enough that the prior has real influence. These conditions describe a large fraction of scientific research.

Exercises

1. (Conceptual) A pharmaceutical researcher sets a prior of \(\text{Normal}(0, 0.1)\) on a drug effect size, arguing that "most drugs have small effects." A colleague objects that this prior is too informative and will bias the posterior toward zero. Under what conditions is the colleague right? Design a simulation to quantify how much the prior influences the posterior as a function of sample size \(n \in \{10, 50, 200, 1000\}\).
2. (Coding) Extend the Hall-Petch model to include an inverse Hall-Petch correction term: \(\sigma_y = \sigma_0 + k_{\text{HP}} \, d^{-1/2} - k_{\text{IHP}} \, d^{-1}\), where \(k_{\text{IHP}} \geq 0\). Implement this in PyMC, run prior predictive checks, fit the model, and compare the posterior predictive distributions of the two models. Does the data prefer the extended model?
3. (Analysis) Take the Hall-Petch model and systematically vary the prior on \(k_{\text{HP}}\): try \(\text{HalfNormal}(10)\), \(\text{HalfNormal}(100)\), \(\text{HalfNormal}(1000)\), and \(\text{HalfNormal}(10000)\). For each prior, compute the posterior mean, 94% HDI width, and prior predictive range. At what point does the prior stop influencing the posterior? Relate this to the concept of "prior sensitivity analysis" in Bayesian workflow.