1. Setup: Ordinary Linear Regression as a Bayesian Model

We consider the normal linear regression model:

$y \mid \beta, \sigma, X \sim N(X\beta,\; \sigma^2 I_n)$

  • $y$: $n \times 1$ vector of outcomes
  • $X$: $n \times k$ matrix of predictors
  • $\beta$: $k \times 1$ vector of regression coefficients
  • $\sigma^2$: error variance
  • $I_n$​: $n \times n$ identity matrix

This is the usual ordinary linear regression with:

  • independent errors,
  • constant variance $\sigma^2$,
  • zero conditional correlations given $(X,\beta,\sigma^2)$.

More complicated variance structures are treated later, but here we stick to this simplest case.

2. “Noninformative” Prior for $(\beta, \sigma^2)$

A common default prior is:

  • Uniform in $(\beta, \log \sigma)$, which is equivalent to $p(\beta, \sigma^2 \mid X) \propto \frac{1}{\sigma^2}.$

This is an improper prior (it does not integrate to 1), but under regular conditions it leads to a proper posterior.

When:

  • the sample size $n$ is large relative to the number of parameters $k$,
  • and the design matrix $X$ is well-behaved,

this noninformative prior often gives reasonable results and matches classical regression estimates.

For small $n$ or large $k$, this prior is more dangerous (likelihood is not sharply peaked), and informative or hierarchical priors become more important.

3. Posterior Distribution: Factorization and Forms

The joint posterior is:

$p(\beta, \sigma^2 \mid y, X) \propto p(\beta, \sigma^2 \mid X)\, p(y \mid \beta, \sigma^2, X)$

and can be factored as:

$p(\beta, \sigma^2 \mid y) = p(\beta \mid \sigma^2, y)\, p(\sigma^2 \mid y).$

(Dependence on $X$ is suppressed in notation.)

3.1 Conditional posterior of $\beta$ given $\sigma^2$

Using standard quadratic completion, the conditional posterior is multivariate normal:

$\beta \mid \sigma^2, y \sim N\big(\hat{\beta},\; V_\beta \sigma^2\big)$

where

$\hat{\beta} = (X^T X)^{-1} X^T y$, $V_\beta = (X^T X)^{-1}.$

So:

  • $\hat{\beta}$​ is the usual least squares estimate,
  • $V_\beta \sigma^2$ is the covariance matrix of $\beta$ given $\sigma^2$.

3.2 Marginal posterior of $\sigma^2$

Integrating out $\beta$ yields:

$\sigma^2 \mid y \sim \text{Inv-}\chi^2(n – k,\; s^2)$

where

$s^2 = \frac{1}{n – k} (y – X\hat{\beta})^T (y – X\hat{\beta})$

i.e., the usual residual variance estimate from classical regression.

3.3 Marginal posterior of $\beta$

If you integrate out $\sigma^2$, then

$\beta \mid y \quad\text{is multivariate Student } t_{n-k}$

with center $\hat{\beta}$​ and scale matrix proportional to $s^2 V_\beta$​.

In practice, for simulation-based inference, it’s usually easier to:

  1. Draw $\sigma^2$ from $\text{Inv-}\chi^2(n-k, s^2)$.
  2. Then draw $\beta \mid \sigma^2$ from $N(\hat{\beta}, V_\beta \sigma^2)$.

4. Connection to Classical Regression

  • The classical point estimates are:
    • $\hat{\beta}$​ (OLS estimate),
    • $s$ (residual standard deviation).
  • Under this noninformative prior, the Bayesian posterior mean and the classical estimates coincide:
    • Posterior mean of $\beta$ is $\hat{\beta}$​.
    • Posterior mode or mean of $\sigma^2$ is closely related to $s^2$.
  • The classical standard error for $\beta$ is obtained by setting $\sigma^2 = s^2$ in the normal posterior covariance formula.

So standard regression output (estimates, standard errors) can be interpreted as Bayesian posterior summaries under the noninformative prior.

5. When Is the Posterior Proper?

Because the prior is improper, we must check that the resulting posterior is proper (integrable).

Conditions:

  1. $n > k$: more observations than regression coefficients.
  2. The rank of $X$ is $k$: columns of $X$ are linearly independent (no perfect multicollinearity).

If both hold, the likelihood is informative enough to convert the improper prior into a proper posterior.

Intuitively:

  • You need enough data points to identify all coefficients.
  • No column of $X$ can be an exact linear combination of the others.

