Probabilistic Inference on Noisy Time Series (PINTS)

Time series models are ubiquitous in science, arising in any situation where researchers seek to understand how a system's behaviour changes over time. A key problem in time series modelling is \emph{inference}; determining properties of the underlying system based on observed time series. For both statistical and mechanistic models, inference involves finding parameter values, or distributions of parameters values, for which model outputs are consistent with observations. A wide variety of inference techniques are available and different approaches are suitable for different classes of problems. This variety presents a challenge for researchers, who may not have the resources or expertise to implement and experiment with these methods. PINTS (Probabilistic Inference on Noisy Time Series - https://github.com/pints-team/pints is an open-source (BSD 3-clause license) Python library that provides researchers with a broad suite of non-linear optimisation and sampling methods. It allows users to wrap a model and data in a transparent and straightforward interface, which can then be used with custom or pre-defined error measures for optimisation, or with likelihood functions for Bayesian inference or maximum-likelihood estimation. Derivative-free optimisation algorithms - which work without harder-to-obtain gradient information - are included, as well as inference algorithms such as adaptive Markov chain Monte Carlo and nested sampling which estimate distributions over parameter values. By making these statistical techniques available in an open and easy-to-use framework, PINTS brings the power of modern statistical techniques to a wider scientific audience.


Introduction
Time series models are common in science, where they are used to describe the dynamics of system behaviours. In many cases, these models are non-linear and impossible to solve analytically, so that the forward problem (predicting the model output for a given set of parameters) is computationally hard. For such models, there is no single method which can reliably solve the inverse problem of estimating parameter values from a noisy time trace. Much like there is a variety of forward models, there is a diversity of approaches for parameter inference. Further, it is often unclear which approach to apply when, meaning that researchers are required to implement a range of methods before successfully fitting their model to data.
PINTS is a software framework that allows users to easily trial and apply different inference methods to their problem. The inference methods supplied by PINTS fall into two broad categories: optimisers, which attempt to find a single best parameter vector, and samplers, which aim to estimate a probability distribution over parameter values that are compatible with observed results. Users are expected to already have a forward model (for example, a simulation) at their disposal, which they make available to PINTS by writing a simple Python wrapper. They then define a Problem (a forward model plus a data set), on which either an ErrorMeasure (for optimisation) or a LogPDF (for optimisation and sampling) is defined. Currently available optimisers include CMA-ES [6], XNES [5], SNES [21], and Particle Swarm Optimisation (PSO) [12]. Sampling methods include Random Walk Metropolis Markov chain Monte Carlo (MCMC) [13,15], adaptive covariance MCMC [10], Population MCMC [8], Differential Evolution [23], DREAM [25], emcee [3], Hamiltonian MCMC [17], and MALA [4]. In addition, ellipsoidal [16] and rejection nested samplers [22] are provided. Convenience plotting methods are provided to quickly visualise the results, as well as diagnostic tools to inspect the validity of the results. An example of an optimisation problem and its solution using PINTS is shown in Figure 1.
PINTS was developed as a community effort by researchers in electrochemistry, cardiac electrophysiology, and statistics, to compare different methods for solving inverse problems in a common framework. It features a clean and transparent object-orientated API that is designed to easily accommodate new error measures, loglikelihoods, optimisers and samplers, allowing users to utilise pre-built components as much as possible, while adding their own code for problem-specific areas. The PINTS team aims for full test coverage, and includes unit testing and extended statistical tests to verify the correct operation of all methods.
Early research using PINTS in electrochemistry has included fitting a differential-algebraic equation (DAE) model of reduction-oxidation to voltammetry measurements of a Polyoxometalates molecule [19], and the design and application of a custom hierarchical statistical model for repeat voltammetry experiments of a Ferricynide process [18].
Optimisation algorithms are implemented in many different software packages, (see, for example, the Python scipy.optimize module), but are often biased towards gradient-based methods, which can perform poorly for many ordinary and partial differential equations used in time series modelling. PINTS therefore focuses on derivative-free optimisers, although we plan to add gradient-based methods for comparison. In contrast with more general-purpose optimisation software, PINTS contains a number of error measures specifically suited to time series models, and adds the ability to use any PINTS log-likelihood class as an error measure in order to perform maximum likelihood estimation (MLE).
Dakota [1] is a widely regarded package for parameter fitting and uncertainty quantification and is most similar to PINTS in that it offers a generic interface to call an (assumed expensive) model, as well as a wide variety of optimisers and samplers. In contrast to the PINTS Python API, Dakota uses either a C++ or file input/output process for communication between user models and the library, and does not provide options for specifying either the error measures or the log-likelihoods of the inverse problem. However, Dakota has some features not yet available in PINTS, such as the option to train a surrogate model, useful for very expensive model evaluations.
Other software packages that enable parameter inference and sampling for ODE models include BioBayes [26], ABC-SysBio [14], SYSBIONS [9], PyMC3 [20], and optimisation procedure with PINTS. Note that the actual simulation code is omitted from the example model wrapper at the top: this is the user-provided part, and can be written in Python or any other language that can interface with Python, allowing computationally heavy forward simulations to be handled entirely outside of PINTS. This image, and the full example code, can also be found in the PINTS repository.
Stan [2]. These packages use either a common model description format (for example, SBML) or their own language (for example, Stan's probabilistic programming language) to specify the model, presenting additional learning hurdles for a user and often restricting the class of models which can be fit. By contrast, PINTS aims to be as general as possible to support a wider variety of models (for example, PDEs). PyMC3 [20] does provide a similar generic model interface to PINTS but, as with the other packages, specialises in one sampling method, whereas PINTS aims to support a wide variety of methods with the assumption that no one sampling method is suitable for all models of interest. BCM [24] offers both a generic interface (via C++) and a wide variety of samplers, but does not supply any likelihood functions and is unfortunately largely undocumented.

Implementation and architecture
PINTS is designed around two core ideas: 1. PINTS should work with a wide range of time series models, and make no demands on how they are implemented other than a minimal input/output interface. 2. It is assumed that model evaluation (simulation) is the most costly step in any optimisation or sampling routine. The decision to use Python fits both these criteria: Python interfaces well with C and C++, which are typically used for high-performance simulation, and any performance hit of using the high-level, easy to read and write language Python is overshadowed by simulation time.

Defining an optimisation or sampling problem
All optimisers operate on a callable ErrorMeasure object that describes a function to minimise or on a callable LogPDF object that describes a probability density function (PDF) to maximise. Similarly, all samplers start from a callable LogPDF, so that the same probability function can be used with both optimisers and samplers. The natural logarithm of the PDF is used for computational efficiency and accuracy, and we allow the probability density to be unnormalised (i.e. its integral does not have to sum to 1). Figure 2 shows how a user-defined model can be wrapped in a PINTS ForwardModel and combined with time points and measured values to create a Problem from which several standard ErrorMeasures and LogPDFs can be created. For inference in a Bayesian context, a LogPosterior class and several LogPrior distributions are provided. If a given LogPDF or ErrorMeasure cannot be constructed from PINTS classes, users can also define their own classes.

Implementation of optimisers and samplers
Most PINTS samplers and optimisers are implemented using a so-called ask-and-tell interface, inspired by the Python implementation of CMA-ES [6] (https://github. com/CMA-ES/pycma). In this framework, the details of solving the forward problem are partitioned away from the rest of the sampling or optimising algorithm. For each iteration the following steps are undertaken (Figure 3): first, the user calls ask() to obtain one or more parameter values from their chosen method -these values are typically vectors generated stochastically conditional on an internal system state; second, the user solves the forward model and generates a score for each parameter vector, for example, an error measure or (unnormalised) posterior probability; third, the user calls tell() to pass the score back to the method, which can then update its internal state and finish the iteration. For example, in Figure 2: An overview of the main PINTS classes used to define error measures (for optimisation) and PDFs (for optimisation or sampling). Users write a wrapper class for their model, making it available to pints, and must provide the experimental data using any Python sequence structure (for example, a list or a NumPy array). With these ingredients, a (single or multi-output) problem can be defined that can then be used with any of the available error measure or likelihood classes. Alternatively, users implement their own ErrorMeasure or LogPDF, which allows for further customisation and for problems other than time series problems to be solved. many MCMC methods each ask() call returns a single proposed sample to be evaluated by the user, and the following tell() then either accepts or rejects this point based on its probability. For optimisers such as CMA-ES, ask() returns a set of points in parameter space, and the scores passed in via tell() are then used to estimate the local gradient, which is used to move towards the estimated optimum. This framework has a number of advantages: since optimisation or sampling can take hours or even days, this allows programs using PINTS to provide regular user feedback and logging (which is not possible when the routine is implemented as a single monolithic function call); allows users with access to CPU clusters or GPU machines to implement their own parallelised evaluation of ErrorMeasures and LogPDFs; lets users implement their own strategies (for example, by dynamically changing hyperparameters) and/or stopping criteria; and, finally, by delineating the sampling or optimisation algorithm's steps from the methods used to solve the forward problem, encourages development of transparent and modular code.

Running optimisation and sampling
For more casual users (whom we expect will be the majority), PINTS includes ' controller' classes (e.g. OptimisationController or MCMCController) that provide a higher level interface (see e.g. Figure  1). When creating a controller, the user specifies the name of the lower-level method to use (e.g. CMAES, or AdaptiveCovarianceMCMC) and passes in an ErrorMeasure or LogPDF. The controller then instantiates and manages this method, as well as providing user-configurable stopping criteria (e.g. maximum number of iterations) and logging to screen and/or the filesystem. An important difference with the lower-level ask-and-tell interface, is that controllers handle function evaluation, and can distribute evaluations over multiple CPU cores using Python multiprocessing. This provides users with out-of-the-box parallelisation, while users looking to implement custom parallelisation can fall back on the ask-and-tell interface. (2) Availability Operating system PINTS uses no functions specific to any operating system (OS) and so can run on any OS that provides Python. Optional parallelisation is provided that uses the Python multiprocessing module, which works best on UNIX-based systems (for example, Linux and OS/X), but runs on Windows with slightly reduced performance.

Programming language
PINTS requires Python 2.7 or higher, or Python 3.4 or higher.

Additional system requirements
PINTS has a minimal disk space footprint (approximately 2MB) and can be run on single-processor devices or headless on multi-processor machines (for example, via ssh).

Dependencies
PINTS uses the NumPy (version 1.8 or higher) and SciPy (version 0.14 or higher, [11]) libraries extensively. The default optimisation method is CMA-ES, for which the cma package (version 2 or higher) is used. The remaining Figure 3: The three steps iterated in an ask-and-tell interface. The stars here represent code specific to the chosen sampling or optimiser method and θ′ is the input parameter vector proposed by ask().
The state(.) of the system varies according to the method but typically holds a set of input parameter vectors and other constant or dynamic variables used by the ask() and tell() steps.