  Derivatives are a key nonparametric functional in wide-ranging applications
where the rate of change of an unknown function is of interest. In the Bayesian
paradigm, Gaussian processes (GPs) are routinely used as a flexible prior for
unknown functions, and are arguably one of the most popular tools in many
areas. However, little is known about the optimal modelling strategy and
theoretical properties when using GPs for derivatives. In this article, we
study a plug-in strategy by differentiating the posterior distribution with GP
priors for derivatives of any order. This practically appealing plug-in GP
method has been previously perceived as suboptimal and degraded, but this is
not necessarily the case. We provide posterior contraction rates for plug-in
GPs and establish that they remarkably adapt to derivative orders. We show that
the posterior measure of the regression function and its derivatives, with the
same choice of hyperparameter that does not depend on the order of derivatives,
converges at the minimax optimal rate up to a logarithmic factor for functions
in certain classes. We analyze a data-driven hyperparameter tuning method based
on empirical Bayes, and show that it satisfies the optimal rate condition while
maintaining computational efficiency. This article to the best of our knowledge
provides the first positive result for plug-in GPs in the context of inferring
derivative functionals, and leads to a practically simple nonparametric
Bayesian method with optimal and adaptive hyperparameter tuning for
simultaneously estimating the regression function and its derivatives.
Simulations show competitive finite sample performance of the plug-in GP
method. A climate change application for analyzing the global sea-level rise is
discussed.