6. Efficient Sampling from the Posterior

To simulate $(\beta, \sigma^2)$:

  1. Compute $\hat{\beta}$​, $V_\beta$​, and $s^2$.
  2. Repeat:
    • Draw $\sigma^2 \sim \text{Inv-}\chi^2(n-k, s^2)$.
    • Draw $\beta \mid \sigma^2 \sim N(\hat{\beta}, V_\beta \sigma^2)$.

Using QR factorization for numerical stability

To compute $\hat{\beta}$ and $V_\beta$​ efficiently and stably:

  1. Compute QR factorization:
    • $X = Q R,$
    • where:
      • $Q$: $n \times k$ with orthonormal columns,
      • $R$: $k \times k$ upper triangular.
  2. Solve $R \hat{\beta} = Q^T y$ to get $\hat{\beta}$​.
  3. Then:
    • $R^{-1}$ is easy to compute.
    • $V_\beta = (X^T X)^{-1} = R^{-1} (R^{-1})^T$,
      so $R^{-1}$ is a Cholesky factor of $V_\beta$​.
  4. To simulate $\beta$:
    • Generate $z \sim N(0, I_k)$,
    • Set $\beta = \hat{\beta} + R^{-1} z \cdot \sigma$.

This is important for large datasets or many predictors, where numerical stability and speed matter.

7. Posterior Predictive Distribution for New Data

Suppose we have new predictor values $\tilde{X}$ and we want to predict new outcomes $\tilde{y}$​.

If $(\beta, \sigma^2)$ were known:

$\tilde{y} \mid \beta, \sigma^2, \tilde{X} \sim N(\tilde{X} \beta,\; \sigma^2 I)$

But $(\beta, \sigma^2)$ are unknown and have posterior uncertainty.

7.1 Posterior predictive simulation

To draw $\tilde{y}$​ from the posterior predictive distribution $p(\tilde{y} \mid y)$:

  1. Draw $(\beta, \sigma^2)$ from the posterior $p(\beta, \sigma^2 \mid y)$.
  2. Draw $\tilde{y} \sim N(\tilde{X} \beta,\; \sigma^2 I)$.

Repeat to approximate the full predictive distribution.

This automatically includes:

  • Irreducible noise (σ²),
  • Parameter uncertainty in $\beta$ and $\sigma^2$.

7.2 Analytic form (given σ², then integrating σ² out)

Given $\sigma^2$:

  • Posterior predictive mean:
    • $E(\tilde{y} \mid \sigma^2, y) = \tilde{X}\hat{\beta}$
  • Posterior predictive variance:
    • $\operatorname{var}(\tilde{y} \mid \sigma^2, y) = \sigma^2 \left(I + \tilde{X} V_\beta \tilde{X}^T\right)$

Two components:

  1. $\sigma^2 I$: new observation noise.
  2. $\sigma^2 \tilde{X} V_\beta \tilde{X}^T$: uncertainty from estimating $\beta$.

Averaging over the posterior of $\sigma^2$ (Inverse-χ²) yields a multivariate t posterior predictive distribution:

$\tilde{y} \mid y \sim t_{n-k}\Big(\text{center}=\tilde{X}\hat{\beta},\; \text{scale}^2 = s^2(I + \tilde{X} V_\beta \tilde{X}^T)\Big)$

In practice, simulation is usually easier than using the closed form.

7.3 When $\tilde{X}$ is not fully observed

If some components of $\tilde{X}$ are missing or random, we must model the distribution of $\tilde{X}$ itself and integrate over that uncertainty. This becomes a missing-data / multivariate modeling problem and is more complex.

8. Model Checking and Robustness

Standard regression diagnostics (from classical statistics) can be interpreted naturally as posterior predictive checks. For example:

  • Plot residuals vs. fitted values or vs. predictors.
  • Look for patterns that violate assumptions (nonlinearity, heteroscedasticity, outliers).

Bayesian advantage:

  • Instead of trying to derive sampling distributions for test statistics analytically, we can:
    1. Simulate replicated data $y^{\text{rep}}$ from the posterior predictive distribution.
    2. Compute any statistic (e.g., residual patterns, correlations, test statistics) on both $y$ and $y^{\text{rep}}$.
    3. Compare the observed statistic to its posterior predictive distribution.

This checks both:

  • Statistical significance (is the pattern rare under the model?),
  • Practical significance (how big is the difference?).