  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 LCDM 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.
