# Gaussian process prediction interval

How can the prediction interval of a Gaussian process be evaluated? I don't know how to estimate this interval though I can find a 95 % confidence interval for the mean line.

• The target of any prediction interval is the value of a particular random variable. Which random variable do you have in mind? – whuber Oct 19 '15 at 16:12
• the variable is y hat for a new x which is not included in the deisgn points – Wis Oct 19 '15 at 18:00
• Do you mean as in stats.stackexchange.com/questions/33433/… (where @gung has provided a detailed answer)? Or perhaps the more general setting as addressed by Rob Hyndman at stats.stackexchange.com/a/9144? – whuber Oct 19 '15 at 21:33
• @raK1 Obviously it's up to you because it's your question but I think DeltalV deserves the bounty for a fantastic answer so if you agree it will need accepting before tomorrow when the bounty runs out. Just a heads up :) cheers! – EHH Dec 13 '16 at 9:14

I will answer your question inside the Bayesian framework. If you specifically need a frequentist solution, you can get one by slightly modifying my answer, but I think it will underestimate the actual uncertainty: you would need a fully frequentist approach, but I don't know how to do that in this specific case.

To briefly recap the Bayesian GPR (Gaussian Process Regression) framework, you assume the model

$$y=f(x|\boldsymbol{\theta})+\epsilon$$

