# Are predictions from Bayesian Gaussian Process Regression normally distributed?

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?