STAT 238 - Bayesian Statistics Lecture Twenty

Spring 2026, UC Berkeley

Bayesian Inference in a High-Dimensional Linear Regression Model¶

We studied the following model in the last lecture.

We have a response variable $y$ and a single covariate $x$. In our example, $y$ denotes weekly earnings and $x$ denotes years of experience. The covariate $x$ takes the values $0, 1, \dots, m$ for some fixed integer $m$.

Our data is $(x_i, y_i), i = 1, \dots, n$. The model is:

where $\relu(u) := \max(u, 0)$, and $\epsilon_i \overset{\text{i.i.d}}{\sim} N(0, \sigma^2)$. Note we did not include $\relu(x_i - m)$ because it always equals 0.

The model (1) can be rewritten in the usual regression way as:

Here $X$ is the $n \times (m+1)$ matrix with columns $1, x, \relu(x-j)$ for $j = 1, \dots, m-1$, where $x$ denotes observed values of the experience variable.

The usual least squares analysis coincides with Bayesian inference using the prior:

for $C \rightarrow \infty$. In Problem 7 of Homework 3, it is proved that the posterior for the above prior is:

where $\hat{\beta} = (X^T X)^{-1} X^T y$ is the least squares estimator. If we marginalize $\sigma$ and write the posterior distribution of $\beta$, we will get a $t$-distribution. The posterior distribution for $\sigma$ above can be written in terms of the Gamma (or chi-squared) distribution (this is problem 7(d) in Homework 3) as:

Because least squares does not give sensible results in this example, we change the prior in (3) to:

Therefore, we are bringing in a new parameter $\tau$ which controls the scale of $\beta_2, \dots, \beta_m$. For now, treat $\tau$ as fixed. So the prior on $\beta$ and $\sigma$ is now:

where $Q$ is the $(m+1) \times (m+1)$ diagonal matrix with diagonal entries $C, C, \tau^2, \dots, \tau^2$.

The likelihood is unchanged:

So the posterior for $\beta, \sigma$ is:

Using the completing the square formula:

where

one can then show that

Note that when $Q \rightarrow \infty$ (i.e., when $C \rightarrow \infty$ and $\tau \rightarrow \infty$), this reverts to$N((X^T X)^{-1} X^T y, \sigma^2 (X^T X)^{-1})$. To get the posterior of $\sigma$ given $\tau$, we simply integrate $\beta$ from the joint posterior to obtain:

Integrating $\beta$ from the joint posterior gives the posterior of $\gamma, \sigma$:

Simplifying by using $\det Q \propto (\tau^2)^{m-1}$ and $Q^{-1} \approx J/\tau^2$ where $J$ is the $(m+1) \times (m+1)$ diagonal matrix with diagonals $0, 0, 1, \dots, 1$, we get

This distribution is not easy to handle because of the presence of the term $\sigma^2/\tau^2$ inside the determinant as well as inside the inverse (in the exponent). One trick to simplify this is to reparametrize and assume $\tau = \sigma \gamma$. With this, the above conditional density becomes

Ignoring terms above which do not depend on $\gamma$, we get

The right hand side above is actually the pdf of an inverse gamma density (see Inverse-gamma distribution). This can be seen by converting it into the density of $1/\sigma^2$:

Thus

Thus when $\gamma$ is fixed, we can perform inference on $\beta, \sigma$ using the above closed form formulae. Since $\gamma$ is also unknown, we can place a prior $p(\gamma)$ on it, and then calculate its posterior.

Finally, we can marginalize $\sigma$ to obtain the posterior of $\gamma$ alone as follows:

Letting

we get

By the change of variable $\sigma = s \sqrt{A}$, we obtain

We usually choose $p(\gamma)$ as:

for two fixed values $\text{low}$ and $\text{high}$. These could be the values of $\gamma$ which lead to the posterior mean (which coincides with ridge regression) leading to underfitting and overfitting respectively.

Inference can be carried out by first taking a grid of $\gamma$ values and computing the above posterior (on the logarithmic scale) at the grid points. This posterior can be used to obtain posterior samples of $\gamma$. For each sample of $\gamma$, we can then sample $\sigma$ using the distribution (18). Given samples from both $\gamma$ and $\sigma$, we can then sample $\beta$ using (17).

Induced Prior for the Regression Function¶

Our regression function is being modeled as:

where $\beta_0, \beta_1, \dots, \beta_m$ are all independent with

This can also be seen as a prior on the regression function $f$ more directly. For every $0 \leq u_1 < \dots < u_k < \infty$, the joint distribution of $(f(u_1), \dots, f(u_k))$ induced by (22) will be multivariate normal. This implies that $(f(u), u \geq 0)$ is a Gaussian Process. The description of the GP can be done in terms of the mean function $\E f(u)$ and the covariance kernel $\text{Cov}(f(u), f(v))$.

The mean function is clearly zero (because each $\E \beta_j = 0$). The covariance kernel is given by:

where $k = \lfloor u \wedge v \rfloor$ (and $u \wedge v := \min(u, v)$).

Suppose now we use the approximation

Then the covariance kernel becomes:

It turns out that the right hand side above is precisely the covariance kernel of:

where $\beta_0, \beta_1 \overset{\text{i.i.d}}{\sim} N(0, C)$ and $I_x$ is Integrated Brownian Motion on $[0, \infty)$. We will revisit this in the next lecture.