Bayesian Linear Regression

Definition

Bayesian linear regression applies Bayes's theorem to the normal linear model $y = X\beta + \varepsilon$ , $\varepsilon \sim \mathcal{N}(0, \sigma^2 I)$ , combining the likelihood with a prior distribution on $(\beta, \sigma)$ to produce a posterior distribution over the regression coefficients. The form of the posterior depends critically on what is assumed about the error variance(s) and on the choice of prior. When a second sample shares the same $\beta$ but has an independent unknown variance, the posterior is a product of two multivariate t-distributions — the "double-t" distribution (Tiao-Zellner 1964).

Key Ideas

Posterior = prior × likelihood: Bayes's theorem combines $p(\beta,\sigma)$ with $\ell(\beta,\sigma|y)$ to give $p(\beta,\sigma|y) \propto p(\beta,\sigma)\ell(\beta,\sigma|y)$ .
Noninformative (Jeffreys-Savage) prior: When little is known a priori, set $p(\beta) \propto 1$ and $p(\log\sigma) \propto 1$ (i.e., $p(\sigma) \propto 1/\sigma$ ). The marginal posterior of $\beta$ is then a multivariate t-distribution.
Conjugate (Normal-Inverse-Gamma) prior: A normal prior on $\beta$ and an inverse-gamma prior on $\sigma^2$ yield a multivariate t marginal posterior with closed-form parameters — the standard textbook result.
Heterogeneous variances → double-t: When two samples share $\beta$ but have independent unknown variances $\sigma_1^2$ and $\sigma_2^2$ , independent noninformative priors on $\log\sigma_1$ and $\log\sigma_2$ yield a product-of-t posterior that does not reduce to a single t-distribution.
Precision-weighted pooling: The leading-term posterior mean under the double-t is $\bar\beta = D^{-1}(M_1\hat\beta_1 + M_2\hat\beta_2)$ where $D = M_1 + M_2$ and $M_i = Z_i/s_i^2$ , the precision-weighted average of the two ordinary least squares (OLS) estimates. This is the Bayesian analogue of the Theil (1963) mixed estimator.

How It Works

Noninformative Prior (Single Sample)

The likelihood for the regression model $y = X\beta + \varepsilon$ with $\varepsilon \sim \mathcal{N}(0,\sigma^2 I)$ is:

$\ell(\beta,\sigma|y) = \left(\frac{1}{\sigma\sqrt{2\pi}}\right)^T \exp\!\left\{-\frac{1}{2\sigma^2}\left[\nu s^2 + Q(\beta,\hat\beta,Z)\right]\right\}$

where $Z = X'X$ , $\hat\beta = Z^{-1}X'y$ , $\nu = T-p$ , $s^2 = (y-X\hat\beta)'(y-X\hat\beta)/\nu$ , and $Q(\beta,\eta,A) \equiv (\beta-\eta)'A(\beta-\eta)$ .

Under $p(\beta) \propto 1$ , $p(\sigma) \propto 1/\sigma$ , the joint posterior is:

$p(\beta,\sigma|y) \propto \sigma^{-(T+1)}\exp\!\left\{-\frac{1}{2\sigma^2}\left[\nu s^2 + Q(\beta,\hat\beta,Z)\right]\right\}$

Integrating out $\sigma$ gives the marginal posterior of $\beta$ :

$p(\beta|y) \propto \left\{1 + \frac{Q(\beta,\hat\beta,Z)}{\nu s^2}\right\}^{-\frac{1}{2}(\nu+p)}$

This is a multivariate t-distribution with $\nu = T-p$ degrees of freedom, mean $\hat\beta$ , and precision matrix $Z/s^2$ . The marginal posterior of any single element $\beta^i$ is univariate t with $T-p$ degrees of freedom.

Conjugate (Normal-Inverse-Gamma) Prior

With $\beta|\sigma^2 \sim \mathcal{N}(\beta_0, \sigma^2 V_0)$ and $\sigma^2 \sim \text{IG}(a_0/2, b_0/2)$ , the posterior is:

$\beta|\sigma^2, y \sim \mathcal{N}(\beta_n, \sigma^2 V_n), \quad \sigma^2|y \sim \text{IG}(a_n/2, b_n/2)$

