STAT 238 - Bayesian Statistics Lecture Fifteen

Spring 2026, UC Berkeley

Linear Regression¶

Our next topic is (multiple) linear regression. Here, we have one response variable $y$ and $m$ covariates $x_1, \dots, x_m$ ($m = 1$ corresponds to simple linear regression). We observe data on $n$ instances or subjects for all these variables: $(y_i, x_{i1}, \dots, x_{im})$ for $i = 1, \dots, n$. The multiple linear regression model (with normal errors) is given by:

The default prior for linear regression is:

for a very large positive $C$. The joint posterior density of $\beta_0, \dots, \beta_m, \sigma$ is then given by

where we use the notation

for the sum of squares.

The posterior over only the coefficient parameters $\beta_0, \dots, \beta_m$ can be obtained by integrating (or marginalizing) the parameter $\sigma$.

where $\hat{\beta}_0, \dots, \hat{\beta}_m$ denote the least squares estimators of $\beta_0, \dots, \beta_m$ (i.e., $(\hat{\beta}_0, \dots, \hat{\beta}_m)$ minimizes $S(\beta_0, \dots, \beta_m)$ over all values of $\beta_0, \dots, \beta_m$).

Our posterior density for $\beta_0, \dots, \beta_m$ is thus:

It turns out that this is a multivariate $t$-density. Here is some basic background on multivariate $t$-densities.

$t$-density¶

The formula for the density corresponding to the $t$-distribution $t_{p}(\mu, \Sigma, \nu)$ is (see (Multivariate t-distribution):

Here:

1. $p$ denotes dimension of the vector $x$ (this is a $p$-variate joint density)

2. $\mu$ is a $p \times 1$ vector called the location

3. $\Sigma$ is a $p \times p$ matrix called the scale matrix

4. $\nu > 0$ denotes the degrees of freedom.

Here is some more information about the $t$-density (5):

1. Connection to the Multivariate Normal Density: The most important term in the formula (5) is $(x - \mu)^T\Sigma^{-1}(x - \mu)$. This exact term also appears in the multivariate normal density. If $X \sim N(\mu, \Sigma)$, then the density of $X$ is given by:

   $$\frac{1}{(2\pi)^{p/2}\sqrt{\det \Sigma}} \exp \left(-\frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu) \right).$$

   This suggests that the $t$-density is closely related to the multivariate normal density. Here is the connection. Suppose $X \sim N_p(\mu, \Sigma)$ and $V \sim \chi^2_{\nu}$ (this is the chi-squared distribution with $\nu$ degrees of freedom) are independent. Then

   $$T := \mu + \frac{X - \mu}{\sqrt{V/\nu}} \sim t_p(\mu, \Sigma, \nu).$$

   (6)

   Thus, in the notation $t_{p}(\mu, \Sigma, \nu)$, $\nu$ denotes degrees of freedom, $p$ denotes dimension, $\mu$ and $\Sigma$ denote the mean vector and covariance matrix of the corresponding normal random vector $X$. For completeness, we include a proof of (6) in Section ??.

2. Individual Components as well as Linear Combinations of Components of $T$ are also $t$-distributed: Suppose $T \sim t_p(\mu, \Sigma, \nu)$ and the components of $T$ are $T_1, \dots, T_p$. Then each individual component $T_j$ is also $t$-distributed. Also every linear combination $a_0 + a_1 T_1 + a_2 T_2 + \dots + a_p T_p$ is also $t$-distributed. To see this, first write

   $$a_0 + a_1 T_1 + \dots + a_p T_p = a_0 + a^T T$$

   where $a$ is the $p \times 1$ vector with components $a_1, \dots, a_p$. Using the formula (6), we can write

   $$a_0 + a^T T = (a_0 + a^T \mu) + \frac{(a_0 + a^TX) - (a_0 + a^T\mu)}{\sqrt{V/\nu}}$$

   Because $a_0 + a^T X \sim N(a_0 + a^T \mu, a^T \Sigma a)$, the same fact (6) applied to this case gives:

   $$a_0 + a^T T \sim t_1(a_0 + a^T \mu, a^T \Sigma a, \nu).$$

   In particular, this implies that for each $j = 1, \dots, p$,

   $$T_j \sim t_1(\mu_j, \Sigma(j, j), \nu)$$

   where $\mu_j$ is the $j$th component of $\mu$ and $\Sigma(j, j)$ is the $(j, j)$th entry of $\Sigma$.

3. When $\nu$ is large, $t$ is very close to normal: This can intuitively be seen by noting that when $\nu$ is large, the term $(x - \mu)^T \Sigma^{-1} (x - \mu)/\nu$ is small so that

   $$1 + \frac{1}{\nu} (x - \mu)^T \Sigma^{-1} (x - \mu) \approx \exp \left(\frac{1}{\nu} (x - \mu)^T \Sigma^{-1} (x - \mu) \right),$$

   where we used the observation that $1 + z \approx e^z$ when $z$ is small. Thus the $t$-density (5) for large $\nu$ becomes approximately:

   $$\exp \left(-\frac{\nu + p}{2 \nu} (x - \mu)^T \Sigma^{-1} (x - \mu) \right) \approx \exp \left(-\frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu) \right)$$

   because $\frac{\nu + p}{\nu} \approx 1$ when $\nu$ is large. This gets us the normal density:

   $$t_{p}(\mu, \Sigma, \nu) \approx N_p(\mu, \Sigma) ~~\text{if $\nu$ is large}.$$

