I read this paper (https://arxiv.org/pdf/1206.2944.pdf) discussing about practical issues of Bayesian optimization and they mentioned that integrating out hyperparameters of Gaussian process using fully-Bayesian treatment would be sensible and give more stable results. The way they do that is to use following equation (Eq. 6 in page 4 of the paper):
$$\hat{a}(x;\{x_n,y_n\})=\int a(x;\{x_n,y_n,\theta\})p(\theta|\{x_n,y_n\})d\theta$$
Does anyone know how this integral is implemented? I know that $a(x;\{x_n,y_n,\theta\})$ can be computed using Gaussian process regression, and in theory $p(\theta|\{x_n,y_n\})$ can be expressed as
$$p(\theta|\{x_n,y_n\})=\frac{p(\{x_n,y_n\}|\theta)p(\theta)}{p(\{x_n,y_n\})}$$
and I expect $p(\{x_n,y_n\}|\theta)$ should be tractable (it is just a multi-variate Gaussian), but what is prior for $p(\theta)$ and how to compute $p(\{x_n,y_n\})$ and more importantly how to combine them to implement the complete integral?
Thanks!
