This paper advocates expansion of the role of Bayesian statistical inference when formally quantifying uncertainty in computer models defined by systems of ordinary or partial differential equations. We adopt the perspective that implicitly defined infinite dimensional functions representing model states are objects to be inferred probabilistically. We develop a general methodology for the probabilistic integration of differential equations via model based updating of a joint prior measure on the space of functions and their temporal and spatial derivatives. This results in a posterior measure over functions reflecting how well they satisfy the system of differential equations and corresponding initial and boundary values. We show how this posterior measure can be naturally incorporated within the Kennedy and O'Hagan framework for uncertainty quantification and provides a fully Bayesian approach to model calibration. By taking this probabilistic viewpoint, the full force of Bayesian inference can be exploited when seeking to coherently quantify and propagate epistemic uncertainty in computer models of complex natural and physical systems. A broad variety of examples are provided to illustrate the potential of this framework for characterising discretization uncertainty, including initial value, delay, and boundary value differential equations, as well as partial differential equations. We also demonstrate our methodology on a large scale system, by modeling discretization uncertainty in the solution of the Navier-Stokes equations of fluid flow, reduced to over 16,000 coupled and stiff ordinary differential equations. Finally, we discuss the wide range of open research themes that follow from the work presented. Keywords: Bayesian Numerical Analysis, Uncertainty Quantification, Gaussian Processes, Differential Equation Models, Uncertainty in Computer Models. In all of the sciences, economics and engineering there is a fundamental reliance on the use of differential equation models to describe complex phenomena concisely, using few but readily interpretable parameters, which we denote by θ. In systems of differential equations, the derivatives with respect to spatial variables, x∈ D ⊂ R d , and temporal variables, t ∈ [a, b] ⊂ R + , are related to the implicitly defined states, u(x, t) ∈ R P , which are hence often analytically intractable. The main challenge of working with differential equation models, from both mathematical and statistical perspectives, is that the solutions are generally not available in closed form. Consequently, numerical solvers are used to approximate a system solution over a discretisation grid. Numerical integration of a system of differential equations yields a deterministic approximation based on subjective choices regarding the order of the numerical method, the specified error tolerance, and the implied discretisation grid, all of which impact the quality of the approximation (Butcher, 2008; Higham, 1996) . These methods implicitly make strong assumptions, for example, that small changes in the step size result in small cumulative changes in the output, and that there exists a unique, well-conditioned solution. Additionally, meaningful inferences from the model are dependent on the assumption that the approximated states contain negligible numerical error. Although the specified error tolerance and associated point-wise errors can be reduced by refining the discretisation grid, there is naturally a trade off between the computational effort involved and the accuracy of the approximated solution, which determines whether the assumption of negligible error is reasonable in * oksana.chkrebtii@gmail.com; address: Department of Statistics, The Ohio State University, 1958 Neil Avenue, 404 Cockins Hall, Columbus, OH 43210-1247 1 arXiv:1306.2365v2 [stat.ME] 28 Apr 2014 practice. Local and global numerical errors are defined point-wise and relate to the asymptotic behaviour of the deterministic approximation of the model. This form of numerical error analysis is not well suited for quantifying the functional uncertainty in the solution for the purpose of fully probabilistic model inference. An additional challenge when working with differential equation models is that in some cases classical solution approximations provide potentially misleading trajectories. Long-time numerical solutions may be globally sensitive to truncation errors introduced at each discretisation point (e.g. Sauer et al., 1997) . Illconditioned models can give rise to another source of uncertainty in the form of solution multiplicity (Ascher et al., 1988; Beyn and Doedel, 1981) . This often cannot be readily verified analytically and clearly poses a problem for classical numerical integration methods that produce a single deterministic solution (Keller, 1968) . In this paper we address the modelling challenges that arise when using the following classes of deterministic ordinary and partial differential equation models, for which analytical solutions are generally not available. Ordinary Differential Equation (ODE) models implicitly represent the derivative of states u(t), with respect to time t through u t (t) := d dt u(t) = f (t, u(t), θ), where we make the standard assumption that f is continuous in the first argument and Lipschitz continuous in the second argument. Inputs and boundary constraints define several model variants. The Initial Value Problem (IVP) models the system states with fixed initial condition u * (a), evolving according to the ODE as follows, u t (t) = f t, u(t), θ , t ∈ [a, b], u(a) = u * (a). The existence of a solution is guaranteed under mild conditions (see for example, Butcher, 2008; Coddington and Levinson, 1955) . Such models may be high dimensional and contain other complexities such as algebraic components, functional inputs, or higher order terms. While IVP models specify a fixed initial condition on the system states, the Mixed Boundary Value Problem (MBVP) may constrain different states at different time points. Typically these constraints are imposed at the ends of the time domain giving the general form for two state mixed boundary value problems, which can be straightforwardly generalised to higher dimensions and extrapolated beyond the final time point b. Whereas a unique IVP solution exists under relatively mild conditions, imposing mixed boundary constraints can result in multiple solutions (e.g. Keller, 1968) introducing severe problems for parameter estimation methods. The Delay Initial Function Problem (DIFP) generalises the initial constraint of IVPs to an initial function, φ(t), thereby relating the derivative of a process to both present and past states at lags τ j ∈ [0, ∞), DIFPs are well suited to describing biological and physical dynamics that take time to propagate through systems. However, they pose challenges to numerical techniques due to potentially large and structured truncation error (Bellen and Zennaro, 2003) . Furthermore, the sensitivity of DIFP solutions to small changes in the initial function, such as those due to interpolation error, may push an otherwise well behaved system into a chaotic regime (Taylor and Campbell, 2007) . Partial Differential Equation (PDE) models represent the derivative of states with respect to multiple arguments, for example time and spatial variables. The main classes of PDE models are based on elliptic, parabolic and hyperbolic equations. Further adding to their complexity, functional boundary constraints and initial conditions make PDE models more challenging, while their underlying theory is less developed compared to ODE models (Polyanin and Zaitsev, 2004  Numerical discretisation error is often characterised in the form of an upper bound, referred to as verification error in the mathematics and engineering literature (Oberkampf and Roy, 2010) . In the many cases where the error is not negligible, for example in large scale Navier-Stokes equations used in modelling fluid flow, an accurate representation of the error propagation resulting from the approximation is required. An illustrative example of quantifying this numerical error in fluid dynamics is provided by Oliver et al. (2014) and a preliminary attempt at a Bayesian formalisation of the uncertainty induced is described by Oliver and Moser (2011) . The approach advocated in this paper is different in nature. We characterise discretisation uncertainty using a probabilistic representation of the solution conditional on fixed initial conditions and solver parameters. Our methodology challenges some of the assumptions used by numerical solvers and yields algorithms robust to many of the inferential issues that occur. Increasing attention is being paid to the challenges associated with quantifying uncertainty in system models, and in particular those based on mathematical descriptions of natural or physical processes using differential equations (Ghanem and Spanos, 2003; Huttunen and Kaipio, 2007; Kaipio et al., 2004; Marzouk and Najm, 2009; Marzouk et al., 2007; Stuart, 2010) . The important question of how uncertainty propagates through complex systems arises in a variety of contexts, and the focus is often on two main problems. The forward problem is concerned with how input uncertainty (i.e. in the parameters and initial conditions) propagates over the numerical trajectory. This type of problem has been investigated in the engineering and applied mathematics literature, making use of sampling based methods, perturbation of initial system states (Mosbach and Turner, 2009), moment closure approaches, and more recently polynomial chaos methods (see Ghanem and Spanos, 2003; Xiu, 2009 , for a complete overview). The inverse problem concerns the uncertainty in the model inputs given measurements with error of the model outputs, and much work in this area has been presented in the statistics literature (Bock, 1983; Brunel, 2008; Calderhead and Girolami, 2011; Calderhead et al., 2009; Campbell and Lele, 2013; Campbell and Steele, 2011; Campbell and Chkrebtii, 2013; Dowd, 2007; Gugushvili and Klaassen, 2012; Ionides et al., 2006; Liang and Wu, 2008; Ramsay et al., 2007; Xue et al., 2010; Xun et al., 2014). Additionally, the inverse problem itself may be ill-conditioned when the data does not lie in the solution space of the differential equation model (Brynjarsdóttir and O'Hagan, 2014; Kennedy and O'Hagan, 2001) . The issue of discretisation uncertainty is intrinsic to both the forward and inverse problems, and it seems natural to consider the unknown solution of the differential equation itself as part of the inferential framework. There are many examples suggesting that a probabilistic functional approach is well suited for modelling solution uncertainty for differential equation models; indeed solving PDEs involves estimation of a field, DIFPs require estimation of an unknown function, and MBVPs may have multiple solutions. All of these problems may be naturally framed in terms of measures over a function space (Stuart, 2010) . A Bayesian functional approach offers a way of defining and updating measures over possible trajectories of a system of differential equations and propagating the resulting uncertainty consistently through the inferential process. The mathematical framework for defining numerical problems on function spaces has been developed through the foundational work of Stuart (2010), Somersalo (2007), and Skilling (1991) , and this is the natural setting that we adopt in our contribution. We develop these ideas further and investigate general approaches to working with classes of differential equations from a Bayesian perspective. The use of probability in numerical analysis has a long history. In the 1950s, Monte Carlo methods were originally developed to estimate analytically intractable high dimensional integrals, and extensions of this methodology have since had an immense impact in a wide variety of areas. Markov chain Monte Carlo (MCMC) methods have allowed the rapid development and application of Bayesian statistical approaches to quantifying uncertainty in scientific and engineering problems. Theoretical analyses of Bayesian statistical approaches to classical numerical analysis problems appear to date back all the way to Poincaré, as eloquently summarised by Diaconis (1988) . More recently, O'Hagan (1992); Skilling (1991) also provided practical motivations for developing and employing Bayesian techniques to characterise the uncertainty that results from computational restrictions in the case of integration, since we cannot evaluate a function at every point in its input space and must therefore account for discretisation uncertainty. Such work is closely related to the idea of computer emulation of large and complex computer models; the use of Gaussian processes to model functions that are computationally very expensive to evaluate has found application in many challenging areas of science, such as climate prediction and geophysics (Conti and O'Hagan, 2010; Kennedy 3 and O'Hagan, 2001; Tokmakian et al., 2012) . Other recent examples of applying Bayesian approaches to numerical integration include Hennig and Hauberg (2014) , who apply Skilling's approach to the calculation of Riemannian statistics, and Kennedy (1998); Osborne et al. (2012) , whose methods allow for a more informed probabilistic approach to the choice of integration points, based on the predicted variance of the surrogate model given limited function evaluations. However, such approaches assume no error on the function evaluation itself, and are therefore not applicable to differential equation models, where the function evaluation gives only an approximate derivative of the solution of the system at the next discrete time point, based on the current predictive distributions over the states. The aim of this paper is to make a case for probabilistically characterising discretisation error when solving general classes of differential equations, and for viewing this epistemic uncertainty as an important and integral part of the overall model inference framework. We first develop a probabilistic differential equation solver that characterises discretisation uncertainty in the solution, by formalising and substantially extending the ideas first presented by Skilling (1991) . The sequential sampling approach suggested is made computationally feasible through recursive Bayesian updating, and it closely follows the sequential construction of a proof of consistency that we present in the Appendix, ensuring the algorithm converges to the exact solution as the time-step tends to zero. In our approach, we treat the sequential model evaluations as auxiliary parameters, over which the solution can then be marginalised, and we further extend the basic method to the cases of mixed boundary value problems exhibiting multiplicity of solutions, delay initial function problems with partially observed initial functions, stiff and chaotic PDEs via spectral projection, as well as providing a framework for directly solving PDEs. We give an example of the scalability of our approach using a Navier-Stokes system involving the solution of over 16,000 coupled stiff ordinary differential equations. Finally, we adopt a flexible and general forward-simulation approach that allows us to embed the formal quantification of discretisation uncertainty for differential equations within the Kennedy and O'Hagan framework, hence defining it as part of the full inferential procedure for inverse problems. We provide code for our Probabilistic Differential Equation Solver (PODES), which allows replication of all results presented in this paper. This is available at http://web.warwick.ac.uk/PODES. We begin in Section 2 by discussing model discrepancy within the Bayesian approach to uncertainty quantification, as first presented by Kennedy and O'Hagan (2001) and subsequently adopted in practice throughout the statistics literature. We discuss how existing numerical approaches to differential equation integration ignore discretisation uncertainty when the solution is not available in closed form. We then develop a general Bayesian framework for modelling this epistemic uncertainty in differential equations, inspired by the work of Skilling (1991) and O'Hagan (1992) . We begin in Section 3 by defining functional priors on the states and derivatives for general ordinary and partial differential equations. Then, in Section 4, we develop a Bayesian approach for updating our prior beliefs given model information obtained sequentially and self-consistently over a given discretisation grid, thus probabilistically characterising the finite dimensional representation of an underlying infinite dimensional solution given model parameters and initial conditions. We show that the resulting posterior trajectory converges in probability to the exact solution under standard assumptions. We then incorporate this methodology into the inverse problem in Section 5, and provide examples in Section 6. Finally, we discuss our proposed approach to Bayesian uncertainty quantification for differential equation models and highlight open areas of research. The Appendix contains proofs and further mathematical and algorithmic details. The importance of formally quantifying sources of uncertainty in computer simulations when mathematically modelling complex natural phenomena, such as the weather, ocean currents, ice sheet flow, and cellular protein transport, is widely acknowledged. These model based simulations inform the reasoning process 4 when, for example, assessing financial risk in deciding on oil field bore configurations, or forming government policy in response to extreme weather events. As such, accounting for all sources of uncertainty and propagating them in a coherent manner throughout the entire inference and decision making process is of great importance. Consider data, y(t), observed at discrete time points t = [t 1 , t 2 , . . . , t T ] and a set of model parameters θ. Using the exact solution of a mathematical model represented by u * (t, θ), a simplified observation model based on some measurement error structure (t) is, In full generality, a nonlinear transformation G of u * (t, θ) may be observed, however for expositional clarity we assume this is an identity. In our setting, u * (t, θ) is the unique function satisfying a general system of differential equations. When performing statistical inference over such models, we are interested in the joint posterior distribution over all unknowns, possibly including the initial and boundary conditions. Throughout the paper, θ will be augmented to contain all the variables of interest. In the landmark paper of Kennedy and O'Hagan (2001) a number of sources of uncertainty were identified when modelling a natural or physical process. Their acknowledgement of uncertainty from incomplete knowledge of the process and the corresponding model inadequacy motivates a probabilistic view of a modelreality mismatch δ(t) also studied in recent papers by Brynjarsdóttir and O'Hagan (2014); Huttunen and Kaipio (2007); Kaipio and Somersalo (2007); and Stuart (2010) . In the engineering literature, this form of error is also know as validation error (see for example, Oberkampf and Roy, 2010). Following Kennedy and O'Hagan by defining δ(t) as a random function drawn from a Gaussian Process (GP), the observational model becomes, Due to the lack of an analytical solution for most nonlinear differential equations, the likelihood p y(t) | u * (t, θ), θ cannot be obtained in closed form. This issue is dealt with throughout the statistics literature by replacing the exact likelihood with a surrogate, p y(t) |û N (t, θ), θ , based on an N -dimensional approximate solution,û N (t, θ), obtained using numerical integration methods (for example, a Runge-Kutta solver with N time steps) whose accuracy has been well studied (see for example, Henrici, 1964). However, there are still many scenarios for which standard approaches for characterising upper bounds on numerical error, and subsequently propagating this uncertainty throughout the rest of the inferential process are unsatisfactory. Limited computation and coarse mesh size are contributors to the numerical error, as described in the comprehensive overview of Oberkampf and Roy (2010) . The seemingly innocuous assumption of negligible numerical integration error can subsequently lead to serious statistical bias and misleading inferences for certain classes of differential equation models, and we illustrate this point in Section 6.1 using a simple example. We may represent this additional uncertainty using the term ζ(t, θ) = u Kennedy and O'Hagan (2001) suggest modelling the computer code using a Gaussian process that is agnostic to the specific form of the underlying mathematical model, which is therefore considered as a "black box". In the discussion of their paper, H. Wynn argues that the sensitivity equations and other underlying mathematical structures that govern the computer simulation could also be included in the overall uncertainty analysis. In the current paper, we "open the black box" by explicitly modelling the solution and associated discretisation uncertainty,û N (t, θ) + ζ(t, θ). This allows the Kennedy and O'Hagan framework to be further enriched by incorporating detailed knowledge of the mathematical model being employed. Noting that the posterior measure over the computer model is defined independently of the observed data suggests a further construction by probabilistically characterising the uncertainty due to system states being defined implicitly by differential equations. This leads to the development of a fully probabilistic scheme for solving differential equations and suggests a powerful new approach to modelling uncertainty in computer models. We model uncertainty in a finite dimensional representation of the infinite dimensional solution through a probability statement on a space of suitably smooth functions. Restricting ourselves to Hilbert spaces for 5 modelling our knowledge of u * (t, θ), we define a Gaussian prior measure on the function space (Stuart, 2010) . We then directly model our knowledge about the solution via the stochastic process u(t, θ), thus replacing (6) with y(t) = u(t, θ) + δ(t) + (t). The contribution in this paper is in our definition and use of u(t, θ), therefore for expositional clarity we focus our attention on the joint posterior measure over differential equation model states, parameters, and associated auxiliary parameters, Ψ, of our probabilistic model of uncertainty, Our Bayesian probabilistic integration framework for differential equation models fully explores the space of trajectories that approximately satisfy model dynamics under theoretical guarantees of consistency. The framework presented in this paper is demonstrated on models described using nonlinear ordinary, delay, mixed boundary value and partial differential equations, and we explicitly address multiplicity of solutions in Section 6.2. For notational simplicity we hereafter omit the dependence of u on θ. When considering uncertainty in the infinite dimensional solution of differential equations, the natural measure for spaces of a wide class of functions is Gaussian (Stuart, 2010). As such we define a Gaussian process (GP) prior measure jointly on the state u and its time derivative u t . We will now discuss the prior construction for ordinary differential equations, and then consider the case for partial differential equations. In the remainder of the paper, we will denote by R λ a deterministic, square integrable kernel function with length-scale λ ∈ (0, ∞), and its integrated version by Q λ (t 1 , t 2 ) = t1 a R λ (s, t 2 )ds. Here a represents the lower boundary of our temporal domain. Our model for the derivative has covariance operator cov( , where α is a prior precision parameter. Therefore the covariance on the state is its integrated version, cov(u(t 1 ), u(t 2 )) = α The cross covariance terms are defined in a similar manner and denoted as RQ(t 1 , t 2 ) and QR(t 1 , t 2 ) respectively. We note that RQ(t 1 , t 2 ) = QR † (t 1 , t 2 ), where † represents the adjoint. We assume a joint Gaussian prior measure on the state and its derivative, The Gaussian conditional prior measures for e.g. p (u t (t) | u(t)) take the standard forms (Stuart, 2010) . The choice of prior means and covariance structure, as well as the impact and choice of auxiliary parameters, Ψ = [α, λ], are discussed in the Appendix (Section 8.2). For exposition, we take the auxiliary parameters Ψ as fixed, although we note that propagating a distribution over Ψ through the probabilistic solver is straightforward, and estimation of Ψ is addressed in Section 6.4. This model straightforwardly generalises to ODE problems of order greater than one, by defining the prior jointly on the state, u, and any required derivatives. Alternatively, higher order ODE problems may be restated in the first order by introducing additional states. We have found that adding states to the model is computationally straightforward, so we adopt this approach in practice, working with the prior in equation (8). PDE problems require a prior measure over multivariate trajectories, as well as modelling single or higher order derivatives with respect to spatial inputs. Therefore the prior specification will depend on the PDE model. Consider as an illustrative example the parabolic heat equation, modelling the heat diffusion over time along a single spatial dimension by, 6 One modelling choice for incorporating the spatial component of the PDE into the covariance is to adopt a product structure. The covariance over time is defined as above, and over space may be defined using a similar construction. Let R µ be a kernel function with length-scale µ. Q µ (z, x 2 )dz, where c denotes a spatial lower boundary (in (9), this is c = 0). The spatial covariance structures are defined similarly to the temporal form in the ODE example as, where β is a spatial prior precision parameter. Cross covariances are defined analogously. The prior construction follows by defining a product structure for space and time, In the next section, we describe in detail the framework for sequentially linking the prior in equation (8) with the ODE model. We describe the corresponding algorithm for updating the prior in equation (10) under the heat equation PDE in the Appendix (Section 8.1). We provide examples illustrating its use experimentally in Section 6. The joint prior on the state and its derivative is defined in equation (8) in terms of mean functions m, m t and a covariance structure determined by the choice of kernel function R λ , which is parameterised by the auxiliary variables α and λ. The choice of the covariance kernel should reflect our assumptions regarding the smoothness of the exact but unknown differential equation solution. In the Appendix (Section 8.3), we provide covariance structures based on two kernels: the infinitely differentiable squared exponential and the non-differentiable uniform kernel. Gaussian process models typically match the kernel for a given application to prior information about smoothness of the underlying function (see, for example, Rasmussen and Williams, 2006) . It is also important to be aware that imposing unrealistically strict smoothness assumptions on the state space by choice of covariance structure may introduce estimation bias if the exact solution is not at least as smooth. Therefore, in cases where the solution smoothness is not known a priori, one can err on the side of caution by using less regular kernels. In the examples reported in this paper we chose to work with stationary derivative covariance structures for simplicity, however we point out that there will be classes of problems where non-stationary kernels may be more appropriate and can likewise be incorporated into the Gaussian process framework. The mean function is chosen based on any prior knowledge about how the system evolves over time. Typically, however, such prior detailed knowledge is unavailable, in which case it is reasonable to use a constant mean function. Additionally, we condition on the known boundary values, enforcing them by choice of the prior mean. For example, for IVPs, we satisfy the boundary constraint exactly by choosing a differentiable prior mean function m such that m(a) := u * (a). Now that the form of the prior measure has been defined, the corresponding posterior is obtained in the following two sections in a sequential manner. This iterative updating closely follows the sequential structure of the proof of consistency, provided in the Appendix (Section 8.6). Here we point out that, although the amount and quality of prior information regarding the true solution will affect the efficiency of our probabilistic solvers, we can still expect convergence as the time step tends to zero, subject to some standard assumptions on the kernel function. Furthermore, in Section 6.1, we illustrate the extent of prior influence on the probabilistic solution for the heat equation PDE. We find that the use of a zero-mean prior on the spatial and temporal derivatives has minimal impact on the solution, even for reasonably rough discretisation grids. We now introduce a sequential framework to characterise the epistemic model uncertainty component in equation (7). We restrict our attention to the ODE initial value problem (1) and consider more general ODE and PDE problems in Section 6. We model discretisation uncertainty by updating the joint prior on the unknown state and derivative iteratively, by evaluating the deterministic ODE model at a finite number N of grid points over the temporal domain, which we denote by f 1:N . Algorithm 1 produces a sample from the joint distribution, [u, Let us consider the problem of modelling uncertainty in the exact but unknown solution of a nonlinear ODE initial value problem. We are given the differential equation model, u t (t) = f (t, u(t), θ), which implicitly defines the derivatives in terms of a Lipschitz continuous function f of the states given θ and a fixed initial value, u * (a). We model our uncertainty regarding the unique explicit state, u, and its derivative, u t , satisfying the initial value problem, via the joint GP prior measure in equation (8). The following procedure sequentially links the prior on the state with the ODE model defined by f . Consider a discretisation grid made up of N time points We begin by fixing the known initial value, u(s 1 ) := u * (a), and computing the exact derivative f 1 := f (s 1 , u(s 1 ), θ) at s 1 , via the deterministic ODE model. We then update our joint GP prior given the computed exact derivative f 1 , obtaining the conditional predictive distribution for the state at the subsequent grid location s 2 , with, This predictive distribution describes our current uncertainty about the solution at time s 2 . We now sample a realisation, u(s 2 ), of the predictive process, and again link our prior to the deterministic ODE model by computing f 2 := f (s 2 , u(s 2 ), θ). In contrast to the first time point s 1 , we can no longer guarantee that the realisation at the second time point, u(s 2 ), and its derivative, u t (s 2 ), exactly satisfy the ODE model, i.e. that u t (s 2 ) = f 2 . Therefore, at time s 2 we explicitly model the mismatch between the ODE evaluation f 2 and the process derivative, u t (s 2 ), as, where the magnitude of the mismatch may be described by the variance, C t (s 2 , s 2 ), of the predictive posterior over the derivative given by, For systems in which we believe this mismatch to be strongly non-Gaussian, we may appropriately modify the above model and associated sampling strategy, as described in the Discussion section. We may now update our joint posterior from the previous iteration by conditioning on the augmented vector f 1:2 := [f 1 , f (s 2 , u(s 2 ), θ)]. We therefore define the matrix Λ 2×2 := diag{0, C t (s 2 , s 2 )} to describe the mismatch between the process derivative and the ODE function evaluations at s 1 and s 2 . The new predictive posterior, As in the previous two steps, we sample the state realisation u(s 3 ) from the above predictive posterior distribution. We then apply the deterministic transformation f , to obtain f 3 := f (s 3 , u(s 3 ), θ), whose mismatch with the realised derivative u t (s 3 ) is again modelled as, We next augment f 1:3 := [f 1:2 , f 3 ]. The corresponding mismatch matrix therefore has a diagonal structure, Λ 3×3 := diag{Λ 2×2 , C t (s 3 , s 3 )}, where the step-ahead predicted covariance in the derivative space is, The diagonal elements of Λ, which are step-ahead predictive derivative variances, are non-decreasing, reflecting our growing uncertainty about how well the realised state obeys the ODE model as we take our sample further and further away from the known initial state. It is also shown in the Appendix (Section 8.7) that its elements tend to zero with the step size but at a much faster rate than the step size. The general scheme can be written according to Algorithm 1, and an illustration of this is provided in Figure 1 . Algorithm 1 Sample from the joint posterior distribution of u and f 1:N for an ODE initial value problem given θ, Ψ, N At time s 1 := a, initialise the derivative f 1 := f s 1 , u(s 1 ), θ for initial state u(s 1 ) := u * (a), and define associated model-derivative mismatch, Λ 1×1 := 0; for n = 1 : N − 1 do Define the predictive state mean and variance, Sample step-ahead realisation u(s n+1 ) from the predictive distribution of the state, Evaluate the ODE model f n+1 := f (s n+1 , u(s n+1 ), θ) for realisation u(s n+1 ) at the subsequent grid point, s n+1 , and augment the vector f 1:n+1 := [f 1:n , f n+1 ]; Define the predictive derivative variance, and augment the matrix We emphasise that integration of the ODE model proceeds probabilistically via Gaussian process integration, without the use of numerical approximations. We obtain a posterior distribution over trajectories 9 Figure 1 : Illustration of Algorithm 1 for generating a sample from the joint distribution of derivative observations and possible trajectories with density p u(t), f 1:N | θ, Ψ . Given two derivative model realisations (red points), we obtain a posterior distribution over the derivative space (top left) and over the state space (bottom middle). A sample is then drawn from the predictive posterior over the states at the next time point s n (bottom middle), and a model realisation is obtained by mapping u(s n ) to the derivative space via the function f (top middle, rightmost red point). Given these three model evaluations, this procedure may be repeated (bottom right, top right) in an analogous manner. 10 u that are governed by the differential equation model over a discrete number of grid points. In the Appendix (Section 8.6), we provide a proof that under certain conditions, the posterior process u(t) obtained via Algorithm 1 tends to the unique solution u * (t) satisfying IVP (1), as follows, where h is the maximum length between consecutive discretisation grid points. We now present a framework to propagate the model discretisation uncertainty characterised in Section 4 through the inverse problem, where we wish to infer model parameters from measurement data. We describe one possible Markov chain Monte Carlo procedure that generates a sample from the joint posterior distribution of the state at the data locations t = [t 1 , · · · , t T ] and the unknown parameters θ conditional on noisy observations y(t) = G u(t), θ + (t) of the states. We assume the parameters θ also include any unknown initial conditions, or auxiliary parameters. Once again, for expositional simplicity, we focus on direct observations of states governed by an ODE initial value problem. However, extension to other ODE and PDE problems is straightforward given a probabilistic solution. We provide an application where states are indirectly observed through the nonlinear transformation of states G in Section 6.4. Algorithm 2 targets the posterior density in equation (7) by forward model proposals via conditional simulation from Algorithm 1. This proposal step avoids the need to explicitly calculate the intractable marginal density p u(t), f 1:N | θ, Ψ, N df 1:N , and can be implemented efficiently as described in the Appendix (Section 8.1). Such partially likelihood-free MCMC implementations (Marjoram et al., 2003) are widely used in the area of inference for stochastic differential equations (see, for example, Golightly and Wilkinson, 2011) for simulating sample paths within the inverse problem. Algorithm 2 Draw K samples from the posterior distribution with density p θ, u(t) | y(t), Ψ Initialise θ and conditionally sample a realisation of the state u(t) via Algorithm 1; for k = 1 : K do Propose θ ∼ q(θ | θ), where q is a proposal density; Sample a probabilistic realisation of the state u (t) conditioned on θ via Algorithm 1; Compute: We now have all the components to take a fully Bayesian approach for quantifying uncertainty on differential equation models of natural and physical systems. In this section we use the framework developed in the previous two sections in applications to a wide range of systems, including ODE and PDE boundary value problems and delay initial function problems. 11 As an illustrative example, we demonstrate the use of our probabilistic framework on the heat equation presented in (9). We model our uncertainty about the solution through the prior (10) defined over time and space using a product covariance structure as described in Section 3. We use our probabilistic framework to integrate the state over time at each of the spatial discretization grid points, employing a uniform covariance kernel for the time component and a squared exponential covariance kernel for the spatial component. We therefore characterise the discretisation uncertainty in both the time and spatial domains by conditioning on PDE model evaluations corresponding to approximate time derivatives at each discretisation grid location, dependent on sequentially sampled realisations from the predictive posterior distribution over the second order spatial derivatives. We describe the full algorithmic construction in the Appendix (Algorithm 3). In the following numerical simulations, we consider dynamics with κ = 1 and initial function u * (x, 0) = sin (xπ) , x ∈ [0, 1]. We firstly consider the probabilistic forward problem using a variety of discretisation grid sizes. In Figure 2 we compare the exact solution with the probabilistic solution on two different grids; a coarse discretisation of 15 points in the spatial domain and 50 points in the temporal domain, and a finer discretisation of 29 points in the spatial domain and 100 points in the temporal domain. We observe from the simulations that as the mesh size becomes finer, the uncertainty in the solution decreases, in agreement with the consistency result (Appendix, Section 8.6), where it is shown that the probabilistic solution should tend to the exact solution as the grid spacing tends to zero. Although this is an illustrative example using a simple toy system, the characterisation of spatial uncertainty is vital for more complex models, where there are computational constraints limiting the number of system evaluations that may be performed. We can see the effect of such discretisation uncertainty by performing posterior inference for the parameter κ given data simulated from an exact solution with κ = 1. Figure 3 shows the posterior distribution over κ obtained by using both a "forward in time, centred in space" (FTCS) finite difference solver and a probabilistic solver under a variety of discretisation grids. The use of a deterministic solver illustrates the problem of inferential bias and overconfident posterior variance that may occur if discretisation uncertainty is not taken into account and too coarse a grid is employed. In this illustrative setting, the use of a probabilistic solver propagates discretisation uncertainty in the solution through to the posterior distribution over the parameters. We obtain parameter estimates that assign positive probability mass to the true value of κ, even when using a coarsely discretised grid, avoiding the problem of overconfident parameter inferences that exclude the true value. MBVPs introduce challenges for many existing numerical solvers, which typically rely on optimisation and the theory of IVPs to estimate the unspecified initial state, u(a) in (2). The optimisation over u(a) is performed until the corresponding IVP solution satisfies the specified boundary condition u * (b) to within a user specified tolerance. For expositional simplicity, we consider a boundary value problem with one constraint located at each boundary of the domain, namely, v(a), u(b) = v * (a), u * (b) . We treat the boundary value, u * (b), as a data point and consider the solution as an inference problem over the unspecified initial value, u(a). The likelihood therefore defines the mismatch between the boundary value u(b), obtained from the realised probabilistic solution given some u(a), with the exact boundary value, u * (b), as follows, where m (u) (b) and C (u) (b, b) are the posterior mean and covariance for state u at time point b obtained from evaluating Algorithm 1. The posterior distribution of the states therefore has density, which exhibits multimodality over the states, as shown on the left side of Figure 4 . While deterministic numerical solvers rely on an ad hoc end point mismatch tolerance, the probabilistic framework naturally Figure 2: We illustrate the probabilistic output of the solution to the heat equation PDE, with κ = 1, integrated between t = 0 and t = 0.25 using two grid sizes; the coarser mesh (shown in blue) consists of 15 spatial discretisation points and 50 time discretisation points, the finer mesh consists of 29 spatial discretisation points and 100 time discretisation points. We show the spatial posterior predictions at three time points; t = 0.02 (top), t = 0.12 (middle) and t = 0.22 (bottom). The exact solution at each time point is represented by the green line. The error bars show the mean and 2 standard deviations for each of the probabilistic solutions calculated using 50 simulations. Figure 3: We illustrate the inverse problem by performing inference over the parameter κ in the heat equation, integrated between t = 0 and t = 0.25. We generate data over a grid of 8 spatial discretisation points and 25 time discretisation points by using the exact solution with κ = 1, then adding noise with standard deviation of 0.005. We firstly use the probabilistic differential equation solver (PODES) using three grid sizes; a coarse mesh consisting of 8 spatial discretisation points and 25 time discretisation points (far left), a finer mesh consisting of 15 spatial discretisation points and 50 time discretisation points (second from left), and a further finer mesh consisting of 29 spatial discretisation points and 100 time discretisation points (second from right). Note the change in scale as the posterior variance decreases with increasing resolution of the discretisation. As an illustrative comparison, we show the posterior distributions using a deterministic forward in time, centred in space (FTCS) integration scheme (far right). If the discretisation is not fine enough, we obtain an overconfident biased posterior that assigns negligible probability mass to the true value of κ. In contrast, use of the exact solution produces a perfectly unbiased posterior, as expected. Figure 4: 3,000 samples were drawn from the posterior probability density (13) and the trajectories are shown for both states (left, above and below). The marginal posterior probability density over the unknown initial condition u( 1 2 ) is also shown (right). defines that tolerance through the predictive distribution of (12). An important consideration for solving MBVPs probabilistically is to sample efficiently from potentially multimodal posteriors over the initial value. We therefore recommend the use of an appropriate MCMC scheme such as parallel tempering (Geyer, 1991) , which can quickly identify and explore disjoint regions of high posterior probability. We provide one such implementation of parallel tempering in the Appendix (Algorithm 5). As a demonstration of probabilistically solving a mixed boundary value problem, we consider a special case the Lane-Emden model, which is used to describe the density u of gaseous spherical objects, such as stars, as a function of its radius t, (Shampine, 2003) . We rewrite the canonical second order ODE as a system of a first order equations with a mixed boundary value, with boundary conditions, u * (1) = √ 3/2 and v * ( 1 2 ) = −288/2197. The unknown initial state u * ( 1 2 ) is assigned a diffuse Gaussian prior with mean 1.5 and standard deviation 2 u * (1) − v * ( 1 2 ) , which reflects the possibility that multiple solutions may be present over a wide range of initial states. In this example, we chose the squared exponential covariance to model what we expect to be a very smooth solution. The discretisation grid consists of 100 equally spaced points. The length-scale is set to twice the discretisation grid step size and the prior precision is set to 1. Figure 4 shows a posterior sample from Equation 13 and identifies two high density regions in the posterior corresponding to distinct trajectories that approximately satisfy model dynamics given our discretisation grid. Figure 4 also shows multimodality through the marginal posterior over the unknown initial state u( 1 2 ). This example illustrates the need for accurately modelling solution uncertainty in a functional manner. A numerical approximation, even with corresponding numerical error bounds, will fail to detect the existence of a second solution, whereas a probabilistic approach allows us to determine the number and location of possible solutions. In this case, a probabilistic approach to solving the differential equation problem addresses model bias caused by multiplicity of solutions. 14 We present the following example as a proof of concept that the probabilistic solver may be applied reliably and straightforwardly to a very high-dimensional dynamical system. The Navier-Stokes system is a fundamental model of fluid dynamics, incorporating laws of conservation of mass, energy and linear momentum, as well as physical properties of an incompressible fluid over some domain given constraints imposed along the boundaries. It is an important component of complex models in oceanography, weather, atmospheric pollution, and glacier movement. Despite its extensive use, the dynamics of Navier-Stokes models are poorly understood even at small time scales, where they can give rise to turbulence. We consider the Navier-Stokes PDE model for the time evolution of 2 components of the velocity, u : D → R 2 , of an incompressible fluid on a torus, D := [0, 2π, ) × [0, 2π], expressed in spherical coordinates. The Navier-Stokes boundary value problem is defined by: where ∆ := is the Laplacian operator such that ∆u = u x2x2 , and is the gradient operator such that ∇u = u x1 + u (15) is not known in closed form. Often, the quantity of interest is the local spinning motion of the incompressible fluid, called vorticity, which we define as, = −∇ × u, where ∇ × u represents the rotational curl defined as the cross product of ∇ and u, with positive vorticity corresponding to clockwise rotation. This variable will be used to better visualise the probabilistic solution of the Navier-Stokes system by reducing the two components of velocity to a one dimensional function. We discretize the Navier-Stokes model (15) over a grid of size 128 in each spatial dimension. Therefore, a pseudo spectral projection in Fourier space yields 16,384 coupled, stiff ODEs with associated constraints. Full details of the pseudo spectral projection are provided to allow full replication of these results and are available on the accompanying website. Figure 6 .3 shows four forward simulated vorticity trajectories (along rows), obtained from two components of velocity governed by the Navier-Stokes equations (15) at four distinct time points (along columns). Slight differences in the state dynamics can be seen at the last time point, where the four trajectories visibly diverge from one another. These differences express the epistemic uncertainty resulting from discretising the exact but unknown infinite dimensional solution. The probabilistic approach for forward simulation of DIFPs is described in Algorithm 6 in the Appendix. In addition to modelling discretisation uncertainty in a structured, functional way, our framework allows us to straightforwardly incorporate uncertainty in the initial function through the forward simulation; indeed, the initial function φ may iteself only be available at a finite number of nodes. The probabilistic approach quantifies the uncertainty associated with the estimation of φ(t) and propagates it recursively through the states. Even when φ(t) is fully specified in advance, the dependence of the current state on previously estimated states impacts the accuracy of the numerical solution using standard solvers, even over the short term. In the following example we account for the uncertainty associated with current and delayed estimates of system states, and fully address these potential sources of bias in the estimation of model parameters. 16 their initial state, returning to the cytoplasm to be used in the next activation cycle. This last stage is not well understood and is proxied in the model by the unknown time delay τ . The model for this mechanism describes changes in 4 reaction states of STAT-5 through the nonlinear DIFP, The initial function components φ (2) (t) = φ (3) (t) = φ (4) (t) are everywhere zero, while the constant initial function φ (1) (t) is unknown. The states for this system cannot be measured directly, but are observed through a nonlinear transformation, and parameterised by the unknown scaling factors k 5 and k 6 . The indirect measurements of the states are assumed contaminated with additive zero-mean Gaussian noise, ε(t), with experimentally determined standard deviations, Our analysis is based on experimental data measured at the locations, t, from Swameye et al. (2003) , which consists of 16 measurements for the first two states of the observation process. Raue et al. (2009) further utilise an additional artificial data point for each of the third and fourth observation process states to deal with lack of parameter identifiability for this system; we therefore adopt this assumption in our analysis. The forcing function, EpoR A : [0, 60] → R + , is not known, but measured at 16 discrete time points t EP O . As per Raue et al. (2009) , we assume that observations are measured without error. We further assume that this function shares the same smoothness as the solution state (piecewise linear first derivative). The full conditional distribution of the forcing function is given by a GP interpolation of the observations. Forward inference for this model will be used within the statistical inverse problem of recovering unknown parameters and first initial state, θ = [k 1 , . . . , k 6 , τ, u (1) (0)], from experimental data. We demonstrate fully probabilistic inference for state trajectories and parameters of the challenging 4 state delay initial function model (16) describing the dynamics of the JAK-STAT cellular signal transduction pathway (Raue et al., 2009 ). There have been several analyses of the JAK-STAT pathway mechanism based on this data (e.g., Campbell and Chkrebtii, 2013; Raue et al., 2009; Schmidl et al., 2003; Swameye et al., 2003) . Despite the variety of modelling assumptions considered by different authors, as well as distinct inference approaches, some interesting common features have been identified that motivate the explicit modelling of discretisation uncertainty for this application. Firstly, the inaccuracy and computational constraints of numerical techniques required to solve the system equations have led some authors to resort to coarse ODE approximations or even indirectly bypassing numerical solution via Generalised Smoothing. This motivates a formal analysis of the structure and propagation of discretisation error through the inverse problem. A further issue is that the model (16) and its variants suffer from model misspecification. The above studies 17 suggest the model is not flexible enough to fit the available data, however it is not clear how much of this misfit is due to model discrepancy, and how much may be attributed to discretisation error. Our analysis proceeds by defining prior distributions on the unknown parameters as follows, We obtained samples from the posterior distribution of the model parameters, θ = [k 1 , . . . , k 6 , τ, u (1) (0)], solution states, u(t, θ), and auxiliary variables, Ψ given the data y(t) using MCMC. In order to construct a Markov chain that efficiently traverses the parameter space of this multimodal posterior distribution, we employed a parallel tempering sampler (Geyer, 1991) with 10 parallel chains along a uniformly spaced temperature profile over the interval [0.5, 1]. Each probabilistic DIFP simulation was generated using Algorithm 6 under an equally spaced discretisation grid of size N = 500. Full algorithmic details are provided in the Appendix (Algorithm 8). We obtained two groups of posterior samples, each of size 50,000. Within chain convergence was assessed by testing for equality of means between disjoint iteration intervals of the chain (Geweke 1992). Between chain convergence was similarly assessed. Additionally we ensured that the acceptance rate fell roughly within the accepted range of 18%-28% for each of the two parameter blocks, and that the total acceptance rate for moves between any two chains remained roughly within 5%-15%. Correlation plots and kernel density estimates with priors for the marginal parameter posteriors are shown in Figure 6 . All parameters with the exception of the prior precision α are identified by the data, including the rate parameter k 2 which appears to be only weakly identified. We observe strong correlations between parameters, consistent with previous studies on this system. For example, there is strong correlation among the scaling parameters k 5 , k 6 and the initial first state u (1) (0). Interestingly, there appears to be a correlation between the probabilistic solver's length-scale λ and the first, third and fourth reaction rates. Furthermore, this correlation has a nonlinear structure, where the highest parameter density region seems to change with length scale implying strong sensitivity to the solver specifications. The length-scale is the probabilistic analogue, under a bounded covariance, of the step number in a numerical method. However, in analyses based on numerical integration, the choice of a numerical technique effectively fixes this parameter at an arbitrary value that is chosen a priori. Our result here suggests that for this problem, the inferred parameter values are highly and nonlinearly dependent on the choice of the numerical method used. We speculate that this effect may become quite serious for more sensitive systems, such that ignoring discretisation uncertainty within the inverse problem may result in inferential bias. A sample from the marginal posterior of state trajectories and the corresponding observation process are shown in Figure 7 . The error bars on the data points show two standard deviations of the measurement error from Swameye et al. (2003) . It is immediately clear that our model, which incorporates discretisation uncertainty in the forward problem, still does not fully capture the dynamics of the observed data. This systematic lack of fit suggests the existence of model discrepancy beyond that described by discretisation uncertainty. This paper has presented a probabilistic formalism to describe the structure of approximate solution uncertainty for general systems of differential equations. Rather than providing a single set of discrete function values that approximately satisfy the constraints imposed by the system of differential equations, our approach yields a probability measure over the space of such infinite dimensional functions. This is a departure : Marginal posterior distribution of the model parameters using probabilistic forward simulation based on a sample of size 50,000, generated using a parallel tempering algorithm with ten chains. Prior probability densities are shown in black. 19 from the existing accepted practice of employing a deterministic numerical integration code to obtain an approximate statistical error model as part of the inference process. This enables adoption of probability calculus for coherently propagating functional uncertainty when solving the statistical inverse problem. We have demonstrated that our probabilistic framework can already be applied to complex and high dimensional systems using the Navier-Stokes example. However, further work is needed to improve computational efficiency to a point where it is comparable with existing implementations of numerical integration codes. In principle, as the probabilistic integration method we have proposed relies on solutions of linear systems, it shares the same algorithmic scaling as most implicit numerical integration codes (e.g. Crank-Nicholson), and it is anticipated that, over time, algorithmic and code development will ensure our methodology attains similar levels of performance, becoming part of the standard toolbox for climate researchers, geoscientists, and engineers. Algorithmic advances and code development will be an exciting and fruitful area of ongoing investigation with high impact. Extension of the methodology to systems that exhibit directional derivative errors is immediate, for example in fluid dynamics applications. Indeed, in Section 4 we have modelled the mismatch between the prior on the derivative and the ODE model as Gaussian. Relaxing this assumption no longer guarantees a closed form representation of the updated prior on the state derivative, however this can be overcome by an additional layer of Monte Carlo sampling within Algorithm 1. Such questions would additionally motivate the development of efficient MCMC algorithms in this context. In this work we have been able to directly exploit the structure of the mathematical model within the inference process in a more informed manner than by treating the intractable system of differential equations and associated numerical solver code as a "black box". We may also exploit the intrinsic manifold structure induced by the probabilistic model posterior to increase efficiency of a variety of MCMC methods (Girolami and Calderhead, 2011). As we have already noted, controlling the auxiliary parameters Ψ associated with smoothness of the space of possible trajectories can also add flexibility to MCMC sampling methods. This feature of our probabilistic method allows one to define intermediate target densities for population MCMC methods such as Smooth Functional Tempering (Campbell and Steele, 2011), thereby increasing efficiency when exploring the complex posterior densities arising in applications described by systems of nonlinear ordinary or partial differential equations. The proposed probabilistic approach takes the form of a linear functional projection, allowing direct computation of sensitivities of the system states with respect to parameters. Local sensitivity analysis can be useful in engineering applications, inference, or in aiding optimisation. This would also allow, for example, our probabilistic framework to be incorporated into a classical nonlinear least squares framework (Bates and Watts, 2007) . Fundamentally, the probabilistic approach we advocate allows one to define a formal tradeoff between uncertainty induced by numerical accuracy and the size of the inverse problem, by choice of discretisation grid. This problem is of deep interest in the uncertainty quantification community (see, for example, Arridge et al., 2006; Kaipio et al., 2004) . We may now quantify and compare the individual contributions to overall uncertainty from discretisation, model misspecification, and Monte Carlo error within the probabilistic framework. Our probabilistic model of discretisation uncertainty may be used to guide mesh refinement and indeed mesh design for complex models, in that a predictive distribution is now available, which forms the basis of experiment design approaches to mesh development. Numerical solutions of ODEs and PDEs for complex models also rely on adaptive mesh selection. The probabilistic formalism presented here can now inform mesh selection in a probabilistic way, by optimising a chosen design criterion. Some work in this direction has already been conducted in Chkrebtii (2013) , where the natural Kullback-Leibler divergence criterion is successfully used to adaptively choose the discretisation grid probabilistically. Such an approach was also suggested by Skilling (1991) using the cross-entropy. Inference and prediction for computer experiments (Sacks et al., 1989 ) relies on numerical solutions of large-scale system models. Currently, numerical uncertainty in the model is largely ignored, although it is informally incorporated through a covariance nugget in the emulator (see, for example, Gramacy and Lee, 2012) . Adopting a probabilistic approach on a large scale will have practical implications in this area by permitting relaxation of the error-free assumption adopted when modelling computer code output, leading to more realistic and flexible emulators. In modelling uncertainty about an exact but unknown solution, the probabilistic integration formalism may be viewed as providing both an estimate of the exact solution and its associated functional error analysis. Indeed, many existing numerical methods can be interpreted in the context of this general probabilistic framework. This suggests the possibility to both generalise existing numerical solvers and to develop new sampling schemes to complement the probabilistic solutions we describe in this contribution. Chaotic systems arise in modelling a large variety of physical systems, including laser cavities, chemical reactions, fluid motion, crystal growth, weather prediction, and earthquake dynamics (Baker and Gollub, 1996, chapter 7). Extreme sensitivity to small perturbations characterises chaotic dynamics, where the effect of discretisation uncertainty becomes an important contribution to global solution uncertainty. Although the long term behaviour of the system is entrenched in the initial states, the presence of numerical discretisation error results in exponential divergence from the exact solution. Consequently, information about initial states rapidly decays as system solution evolves in time. This insight, demonstrated and explained by (Berliner, 1991) , is showcased in the following example. We consider the classical Lorenz initial value problem (Lorenz, 1963) , a deceptively simple three-state ODE model of convective fluid motion induced by a temperature difference between an upper and lower surface. A sample from the probabilistic posterior of this system is shown in Figure 8 given model parameters in the chaotic regime and a fixed initial state. There is a short time window within which there is negligible uncertainty in the solution, but the accumulation of uncertainty quickly results in divergent, yet highly structured, solutions as the flow is restricted around a chaotic attractor. This example highlights the need for a functional model of discretization uncertainty to replace numerical pointwise error bounds. Within the inverse problem, a probabilistic approach also allows us to quantify the loss of model information as we move away from the initial state. The study of chaotic systems is clearly a very important area of research, and it is anticipated that the probabilistic approach proposed may help to accelerate progress. We now return to an important discussion regarding the class of problems for which probabilistic integration is well suited. As we have seen with inference for the JAK-STAT model, even relatively small discretisation uncertainty can become amplified through the nonlinear forward model and may introduce bias in the estimation of parameters to which the trajectory is highly sensitive. Nevertheless, it may be conjectured that inference for low dimensional, stable systems, without any topological restrictions on the solution may not suffer from significant discretisation effects. However, the benefits of our formulation undoubtedly lie at the frontier of research in uncertainty quantification, dealing with massive nonlinear models, exhibiting complex or even chaotic dynamics, and strong spatial-geometric effects (e.g. subsurface flow models). Indeed, solving such systems approximately is a problem that is still at the edge of current research in numerical analysis. We suggest that a probabilistic approach would provide an important additional set of tools for the study of such systems. In conclusion, we hope that our paper will encourage much more research at this exciting interface between mathematical analysis and statistical methodology.
A similar approach is followed in REF , yielding a nonparametric posterior rather than a Gaussian approximation.
Bayesian Solution Uncertainty Quantification for Differential Equations
arXiv: Methodology
Mathematics