Uncertainty in Predictions via Bayesian Inference¶
Preface¶
From: Bayesian Linear Regression, we've explained how to obtain credible intervals for \(p(\mathbf{w} | \mathcal{D})\).
Our goal now is to extend this idea further and obtain credible intervals for our prediction function, \(f(\mathbf{x})\).
Definition¶
Given \(p(\mathbf{w}| \mathcal{D})\), we would like to reason about \(p(f(\mathbf{x})|\mathcal{D})\) where \(f(\mathbf{x}) = \mathbf{xw}\) for any \(\mathbf{x}\). Notice this isn't about the stochasticity in \(y\): it isn't about computing \(p(y | \mathbf{x})\).
This is about the uncertainty in our prediction \(f(\mathbf{x}) \approx \mathbb{E}[Y | \mathbf{x}]\).
What is this uncertainty that we're talking about?¶
In Bayesian Linear Regression, we want to find the best linear function in our set of linear functions \(\mathcal{F}\). As we get more samples \(n\), our predictor function will reach the best linear function in \(\mathcal{F}\), in terms of minimizing the squared errors across the entire space.
For any given point \(\mathbf{x}\), we have the prediction \(f(\mathbf{x})\), and then an actual outcome \(y\) is revealed, for which we reason about the three sources of error.
Let \(f_{\text{best}} \in \mathcal{F}\) be the best function in our linear function class, \(f^{*} = \mathbb{E}[Y | \mathbf{x}]\) the true function, which may not be in \(\mathcal{F}\). Then
Our credible interval is only reasoning about the uncertainty in the first term: \(f(\mathbf{x}) - f_{\text{best}}(\mathbf{x})\). We can map this term to the variance component of the reducible error of the generalization error.
The irreducible error is due to the variance in \(Y\) given \(x\): due to \(\sigma_{y}^2\), which is the inherent stochasticity in \(y\) that we simply can't explain given our current information (there may be latent variables that affect \(y\) that we don't account for, or simply just stochasticity). We usually map this stochasticity as: \(f(\mathbf{x}) + \epsilon\) where \(\epsilon \sim \mathcal{N}(0, \sigma_{y}^2)\)
This results in us having:
To summarize, we are talking about the uncertainty in \(f(\mathbf{x})\) due to insufficient data.
Calculating and Reasoning about this uncertainty¶
We use the result that a linear weighting of Gaussian variables results in a Gaussian variable: Let \(\mathbf{v} \in \mathbb{R}^{d}\) be a multivariate Gaussian vector with mean \(\mathbf{\mu} \in \mathbb{R}^{d}\) and covariance \(\boldsymbol{\Sigma} \in \mathbb{R}^{d}\). Let \(\mathbf{b} \in \mathbb{R}^d\) be any vector of coefficients. Then \(z = \mathbf{b}^{\top}\mathbf{v}\) is Gaussian RV, with mean \(\mu = \mathbf{b}^{\top}\boldsymbol{\mu}\) and variance \(\sigma^2 = \mathbf{b}^{\top}\boldsymbol{\Sigma}\mathbf{b}\).
From the result above, if we know the variance \(\sigma_{y}^2\) is known for the likelihood \(p(y | \mathbf{x})\), then \(p(\mathbf{w} | \mathcal{D})\) is a multivariate Gaussian distribution. This is due to our initial assumptions on the model, where \(f(\mathbf{x}) = \mathbf{w^{\top}x} + \epsilon \sim \mathcal{N}(0, \sigma_{y}^2)\). Thus the variance here is fixed.
The random variable \(f(\mathbf{x}) = \mathbf{x^{\top}w}\) is a linear weighting of the multivariate Gaussian RV \(\mathbf{w}\). If \(p(\mathbf{w}) | \mathcal{D}\) is a \(\mathcal{N}(\boldsymbol{\mu}_{n}, \boldsymbol{\Sigma}_{n})\), then we know that:
The 95% credible interval for our prediction \(f(\mathbf{x})\) is \([\mu_{x} - 1.96\sigma_{x}, \mu_{x} + 1.96 \sigma_{x}]\)
If we didn't know \(\sigma_{y}^2\), we must treat it as a random variable and marginalize it out, and through some calculations, the integral we get leads to a Student-t distribution.
Most of the time, we dont know \(\sigma_{y}^2\) though, and a Student-t distribution does not have the same nice properties as a Gaussian. Instead, when we want to use Bayesian linear regression, we assume \(\sigma_{y}\) to be the mode given by our Normal Inverse Gamma distribution: \(\sigma_{y}^2 = \frac{b_{n}}{a_{n}- 1}\).
Given a specific \(\sigma_{y}\), the distribution over \(\mathbf{w}\) is Gaussian, and so we can use the formula above once more. This is simply a property of the NIG, which is explained in Bayesian Linear Regression.
Procedure¶
- Pick the hyperparameters \(\lambda > 0\) for the prior over weights, and \(a_{0}\) and \(b_{0}\) for the prior over \(\sigma_{y}^{2}\). (ie. \(\lambda = 0.01, \ a_{0}=b_{0}=1\)).
- Estimate the NIG \(p(\mathbf{w}, \sigma_{y}^{2} \ | \ \mathcal{D})\) using the formula in Bayesian Linear Regression.
- Get Gaussian posterior \(p(\mathbf{w} | \sigma_{y}^2, \mathcal{D})\) by selecting \(\sigma_{y}^{2} = \frac{b_{n}}{a_{n}-1}\) and computing \(\boldsymbol{\Sigma}_{n} = \sigma_{y}^2 \boldsymbol{\Lambda}^{-1}\) for \(\boldsymbol{\Lambda} = \sum_{i=1}^n \mathbf{x}_{i}^{\top}\mathbf{x}_{i} + \lambda \mathbf{I}\) and \(\boldsymbol{\mu}_{n} = \Lambda^{-1}\sum_{i=1}^n \mathbf{x}_{i}^{\top}y_{i}\) to get \(\mathcal{N}(\boldsymbol{\mu}_{n}, \boldsymbol{\Sigma}_{n})\)
- For any input \(\mathbf{x}\), you can obtain \(p(f(\mathbf{x})| \mathcal{D})\) as a Gaussian with \(\mathcal{N}(\boldsymbol{x\mu}_{n}, \mathbf{x}\boldsymbol{\Sigma}_{n}\mathbf{x}^{\top})\)
We can constantly update the NIG posterior with new data, improving our estimates of \(\mathbf{w}\) and \(\sigma_{y}^2\).
Bayesian Linear Regression therefore elegantly facilitates incorporating new data. Notice that this is one of the reasons we maintain \(\boldsymbol{\Lambda}_{n}\) instead of \(\boldsymbol{\Sigma}_{n}\), because we cannot easily update the inverted matrix with new data.
As we get more data, the variance of the NIG shrinks, as does the variance of \(p(f(\mathbf{x}) | \mathcal{D})\). Notice that \(\boldsymbol{\Sigma}_{n}\) shrinks as \(n\) gets larger because \(\boldsymbol{\Lambda}_{n} = \sum_{i=1}^n \mathbf{x}_{i}^{\top}\mathbf{x}_{i} + \lambda \mathbf{I}\) grows with \(n\). As a consequence, the variance \(\sigma_{x}^2 = \mathbf{x}\boldsymbol{\Sigma}_{n}\mathbf{x}^{\top}\) shrinks as \(n\) grows larger.
This indicates a reduction in uncertainty in our predictions as we get more data, in terms of identifying the best function in our linear function class \(\mathcal{F}\). If \(f^{*}\) is in \(\mathcal{F}\), then the mean converges to \(f^*(\mathbf{x})\) and the interval shrinks to zero around this \(f^*(\mathbf{x})\).
References¶
Generalization Error Student-t Distribution Bayesian Linear Regression Normal-Inverse Gamma Distribution