Shaw [@Shaw:1988]) the randomized setting also establishes a correspondence between Numerical Analysis and Bayesian Inference providing a natural framework for the statistical description of numerical errors (in which confidence intervals can be derived from posterior distributions). Furthermore [@PalastiRenyi1956; @Diaconis:1988], classical min max numerical quadrature rules can be formulated as solutions of Bayesian inference problems with carefully chosen priors [@Diaconis:1988] and, as shown by O’Hagan [@Hagan:1991; @Hagan:1992], this correspondence can be exploited to discover new and useful numerical quadratures rules. As envisioned by Skilling [@Skilling1992], by placing a (carefully chosen) probability distribution on the solution space of an ODE and conditioning on quadrature points, one obtains a posterior distribution on the solution whose mean may coincide with classical numerical integrators such as Runge-Kutta methods [@schober2014nips]. As shown in [@ChkrebtiiCampbell2015] the statistical approach is particularly well suited for chaotic dynamical systems for which deterministic worst case error bounds may provide little information. While in [@Skilling1992; @schober2014nips; @ChkrebtiiCampbell2015] the probability distribution is directly placed on the solutions space, for PDEs [@Owhadi:2014] argues that the prior distribution must be placed on source terms (or on the image space of an integro-differential operator) and propagated/filtered through the inverse operator to reflect the structure of the 