← All examples

Linear & Regression Models

The Gaussian linear model as a starting point, extended in four directions. On the error: a robust Student-t model that downweights outliers via a scale mixture, and a mixture-of-normals regression whose flexible disturbance captures skew and fat tails a single Normal cannot. On the parameters: a linear model with inequality constraints, drawing from a truncated posterior where economic sign restrictions bind and rejection sampling fails. Across equations: multivariate regression (shared design, exact conjugate draws) and seemingly-unrelated regression (different regressors, cross-equation error correlation). And across units: hierarchical linear and hierarchical SUR models that pool coefficients toward a common mean, plus a mixture-heterogeneity variant that asks whether that population is one continuous distribution or discrete segments. Each is built from scratch; most are then cross-checked against a second engine — PyMC or NumPyro — and against R (bayesm and related packages).

Three words that recur. An elasticity is a percentage response to a percentage change — what a regression coefficient becomes when both sides are measured in logs. An own-price elasticity of −2 says that a 1% price rise costs about 2% of unit sales. Below −1 is called elastic (raising the price loses revenue on net); between −1 and 0, inelastic. A credible interval (CrI) is the Bayesian counterpart of the confidence interval, and it means what people usually, and wrongly, take a confidence interval to mean: a 95% CrI is a range the parameter lies in with 95% posterior probability. An interval straddling zero says the data do not settle even the sign. A parameter is identified when the data could in principle distinguish it from the alternatives. If a shop never changes its price, no quantity of sales records from that shop can reveal how its customers would react to a price change — the coefficient is unidentified, and that is a statement about the information in the data, not about the software. It is exactly the condition that pooling repairs.

What a Gibbs sampler is. Nearly every sampler built from scratch in this section is a Gibbs sampler, described below as 4-block or 6-block, so it is worth saying what a block is. The obstacle in a Bayesian model is that the joint posterior over all the unknowns at once is intractable. But very often each unknown on its own, with everything else held fixed, has a posterior you can write down and draw from exactly. Gibbs sampling lives off that asymmetry: hold the variances fixed and draw the coefficients from their conditional; hold the coefficients fixed and draw the variances from theirs; keep rotating. One pass through every block is a sweep, and the striking result is that the sequence of sweeps converges to the joint posterior even though no single step ever evaluates it. A block is simply a group of unknowns updated together, and fewer, larger blocks generally give a faster-mixing chain — which is why the page headings count them. The price is the one NUTS does not pay: every conditional must be available in closed form, which is what conjugacy delivers and what a single non-conjugate prior takes away.

What pooling buys. Anyone holding data on many similar units — 88 supermarkets, 83 stores, 30 rats — faces a choice between two bad options. Fit each unit separately and every estimate is built from that unit's own thin, noisy slice of data. Fit them all together as though identical and any real differences between them are erased. A hierarchical (multilevel) model does a third thing: it treats the unit-level coefficients as themselves drawn from a common population distribution, and estimates that distribution at the same time as the coefficients. The consequence is shrinkage, or partial pooling — each unit's estimate is pulled toward the population mean by an amount set by how little its own data has to say. A shop whose price barely moves gets an estimate close to the pooled one, because its data cannot argue otherwise; a shop with genuinely variable prices keeps its own answer. The cheese panel makes this concrete. Fitted store by store, the log-price elasticities scatter from −12.65 to +2.86 with a standard deviation of 1.77, and several come out positive — meaning those shops supposedly sell more cheese when they charge more, which is not a finding but an artefact of too little price variation per store. Partial pooling cuts that spread to 0.76, and it does so by removing noise rather than signal.

Where the data comes from. Five of the pages below run on retail scanner panels, and since the same three datasets recur it is worth saying once what they are. A scanner panel is supermarket checkout data, aggregated into a row per store per week per brand: how many units sold, the shelf price, and whether the item was being promoted. The three used here are the IRI cheese panel (88 retailers, 5,555 store-weeks), and two from Dominick’s Finer Foods in Chicago — canned tuna (7 brands, 338 weeks) and orange juice (11 brands across 83 stores, roughly 117 weeks each). Economists reach for this kind of data for two reasons. Prices move for reasons that have little to do with the shopper standing in the aisle — promotion calendars, wholesale costs, competitive responses — which makes the variation closer to an experiment than most observational data gets. And the panel shape, many stores each watched over many weeks, is precisely the structure a hierarchical model is built for. The questions asked of them are ordinary commercial ones: if we raise this price, how much do we sell? Is the answer the same in every store? Do the brands on this shelf take sales from each other or move together?

The tooling, introduced once. This is where the probabilistic programming languages first appear, and they recur through every section that follows, so it is worth saying what they are. PyMC and NumPyro are languages for writing down a model rather than an algorithm: you declare the priors and the likelihood, and the library works out how to sample from the posterior. That is the whole appeal — the from-scratch samplers here are hand-derived for one specific model and break the moment it changes, while a declared model can be edited freely. NumPyro is the same idea built on JAX, which compiles the model and runs it on the GPU, so it is generally the faster of the two on large problems.

What NUTS does. The sampler both libraries reach for by default is the No-U-Turn Sampler, and the reason it displaced the Gibbs and Metropolis samplers written from scratch in these notebooks is worth understanding. Metropolis proposes a random step and accepts or rejects it — in high dimensions almost every guess is bad, so it crawls. Gibbs updates one parameter at a time from its exact conditional, which is efficient but requires those conditionals to be available in closed form, which is exactly what conjugacy buys and what breaks the moment a model steps outside it. Hamiltonian Monte Carlo instead uses the gradient of the posterior to travel along its contours, taking long, informed steps rather than blind ones. It has one awkward tuning knob — how far to travel before stopping — and NUTS removes it by doubling the trajectory until it starts to double back on itself, hence the name. The practical effect is that NUTS samples models with hundreds of correlated parameters that Metropolis could not finish in a lifetime, at the price of needing a differentiable model, which is why discrete parameters have to be integrated out or handled another way.

Reading the diagnostics. Three numbers appear beside almost every posterior in the collection. R-hat compares the variance between chains against the variance within them: chains started in different places should end up describing the same distribution, so R-hat near 1.00 is the pass and anything above about 1.01 means they have not converged on a common answer. ESS, the effective sample size, is how many independent draws the correlated chain is worth — 4,000 samples from a sticky chain can carry the information of 200, and it is the ESS, not the raw count, that sets the precision of a posterior mean. Divergences are unique to Hamiltonian methods: they flag places where the trajectory simulation broke down, usually in a sharply curved region of the posterior, and even a handful means the sampler may be systematically failing to visit part of the distribution — a bias, not just noise. Zero divergences, R-hat at 1.00 and ESS in the hundreds or thousands is the standard applied throughout.