It turns out that (4) is a special case of (5) for some $p, \mu, \Sigma, \nu$. To see this, we need to first rewrite (4) using matrix notation which we do in the next section.

Matrix Notation for Multiple Linear Regression¶

This notation is used not just to write formulae for linear regression, but also in code. For example, the OLS function in statsmodels uses the syntax sm.OLS(y, X).fit() to fit the linear regression model, where $y$ ($n \times 1$ vector) and $X$ ($n \times (m+1)$ matrix) are defined above.

With this notation, one can write the sum of squares $S(\beta_0, \dots, \beta_m)$ as:

There are two important facts about $S(\beta)$:

1. Fact 1: the least squares estimator $\hat{\beta}$ is given by the formula:

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

   (7)

   The proof of (7) is as follows. The gradient of $S(\beta)$ is given by

   $$\begin{align*} \nabla S(\beta) &= \nabla \left[ \|y - X \beta\|^2 \right] \\ &= \nabla \left[ (y - X\beta)^T (y - X \beta) \right] \\ &= \nabla \left[ y^T y - \beta^T X^T y - y^T X \beta + \beta^T X^T X \beta\right] = 2 X^T y - 2 X^T X \beta. \end{align*}$$

   (8)

   Because $\hat{\beta}$ minimizes $S(\beta)$, the gradient should equal zero when $\beta = \hat{\beta}$, and this leads to

   $$\begin{align} X^T (y - X \hat{\beta}) = 0 \implies X^T X \hat{\beta} = X^T y \implies \hat{\beta} = (X^T X)^{-1} X^T y. \end{align}$$

   (9)

2. Fact 2: The following Pythagorean identity holds:

   $$\begin{align} S(\beta) = S(\hat{\beta}) + \|X \beta - X \hat{\beta}\|^2 = S(\hat{\beta}) + (\beta - \hat{\beta})^T X^T X (\beta - \hat{\beta}). \end{align}$$

   (10)

   To prove (10), write

   $$\begin{align*} S(\beta) &= \|y - X \beta\|^2 \\ &= \|y - X\hat{\beta} + X \hat{\beta} - X \beta\|^2 \\ &= \|y - X \hat{\beta}\|^2 + \|X \hat{\beta} - X \beta\|^2 + 2 \left<y - X \hat{\beta}, X \hat{\beta} - X \beta \right>. \end{align*}$$

   (11)

   The cross product is zero (leading to (10)) because:

   $$\begin{align*} \left<y - X \hat{\beta}, X \hat{\beta} - X \beta \right> &= \left(X \hat{\beta} - X \beta \right)^T \left(y - X \hat{\beta} \right) \\ &= \left( \hat{\beta} - \beta \right)^T X^T (y - X \hat{\beta}) = \left( \hat{\beta} - \beta \right)^T \left(X^T y - X^T X \hat{\beta} \right) = 0 \end{align*}$$

   (12)

   where we used (9).

Using (10), we can write the posterior density (4) as

The above formula is a special case of (5) with

or equivalently

We thus have

With the posterior density (14), one can do uncertainty quantification about the parameters $\beta_0, \beta_1, \dots, \beta_m$. One can generate multiple samples from $t_{m+1}(\hat{\beta}, (S(\hat{\beta})/(n-m-1)) (X^T X)^{-1}, n-m-1)$ and plot the resulting fitted values to visualize the uncertainty in the coefficients.

In (14), the quantity $S(\hat{\beta})/(n-m-1)$ is the frequentist unbiased estimator for $\sigma^2$, so we denote it by $\hat{\sigma}^2$:

$\hat{\sigma}$ can also be justified as a Bayesian estimator of $\sigma$ (this will be a question in Homework three). The terminology Residual Standard Error is sometimes used for $\hat{\sigma}$.

With the notation for $\hat{\sigma}$, the posterior (14) becomes:

By one of the facts mentioned about the $t$-distribution, the posterior of each individual $\beta_j$ is also $t$:

where $(X^T X)^{j+1, j+1}$ is the $(j+1)^{th}$ diagonal entry of $(X^T X)^{-1}$ (note that we are using the $(j+1)$th diagonal entry of $X^TX$ because $\beta_j$ is the $(j+1)$th component of $\beta$). Writing this density out, we have

which implies that

This can be used to obtain uncertainty intervals for $\beta_j$. If $t_{n-m-1, \alpha/2}$ is the point beyond which the $t$-distribution (with $n-m-1$ degrees of freedom) assigns probability $\alpha/2$, then

which is same as:

This interval:

is called the $100(1-\alpha)\%$ Bayesian Credible interval for $\beta_j$. It exactly coincides with the frequentist $100(1 - \alpha)\%$ confidence interval for $\beta_j$.

The degrees of freedom corresponding to the $t$-density in (14) is $n - m - 1$ where $n$ is the number of observations, and $m$ is the number of covariates. Thus if $n-m-1$ is large, then the posterior distribution (which is actually $t$) is approximately normal:

In other words, when $n-m-1$ is large, the $t$-density (14) is approximately equal to the $N_{m+1}(\hat{\beta}, \hat{\sigma}^2 (X^T X)^{-1})$. Further, when $n-m-1$ is large, the distribution (16) will be close to the normal distribution $N(\hat{\beta}_j, \hat{\sigma}^2 (X^T X)^{j+1, j+1})$. The quantity $\hat{\sigma} \sqrt{(X^T X)^{j+1, j+1}}$ is known as the standard error corresponding to $\beta_j$.

Proof of (6)¶