2.2.1 Approximating distributions

At the most fundamental level, sampling is about approximating distributions. Say we want to know the distribution of the height of people in New York. We could go out and measure everyone, but that would involve measuring the height of 8.4 million people! While this would give us our most accurate answer, it’s also a deeply impractical approach.

Instead, we can sample from the population. This gives us a basic example of Monte Carlo sampling, where we use random sampling to provide data from which we can approximate a distribution. For example, given a database of New York residents, we could select – at random – a sub-population of residents, and use this to approximate the height distribution of all residents. With random sampling – and any sampling, for that matter – the accuracy of the approximation is dependent on the size of the sub-population. What we’re looking to achieve is a statistically significant sub-sample, such that we can be confident in our approximation.

To get a better imdivssion of this, we’ll simulate the problem by generating 100,000 data points from a truncated normal distribution, to approximate the kind of height distribution we may see for a population of 100,000 people. Say we draw 10 samples, at random, from our population. Here’s what our distribution would look like (on the right) compared with the true distribution (on the left):

Figure 2.3: Plot of true distribution (left) versus sample distribution (right)

As we can see, this isn’t a great representation of the true distribution: what we see here is closer to a triangular distribution than a truncated normal. If we were to infer something about the population’s height based on this distribution alone, we’d arrive at a number of inaccurate conclusions, such as missing the truncation above 200 cm, and the tail on the left of the distribution.

We can get a better imdivssion by increasing our sample size – let’s try drawing 100 samples:

Figure 2.4: Plot of true distribution (left) versus sample distribution (right).

Things are starting to look better: we’re starting to see some of the tail on the left as well as the truncation toward 200 cm. However, this sample has sampled more from some regions than others, leading to misrepresentation: our mean has been pulled down, and we’re seeing two distinct peaks, rather than the single peak we see in the true distribution. Let’s increase our sample size by a further order of magnitude, scaling up to 1,000 samples:

Figure 2.5: Plot of true distribution (left) versus sample distribution (right)

This is looking much better – with a sample set of only one hundredth the size of our true population, we now see a distribution that closely matches our true distribution. This example demonstrates how, through random sampling, we can approximate the true distribution using a significantly smaller pool of observations. But that pool still has to have enough information to allow us to arrive at a good approximation of the true distribution: too few samples and our subset will be statistically insufficient, leading to poor approximation of the underlying distribution.

But simple random sampling isn’t the most practical method for approximating distributions. To achieve this, we turn to probabilistic inference. Given a model, probabilistic inference provides a way to find the model parameters that best describe our data. To do so, we need to first define the type of model – this is our prior. For our example, we’ll use a truncated Gaussian: the idea here being, using our intuition, it’s reasonable to assume people’s height follows a normal distribution, but that very few people are above, say, 6’5.” So, we’ll specify a truncated Gaussian distribution with an upper limit of 205 cm, or just over 6’5.” As it’s a Gaussian distribution, in other words, 𝒩(μ,σ), our model parameters are 𝜃 = {μ,σ} – with the additional constraint that our distribution has an upper limit of b = 205.

This brings us to a fundamental class of algorithms: Markov Chain Monte Carlo, or MCMC methods. Like simple random sampling, these allow us to build a picture of the true underlying distribution, but they do so sequentially, whereby each sample is dependent on the sample before it. This sequential dependence is known as the Markov property, thus the Markov chain component of the name. This sequential approach accounts for the probabilistic dependence between samples and allows us to better approximate the probability density.

MCMC achieves this through sequential random sampling. Just as with the random sampling we’re familiar with, MCMC randomly samples from our distribution. But, unlike simple random sampling, MCMC considers pairs of samples: some previous sample x_t−1 and some current sample x_t. For each pair of samples, we have some criteria that specifies whether or not we keep the sample (this varies depending on the particular flavor of MCMC). If the new value meets this criteria, say if x_t is ”preferential to” our previous value x_t−1, then the sample is added to the chain and becomes x_t for the next round. If the sample doesn’t meet the criteria, we stick with the current x_t for the next round. We repeat this over a (usually large) number of iterations, and in the end we should arrive at a good approximation of our distribution.

