Random time step probabilistic methods for uncertainty quantification in chaotic and geometric numerical integration

A novel probabilistic numerical method for quantifying the uncertainty induced by the time integration of ordinary differential equations (ODEs) is introduced. Departing from the classical strategy to randomise ODE solvers by adding a random forcing term, we show that a probability measure over the numerical solution of ODEs can be obtained by introducing suitable random time steps in a classical time integrator. This intrinsic randomisation allows for the conservation of geometric properties of the underlying deterministic integrator such as mass conservation, symplecticity or conservation of first integrals. Weak and mean square convergence analysis is derived. We also analyse the convergence of the Monte Carlo estimator for the proposed random time step method and show that the measure obtained with repeated sampling converges in the mean square sense independently of the number of samples. Numerical examples including chaotic Hamiltonian systems, chemical reactions and Bayesian inferential problems illustrate the accuracy, robustness and versatility of our probabilistic numerical method.     

Laplace based approximate posterior inference for differential equation models

Ordinary differential equations are arguably the most popular and useful mathematical tool for describing physical and biological processes in the real world. Often, these physical and biological processes are observed with errors, in which case the most natural way to model such data is via regression where the mean function is defined by an ordinary differential equation believed to provide an understanding of the underlying process. These regression based dynamical models are called differential equation models. Parameter inference from differential equation models poses computational challenges mainly due to the fact that analytic solutions to most differential equations are not available. In this paper, we propose an approximation method for obtaining the posterior distribution of parameters in differential equation models. The approximation is done in two steps. In the first step, the solution of a differential equation is approximated by the general one-step method which is a class of numerical numerical methods for ordinary differential equations including the Euler and the Runge-Kutta procedures; in the second step, nuisance parameters are marginalized using Laplace approximation. The proposed Laplace approximated posterior gives a computationally fast alternative to the full Bayesian computational scheme (such as Makov Chain Monte Carlo) and produces more accurate and stable estimators than the popular smoothing methods (called collocation methods) based on frequentist procedures. For a theoretical support of the proposed method, we prove that the Laplace approximated posterior converges to the actual posterior under certain conditions and analyze the relation between the order of numerical error and its Laplace approximation. The proposed method is tested on simulated data sets and compared with the other existing methods.     

Parameter estimation of complex mixed models based on meta-model approach

Complex biological processes are usually experimented along time among a collection of individuals, longitudinal data are then available. The statistical challenge is to better understand the underlying biological mechanisms. A standard statistical approach is mixed-effects model where the regression function is highly-developed to describe precisely the biological processes (solutions of multi-dimensional ordinary differential equations or of partial differential equation). A classical estimation method relies on coupling a stochastic version of the EM algorithm with a Monte Carlo Markov Chain algorithm. This algorithm requires many evaluations of the regression function. This is clearly prohibitive when the solution is numerically approximated with a time-consuming solver. In this paper a meta-model relying on a Gaussian process emulator is proposed to approximate the regression function, that leads to what is called a mixed meta-model. The uncertainty of the meta-model approximation can be incorporated in the model. A control on the distance between the maximum likelihood estimates of the mixed meta-model and the maximum likelihood estimates of the exact mixed model is guaranteed. Eventually, numerical simulations are performed to illustrate the efficiency of this approach.     

Statistical analysis of error for fourth-order ordinary differential equation solvers

We develop an autoregressive integrated moving average (ARIMA) model to study the statistical behavior of the numerical error generated from three fourth-order ordinary differential equation solvers: Milne's method, Adams–Bashforth method and a new method that randomly switches between the Milne and Adams–Bashforth methods. With the actual error data based on three differential equations, we desire to identify an ARIMA model for each data series. Results show that some of the data series can be described by ARIMA models but others cannot. Based on the mathematical form of the numerical error, other statistical models should be investigated in the future. Finally, we assess the multivariate normality of the sample mean error generated by the switching method.     

Estimating stellar-oscillation-related parameters and their uncertainties with the moment method

