  Gaussian processes (GPs) are frequently used in machine learning and
statistics to construct powerful models. However, when employing GPs in
practice, important considerations must be made, regarding the high
computational burden, approximation of the posterior, choice of the covariance
function and inference of its hyperparmeters. To address these issues, Hensman
et al. (2015) combine variationally sparse GPs with Markov chain Monte Carlo
(MCMC) to derive a scalable, flexible and general framework for GP models.
Nevertheless, the resulting approach requires intractable likelihood
evaluations for many observation models. To bypass this problem, we propose a
pseudo-marginal (PM) scheme that offers asymptotically exact inference as well
as computational gains through doubly stochastic estimators for the intractable
likelihood and large datasets. In complex models, the advantages of the PM
scheme are particularly evident, and we demonstrate this on a two-level GP
regression model with a nonparametric covariance function to capture
non-stationarity.
