I had a quick read over breakfast of a short paper by Busti, Clarkson, & Seikel (2014) in which they examine the discrepancy between local measurements of the Hubble constant and the corresponding estimate from *Planck*. To this end the authors use a Gaussian process (GP) regression to fit 19 estimates of H(z) at a range of different z, and from their best-fit GP posterior extrapolate a parsimonious estimate for H(0) (which they then compare against *Planck* and various local measurements). The motivation for the authors’ strategy of GP regression—rather than, say, fitting one of the popular parametric cosmological models—is supposedly to offer new insight from a *Bayesian**,** *** non-parametric** computation. However, to my mind the work presented is neither convincingly Bayesian nor convincingly non-parameteric.

Perhaps my biggest grumble is that the authors present as their best prediction for H(0) only the posterior predictive under an empirical Bayes solution: i.e., with maximum likelihood estimation of the two hyper-parameters of their favoured covariance function (the squared exponential)—as opposed to the fully Bayesian solution of marginalization over the posteriors for these two hyper-parameters (cf. e.g. Neal 1998). No doubt this step would give much less impressive bounds on H(0), but should really be done (in my opinion) since the uncertainty in the hyper-parameters is a key uncertainty component of the regression model and for this small dataset would be extremely fast to compute. Likewise, I was disappointed that although the results for a number of alternative covariance functions (in the Matern class) were presented—one of which gives the opposite conclusion to their main analysis—there was no attempt to weigh the plausibility of each via the fully marginalised likelihood (i.e., Bayesian model selection).

Another problem of the analysis I feel is the lack of motivation for the GP prior with covariance function depending on the separation of points in redshift. Of course, it can become a circular argument to attempt to convert the observed z’s to lookback distances (as these depend on H(z) itself) but I would have liked to see some exploration of alternative distance metrics (as transforms of z) in the covariance function. A similar idea would be to use a GP (potentially with look-back distance-based covariance function) as error model for extra-variance in the observed data against the LambaCDM fit, as a means to give further insight into the likely value of H(0) under this benchmark cosmological model.

Finally, although it should not be a significant problem for the present study (as the extrapolation is over a small redshift interval) one needs to be extremely cautious when using (stationary) Gaussian processes for extrapolation beyond the range of available data as the GP posterior has the property of ultimately (at long distance) reverting to its prior mean.