The result is an efficient sampling method that is able to closely approximate the true parameters of our distribution. Let’s see how this applies to our height distribution example. Using MCMC with just 10 samples, we arrive at the following approximation:

Figure 2.6: Plot of true distribution (left) versus approximate distribution via MCMC (right)

Not bad for ten samples – certainly far better than the triangular distribution we arrived at with simple random sampling. Let’s see how we do with 100:

Figure 2.7: Plot of true distribution (left) versus approximate distribution via MCMC (right)

This is looking pretty excellent – in fact, we’re able to obtain a better approximation of our distribution with 100 MCMC samples than we are with 1,000 simple random samples. If we continue to larger numbers of samples, we’ll arrive at closer and closer approximations of our true distribution. But our simple example doesn’t fully capture the power of MCMC: MCMC’s true advantage comes from being able to approximate high-dimensional distributions, and has made it an invaluable technique for approximating intractable high-dimensional integrals in a variety of domains.

In this book, we’re interested in how we can estimate the probability distribution of the parameters of machine learning models – this allows us to estimate the uncertainty associated with our predictions. In the next section, we’ll take a look at how we do this practically by applying sampling to Bayesian linear regression.

2.2.2 Implementing probabilistic inference with Bayesian linear regression

In typical linear regression, we want to predict some output ŷ from some input x using a linear function f(x), such that ŷ = βx + ξ. With Bayesian linear regression, we do this probabilistically, introducing another parameter, σ², such that our regression equation becomes:

That is, ŷ follows a Gaussian distribution.

Here, we see our familiar bias term ξ and intercept β, and introduce a variance parameter σ². To fit our model, we need to define a prior over these parameters – just as we did for our MCMC example in the last section. We’ll define these priors as:

Note that equation 2.15 denotes the half-normal of a Gaussian distribution (the positive half of a zero-mean Gaussian, as standard deviation cannot be negative). We’ll refer to our model parameters as 𝜃 = β,ξ,σ², and we’ll use sampling to find the parameters that maximise the likelihood of these given our data, in other words, the conditional probability of our parameters given our data D: P(𝜃|D).

There are a variety of MCMC sampling approaches we could use to find our model parameters. A common approach is to use the Metropolis-Hastings algorithm. Metropolis-Hastings is particularly useful for sampling from intractable distributions. It does so through the use of a proposal distribution, Q(𝜃′|𝜃), which is proportional to, but not exactly equal to, our true distribution. This means that, for example, if some value x₁ is twice as likely as some other value x₂ in our true distribution, this will be true of our proposal distribution too. Because we’re interested in the probability of observations, we don’t need to know what the exact value would be in our true distribution – we just need to know that, proportionally, our proposal distribution is equivalent to our true distribution.

Here are the key steps of Metropolis-Hastings for our Bayesian linear regression.

First, we initialize with an arbitrary point 𝜃 sampled from our parameter space, according to the priors for each of our parameters. Using a Gaussian distribution centered on our first set of parameters 𝜃, select a new point 𝜃′. Then, for each iteration t ∈ T, do the following:

1. Calculate the acceptance criteria, defined as:

   

2. Generate a random number from a uniform distribution 𝜖 ∈ [0,1]. If 𝜖 <= α, accept the new candidate parameters – adding these to the chain, assigning 𝜃 = 𝜃′. If 𝜖 > α, keep the current 𝜃 and draw a new value.

This acceptance criteria means that, if our new set of parameters have a higher likelihood than our last set of parameters, we’ll see α > 1, in which case α < 𝜖. This means that, when we sample parameters that are more likely given our data, we’ll always accept these parameters. If, on the other hand, α < 1, there’s a chance we’ll reject the parameters, but we may also accept them – allowing us to explore regions of lower likelihood.

