STAT 238 - Bayesian Statistics Lecture Twenty Three

Spring 2026, UC Berkeley

Last lecture: interpolation with Gaussian processes¶

The interpolation problem is: Suppose we are given values of $f$ at points $x_1, \dots, x_n$ in the domain $\Omega$ of $f$. Let $x \in \Omega$ be a new point (i.e., distinct from $x_1, \dots, x_n$). What can we say about $f(x)$?

In the last lecture, we saw how to use Gaussian processes to solve this problem. We model $\{f(x), x \in \Omega\}$ as a Gaussian process with mean zero and covariance function or kernel given by $K(x, x')$ i.e.,

Under this modeling assumption, the answer to the interpolation question is given by the conditional distribution of $f(x)$ given $f(x_1), \dots, f(x_n)$ which is calculated as follows. Note that:

Using the notation $K = (K(x_i, x_j))_{n \times n}$ and $\mathbf{k} = (K(x_i, x))_{n \times 1}$, we can write the conditional distribution of $f(x)$ given $f(x_1), \dots, f(x_n)$ as

Thus the posterior mean (or mode) estimate of $f(x)$ is given by

Regression with Gaussian Processes¶

Today, we will see how to perform regression using Gaussian processes. The key difference between regression and interpolation is that, in regression, the values $f(x_1), \dots, f(x_n)$ are not observed exactly; instead, they are observed with noise.

More precisely, we are given observations $y_1, \dots, y_n$ modeled as

The parameter $\sigma$ controls the level of noise, i.e., how much each observation $y_i$ deviates from the true value $f(x_i)$.

As in the interpolation setting, our goal is to estimate $f(x)$ at a test point $x$. This test point may be different from the observed inputs $x_1, \dots, x_n$, or it may coincide with one of them.

Importantly, because the observations are noisy, it is meaningful to estimate $f(x_i)$ even at observed points. In fact, by combining information from all observations, we can often obtain a better estimate of $f(x_i)$ than the raw observation $y_i$.

Another key difference from interpolation is that the input points $x_1, \dots, x_n$ need not necessarily be distinct. Since each observation contains noise, having repeated measurements at the same input can help improve the overall estimate.

The solution to the regression problem is very similar to that of the interpolation problem with only one difference. The goal is to estimate $f(x)$ at the test point $x$ given the data $(x_1, y_1), \dots, (x_n, y_n)$. We will use the conditional distribution of $f(x)$ given $y_1, \dots, y_n$ (we will assume that $x_1, \dots, x_n, x$ are deterministic). To calculate this conditional distribution, first note that the marginal distribution of $(y_1, \dots, y_n, f(x_*))$ is given by

Using the notation $K = (K(x_i, x_j))_{n \times n}$ and $\mathbf{k} = (K(x_i, x))_{n \times 1}$, we can write the conditional distribution of $f(x)$ given $y_1, \dots, y_n$ as

Thus the posterior mean (or mode) estimate of $f(x)$ is given by

where $y$ is the $n \times 1$ vector with entries $y_1, \dots, y_n$.

So the only difference between (6) and (3) is the presence of the additional $\sigma^2 I_n$ term in case of regression.

Often, we would need to estimate $\sigma$ from the observed data (and also additional hyperparameters present in the kernel $K$). For these, the marginal likelihood of $y_1, \dots, y_n$ is important. This is simply the multivariate normal distribution with mean vector 0 and covariance $K + \sigma^2 I_n$.

To illustrate the calculations, let us take the special case of the Integrated Brownian Motion prior.

Calculations for the IBM prior¶

We take the prior

where $\beta_0, \beta_1$ are i.i.d $N(0, C)$ and $I(x)$ is integrated Brownian Motion.

This is a Gaussian process prior with mean zero and covariance kernel:

where $K_I(u, v)$ is the kernel corresponding to IBM, which is:

Note that the kernel has the unknown parameter $\tau$ (which we write as $\tau = \gamma \sigma$). The constant $C$ is assumed to be large and it is not to be estimated (ideally we want to take $C = +\infty$).

Given data $(x_1, y_1), \dots, (x_n, y_n)$, what is the posterior of $f(x)$? First let us assume that $\tau, \sigma$ are given. The estimate is simply (6). We simplify this expression below. Let $X$ denote the $n \times 2$ matrix with columns 1 and $x_i$ (this is the usual $X$ matrix in the simple linear regression of $y$ on $x$ based on the data $(x_1, y_1), \dots, (x_n, y_n)$).

Note that

where

and $(1, x)$ denotes the row vector with entries 1 and $x$.

Further the $n \times n$ matrix $K$ has the $(i, j)$-th entry:

so that

where $K_I$ is the $n \times n$ matrix with $(i, j)$-th entry $K_I(x_i, x_j)$.

As a result,

This expression depends on the large constant $C$. Direct computation with a large $C$ might make it unstable. It is therefore natural to compute the limit as $C \rightarrow \infty$. Using the Sherman-Morrison-Woodbury identity,

where

Further

We thus get

Write $\tau = \gamma \sigma$. In Fact 1, it is proved that, as $C \rightarrow \infty$ the above converges to:

where

This expression for $\widehat{f(x)}$, which only depends on $\gamma = \tau/\sigma$ and not on $\tau$ and $\sigma$ individually, can be used for computation.

Estimation of $\gamma$ and $\sigma$¶

The hyperparameters $\gamma$ and $\sigma$ also need to be estimated from the observed data $(x_1, y_1), \dots, (x_n, y_n)$ (recall that $\tau = \gamma \sigma$). For this, the marginal likelihood of the data given $\sigma, \gamma$ is important. This is calculated using:

where $K$ is given by (14). In other words,

For $(K + \sigma^2 I_n)^{-1}$, we use (16) and let $C \rightarrow \infty$ to get

Further, using the matrix determinant lemma, we get

So the marginal likelihood of $y$ given $\tau, \sigma$ becomes:

with $\Sigma$ defined in (17).

We now take $\tau = \gamma \sigma$ so that

This gives

Combining this with the prior $\log \sigma \mid \gamma \sim \text{uniform}(-\infty, \infty)$ gives

Finally if we take a prior $p(\gamma)$ for $\gamma$, then the posterior of $\gamma$ becomes:

Proof of the $C \rightarrow \infty$ fact¶