with $V_n^{-1} = V_0^{-1} + X'X$ , $\beta_n = V_n(V_0^{-1}\beta_0 + X'y)$ , $a_n = a_0 + T$ , $b_n = b_0 + (y-X\beta_n)'(y-X\beta_n) + (\beta_n-\beta_0)'V_0^{-1}(\beta_n-\beta_0)$ . Marginalizing $\sigma^2$ : $\beta|y \sim t_{a_n}(\beta_n, b_n V_n / a_n)$ .

Double-t Distribution (Heterogeneous Variances)

When two samples $(y_1, X_1, T_1)$ and $(y_2, X_2, T_2)$ share the same $\beta$ but have independent unknown variances $\sigma_1^2$ and $\sigma_2^2$ , independent log-uniform priors on $\sigma_1$ and $\sigma_2$ yield:

$p(\beta|y_1,y_2) = k^{-1} f(\beta), \quad f(\beta) = \left\{1 + \frac{Q_1}{\nu_1}\right\}^{-\frac{1}{2}(\nu_1+p)} \left\{1 + \frac{Q_2}{\nu_2}\right\}^{-\frac{1}{2}(\nu_2+p)}$

where $Q_i = Q(\beta,\hat\beta_i,M_i)$ with $M_i = Z_i/s_i^2$ , and $k = \int f(\beta)\,d\beta$ is a p-dimensional integral without a closed form. This is the multivariate double-t distribution (Tiao-Zellner 1964, eq. 2.12) — the multivariate extension of the Behrens-Fisher problem.

Asymptotic expansion. Expanding each t-factor as $\{1 + Q_i/\nu_i\}^{-\frac{1}{2}(\nu_i+p)} = e^{-\frac{1}{2}Q_i}\sum_{j=0}^\infty p_j\nu_i^{-j}$ and multiplying, the posterior becomes:

$p(\beta|y_1,y_2) = W^{-1}\frac{|D|^{\frac{1}{2}}}{(2\pi)^{\frac{1}{2}p}}\exp\!\left\{-\tfrac{1}{2}Q(\beta,\bar\beta,D)\right\}\sum_{i=0}^\infty\sum_{j=0}^\infty d_{ij}\nu_1^{-i}\nu_2^{-j}$

with $D = M_1 + M_2$ , $\bar\beta = D^{-1}(M_1\hat\beta_1 + M_2\hat\beta_2)$ . The leading term ( $d_{00}=1$ ) is $\mathcal{N}(\bar\beta, D^{-1})$ ; higher-order terms give finite-sample corrections via mixed cumulants of $Q_1$ and $Q_2$ . As $\nu_1,\nu_2\to\infty$ , the posterior converges to this normal exactly.

Behrens-Fisher connection. For $p=1$ , the univariate double-t agrees with Fisher's (1961) fiducial solution to the problem of combining two normal samples with unequal unknown variances.

Theil mixed estimator as a limiting case. If $\sigma_1^2$ is treated as known and $\nu_1\to\infty$ , the first-sample posterior collapses to a degenerate normal and the posterior mean converges to Theil's (1963) mixed estimator $\bar\beta = D^{-1}(M_1\hat\beta_1 + M_2\hat\beta_2)$ . The Bayesian framework thus provides the full posterior distribution around what Theil obtained as a point estimate.

Why It Matters

The noninformative-prior result establishes the multivariate t-distribution as the fundamental posterior for regression coefficients when variance is unknown, replacing classical interval estimates with a complete probabilistic description. The double-t generalizes this to heterogeneous-variance settings that arise naturally whenever data come from sources with different noise levels — manufacturing process comparisons, cross-country regressions with country-specific volatility, or panel data with heteroskedastic equations. The asymptotic expansion technique developed by Tiao-Zellner (1964) anticipates modern Laplace approximation methods; the underlying intractability they identified is now resolved routinely by Markov Chain Monte Carlo (MCMC).

Open Questions

For large p, the double-t normalizing constant is still challenging outside MCMC — variational and Laplace approximations may not be accurate.
The double-t assumes independence of $\sigma_1^2$ and $\sigma_2^2$ ; correlated scale parameters (as in full Seemingly Unrelated Regression, SUR) require the full Zellner (1971) / Chib-Greenberg (1995) treatment.