Summary.  The moment method is a well-known astronomical mode identification technique in asteroseismology which uses a time series of the first three moments of a spectral line to estimate the discrete oscillation mode parameters l and m . The method, in contrast with many other mode identification techniques, also provides estimates of other important continuous parameters such as the inclination angle α and the rotational velocity v _e. We developed a statistical formalism for the moment method based on so-called generalized estimating equations. This formalism allows an estimation of the uncertainty of the continuous parameters, taking into account that the different moments of a line profile are correlated and that the uncertainty of the observed moments also depends on the model parameters. Furthermore, we set up a procedure to take into account the mode uncertainty, i.e. the fact that often several modes ( l ,  m ) can adequately describe the data. We also introduce a new lack-of-fit function which works at least as well as a previous discriminant function, and which in addition allows us to identify the sign of the azimuthal order m . We applied our method to star HD181558 by using several numerical methods, from which we learned that numerically solving the estimating equations is an intensive task. We report on the numerical results, from which we gain insight in the statistical uncertainties of the physical parameters that are involved in the moment method.     

Strong convergence rates of probabilistic integrators for ordinary differential equations

Probabilistic integration of a continuous dynamical system is a way of systematically introducing discretisation error, at scales no larger than errors introduced by standard numerical discretisation, in order to enable thorough exploration of possible responses of the system to inputs. It is thus a potentially useful approach in a number of applications such as forward uncertainty quantification, inverse problems, and data assimilation. We extend the convergence analysis of probabilistic integrators for deterministic ordinary differential equations, as proposed by Conrad et al. (Stat Comput 27(4):1065–1082, 2017. https://doi.org/10.1007/s11222-016-9671-0), to establish mean-square convergence in the uniform norm on discrete- or continuous-time solutions under relaxed regularity assumptions on the driving vector fields and their induced flows. Specifically, we show that randomised high-order integrators for globally Lipschitz flows and randomised Euler integrators for dissipative vector fields with polynomially bounded local Lipschitz constants all have the same mean-square convergence rate as their deterministic counterparts, provided that the variance of the integration noise is not of higher order than the corresponding deterministic integrator. These and similar results are proven for probabilistic integrators where the random perturbations may be state-dependent, non-Gaussian, or non-centred random variables.

     

Smooth functional tempering for nonlinear differential equation models

Differential equations are used in modeling diverse system behaviors in a wide variety of sciences. Methods for estimating the differential equation parameters traditionally depend on the inclusion of initial system states and numerically solving the equations. This paper presents Smooth Functional Tempering, a new population Markov Chain Monte Carlo approach for posterior estimation of parameters. The proposed method borrows insights from parallel tempering and model based smoothing to define a sequence of approximations to the posterior. The tempered approximations depend on relaxations of the solution to the differential equation model, reducing the need for estimating the initial system states and obtaining a numerical differential equation solution. Rather than tempering via approximations to the posterior that are more heavily rooted in the prior, this new method tempers towards data features. Using our proposed approach, we observed faster convergence and robustness to both initial values and prior distributions that do not reflect the features of the data. Two variations of the method are proposed and their performance is examined through simulation studies and a real application to the chemical reaction dynamics of producing nylon.     

Stochastic models for greenhouse gas emission rate estimation from hydroelectric reservoirs: a Bayesian hierarchical approach

Herein, we propose a fully Bayesian approach to the greenhouse gas emission problem. The goal of this work is to estimate the emission rate of polluting gases from the area flooded by hydroelectric reservoirs. We present models for gas concentration evolution in two ways: first, by proposing them from ordinary differential equation solutions and, second, by using stochastic differential equations with a discretization scheme. Finally, we present techniques to estimate the emission rate for the entire reservoir. In order to carry out the inference, we use the Bayesian framework with Monte Carlo via Markov Chain methods. Discretization schemes over continuous differential equations are used when necessary. These models applied to greenhouse gas emission and Bayesian inference for this purpose are completely new in statistical literature, as far as we know, and contribute to estimate the amount of polluting gases released from hydroelectric reservoirs in Brazil. The proposed models are applied in a real data set and results are presented.     

Confidence intervals for finite difference solutions

Although applications of Bayesian analysis for numerical quadrature problems have been considered before, it is only very recently that statisticians have focused on the connections between statistics and numerical analysis of differential equations. In line with this very recent trend, we show how certain commonly used finite difference schemes for numerical solutions of ordinary and partial differential equations can be considered in a regression setting. Focusing on this regression framework, we apply a simple Bayesian strategy to obtain confidence intervals for the finite difference solutions. We apply this framework on several examples to show how the confidence intervals are related to truncation error and illustrate the utility of the confidence intervals for the examples considered.     

Estimable functions in the nonlinear model

When the method of least squares is used to estimate the parameters in a general model and the generated system of normal equations is linearly dependent, the estimate of the vector of parameters which satisfies the criterion is not unique. However, there exist certain functions of the estimated vector of parameters which are invariant to the least squares solution obtained from the normal equations. We define those invariant functions to be estimable, and present a technique to determine the functions of the parameters which are estimable for the general model. The method results in solving either a linear first order partial differential equation or a system of linear first order partial differential equations corresponding, respectively, to a single or multiple dependency between columns of the Jacobian matrix of the mean of the model. The usual results concerning estimability for linear models are a special case of the general results developed.     

Stochastic differential equation systems for an SIS epidemic model with vaccination and immigration

Two Itô stochastic differential equation (SDE) systems are constructed for a Susceptible-Infected-Susceptible epidemic model with temporary vaccination. A constant number of new members enter the population and total size of the population is variable. Some conditions for disease extinction in the stochastic models are established and compared with conditions in deterministic one. It is shown that the two stochastic models are equivalent in the sense that their solutions come from same distribution. In addition, the SDE models are simulated and the equivalence of the two stochastic models is confirmed by numerical examples. The probability distribution for extinction is also obtained numerically, provided there exists a probability for disease persistence whereas the expected duration of epidemic is acquired when extinction occurs with probability 1.     

Non-identifiable parametric probability models and reparametrization

Identifiability is a primary assumption in virtually all classical statistical theory. However, such an assumption may be violated in a variety of statistical models. We consider parametric models where the assumption of identifiability is violated, but otherwise satisfy standard assumptions. We propose an analytic method for constructing new parameters under which the model will be at least locally identifiable. This method is based on solving a system of linear partial differential equations involving the Fisher information matrix. Some consequences and valid inference procedures under non-identifiability have been discussed. The method of reparametrization is illustrated with an example.     

Importance sampling for Kolmogorov backward equations

The solution of the Kolmogorov backward equation is expressed as a functional integral by means of the Feynman–Kac formula. The expectation value is approximated as a mean over trajectories. In order to reduce the variance of the estimate, importance sampling is utilized. From the optimal importance density, a modified drift function is derived which is used to simulate optimal trajectories from an Itô equation. The method is applied to option pricing and the simulation of transition densities and likelihoods for diffusion processes. The results are compared to known exact solutions and results obtained by numerical integration of the path integral using Euler transition kernels. The importance sampling leads to strong variance reduction, even if the unknown solution appearing in the drift is replaced by known reference solutions. In models with low-dimensional state space, the numerical integration method is more efficient, but in higher dimensions it soon becomes infeasible, whereas the Monte Carlo method still works.     

Variance-based sensitivity analysis with the uncertainties of the input variables and their distribution parameters

For the structural systems with both the uncertainties of input variables and their distribution parameters, three sensitivity indices are proposed to measure the influence of input variables, distribution parameters and their interactive effects. With those sensitivity indices, analysts can make a decision that whether it is worth to accumulate data of one distribution parameter to reduce its uncertainty. Due to the large computational cost, the analytical solutions are derived for quadratic polynomial output responses. Whereas for the complex models, state dependent parameter (SDP) method is utilized to solve the proposed sensitivity indices efficiently.     

A Nonstationary Statistical Model for Computationally Intensive Numerical Ordinary Differential Systems

ABSTRACT

Computer models depending on unknown inputs are used in computer experiments in order to study the input-output relationship. Some computer models are computationally intensive and only a few computer runs can be made. A nonstationary statistical model that is used as a faster-running surrogate for computationally intensive numerical solvers of ordinary differential systems is proposed in this article. This statistical model reflects the dynamics of the system and, as a population dynamics example will show, it can be more accurate than a commonly used black box statistical model.     

Method for simulating non-linear stochastic differential equations in ℝ1

Very few specific stochastic differential equations have explicitly known solutions. The most common procedure to obtain a simulated path of a solution is based on a discretization of the stochastic differential equations. However, there are some cases where the discrete-time discretization cannot be used. In this article, we propose a new method to simulate the solution of a non-linear stochastic differential equation, which, in principle, is exempt from error of simulation and can be widely applied including in cases where the discrete-time discretization cannot be used.     

Parameter estimation for semiparametric ordinary differential equation models

Abstract

We propose a new class of two-stage parameter estimation methods for semiparametric ordinary differential equation (ODE) models. In the first stage, state variables are estimated using a penalized spline approach; In the second stage, form of numerical discretization algorithms for an ODE solver is used to formulate estimating equations. Estimated state variables from the first stage are used to obtain more data points for the second stage. Asymptotic properties for the proposed estimators are established. Simulation studies show that the method performs well, especially for small sample. Real life use of the method is illustrated using Influenza specific cell-trafficking study.     

Biases and Uncertainty in Climate Projections

Abstract. We study statistical procedures to quantify uncertainty in multivariate climate projections based on several deterministic climate models. We introduce two different assumptions – called constant bias and constant relation respectively – for extrapolating the substantial additive and multiplicative biases present during the control period to the scenario period. There are also strong indications that the biases in the scenario period are different from the extrapolations from the control period. Including such changes in the statistical models leads to an identifiability problem that we solve in a frequentist analysis using a zero sum side condition and in a Bayesian analysis using informative priors. The Bayesian analysis provides estimates of the uncertainty in the parameter estimates and takes this uncertainty into account for the predictive distributions. We illustrate the method by analysing projections of seasonal temperature and precipitation in the Alpine region from five regional climate models in the PRUDENCE project.     

Sieve Estimation of Constant and Time-Varying Coefficients in Nonlinear Ordinary Differential Equation Models by Considering Both Numerical Error and Measurement Error

This article considers estimation of constant and time-varying coefficients in nonlinear ordinary differential equation (ODE) models where analytic closed-form solutions are not available. The numerical solution-based nonlinear least squares (NLS) estimator is investigated in this study. A numerical algorithm such as the Runge-Kutta method is used to approximate the ODE solution. The asymptotic properties are established for the proposed estimators considering both numerical error and measurement error. The B-spline is used to approximate the time-varying coefficients, and the corresponding asymptotic theories in this case are investigated under the framework of the sieve approach. Our results show that if the maximum step size of the p-order numerical algorithm goes to zero at a rate faster than n(-1/(p∧4)), the numerical error is negligible compared to the measurement error. This result provides a theoretical guidance in selection of the step size for numerical evaluations of ODEs. Moreover, we have shown that the numerical solution-based NLS estimator and the sieve NLS estimator are strongly consistent. The sieve estimator of constant parameters is asymptotically normal with the same asymptotic co-variance as that of the case where the true ODE solution is exactly known, while the estimator of the time-varying parameter has the optimal convergence rate under some regularity conditions. The theoretical results are also developed for the case when the step size of the ODE numerical solver does not go to zero fast enough or the numerical error is comparable to the measurement error. We illustrate our approach with both simulation studies and clinical data on HIV viral dynamics.