These mechanics of Metropolis-Hastings result in samples that can be used to compute high-quality approximations of our posterior distribution. Practically, Metropolis-Hastings (and MCMC methods more generally) requires a burn-in phase – an initial phase of sampling used to escape regions of low density, which are typically encountered given the arbitrary initialization.

Let’s apply this to a simple problem: we’ll generate some data for the function y = x² + 5 + η, where η is a noise parameter distributed according to η ≈𝒩(0,5). Using Metropolis-Hastings to fit our Bayesian linear regressor, we get the following fit using the points sampled from our function (represented by the crosses):

Figure 2.8: Bayesian linear regression on generated data with low variance

We see that our model fits the data in the same way we would expect for standard linear regression. However, unlike standard linear regression, our model produces predictive uncertainty: this is represented by the shaded region. This predictive uncertainty gives an imdivssion of how much our underlying data varies; this makes this model much more useful than a standard linear regression, as now we can get an imdivssion of the sdivad of our data, as well as the general trend. We can see how this varies if we generate new data and fit again, this time increasing the sdivad of the data by modifying our noise distribution to η ≈𝒩(0,20):

Figure 2.9: Bayesian linear regression on generated data with high variance

We see that our predictive uncertainty has increased proportionally to the sdivad of the data. This is an important property in uncertainty-aware methods: when we have small uncertainty, we know our prediction fits the data well, whereas when we have large uncertainty, we know to treat our prediction with caution, as it indicates the model isn’t fitting this region particularly well. We’ll see a better example of this in the next section, which will go on to demonstrate how regions of more or less data contribute to our model uncertainty estimates.

Here, we see that our predictions fit our data pretty well. In addition, we see that σ² varies according to the availability of data in different regions. What we’re seeing here is a great example of a very important concept, well calibrated uncertainty – also termed high-quality uncertainty. This refers to the fact that, in regions where our Predictions are inaccurate, our uncertainty is also high. Our uncertainty estimates are poorly calibrated if we’re very confident in regions with inaccurate predictions, or very uncertain in regions with accurate predictions. As it’s well-calibrated, sampling is often used as a benchmark for uncertainty quantification.

Unfortunately, while sampling is effective for many applications, the need to obtain many samples for each parameter means that it quickly becomes computationally prohibitive for high dimensions of parameters. For example, if we wanted to start sampling parameters for complex, non-linear relationships (such as sampling the weights of a neural network), sampling would no longer be practical. Despite this, it’s still useful in some cases, and later we’ll see how various BDL methods make use of sampling.

In the next section, we’ll explore the Gaussian process – another fundamental method for Bayesian inference, and a method that does not suffer from the same computational overheads as sampling.

2.3.1 Defining our prior beliefs with kernels

GP kernels describe the prior beliefs we have about our data, and so you’ll often see them referred to as GP priors. In the same way that the prior in equation 2.3 tells us something about the probability of the outcome of our two dice rolls, the GP prior tells us something important about the relationship we expect from our data.

While there are advanced methods for inferring a prior from our data, they are beyond the scope of this book. We will instead focus on more traditional uses of GPs, for which we select a prior using our knowledge of the data we’re working with.

In the literature and any implementations you encounter, you’ll see that the GP prior is often referred to as the kernel or covariance function (just as we have here). These three terms are all interchangeable, but for consistency with other work, we will henceforth refer to this as the kernel. Kernels simply provide a means of calculating a distance between two data points, and are exdivssed as k(x,x′), where x and x′ are data points, and k() represents the function of the kernel. While the kernel can take on many forms, there are a small number of fundamental kernels that are used in a large proportion of GP applications.

Perhaps the most commonly encountered kernel is the squared exponential or radial basis function (RBF) kernel. This kernel takes the form:

