  This paper focusses on the formulation of numerical integration as an
inferential task. To date, research effort has largely focussed on the
development of Bayesian cubature, whose distributional output provides
uncertainty quantification for the integral. However, the point estimators
associated to Bayesian cubature can be inaccurate and acutely sensitive to the
prior when the domain is high-dimensional. To address these drawbacks we
introduce Bayes-Sard cubature, a probabilistic framework that combines the
flexibility of Bayesian cubature with the robustness of classical cubatures
which are well-established. This is achieved by considering a Gaussian process
model for the integrand whose mean is a parametric regression model, with an
improper flat prior on each regression coefficient. The features in the
regression model consist of test functions which are guaranteed to be exactly
integrated, with remaining degrees of freedom afforded to the non-parametric
part. The asymptotic convergence of the Bayes-Sard cubature method is
established and the theoretical results are numerically verified. In
particular, we report two orders of magnitude reduction in error compared to
Bayesian cubature in the context of a high-dimensional financial integral.
