How to use GP: effects of the mean function and hyperparameter selection on Gaussian process regression

Gaussian processes have been widely used in cosmology to reconstruct cosmological quantities in a model-independent way. However, the validity of the adopted mean function and hyperparameters, and the dependence of the results on the choice have not been well explored. In this paper, we study the effects of the underlying mean function and the hyperparameter selection on the reconstruction of the distance moduli from type Ia supernovae. We show that the choice of an arbitrary mean function affects the reconstruction: a zero mean function leads to unphysical distance moduli and the best-fit Lambda CDM to biased reconstructions. We propose to marginalize over a family of mean functions and over the hyperparameters to effectively remove their impact on the reconstructions. We further explore the validity and consistency of the results considering different kernel functions and show that our method is unbiased.