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).
Bayesian Linear Model with Inequality Constraints
A Python implementation of Geweke's (1995) Bayesian linear regression with
linear inequality constraints a ≤ Dβ ≤ w. Two key
algorithms: the GHK probability simulator — which estimates the
posterior probability that all constraints hold — and a
component-wise constrained Gibbs sampler that draws from the
truncated posterior regardless of how small that probability is. Applied to the
UCI Automobile Imports dataset (195 cars): economic sign restrictions on engine
size, horsepower, and fuel economy are imposed; bore and stroke OLS estimates
violate the constraints due to multicollinearity, making
p2|1 ≈ 4 × 10−4 and rejection
sampling completely impractical. Includes GHK vs. analytical validation, a
hard-constraint stress test, and a 95% credible-interval coverage simulation.
Cross-checked against HMC (NumPyro NUTS) — which agrees on sign constraints but
struggles on a tight polytope where the truncated-Gibbs sails through — and R
(TruncatedNormal, tmvtnorm, bain).
Bayesian Student-t Linear Model
A Python implementation of Geweke's (1993) Bayesian linear model with
i.i.d. Student-t errors, estimated via a four-block Gibbs sampler.
The scale-mixture-of-normals representation augments each observation with a
latent precision weight λi: outliers receive
low weights, limiting their leverage on β automatically.
The degrees-of-freedom parameter ν is drawn as a discrete categorical variable
on a grid from 2 to 40, avoiding numerical instability near the boundary.
Applied to synthetic contaminated data (recovers ν = 4) and the Nelson-Plosser
(1982) 14 annual U.S. macroeconomic series (1909–1970): low posterior ν for
Real GNP, Industrial Production, and Employment identifies Depression- and
WWII-era shocks as heavy-tailed outliers rather than evidence of a stochastic
trend. Cross-checked against PyMC (continuous ν via NUTS) and R
(hett::tlm Student-t ML, MASS::rlm robust regression, OLS baseline).
Bayesian Mixture-of-Normals Regression
A Python implementation of a linear regression with a flexible
mixture-of-normals error, estimated via a from-scratch 4-block Gibbs
sampler (component indicators → GLS β → Normal-Inverse-Gamma components → Dirichlet
weights). The continuous sibling of Geweke & Keane (1997): rather than the single-Normal
disturbance of ordinary regression, the error follows
εt ~ Σj pj N(μj, σj²),
which can approximate any density — capturing the skew and fat tails a single Normal cannot.
The mixture is recentred to mean 0 each sweep so β stays the usual mean function and the
mixture describes only the disturbance shape; models are compared by WAIC. Applied to PSID
log-earnings (Geweke 2005; OLS residual kurtosis ≈ 23): the mixture beats the Normal by
≈1,580 WAIC, resolving the error into an ≈86% majority component plus a wide ≈14%
low-earnings left tail. Cross-checked in R (flexmix BIC, bayesm's
rnmixGibbs) and complemented by a PyMC notebook that scores the Normal-vs-mixture
question by Bayes factor (marginal likelihood via SMC), confirming the WAIC
verdict while flagging its prior sensitivity. Distinct from the coefficient-mixture
examples, which flex heterogeneity across units rather than the error in one equation.
Bayesian Multivariate Regression
Conjugate Normal–Inverse-Wishart posterior for the multivariate regression model
(shared design matrix across equations). Because the posterior is exact,
draws are IID — no MCMC, no burn-in. Validated on simulated data, then
applied to a tuna demand system (7 brands): B is an
elasticity matrix and Σ the cross-brand error correlation. Three engines
agree — the from-scratch mvr_conjugate.py, R's bayesm rmultireg, and a PyMC/NUTS fit with an LKJ prior on Σ — recovering the same
elasticities and cross-equation correlation.
Bayesian Seemingly Unrelated Regression
A Python implementation of the Bayesian SUR Gibbs sampler, replicating R's
bayesm rsurGibbs. Unlike multivariate regression, each equation
can have different regressors — cross-equation correlation in the error covariance
Ω is what makes joint estimation more efficient than OLS.
Applied to Dominick's Finer Foods tuna scanner data (338 weeks, 7 brands):
own-price and cross-price elasticity estimation. Includes a hierarchical extension
(Chib & Greenberg 1995) that pools brand-level coefficients toward a common mean.
Bayesian Hierarchical Linear Regression
A Python implementation of the Bayesian hierarchical linear model Gibbs sampler,
replicating R's bayesm rhierLinearModel. Estimated on the
IRI cheese scanner panel (88 retailers, 5,555 store-weeks): pooled shrinkage
reduces the cross-store standard deviation of log-price elasticities by 58%
relative to OLS. Includes unit-covariate extension (large vs small markets)
and full R/Python cross-validation.
Bayesian Hierarchical Linear Model — Mixture Heterogeneity
An extension of the hierarchical linear model that replaces the single-normal coefficient
distribution with a mixture of normals
(βi | si=k ~ N(μk, Σk)),
replicating bayesm's rhierLinearMixture (Rossi, Allenby & McCulloch
Ch. 5) via a from-scratch 6-block conjugate Gibbs sampler. Two applications probe whether the
mixture is a flexible heterogeneity density or a segmentation device. On the IRI
cheese panel (88 retailers) a 2-component mixture is only mildly favored by
WAIC, no store commits to a component, and three engines find different modes — the
coefficient distribution is continuous and weakly non-normal, not discrete segments.
On the canonical Rats growth data (Gelfand et al. 1990; 30 rats × 5 ages) the
K=1 BiRats fit (intercept ≈242 g, slope ≈6.2 g/day, positive intercept–slope correlation,
partial-pooling shrinkage) is confirmed by four independent engines — from-scratch Gibbs,
PyMC, bayesm, and lme4 — and the mixture finds no growth types. The
coefficient-mixture sibling of the multinomial-logit
mixture.
Bayesian Hierarchical SUR
A Python implementation of a 4-block Bayesian Hierarchical SUR Gibbs sampler, combining cross-equation error correlation (Ω) with a random-effects hierarchy across stores. Applied to Dominick's Finer Foods orange juice scanner data (83 stores, 11 brands, ≈117 weeks): population-level price elasticities, store-level heterogeneity, store-demographic extensions, and unrestricted cross-price elasticities.
View example →