This introduces us to a couple of common kernel parameters: l and σ². The output variance parameter σ² is simply a scaling factor, used to control the distance of the function from its mean. The length scale parameter l controls the smoothness of the function – in other words, how much your function is expected to vary across particular dimensions. This parameter can either be a scalar that is applied to all input dimensions, or a vector with a different scalar value for each input dimension. The latter is often achieved using Automatic Relevance Determination, or ARD, which identifies the relevant values in the input space.

GPs make predictions via a covariance matrix based on the kernel – essentially comparing a new data point to previously observed data points. However, just as with all ML models, GPs need to be trained, and this is where the length scale comes in. The length scale forms the parameters of our GP, and through the training process it learns the optimal value(s) for the length scale(s). This is typically done using a nonlinear optimizer, such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimizer. Many optimizers can be used, including optimizers you may be familiar with for deep learning, such as stochastic gradient descent and its variants.

Let’s take a look at how different kernels affect GP predictions. We’ll start with a straightforward example – a simple sine wave:

Figure 2.13: Plot of sine wave with four sampled points

We can see the function illustrated here, as well as some points sampled from this function. Now, let’s fit a GP with a periodic kernel to the data. The periodic kernel is defined as:

Here, we see a new parameter: p. This is simply the period of the periodic function. Setting p = 1 and applying a GP with a periodic kernel to the preceding example, we get the following:

Figure 2.14: Plot of posterior predictions from a periodic kernel with p = 1

This looks pretty noisy, but you should be able to see that there is clear periodicity in the functions produced by the posterior. It’s noisy for a couple of reasons: a lack of data, and a poor prior. If we’re limited on data, we can try to fix the problem by improving our prior. In this case, we can use our knowledge of the periodicity of the function to improve our prior by setting p = 6:

Figure 2.15: Plot of posterior predictions from a periodic kernel with p = 6

We see that this fits the data pretty well: we’re still uncertain in regions for which we have little data, but the periodicity of our posterior now looks sensible. This is possible because we’re using an informative prior; that is, a prior that incorporates information that describes the data well. This prior is composed of two key components:

• Our periodic kernel

• Our knowledge about the periodicity of the function

We can see how important this is if we modify our GP to use an RBF kernel:

Figure 2.16: Plot of posterior predictions from an RBF kernel

With an RBF kernel, we see that things are looking pretty chaotic again: because we have limited data and a poor prior, we’re unable to appropriately constrain the space of possible functions to fit our true function. In the ideal case, we’d fix this by using a more appropriate prior, as we saw in Figure 2.15 – but this isn’t always possible. Another solution is to sample more data. Sticking with our RBF kernel, we sample 10 data points from our function and re-train our GP:

Figure 2.17: Plot of posterior predictions from an RBF kernel, trained on 10 observations

This is looking much better – but what if we have more data and an informative prior?

Figure 2.18: Plot of posterior predictions from a periodic kernel with p = 6, trained on 10 observations

The posterior now fits our true function very closely. Because we don’t have infinite data, there are still some areas of uncertainty, but the uncertainty is relatively small.

Now that we’ve seen some of the core principles in action, let’s return to our example from Figures 2.10-2.12. Here’s a quick reminder of our target function, our posterior samples, and the linear interpolation we saw earlier:

Figure 2.19: Plot illustrating the difference between linear interpolation and the true function

Now that we’ve got some idea of how a GP will affect our predictive posterior, it’s easy to see that linear interpolation falls very short of what we achieve with a GP. To illustrate this more clearly, let’s take a look at what the GP prediction would be for this function given the three samples:

Figure 2.20: Plot illustrating the difference between GP predictions and the true function

Here, the dotted lines are our mean (μ) predictions from the GP, and the shaded area is the uncertainty associated with those predictions – the standard deviation (σ) around the mean. Let’s contrast what we see in Figure 2.20 with Figure 2.19. The differences may seem subtle at first, but we can clearly see that this is no longer a straightforward linear interpolation: the predicted values from the GP are being ”pulled” toward our actual function values. As with our earlier sine wave examples, the behavior of the GP predictions are affected by two key factors: the prior (or kernel) and the data.

