Are predictions from Bayesian Gaussian Process Regression normally distributed?

Question

This is not directly related to my other question, though the topic is the same. It's also most probably a very trivial question, but bear with me :) I was discussing with a coworker on the use of Gaussian Process Regression, and he made two assertions I don't agree with:

1. GPR can only be used to model a response when the predictors are normally distributed.
2. the response of a GPR model is always normally distributed.

I believe that the first assertion is false (actually, GPR makes no assumptions at all on the joint distribution of the predictors), while the second is only true if the hyperparameters are fixed. However, if we follow a fully Bayesian approach, and we derive the posterior probability distribution of the hyperparameters, then the posterior predictive distribution is no more normally distributed: it's only the distribution of the response, conditional on the hyperparameters and the observations, which is normally distributed. In formulas:

$$y=f(\mathbf{x})+\epsilon, \quad \epsilon\sim N(0,\sigma^2_{noise})$$

and assume a GP prior on $f(\mathbf{x})$. Let $\{(\mathbf{x_1},y_1,)\dots,(\mathbf{x_d},y_d,)\}$ be a set of observations, then the posterior probability distribution of the hyperparameters is

$$p(\boldsymbol{\theta}|\mathbf{y})\propto p(\mathbf{y}|\boldsymbol{\theta})p(\boldsymbol{\theta})$$

Now, the distribution of a new response vector $\mathbf{y^*}$, conditional on the hyperparameters and on the observations, i.e., $p(\mathbf{y^*}|\boldsymbol{\theta},\mathbf{y})$, is normally distributed (right?). However, the posterior predictive distribution is

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

In the integral, only the term $p(\mathbf{y^*}|\boldsymbol{\theta},\mathbf{y})$ is a (multivariate) Normal pdf. $p(\mathbf{y}|\boldsymbol{\theta})$ and $p(\boldsymbol{\theta})$ may have whatever distribution we consider appropriate to model the statistical problem at hand. There's no reason to think that the integral w.r.t. $\boldsymbol{\theta}$ of the product of these three distributions is normally distributed, thus we cannot say that the vector $\mathbf{y^*}|\mathbf{y}$ is normally distributed. Is this correct?

Answer 1

1. GPR does not make any statistical assumptions about the predictors. They don't even have to be numbers! All you need is a prior mean function and a covariance function, which can also be defined for non-numeric data (discrete joint unions, strings, sets, etc.).
2. This is sort of true or assumed when people talk about GPR, because its most interesting aspect is that it allows for exact inference: it essentially just boils down to linear algebra. The moment you introduce more flexibility, e.g. non-Gaussian noise, priors over hyperparameters, you lose this important property and have to resort to approximate inference. That said, even then there typically are computational advantages when using GPR-based models.