where $f(x|\boldsymbol{\theta})\sim \mathcal{GP}(\mu(x|\boldsymbol{\theta}),k(x,x'|\boldsymbol{\theta}))$, i.e., the latent variables or function values are distributed as a Gaussian Process, conditionally on the hyperparameters $\boldsymbol{\theta}$, and $\epsilon\sim\mathcal{N}(0,\sigma^2)$ is the usual iid Gaussian noise.

Actually, $\sigma^2$ is an hyperparameter too, so it really belongs in $\boldsymbol{\theta}$, but I wanted to underscore that GPR usually assumes a trivial covariance structure for the noise.

The posterior predictive distribution of $y_*$ at a new point $x_*$, conditionally on data $\{(x_1,y_1,)\dots,(x_d,y_d)\}=(\mathbf{x},\mathbf{y})$ and on hyperparameters $\boldsymbol{\theta}$, is $p(y_*|\boldsymbol{\theta},\mathbf{y})$. Now, assume that the mean function of the Gaussian Process is zero: the general case can also be dealt with, but let's try and keep things simple. Then, using the usual GPR machinery, we get

$$p(f_*|\boldsymbol{\theta},\mathbf{y}) = \mathcal{N}(\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{y},k(x_*,x_*)-\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{k}_*)$$

where

$$K=\pmatrix{k(x_{1},x_{1};\boldsymbol{\theta})& & k(x_{1},x_{d};\boldsymbol{\theta}) \\ & \ddots & \\k(x_{1},x_{d};\boldsymbol{\theta})& & k(x_{d},x_{d};\boldsymbol{\theta}) }$$

$$\mathbf{k}_*=\pmatrix{k(x_*,x_{1};\boldsymbol{\theta}) \\ \vdots \\k(x_*,x_{d};\boldsymbol{\theta})}$$

i.e., conditional on observed data and hyperparameters, the distribution of the latent variable at a new point is still Gaussian, with the mean and standard deviation shown above.

However, we are interested in the distribution of a new observation $y_*$, not a new latent variable. This is easy because in our model the noise is additive, independent of all other variables and is normally distributed with zero mean and variance $\sigma^2$, thus we only need to add the noise variance:

$$p(y_*|\boldsymbol{\theta},\mathbf{y}) = \mathcal{N}(\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{y},k(x_*,x_*)-\mathbf{k}_*^T(K+\sigma^2I)^{-1}\mathbf{k}_*+\sigma^2)$$

Note that I'm considering a single new observation $y_*$, so the distribution $p(y_*|\boldsymbol{\theta},\mathbf{y})$ is just a univariate Gaussian, and the variance is actually a variance and not a variance-covariance matrix.

To actually use this expression, you need values for the hyperparameters, which are not known. There are 2 ways out of this:

1. (the most common solution) the hyperparameters are estimated by MLE, or MAP, and the above expression is used. This approach completely neglects the uncertainty in the estimation of the hyperparameters, thus it doesn't seem very safe.
2. in a fully Bayesian approach, you are not really interested in $p(y_*|\boldsymbol{\theta},\mathbf{y})$, but in the predictive distribution of $y_*$ given $\mathbf{y}$, which is obtained by $p(y_*|\boldsymbol{\theta},\mathbf{y})$ after integrating out the hyperparameters:

$$p(y_*|\mathbf{y})=\int{p(y_*,\boldsymbol{\theta}|\mathbf{y})}d\boldsymbol{\theta}=\int{p(y_*|\boldsymbol{\theta},\mathbf{y})p(\boldsymbol{\theta}|\mathbf{y})}d\boldsymbol{\theta}$$

There are two problems here: given a prior distribution for the hyperparameters $p(\boldsymbol{\theta})$, then the posterior distribution $p(\boldsymbol{\theta}|\mathbf{y})$, which appears in the integral, is not known but must be derived using the Bayes' theorem, which for most hyperpriors means having to run a MCMC. Thus we don't have an explicit expression for $p(\boldsymbol{\theta}|\mathbf{y})$, but only samples from the MCMC. And even if we had an expression for $p(\boldsymbol{\theta}|\mathbf{y})$, then the integral giving $p(y_*|\mathbf{y})$ would still be impossible to evaluate in a closed form in most cases. The solution is an hierarchical Bayes simulation: for each sample $\hat{\boldsymbol{\theta}}_i$ obtained from $p(\boldsymbol{\theta}|\mathbf{y})$ with the MCMC, you draw a sample $y^*_i$ from $p(y_*|\hat{\boldsymbol{\theta}}_i,\mathbf{y})$. Use these $m$ samples $y^*_i$ to estimate an HPD interval for $y_*$, and there you are.

From an intuitive point of view, the second solution draws samples from a distribution where the hyperparameters are "not fixed", but allowed to vary randomly according to their posterior distribution $p(\boldsymbol{\theta}|\mathbf{y})$. Thus the prediction interval obtained in the second case takes into account the uncertainty due to our lack of knowledge about the hyperparameters.

• Sorry to nitpick but could you please explain briefly in the answer how the predictive formula is obtained, as you did in the comments to the other answer, this will make the answer self contained and useful for anyone looking for it in the future. Great discussion of hyperparameter issues btw! Thanks :) – EHH Dec 8 '16 at 11:09
• @EHH no problem, I added the bit you were referring to. – DeltaIV Dec 8 '16 at 12:59
• @DeltalV Awesome, this has really helped me clear some stuff up that I was wondering about! Thanks! – EHH Dec 8 '16 at 17:07
• @Mathews24 that's another question and CV policy is one question per post. Search the site to see if such a question has already been asked, otherwise feel free to ask a new question yourself. – DeltaIV Nov 26 '18 at 7:09
• @DeltaIV You state "we are interested in the distribution of a new observation $y_∗$". I suppose it is context-dependent, but are there any general rules for when one is interested in $y_*$ versus $f_*$? What is the physical interpretation for $f_*$? For example, if $y$ correspond to measurements from a diagnostic with some error $\sigma$, would not being able to model $f_*$ actually give us the model for the underlying physical phenomenon? Why would we prefer $y_*$ in this instance? – Mathews24 May 20 at 21:36

If you're referring to Bayesian regression with Gaussian likelihood, the posterior distribution of a Gaussian process is Gaussian: $$p(f(x)\mid X_n,Y_n) = \mathcal{N}\big(\mu_n(x),\sigma_n^2(x)\big)\,,$$ where $X_n$ are the data locations and $Y_n$ are the data values, and $\mu_n$ and $\sigma_n^2$ computed with Bayesian inference : $$\mu_n(x) = \mathbf{k}_n(x)^\top C_n^{-1}Y_n \text{ and } \sigma_n^2(x) = k(x,x) - \mathbf{k}_n(x)^\top C_n^{-1} \mathbf{k}_n(x)\,,$$ where $\mathbf{k}_n(x) = [k(x_t, x)]_{x_t \in X_n}$ is the kernel vector between $x$ and $X_n$, and $C_n = K_n + \eta^2 I$ with $\eta^2$ the standard deviation of the observation noise and $K_n=[k(x_t,x_{t'})]_{x_t,x_{t'} \in X_n}$ the kernel matrix (see the second chapter of Rasmussen and Williams' book).

Therefore a ~95% confidence interval for $x$ is simply $\mu_n(x) \pm 2 \sigma_n(x)$.

• It sounds like you are assuming (a) the OP wants a Bayesian prediction interval and (b) is willing to adopt a Normal conjugate prior for the process. – whuber Oct 19 '15 at 19:03
• Indeed. I include those assumption in my answer. – Emile Oct 19 '15 at 19:15
• How can U find σ2n(x) ? do you know the equation for it . Thank you – Wis Oct 20 '15 at 21:58
• This answer is absolutely wrong. The OP said that they know how to get the confidence intervals, but want to get the prediction intervals which is the 95% probability intervals for future observation $\mathbf{with ~noise}$. I recommend that @raK1 remove this as being the accepted answer. – EHH Dec 5 '16 at 9:54
• @DeltaIV I believe that future measurements variance being the sum of the function plus noise variance is the answer the OP was looking for (I thought it might be this but hadn't found a reference for it and so was hoping for someone who was sure to answer) so you should maybe write this as an answer to the question. About your other point I completely agree that assuming hyperparameters are known underestimates uncertainty, however pragmatically this is what is done when applying GPs in practise so I was answering in that context. – EHH Dec 7 '16 at 11:10