But there’s another crucial detail illustrated in Figure 2.20: the predictive uncertainties from our GP. We see that, unlike many typical ML models, a GP gives us uncertainties associated with its predictions. This means we can make better decisions about what we do with the model’s predictions – having this information will help us to ensure that our systems are more robust. For example, if the uncertainty is too great, we can fall back to a manual system. We can even keep track of data points with high predictive uncertainty so that we can continuously refine our models.

We can see how this refinement affects our predictions by adding a few more observations – just as we did in the earlier examples:

Figure 2.21: Plot illustrating the difference between GP predictions and the true function, trained on 5 observations

Figure 2.21 illustrates how our uncertainty changes over regions with different numbers of observations. We see here that between x = 3 and x = 4 our uncertainty is quite high. This makes a lot of sense, as we can also see that our GP’s mean predictions deviate significantly from the true function values. Conversely, if we look at the region between x = 0.5 and x = 2, we can see that our GP’s predictions follow the true function fairly closely, and our model is also more confident about these predictions, as we can see from the smaller interval of uncertainty in this region.

What we’re seeing here is a great example of a very important concept: well calibrated uncertainty – also termed high-quality uncertainty. This refers to the fact that, in regions where our predictions are inaccurate, our uncertainty is also high. Our uncertainty estimates are poorly calibrated if we’re very confident in regions with inaccurate predictions, or very uncertain in regions with accurate predictions.

GPs are what we can term a well principled method – this means that they have solid mathematical foundations, and thus come with strong theoretical guarantees. One of these guarantees is that they are well calibrated, and this is what makes GPs so popular: if we use GPs, we know we can rely on their uncertainty estimates.

Unfortunately, however, GPs are not without their shortcomings – we’ll learn more about these in the following section.

2.3.2 Limitations of Gaussian processes

Given the fact that GPs are well-principled and capable of producing high-quality uncertainty estimates, you’d be forgiven for thinking they’re the perfect uncertainty-aware ML model. GPs struggle in a few key situations:

• High-dimensional data

• Large amounts of data

• Highly complex data

The first two points here are largely down to the inability of GPs to scale well. To understand this, we just need to look at the training and inference procedures for GPs. While it’s beyond the scope of this book to cover this in detail, the key point here is in the matrix operations required for GP training.

During training, it is necessary to invert a D × D matrix, where D is the dimensionality of our data. Because of this, GP training quickly becomes computationally prohibitive. This can be somewhat alleviated through the use of Cholesky deomposition, rather than direct matrix inversion. As well as being more computationally efficient, Cholesky decomposition is also more numerically stable. Unfortunately, Cholesky decomposition also has its weaknesses: computationally, its complexity is O(n³). This means that, as the size of our dataset increases, GP training becomes more and more expensive.

But it’s not only training that’s affected: because we need to compute the covariance between a new data point and all observed data points at inference, GPs have a O(n²) computational complexity at inference.

As well as the computational cost, GPs aren’t light in memory: because we need to store our covariance matrix K, GPs have a O(n²) memory complexity. Thus, in the case of large datasets, even if we have the compute resources necessary to train them, it may not be practical to use them in real-world applications due to their memory requirements.

The last point in our list concerns the complexity of data. As you are probably aware – and as we’ll touch on in Chapter 3, Fundamentals of Deep Learning – one of the major advantages of DNNs is their ability to process complex, high-dimensional data through layers of non-linear transformations. While GPs are powerful, they’re also relatively simple models, and they’re not able to learn the kinds of powerful feature representations that are possible with DNNs.

All of these factors mean that, while GPs are an excellent choice for relatively low-dimensional data and reasonably small datasets, they aren’t practical for many of the complex problems we face in ML. And so, we turn to BDL methods: methods that have the flexibility and scalability of deep learning, while also producing model uncertainty estimates.