cov_mat = outer(x, x, sqexpcov)8 Gaussian Processes
“Uncertainty is the only certainty there is, and knowing how to live with insecurity is the only security.”—John Allen Paulos
A Gaussian Process (GP) defines a probability distribution over functions rather than estimating parameters for a single function. Given data, a GP returns not just a predicted curve but a full posterior distribution over curves, quantifying uncertainty at every point. GPs are used extensively in machine learning for tasks ranging from robotic control to geospatial analysis and Bayesian optimization.
Formally, a GP is a collection of random variables, any finite number of which have a joint Gaussian distribution. A finite collection of \(n\) points drawn from a GP is completely specified by its \(n\)-dimensional mean vector \(\boldsymbol{\mu}\) and covariance matrix \(\boldsymbol{\Sigma}\). We assume the process is indexed by a real variable \(x\in\mathbb{R}\) (e.g., time or space) and has real-valued outputs. The GP is defined by: 1. Mean function \(m(x) = \E{f(x)}\): The expected value of the function at point \(x\). 2. Covariance (Kernel) function \(k(x, x') = \E{(f(x) - m(x))(f(x') - m(x'))}\): A measure of similarity between values at \(x\) and \(x'\).
We denote this as: \[ f(x) \sim \mathcal{GP}(m(x), k(x, x')). \]
Intuitively, the kernel function determines the “shape” and “smoothness” of the functions we expect to see. If \(x\) and \(x'\) are close, \(k(x, x')\) should be high, implying \(f(x)\) and \(f(x')\) are likely similar.
In practice, we often assume a zero mean, \(m(x)=0\), and focus on the covariance kernel. The choice of kernel encodes our prior beliefs about the data. The most common choice is the Squared Exponential (SE) kernel (also known as the Radial Basis Function or RBF):
\[ k(x, x') = \sigma^2 \exp\left(-\frac{(x - x')^2}{2l^2}\right) \]
Here, we have two hyperparameters:
- \(\sigma^2\) (Signal Variance): Controls the vertical amplitude of the function.
- \(l\) (Length Scale): Controls the horizontal “wiggliness.” A large \(l\) implies the function changes slowly (smooth), while a small \(l\) allows for rapid variations.
Observe that \(k(x,x) = \sigma^2\) and \(k(x,x') \rightarrow 0\) as the distance \(|x-x'| \rightarrow \infty\).
We can illustrate a GP with a simulated example. First generate a sequence of 100 input points (indices)
and then define the mean function and the covariance function
The covariance function depends only on the distance between two points, not on their absolute values. The squared exponential kernel is infinitely differentiable, which means that the GP is a very smooth function. The squared exponential kernel is also called the radial basis function (RBF) kernel. The covariance matrix is then defined as
and we can generate a sample from the GP using the mvrnorm function from the MASS package and plot a sample.
Figure 8.1 displays 100 values of a function \(f(x)\) drawn from a GP with zero mean and a squared-exponential kernel at inputs \(x=(0,0.1,0.2,\ldots,10)\). The realisation is smooth, with most values lying between -2 and 2. Because each diagonal element of the covariance matrix equals \(\sigma^2=1\), the marginal variance is one. By properties of the normal distribution, approximately 95 percent of the points of \(Y\) should therefore fall within 1.96 standard deviations of the mean. The mild oscillations arise because values with neighbouring indices are highly correlated.
We can generate a few more samples from the same GP and plot them together
Each finite sample path differs from the next, yet all share a similar range, a comparable number of bumps, and overall smoothness. That’s what it means to have function realizations under a GP prior: \(Y = f(x) \sim \mathcal{GP}(0, k(x, x'))\).
Simulating from the prior shows us the richness of possible functions we can model. However, our goal is not just to generate random curves, but to learn. We want to constrain these possibilities using actual observations.
8.1 Making Predictions with Gaussian Processes
Suppose our observed data consist of \(n\) inputs \(\mathbf{X}=(x_1,\ldots,x_n)^\top\) and outputs \(\mathbf{y}=(y_1,\ldots,y_n)^\top\). We assume these are a realization of a GP. Our goal is to predict the outputs \(\mathbf{y}_*\) at new inputs \(\mathbf{X}_*\).
By definition of a GP, the joint distribution of the observed data \(\mathbf{y}\) and the predictions \(\mathbf{y}_*\) is a multivariate Gaussian:
\[ \begin{bmatrix} \mathbf{y} \\ \mathbf{y}_* \end{bmatrix} \sim \mathcal{N} \left ( \begin{bmatrix} \boldsymbol{\mu} \\ \boldsymbol{\mu}_* \end{bmatrix}, \begin{bmatrix} \mathbf{K} & \mathbf{K}_* \\ \mathbf{K}_*^\top & \mathbf{K}_{**} \end{bmatrix} \right ) \]
where:
- \(\mathbf{K} = k(\mathbf{X}, \mathbf{X})\in \mathbb{R}^{n\times n}\) is the covariance of the training data.
- \(\mathbf{K}_* = k(\mathbf{X}, \mathbf{X}_*)\in \mathbb{R}^{n\times q}\) is the covariance between training and test data.
- \(\mathbf{K}_{**} = k(\mathbf{X}_*, \mathbf{X}_*) \in \mathbb{R}^{q\times q}\) is the covariance of the test data.
- \(\boldsymbol{\mu} = \E{\mathbf{y}}\) and \(\boldsymbol{\mu}_* = \E{\mathbf{y}_*}\).
The beauty of Gaussians is that conditioning on observed data is closed-form. Assuming noiseless observations (that is, \(\mathbf{y} = f(\mathbf{X})\) exactly, with no measurement error), the conditional distribution of \(\mathbf{y}_*\) given \(\mathbf{y}\) is:
\[ \mathbf{y}_* \mid \mathbf{y}, \mathbf{X}, \mathbf{X}_* \sim \mathcal{N}(\boldsymbol{\mu}_{\mathrm{post}}, \boldsymbol{\Sigma}_{\mathrm{post}}) \]
The posterior mean \(\boldsymbol{\mu}_{\mathrm{post}}\) serves as our prediction, and the posterior covariance \(\boldsymbol{\Sigma}_{\mathrm{post}}\) quantifies our uncertainty:
\[ \boldsymbol{\mu}_{\mathrm{post}} = \boldsymbol{\mu}_* + \mathbf{K}_*^\top\mathbf{K}^{-1} (\mathbf{y} - \boldsymbol{\mu}) \tag{8.1}\]
\[ \boldsymbol{\Sigma}_{\mathrm{post}} = \mathbf{K}_{**} - \mathbf{K}_*^\top \mathbf{K}^{-1} \mathbf{K}_* \tag{8.2}\]
Equation 8.1 and Equation 8.2 are standard properties of the multivariate normal distribution (see Chapter 3 and Appendix A). Intuitively, the posterior variance \(\boldsymbol{\Sigma}_{\mathrm{post}}\) is equal to the prior variance \(\mathbf{K}_{**}\) minus a term representing the information gained from the observed data. This structure reflects a fundamental Bayesian principle: data acts to reduce our prior uncertainty (represented by the second term in Equation 8.2).
Example 8.1 (Gaussian Process for \(\sin\) Function) We can use the GP to make predictions about the output values at new inputs \(x_*\). We use \(x\) in the [0,\(2\pi\)] range and \(y\) to be the \(y = \sin(x)\). We start by simulating the observed \(x\)-\(y\) pairs.
The additive term diag(eps, n) corresponds to adding \(\epsilon \mathbf{I}\), which stabilizes the matrix inversion by ensuring strict positive definiteness; in machine learning practice, this is known as ‘jitter’. Now we implement a function that calculates the mean and covariance of the posterior distribution of \(y_*\) given \(Y\).
Now we generate a new set of inputs \(x_*\) and calculate the covariance matrices \(K_*\) and \(K_{**}\).
Notice that we did not add \(\epsilon I\) to \(K_*\) = KX matrix, but we add it to \(K_{**}\) = KXX to ensure the posterior covariance is invertible. We do not add it to \(K_*\) (KX) as it represents cross-covariance, which does not need to be positive definite. Now we can calculate the mean and covariance of the posterior distribution of \(y_*\) given \(Y\).
Now, we can generate a sample from the posterior distribution over \(y_*\), given \(Y\)
Using our convenience function plot_gp we can plot the posterior distribution over \(y_*\), given \(Y\).
Example 8.2 (Gaussian Process for Simulated Data Using MLE) In the previous example, we assumed fixed values for the hyperparameters: \(\sigma^2 = 1\) and \(l^2 = 1\) (i.e., \(2l^2 = 2\)). In real applications, we don’t know these values; we must estimate them from the data.
We use Maximum Likelihood Estimation (MLE) to find the parameters that maximize the probability of observing our data. This parallels the likelihood-to-loss framing in Chapter 11: we write down a (marginal) likelihood for the data and then optimize it, often via its log-likelihood. This section relies on basic matrix operations (inverse, determinant); Appendix A provides a short refresher. If you have not seen gradient-based optimization, the core intuition appears in Chapter 20; here we use it only as a practical tool to fit GP hyperparameters.
The marginal likelihood of the data \(\mathbf{y}\) (integrating out the function values \(f\)) is:
\[ p(\mathbf{y} \mid \mathbf{X}, \sigma, l) = \frac{1}{(2\pi)^{n/2} |\mathbf{K}|^{1/2}} \exp \left ( -\frac{1}{2} \mathbf{y}^\top \mathbf{K}^{-1} \mathbf{y} \right ) \]
where \(\mathbf{K}\) is the covariance matrix computed with hyperparameters \(\sigma\) and \(l\). For numerical stability, we typically maximize the log-likelihood:
\[ \log p(\mathbf{y} \mid \mathbf{X}, \sigma, l) = -\frac{1}{2} \log |\mathbf{K}| - \frac{1}{2} \mathbf{y}^\top \mathbf{K}^{-1} \mathbf{y} - \frac{n}{2} \log 2\pi. \]
This equation encapsulates Occam’s Razor. The term \(-\frac{1}{2} \mathbf{y}^\top \mathbf{K}^{-1} \mathbf{y}\) rewards the model for fitting the data well. The term \(-\frac{1}{2} \log |\mathbf{K}|\) penalizes model complexity; a more flexible kernel (e.g., smaller length scale) leads to a simpler determinant term that exacts a cost. MLE automatically balances these two competing objectives to prevent overfitting.
We can implement a function that calculates the log-likelihood of the data given the hyperparameters \(\sigma\) and \(l\) and use optim function to find the maximum of the log-likelihood function.
## 1.5 2.4The optim function returns the hyperparameters that maximize the log-likelihood function. We can now use those hyperparameters to make predictions about the output values at new inputs \(x_*\).
We can see that our uncertainty is much narrower: the posterior distribution is considerably tighter. This is because we used the observed data to estimate the hyperparameters. We can also see that the posterior mean is closer to the true function \(y = \sin(x)\). Although our initial guess of \(\sigma^2 = 1\) and \(2l^2 = 2\) was not too far off, the model fits the data much better when we use the estimated hyperparameters.
The default optim function uses numerical approximations for derivatives. While convenient, this can be slow and less precise. For GPs, we can calculate the analytical gradients of the log-likelihood with respect to the hyperparameters, significantly speeding up optimization.
Now we can implement a function that calculates the derivative of the log-likelihood function with respect to \(\sigma\) and \(l\).
Now we can use the optim function to find the maximum of the log-likelihood function and provide the derivative function we just implemented.
par1 = optim(c(1,1), fn=loglik, gr=gnlg ,X=X, Y=Y,method="BFGS")$par
l = par1[2]; sigma = par1[1]
print(par1)
## 1.5 2.4The result is the same compared to when we called optim without the derivative function. Even execution time is the same for our small problem. However, at larger scale, the derivative-based optimization algorithm will be much faster.
Furthermore, instead of coding our own derivative functions, we can use an existing package, such as the laGP package, developed by Bobby Gramacy to estimate the hyperparameters. The laGP package uses the same optimization algorithm we used above, but it also provides better selection of the covariance functions and implements approximate GP inference algorithms for large scale problems, when \(n\) becomes large, and inversion of the covariance matrix \(K\) is prohibitively expensive.
library(laGP)
gp = newGP(X, Y, 1, 0, dK = TRUE)
res = mleGP(gp, tmax=20)
l.laGP = sqrt(res$d/2)
print(l.laGP)
## 2.4In the newGP function defines a Gaussian process with square exponential covariance function and assumes \(\sigma^2 = 1\), then mleGP function uses optimization algorithm to maximize the log-likelihood and returns the estimated hyperparameters d = \(2l^2\), we can see that the length scale is close to the one we estimated above. We will use the predplot convenience function to calculate the predictions and plot the data vs fit.
We can see that there is visually no difference between the two fits. Thus, it seems irrelevant whether we keep sigma fixed \(\sigma=1\) or estimate it using MLE. However, in other applications when uncertainty is larger, the choice of \(\sigma\) is important when we use GP for regression and classification tasks. Even for our example, if we ask our model to extrapolate
Extrapolation: Posterior distribution over \(y_*\), given \(Y\).
We can see that outside of the range of the observed data, the model with \(\sigma=1\) is more confident in its predictions.
Now, instead of using GP to fit a known function (\(\sin\)), we will apply it to a real-world dataset. We will use the motorcycle accident dataset from the MASS package. The dataset contains accelerator readings taken through time in a simulated experiment on the efficacy of crash helmets. These data are non-linear and exhibit varying curvature, making them an excellent candidate for GP regression where linear models would fail.
Example 8.3 (Gaussian Process for Motorcycle Accident Data) We first estimate the length scale parameter \(l\) using the laGP package.
Now we plot the data and the fit using the estimated length scale parameter \(l\).
We can see that our model is more confident for time values between 10 and 30. The confidence interval is wider for time values between 0 and 10 and between 30 and 60, and less confident at the end close to the 60 mark. For some reason, the acceleration values were not measured evenly. If we look at the histogram of time values, we can see that there are more data points in the middle of the time range.
The widening of the confidence intervals in regions with fewer data points (e.g., between 30 and 40) is a natural property of the GP; where data are sparse, the model reverts to the prior covariance, resulting in higher uncertainty.
In summary, Gaussian Processes provide a robust and flexible framework for modeling functions where uncertainty is key. By defining a prior over functions and updating it with data, we obtain a posterior distribution that captures both predictions and the confidence in those predictions. The key features of GPs are:
- Non-parametric: GPs can model functions of arbitrary complexity without a fixed number of parameters.
- Data-efficient: They work well with small datasets and provide uncertainty estimates that are valuable for decision-making.
- Versatile: Through the choice of the kernel function, GPs can capture various structures (smoothness, periodicity, etc.) and are used in fields ranging from environmental modeling to hyperparameter optimization in deep learning (“Bayesian Optimization”).
Exercises
Pen-and-Pencil
Exercise 8.1 (Kernel Hyperparameter Interpretation) For the squared exponential kernel \(k(x, x') = \sigma^2 \exp\left(-\frac{(x-x')^2}{2\ell^2}\right)\):
- What happens to the GP when \(\ell \to 0\)?
- What happens when \(\ell \to \infty\)?
- If you double the length scale, by what factor does the “effective range” of correlation increase?
- A meteorologist uses a GP to model temperature as a function of location. She finds \(\ell = 50\) km. Interpret this value.
Exercise 8.2 (True/False: Gaussian Processes)
A Gaussian Process defines a probability distribution over functions.
For a GP, any finite collection of function values has a joint Gaussian distribution.
The posterior variance of a GP at a test point depends on the observed function values \(y\).
The squared exponential kernel produces smooth (infinitely differentiable) sample functions.
Exercise 8.3 (GP Posterior Predictive by Hand) Consider a noiseless GP with mean function \(m(x) = 0\) and squared exponential kernel \(k(x,x') = \sigma^2\exp\left(-\frac{(x-x')^2}{2\ell^2}\right)\) with \(\sigma^2 = 1\) and \(\ell = 1\). You observe two training points:
| x | y |
|---|---|
| 0 | 1 |
| 1 | 2 |
- Build the \(2\times 2\) training covariance matrix \(\mathbf{K}\) by evaluating the kernel at \(x = 0, 1\).
- Invert \(\mathbf{K}\) by hand, using the explicit \(2\times2\) inverse formula \(\begin{pmatrix} a & b \\ b & d\end{pmatrix}^{-1} = \frac{1}{ad-b^2}\begin{pmatrix} d & -b \\ -b & a \end{pmatrix}\).
- Using Equation 8.1 and Equation 8.2, compute the posterior mean \(\mu_*\) and variance \(\Sigma_*\) at the test point \(x_* = 0.5\).
- Without redoing the arithmetic, explain why \(\mu_*\) must lie strictly between \(y_1 = 1\) and \(y_2 = 2\), and why \(\Sigma_*\) must be strictly less than \(\sigma^2 = 1\).
Exercise 8.4 (GP Posterior Mean as Kernel Ridge Regression) Consider a GP with zero mean and kernel \(k(\cdot,\cdot)\), observed with i.i.d. noise \(y_i = f(x_i) + \epsilon_i\), \(\epsilon_i \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0,\sigma_n^2)\).
- Using Equation 8.1 and Equation 8.2 as a template, argue that the noisy-GP posterior mean at a test point \(x_*\) is \[ \mu_* = \mathbf{K}_*^\top(\mathbf{K}+\sigma_n^2\mathbf{I})^{-1}\mathbf y. \]
- Kernel ridge regression fits \(f(x) = \sum_{i=1}^n \alpha_i k(x,x_i)\) by minimizing the penalized least-squares objective \[ J(\boldsymbol\alpha) = \|\mathbf y - \mathbf K\boldsymbol\alpha\|^2 + \lambda\, \boldsymbol\alpha^\top \mathbf K \boldsymbol\alpha, \] where the penalty term \(\boldsymbol\alpha^\top \mathbf K \boldsymbol\alpha\) is the squared norm of \(f\) in the reproducing kernel Hilbert space associated with \(k\). Differentiate \(J\) with respect to \(\boldsymbol\alpha\) and solve for the minimizer \(\hat{\boldsymbol\alpha}\).
- Substitute \(\hat{\boldsymbol\alpha}\) into \(f(x_*) = \sum_i \alpha_i k(x_*,x_i)\) and compare the result to part (a). For what value of \(\lambda\) are the two predictions identical?
- The two methods produce the same point prediction under that correspondence. What does the GP framework provide that ridge regression does not?
Exercise 8.5 (Gradient of the GP Log-Likelihood and a Closed Form for the Signal Variance) The chapter states the gradient of the GP log marginal likelihood (in the collapsible “Mathematical Details: Gradient Derivation” box) but leaves the derivation to the reader. Fill it in, then use it to obtain a closed form for the signal variance. Work throughout with the noiseless log marginal likelihood \[ \mathcal{L}(\theta) = \log p(\mathbf y \mid \mathbf X, \theta) = -\frac{1}{2}\log|\mathbf K| - \frac{1}{2}\mathbf y^\top \mathbf K^{-1}\mathbf y - \frac{n}{2}\log 2\pi, \] where \(\mathbf K\) is built from the squared exponential kernel \(k(x,x') = \sigma^2\exp\left(-(x-x')^2/(2\ell^2)\right)\) with hyperparameters \(\theta\in\{\sigma,\ell\}\). You may use the two matrix calculus identities stated in the chapter, \[ \frac{\partial \mathbf K^{-1}}{\partial \theta} = -\mathbf K^{-1}\frac{\partial \mathbf K}{\partial \theta}\mathbf K^{-1}, \qquad \frac{\partial \log|\mathbf K|}{\partial \theta} = \mathrm{tr}\left(\mathbf K^{-1}\frac{\partial \mathbf K}{\partial \theta}\right). \]
Differentiating \(\mathcal{L}\) term by term, and being explicit about how the first identity enters when you differentiate the quadratic term \(\mathbf y^\top \mathbf K^{-1}\mathbf y\), derive \[ \frac{\partial \mathcal{L}}{\partial \theta} = -\frac{1}{2}\mathrm{tr}\left(\mathbf K^{-1}\frac{\partial \mathbf K}{\partial \theta}\right) + \frac{1}{2}\mathbf y^\top \mathbf K^{-1}\frac{\partial \mathbf K}{\partial \theta}\mathbf K^{-1}\mathbf y. \] Then, writing \(\boldsymbol\alpha = \mathbf K^{-1}\mathbf y\), show that this collapses to the single trace \[ \frac{\partial \mathcal{L}}{\partial \theta} = \frac{1}{2}\mathrm{tr}\left[\left(\boldsymbol\alpha\boldsymbol\alpha^\top - \mathbf K^{-1}\right)\frac{\partial \mathbf K}{\partial \theta}\right]. \]
Positive hyperparameters are usually optimized in an unconstrained space by taking logarithms, as the held-out motorcycle exercise Exercise 8.10 does for its three hyperparameters. Take \(\sigma = e^{u}\) and optimize over \(u = \log\sigma\). Using the chain rule together with the chapter identity \(\partial \mathbf K/\partial\sigma = (2/\sigma)\mathbf K\), show that \(\partial \mathbf K/\partial u = 2\mathbf K\) (the explicit dependence on \(\sigma\) cancels), and hence that \[ \frac{\partial \mathcal{L}}{\partial u} = \mathbf y^\top \mathbf K^{-1}\mathbf y - n. \]
Write \(\mathbf K = \sigma^2 \mathbf R\), where \(\mathbf R\) is the correlation matrix (the same kernel with \(\sigma^2\) set to \(1\)), so that \(\mathbf K^{-1} = \mathbf R^{-1}/\sigma^2\). Setting the derivative from part (b) to zero, show that with \(\ell\) held fixed, the log-likelihood is maximized in closed form at \[ \hat\sigma^2 = \frac{1}{n}\mathbf y^\top \mathbf R^{-1}\mathbf y, \] and explain why this stationary point is the unique maximum, so the signal variance never needs numerical optimization once \(\ell\) is fixed.
Apply the closed form to the two-point dataset of Exercise 8.3 (\(x = (0,1)\), \(y = (1,2)\), \(\ell = 1\)). Reusing the matrix \(\mathbf K^{-1}\) computed there (where \(\sigma = 1\), so \(\mathbf R = \mathbf K\)), compute \(\mathbf y^\top \mathbf R^{-1}\mathbf y\) and \(\hat\sigma^2\). Is \(\hat\sigma\) larger or smaller than the value \(\sigma = 1\) used in that exercise? Use your answer to predict the sign of \(\partial \mathcal{L}/\partial u\) at \(\sigma = 1, \ell = 1\), then confirm it against the formula in part (b).
Computing
Exercise 8.6 (GP Prior Sampling) Consider a Gaussian Process with mean function \(m(x) = 0\) and squared exponential kernel: \[ k(x, x') = \sigma^2 \exp\left(-\frac{(x-x')^2}{2\ell^2}\right) \]
- For \(\sigma^2 = 1\) and \(\ell = 1\), calculate the covariance matrix \(K\) for points \(x = (0, 1, 2)^T\).
- Sample 5 functions from this GP prior at points \(x = (0, 0.5, 1, 1.5, 2, 2.5, 3)\).
- How does changing \(\ell\) to 0.5 affect the sampled functions?
- How does changing \(\sigma^2\) to 4 affect the sampled functions?
Exercise 8.7 (GP Posterior Prediction) You observe the following data points: | x | y | |—|—| | 0 | 1.0 | | 2 | 0.5 | | 4 | 2.0 |
Assume a GP prior with \(m(x) = 0\), squared exponential kernel with \(\sigma^2 = 1\), \(\ell = 1\), and observation noise \(\sigma_n^2 = 0.1\). Because real observations include noise, this exercise uses the noisy-GP posterior obtained from Equation 8.1 and Equation 8.2 by replacing \(\mathbf{K}\) with \(\mathbf{K} + \sigma_n^2\mathbf{I}\) wherever it appears, equivalent to observing \(y_i = f(x_i) + \epsilon_i\) with \(\epsilon_i \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0,\sigma_n^2)\).
- Write the posterior mean formula \(\mu_*\) for predictions at new points.
- Calculate the posterior mean and variance at \(x_* = 1\).
- Calculate the posterior mean and variance at \(x_* = 3\).
- Plot the posterior mean and 95% credible interval for \(x \in [0, 5]\).
Exercise 8.8 (Log Marginal Likelihood) Using the data from the GP Posterior Prediction exercise above (\(X = (0, 2, 4)\), \(y = (1.0, 0.5, 2.0)\)), the log marginal likelihood for a (noiseless) GP is: \[ \log p(y \mid X, \theta) = -\frac{1}{2}y^T K^{-1} y - \frac{1}{2}\log|K| - \frac{n}{2}\log(2\pi) \] where \(\theta = (\sigma^2, \ell)\) are the kernel hyperparameters. This uses the noiseless log marginal likelihood matching Equation 8.1 and Equation 8.2 as derived in the chapter, which generalizes to the noisy-GP posterior of the previous exercise by replacing \(\mathbf{K}\) with \(\mathbf{K}+\sigma_n^2\mathbf{I}\) throughout and extending \(\theta\) to \((\sigma^2,\ell,\sigma_n^2)\).
- Calculate the log marginal likelihood for \((\sigma^2 = 1, \ell = 1)\).
- Calculate it for \((\sigma^2 = 1, \ell = 0.5)\).
- Which hyperparameters are better supported by the data?
- Interpret each term in the log marginal likelihood formula.
Exercise 8.9 (RBF versus Matern: Sample-Path Roughness) The squared exponential kernel is infinitely differentiable, which forces every function sampled from a GP with this kernel to be extremely smooth. A commonly used alternative is the Matern family of kernels, which controls smoothness directly through a parameter. The Matern-5/2 kernel, with length scale \(\ell\) and signal variance \(\sigma^2\), is \[ k_{5/2}(x,x') = \sigma^2\left(1+\frac{\sqrt5\,r}{\ell}+\frac{5r^2}{3\ell^2}\right)\exp\left(-\frac{\sqrt5\,r}{\ell}\right), \qquad r = |x-x'|, \] and produces sample paths that are twice, but not infinitely, differentiable: rougher than the squared exponential kernel’s paths at small scales.
- On a grid of \(n=200\) points over \([0,10]\), build the covariance matrices for the squared exponential kernel and the Matern-5/2 kernel above, both with \(\sigma^2=1\), \(\ell=1\).
- Using the same underlying standard normal draws for both kernels, so that any difference between the two sets of sample paths is due to the kernel alone and not to randomness, draw 5 sample paths from each GP prior and plot them.
- Compute a roughness statistic, such as the sum of squared successive differences \(\sum_i(f(x_{i+1})-f(x_i))^2\) or the number of local extrema, for each path, and report the average over the 5 paths for each kernel.
- Which kernel produces rougher paths, and does the direction of the difference match what the differentiability of the two kernels predicts?
Exercise 8.10 (GP Regression on the Motorcycle Data) The chapter’s motorcycle accident example (Example 8.3) fits a GP to the full MASS::mcycle dataset using the laGP package. Here you will implement the GP posterior directly with matrix algebra, as in Equation 8.1 and Equation 8.2, and use it to fill in a deliberately held-out gap in the data.
- Load
MASS::mcycleand remove all observations withtimesin \((25, 30]\); call the remaining data the training set and the removed points the held-out test set. - Using a squared exponential kernel, estimate the signal variance \(\sigma^2\), length scale \(\ell\), and noise variance \(\sigma_n^2\) by maximizing the log marginal likelihood on the training set only, following the same MLE approach as Example 8.2 but now including \(\sigma_n^2\) as in the noisy-GP posterior of Exercise 8.4 (replace \(\mathbf{K}\) with \(\mathbf{K}+\sigma_n^2\mathbf{I}\) in Equation 8.1 and Equation 8.2). It is numerically convenient to standardize
accel(subtract the training mean, divide by the training standard deviation) before fitting, and to transform predictions back to the original scale afterward. - Predict the posterior mean and standard deviation at the held-out times, and plot them (mean \(\pm\) 2 standard deviations) against the true held-out observations.
- Compare the root-mean-squared prediction error of the GP on the held-out points to that of a cubic polynomial regression fit on the same training data.