key: cord-0147334-hu8uio07
authors: Sherlock, Chris
title: Direct statistical inference for finite Markov jump processes via the matrix exponential
date: 2018-09-19
journal: nan
DOI: nan
sha: 7e535535eca890aa4f3f8c4b989e121e208818eb
doc_id: 147334
cord_uid: hu8uio07

Given noisy, partial observations of a time-homogeneous, finite-statespace Markov chain, conceptually simple, direct statistical inference is available, in theory, via its rate matrix, or infinitesimal generator, $mathsf{Q}$, since $exp (mathsf{Q}t)$ is the transition matrix over time $t$. However, perhaps because of inadequate tools for matrix exponentiation in programming languages commonly used amongst statisticians or a belief that the necessary calculations are prohibitively expensive, statistical inference for continuous-time Markov chains with a large but finite state space is typically conducted via particle MCMC or other relatively complex inference schemes. When, as in many applications $mathsf{Q}$ arises from a reaction network, it is usually sparse. We describe variations on known algorithms which allow fast, robust and accurate evaluation of the product of a non-negative vector with the exponential of a large, sparse rate matrix. Our implementation uses relatively recently developed, efficient, linear algebra tools that take advantage of such sparsity. We demonstrate the straightforward statistical application of the key algorithm on a model for the mixing of two alleles in a population and on the Susceptible-Infectious-Removed epidemic model.

A reaction network is a stochastic model for the joint evolution of one or more populations of species. These species may be chemical or biological species (e.g. Wilkinson, 2012) , animal species (e.g. Drovandi and McCutchan, 2016) , interacting groups of individuals at various stages of a disease (e.g. Andersson and Britton, 2000) , or counts of sub-populations of alleles (e.g. Moran, 1958) , for example. The state of the system is encapsulated by the number of each species that is present, and the system evolves via a set of reactions: Poisson processes whose rates depend on the current state.

Typically, partial and/or noisy observations of the state are available at a set of time points, and statistical interest lies in inference on the unknown rate parameters, the filtering estimate of the state of the system after the latest observation or prediction of the future evolution of the system. The usual method of choice for exact inference on discretely observed Markov jump processes (MJPs) on a finite or countably infinite state space is Bayesian inference via particle Markov chain Monte Carlo (particle MCMC, Andrieu et al., 2010) using a bootstrap particle filter (e.g. Andrieu et al., 2009; Golightly and Wilkinson, 2011; Wilkinson, 2012; McKinley et al., 2014; Owen et al., 2015; Koblents and Miguez, 2015) . Other MCMC and SMC-based techniques are available (e.g. Kypraios et al., 2017) , and, a further latent-variable-based MCMC method when the statespace is finite (Rao and Teh, 2013) .

Particle MCMC and SMC, however, are relatively complex algorithms, even more so when a bootstrap particle filter (simulation from the process itself) is not suitable and a bridge simulator is necessary, such as when observation noise is small or when there is considerable variability in the state from one observation to the next (Golightly and Wilkinson, 2015; Golightly and Sherlock, 2019) . In cases where the number of states, d, is finite, direct exact likelihood-based inference is available via the exponential of the infinitesimal generator for the continuous-time Markov chain, or rate matrix, Q. Whilst such inference is conceptually straightforward, it has been avoided in practice except in cases where the number of states is very small (e.g. Amoros et al., 2019) . Matrix exponentiation has a computational cost of O(d 3 ), which, together with a lack of suitable tools in R, could explain the lack of uptake of this method. However, conceptually simple statistical inference via the matrix exponential is entirely practical in many cases even when the number of states is in the thousands or higher, for three main reasons:

1. Matrix exponentials themselves are never needed; only the product of a vector and a matrix exponential is ever required.

2. The matrices to be exponentiated are infinitesimal generators and, as such, have a special structure; furthermore, the vector that pre-multiplies the matrix exponential is non-negative.

3. The matrices to be exponentiated are sparse; tools for basic operations with large, sparse matrices in C++ and interfacing the resulting code with R have recently become widely available (Eddelbuettel and Sanderson, 2014; Sanderson and Curtin, 2018) .

The sparsity of Q arises because the number of possible 'next' states given the current state is bounded by the number of reactions, which is typically small. This article describes matrix exponential algorithms suitable for statistical application in many cases, and demonstrates their use for inference, filtering and prediction. Associated code provides easy-to-use R interfaces to C++ implementations of the algorithms, which are typically simpler and often faster than more generally applicable algorithms for matrix exponentiation.

Section 1.1 describes the Susceptible-Infectious-Removed (SIR) model for the evolution of an infectious disease and the Moran model for the mixing of two alleles in a population, then briefly mentions many more such models where the statespace is finite, and a few where it is countably infinite. The two main examples will be used to benchmark and illustrate the techniques in this article. As well as being directly of use for models with finite state spaces, exponentials of finite rate matrices can also be used to perform inference on Markov jump processes with a countably infinite statespace; see Georgoulas et al. (2017) and Sherlock and Golightly (2019) . The latter uses the uniformisation and scaling and squaring algorithms as described in this article, while the former uses the less efficient but more general algorithm of Al-Mohy and Higham (2011) (see Section 3).

Section 2 of this article presents the likelihood for discretely and partially observed data on a finite-statespace continuous-time Markov chain and presents two 'tricks' specific to epidemic models, that allow for a massive reduction in the size of the generators that are needed compared with the size of the statespace. Section 3 describes the Matrix exponential algorithms and Section 4 benchmarks some of the algorithms and demonstrates their use for inference, filtering and prediction. The article concludes in Section 5 with a discussion.

Both by way of motivation and because we shall use them later to illustrate our method, we now present two examples of continuous-time Markov processes, where a finite, sparse rate matrix contains all of the information about the dynamics.

For each Markov process, the set of possible states can be placed in one-to-one correspondance with a subset of the non-negative integers {1, . . . , d}. The off-diagonal elements of the rate matrix, Q, are all non-negative, and the ith diagonal element is Q ii = − d j=1,j =i Q i,j . A chain that is currently in state i leaves this state upon the first event of a Poisson process with a rate of −Q i,i ; the state to which it transitions is j with a probability of Q i,j /(−Q i,i ). Whilst the rate matrix, Q, is a natural description of the process, the likelihood for typical observation regimes involves the transition matrix, exp(Qt), the (i, j)th element of which is exactly P (X t = j|X 0 = i).

Example 1. The SIR model for epidemics. The SIR model for a disease epidemic has 3 species: those who are susceptible to the epidemic, S, those both infected and infectious, I, and those who have recovered from the epidemic and play no further part in the dynamics, R. The non-negative counts of each species are denoted by S, I, and R. For relatively short epidemics the population, n pop , is assumed to be fixed, and so the state of the Markov chain, represented by (S, I), is subject to the additional constraint of S + I ≤ n pop . The two possible reactions and their associated rates are:

Example 2. The Moran model for allele frequency descibes the time evolution of the frequency of two alleles, A 1 and A 2 in a population with a fixed size of n pop . Individuals with allele A 1 reproduce at a rate of α, and those with A 2 reproduce at a rate of β. When an individual dies it is replaced by the offspring of a parent chosen uniformly at random from the whole population (including the individual that dies). The allele that the parent passes to the offspring usually matches its own, however as it is passed down an allele may mutate; allele A 1 switching to A 2 with a probability of u and A 2 switching to A 1 with a probability of v. Let A 1 and A 2 represent individuals with alleles A 1 and A 2 respectively and let N be the number of individuals with allele A 1 . The two reactions are

Setting f N = N/n pop , the corresponding infinitesimal rates are

where the unit of time is the expectation of the exponentially distributed time for an individual to die and be replaced.

The many other examples of interest include the SIS and SEIR models for epidemics (e.g. Andersson and Britton, 2000) , dimerisation and the Michaelis-Menten reaction kinetics (e.g. Wilkinson, 2012) . Further examples but with an infinite statespace include the Schlögel model (e.g. Vellela and Qian, 2009) , the Lotka-Volterra predator-prey model (e.g. Wilkinson, 2012; Drovandi and McCutchan, 2016) and models for the autoregulation of the production of a protein (e.g. Wilkinson, 2012) , all of which are tackled using matrix exponentials in Sherlock and Golightly (2019) .

Denote the statespace of the Markov chain {X t } t≥0 by X = {x (k) } d k=1 . Let the prior mass function across states be ν(x|θ), the infinitesimal generator be Q(θ), and suppose there are observations y 0 , y 1 , . . . , y n at times t 0 , t 1 , . . . , t n , where Y i |(X i = x i ) has a mass function of p(y i |x i , θ), i = 0, . . . , n.

For any continuous-time Markov chain {X t } t≥0 with an infinitesimal generator, or rate matrix of Q, the (x, x )th element of exp(Qt) gives the transition probability (e.g. Norris, 1997) :

x,x , where here and elsewhere we abuse notation by identifying the state x (i) ∈ X with the corresponding index i ∈ {1, . . . , d}.

Defining the diagonal likelihood matrix to be L j (θ) = diag(p(y j |x (1) , θ), . . . , p(y j |x (d) , θ)) and ∆ j = t j − t j−1 , j = 1, . . . , n, the likelihood for the observations is then

where 1 is the d-vector of ones. Similarly, the filtering distribution after observation y m is

Consider the required multiplication from left to right: since the likelihood vectors L j (θ) are diagonal, pre-multiplication by a d-vector is an O(d) operation. Pre-multiplcation of the exponential of a sparse matrix by a d-vector via the uniformisation algorithm is also O(d) (see Section 3.1), so the entire likelihood calculation is O(d). In the case of certain epidemic models d itself can be much smaller than might naively be assumed. Ho et al. (2018) points out that the evolution of the Markov chain in an SIR model from one observation time to the next can be described entirely by the number of new infections and the number of new removals, B I and B R , neither of which can be negative. Consider the case of exact observations and suppose, for example, that in a population of size n pop = 500, x a = (S a , I a , R a ) = (485, 2, 13) and for some t > 0, x a+t = (470, 3, 27). Then b R = R a+t − R a = 14 and b I = S a − S a+t = 15. The size of the statespace for evolution between time a and time a + t, X a+t a , is then reduced from the size of the full statespace, (n pop + 1)(n pop + 2)/2 = 125751 to (b I + 1)(b R + 1) = 240.

In Ho et al. (2018) , a recursive formula for the Laplace transform of the transition probability to a given new state in terms of transition probabilities for old states then permits estimation of the transition vector from a known initial starting point in O(d) operations, where d is the dimension of the statespace actually required.

We may use the same statespace formulation as Ho et al. (2018) , provided we include an additional coffin state, C, with Q C,x = 0 for all x ∈ X a+t a ∪ C. Any births that would leave the statespace (and hence contradict the observation at time a + t) instead go to C. We also provide a further reduction in the size of the statespace, by a factor of up to one half. The current number of infections can never be negative, so, throughout the time interval [a, a + t], b R ≤ I a + b I . In the example above, this reduces the statespace size still further, from 240 to 162.

Although we do not examine it here, a similar reduction of the statespace but for the SEIR model, to (b E , b I , b R ), is described in Ho et al. (2018) . A further reduction, by a factor of up to 6, is possible by observing that neither

In epidemics, typically only removals are observed. Thus, to take advantage of the reduced statespace formulation, latent variables representing the number infected at each observation time must be introduced. This makes Bayesian inference via MCMC feasible for relatively large populations. For example, Ho et al. (2018) performs inference for the SIR model using data from the Ebola outbreak in regions of Guinea, and in Section 4.3 we perform inference on data from the 2013 Measles outbreak in Swansea, Wales.

The exponential of a d×d square matrix, M is defined via its infinite series:

As might be anticipated from the definition, for a d × d matrix, algorithms for evaluating exp (M) take O(d 3 ) operations (see Moler and Van Loan, 2003 , for a review of many such methods). However, for a d-vector, v, the product exp (M t) v is the solution to the initial value problem w(0) = v, dw/dt = Mw, and is the key component of the solution to more complex differential equations such as dw/dt = Mw + Bu(t). For this reason the numerical evaluation of the action of a matrix exponential on a vector has received considerable attention of itself (e.g. Gallopoulos and Saad, 1992; Saad, 1992; Sidje, 1998; Al-Mohy and Higham, 2011) . With double-precision arithmetic, real numbers are stored to an accuracy of approximately 10 −16 . Thus, evaluation of the exponential of a large negative number via its Taylor series is prone to potentially enormous round-off errors due to the almost cancellation of successive large positive and negative terms; a similar problem can affect the exponentiation of a matrix. Such issues are typically circumvented via the identity

applied for a sufficiently large integer K, and evaluated via K successive evaluations of product of exp(M/k) and a vector. The calculation on the right of (4) typically involves many more numerical operations than the direct calculation on the right of (3), so K should be the smallest integer that leads to the required precision by mitigating sufficiently against the cancellation of large positive and negative terms. This minimises both the accumulation of rounding errors and the total compute time given the required accuracy.

One common technique for such multiplication, exemplified in the popular Expokit FORTRAN routines (Sidje, 1998) , estimates e M/K v via its projection on to the Krylov subspace of Span{v, Mv, . . . , M n−1 v}, where n << d. A second method is provided in Al-Mohy and Higham (2011), where the key contributions lie in the method for choosing K and for choosing a suitable truncation point for the infinite series, as well as a means of truncating each series early depending on the behaviour of recent terms.

Both of the above algorithms use the fact that M is sparse and that only the action of exp(M) on a vector is required, but neither uses the special structure of the problem of interest to us: we require ν exp(Qt) where Q is a rate matrix and ν is a non-negative vector. Since Qt is also a rate matrix, we henceforth set t = 1 without loss of generality.

P is a Markov transition matrix, and the key observation is that

Firstly, P has no negative entries so cancellation of terms with alternating signs is no longer a concern. Secondly, exp Q can be interpreted as a mixture over a Poisson(ρ) random variable I, of I transitions of the discrete-time Markov chain with a transition matrix of P.

The next two subsections detail variations on two existing algorithms that utilise this special structure: the uniformisation algorithm and a variation on the scaling and squaring algorithm. For sparse rate matrices, the uniformisation algorithm has a cost of O(ρd), whereas the scaling and squaring algorithm has a cost of O(d 3 log ρ). Thus, the uniformisation algorithm is preferred when ρ is small, and scaling and squaring when ρ is large but d is relatively small. We now describe the two algorithms in detail.

In many statistical applications, the most appropriate algorithm for calculating µ := ν exp Q is the uniformisation algorithm (e.g. Reibman and Trivedi, 1988; Sidje and Stewart, 1999) . This estimates µ by truncating a single series none of whose terms can be negative, rather than truncating multiple series where terms may change sign as in Al-Mohy and Higham (2011). Given an > 0, the algorithm calculates an approximation, µ, to µ by picking a truncation point for the infinite series, such that, if ν were a probability vector, the (guaranteed to be non-negative) amount of true missing probability over all of the d dimensions is controlled:

where µ * is the probability vector that would be calculated if there were no rounding errors, and the only errors were due to the truncation of the infinite series. Typically we aim for to be similar to the machine's precision. We control the absolute truncation error and note that with any truncation of the power series, it is impossible to obtain general control of the relative error in a given component of µ, | µ i /µ i − 1|. Consider, for example, a Moran process (Example 2), where Q is tridiagonal. Then Q k is also banded, with a band width of 2k + 1. For any given m max , and ν = (1, 0, 0, . . . ), set d > m max + 1. The truncated approximation to e Q gives a transition probability of 0 for all states above m max + 1, yet, in truth there is a non-zero probability of such a transition.

From (6),

So the absolute relative error, or (when ν is a probability vector) missing probability mass, due to truncation is

the smallest m required to achieve an error of at most , or, essentially, the quantile function for a Poisson(ρ) random variable, evaluated at 1 − . Chebyshev's inequality applied to

computational cost given earlier in this section.

In many programming languages, standard functions are available to evaluate m (ρ). However, for example, in R we find i.e., an inability to calculate m (ρ) correctly given the small values that we require; the underlying functions are also callable from C++ and lead to the same error. In Appendix A we provide sharp bounds on m (ρ), and this leads to an accurate methodology for its exact calculation, producing the same (correct) answers as the C++ boost library (which we have not been able to use with RCpp) and up to twice as quickly.

The uniformisation algorithm is presented as Algorithm 3.1. For large values of ρ, although there is no problem with large positive and negative terms cancelling, it is possible that the partial sum k i=0 ρ i i! might exceed the largest floating point number storable on the machine. We circumvent this problem by occasionally renormalising the vector partial sum when the most recent contribution is large, and compensating for this at the end; see lines 5, 12 and 14.

One of the simplest, yet most robust methods for exponentiating any square matrix is the scaling and squaring algorithm (e.g. Moler and Van Loan, 2003) . When the square matrix is an infinitesimal generator, this method can be made even more robust using the reformulation in (6). Furthermore, when not exp Q but ν exp Q is required, some further computational savings can be obtained.

The basic scaling and squaring algorithm takes advantage of the identity

where for any integer s, a square matrix is raised to the power of 2 s by squaring it s times. We set M = Q + ρI = ρP from (5).

Algorithm 1 Uniformisation algorithm for ν e Q with a missing mass of at most .

is approximated via the uniformisation algorithm applied to a matrix (e.g. Ross, 1996) :

This quantity is then squared s times. A linear search between tight upper and lower bounds gives the s which minimises the total number of matrix multiplications, m (ρ/2 s ) + s; we set s ← s − min(2, s) as a compromise for the fact that when Q is sparse, the individual matrix multiplications used to evaluate exp(M small ) are cheaper than those used to square it, since exp(M small ) is dense.

When evaluating ν exp(Q) = exp(−ρ)ν exp(M) via scaling and squaring with s > 0 it is never most efficient to first evaluate exp(M). Let s 1 and s 2 be two integers such that s 1 + s 2 = s. Then

with 2 s 2 matrix vector products. The cost of s 1 matrix squares and 2 s 2 vector-matrix products (where the matrix is dense) is s 1 d 3 + 2 s 2 d 2 . We round the minimiser down to the nearest integer for simplicity, setting

Even with d = 2 this gives s 2 = min(s, 1).

We now describe two optional extensions: renormalisation, which improves the accuracy of any matrix exponentiation algorithm used on a rate matrix, and two-tailed truncation, which is unique to the uniformisation algorithm and allows a small computational saving.

Since a := d i=1 µ i = d i=1 ν i there is no need to keep track of the logarithmic offset (c in Algorithm 3.1). Instead the final vector (v sum in Algorithm 3.1) is renormalised at the end so that its components sum to a.

Algorithm 2 Scaling and squaring algorithm for ν e Q with a missing mass of at most .

Find s via linear search; ρ small ← ρ/2 s ; find m (ρ small ); find (s 1 , s 2 ) via (7). 3: M small ← (Q + ρI)/2 s . 4: ν pro ← ν. 5: A pro ← M small ; A sum ← I + M small 6: f ← 2. 7: for j from 2 to m do 8:

A sum ← A sum + A pro .

10:

f ← f + 1.

11: A sum ← e −ρ small A sum 12: for j from 1 to s 1 do

A sum ← A sum × A sum . 14: for j from 1 to 2 s 2 do 15:

Two-tailed truncation (e.g. Reibman and Trivedi, 1988 ) permits a small reduction in the computational cost of the uniformisation algorithm with no loss of accuracy. When ρ is moderate or large, the total mass of probability from the initial value of v sum and the early values accumulated into v sum (Steps 6 and 10 of Algorithm 3.1) is negligible (has a relative value smaller than /2, say) compared with the sum of the later values. In such cases v sum may be initialised to 0 and step 10 omitted for values of j beneath some m lo . Proposition 1 below shows that if m is chosen such that P (Po(ρ) > m) ≤ /2 then setting m lo := max(0, 2 ρ − 0.5 − m) ensures that the missing probability mass is no more than . For large ρ, m − m lo = O( √ ρ), so with two-tailed truncation the cumulative cost of Step 10 dwindles compared with the other O(d) costs, which are repeated O(ρ) times.

Proposition 1. Given ρ > 0, let p n = e −ρ ρ n /n! = P (Poisson(ρ) = n), and let c = ρ−1/2 . Then for a ≤ c − 1,

Proof. For any integer b, and 1 ≤ i ≤ b,

Our C++ implementation uses the recent basic sparse matrix functionality in the C++ Armadillo library Curtin, 2016, 2018) to calculate ν exp Q, where ν is non-negative and Q is a large, sparse rate matrix. Direct function calls from the R programming language are enabled through RcppArmadillo (Eddelbuettel and Sanderson, 2014) . For completeness, the functions can also be called with dense rate matrices. The functions are collected into the expQ package which is downloadable from https: //github.com/ChrisGSherlock/expQ and are briefly outlined in Appendix B.

The speed of a vector multiplication by a sparse-matrix depends on the implementation of the sparse matrix algorithm. In R (R Core Team, 2018) and in C++ Armadillo, sparse matrices are stored in column-major order. Hence pre-multiplication of the sparse matrix by a vector, ν Q, is much quicker than post multiplcation, Qν. In other languages, such as Matlab, sparse matrices are stored in row-major order and post-multiplication is the quicker operation, so Q should be stored and used, rather than Q.

In Al-Mohy and Higham (2011) their new algorithm (henceforth referred to as AMH) is compared across many examples against state-of-the-art competitors, including, in particular, the expokit function expv (Sidje, 1998) . In most of the experiments AMH is found to give comparable or superior accuracy together with superior computational speed. Given these existing comparisons and that the superiority of the uniformisation algorithm over the algorithm of Al-Mohy and Higham (2011) (for rate matrices) is not the main thrust of this paper, we perform a short comparison of accuracy and speed for two different likelihood calculations for an SIR model fitted to data from the Eyam plague. We compare our implementation of the uniformisation algorithm, the algorithm of AMH, the expAtv function which is from the R package expm and uses the method of Sidje (1998) , and the bespoke algorithm for epidemic processes in Ho et al. (2018) . Since it would be unfair to compare the clock-speeds for the Matlab code for AMH directly with those of our RCppArmadillo implementation, we compare the number of sparse vector-matrix multiplications that are required.

The highest accuracy available in C++ using sparse matrices and the armadillo linear algebra library is double precision, which we used throughout in our implementation of both of our algorithms. For the uniformisation and scaling and squaring algorithms we used = 10 −15 , and for AMH we used the double-precision option. For expAtv and for Ho et al. (2018) we use the default package setting.

To examine the speed and accuracy of the algorithm we consider the collection (see the first three rows of Table 1 ) of (S, I) (susceptible and infected) values, which originated in Raggett (1982) and were used in Ho et al. (2018), for the Eyam plague. We set the parameters to their maximum-likelihood estimates, (β, γ) = (0.0196, 3.204) and consider the likelihood for the data in Table 1 . In addition, to mimic the size of potential changes between observation times and the size of the elements of the rate matrix from a larger population, we also evaluated the likelihood for the jump directly from the data at time 0 to the data at time 4. The final three rows of Table 1 refer to the rate matrix for the transition between consecutive observations and provide the dimension the matrix first using the reformulation of Ho et al. (2018) and then applying the improvement described in Section 2.2; the final row is the absolute value of the largest entry of Q, ρ. The rate matrix for the single jump between times 0 and 4 had d HCS = 30789, d = 16082 and ρ ≈ 3439.5. The full statespace has a size of 34453. Thus, for large changes, the main reduction in size arises from the improvement in Section 2.2, but for small jumps this provides a smaller relative reduction compared with that in Ho et al. (2018) .

For the uniformisation and scaling and squaring algorithm, with = 10 −15 , the algorithm of Ho et al. (2018) and the expAtv function from the R package expm package (Sidje, 1998) we found the CPU time for 1000 estimations of the likelihood (20 estimates for the likelihood for the jump from t = 0 to t = 4). We also recorded the error in the evaluation of the log likelihood. Since for uniformisation, using renormalisation and twotailed truncation together produced the fastest and most accurate evaluations, we only considered this combination. Given that the true likelihood is not known, the error using uniformisation, from scaling and squaring and from Al-Mohy and Higham (2011) were approximately bounded by examining their discrepancy from each other. The results are presented in Table 2 .

Scaling and squaring is extremely slow in these high-dimensional scenarios; however, Sherlock and Golightly (2019) provides a bistable example, the Schlögel model, where d ≈ 100 − 200 but ρ > 10 5 , and the scaling and squaring algorithm outperforms uniformisation by orders of magnitude.

Since m = O(ρ) the choice of tolerance, , typically has only a small effect on the speed of the uniformisation algorithm. For the full likelihood evaluation, uniformisation is over twice as fast as the algorithm of Ho et al. (2018) and approximately thirty times as fast as expAtv, and is more accurate than either; it is also over twice as fast as the algorithm of Al-Mohy and Higham (2011), although both are very accurate.

For the single large jump between observations, we see the same pattern in terms of ac- curacy. There is a gain in efficiency by using two-tailed-truncation because ρ is larger (m lo = 3081 and m = 3797), but despite this, the method of Ho et al. (2018) is now more efficient than uniformisation, although considerably less accurate than it. Again, expAtv is over twenty times slower than uniformisation and less accurate, and AMH is over three times slower than uniformisation.

We now consider the Moran model, which has four unknown parameters: (α, β, u, v) and n pop = 1000. Setting (α, β, u, v) = (1, 0.3, 0.2, 0.1), we simulate a path of the process for T = 10000 time units. We then sample 51 observations at times 0, 200, 400, . . . , 10000, by taking the value of the process at each of these times and adding independent noise with a distribution of Bin(800, 0.5) − 400.

We then perform inference on θ = (log α, log β, log[u/(1 − u)], log[v/(1 − v)]) by maximising the likelihood based on all the data and, separately, based on the data up to T = 5000. In each of these two data scenarios we find the filtering distribution, P X T |y 0:T , θ , at time T via (2); finally we predict forward from T in steps of 200 for a further time of T pred = 5000 by repeatedly multiplying the current distribution vector by exp (200Q( θ) ). The true values, observations and filtering and prediction distributions are shown in Figure 1 . The whole process of inference and prediction took less than two minutes on a single i7-3770 CPU running at 3.40GHz. Further, after defining Q, only 10 lines of R code are required to calculate the log-likelihood, and fewer than this to produce the filtering distribution (see Appendix D). x 

The largest measles outbreak in the United Kingdom between 2011 and 2019 centred around Swansea, Wales and occurred between November 2012 and July 2013. Of the 1219 cases in mid-and west-Wales, 664 occurred in the Swansea Local Authority (LA) area, 243 in the nearby Neath and Port Talbot LA and fewer than 80 occurred in any of the other individual LA areas in West or South Wales (http://www.wales.nhs.uk/sitesplus/888/page/ 66389, accessed February 10th 2020). A reduction in uptake of the MMR (Measles, Mumps, Rubella) vaccine has been blamed (e.g. Jakab and Salisbury, 2013) for this, with particularly low rates reported in Swansea (https://en.wikipedia.org/wiki/2013_Swansea_ measles_epidemic, accessed February 10th 2020).

The basic reproduction number, R 0 , is the expected number of secondary infections in a susceptible population that arise directly from the primary infection of a single individual. For the SIR model described in Section 1.1, R 0 = β/γ. For measles, R 0 is often reported as between 14 and 18 (e.g. Anderson and May, 1982) , which fits with the World Health Organisation (WHO) recommendation of vaccination level of at least 93 − 95% (WHO, 2009).

We fit the SIR model to the notification data for the Swansea LA provided in Table 3 so as to estimate the overall R 0 for the partially vaccinated population in Swansea and to demonstrate inference on the unknown number of infectious individuals at each observation time. In fitting the model we are making several assumptions and simplifications, including the following. Firstly, we are ignoring infections from Swansea to other LAs and from these LAs to Swansea; since most of the infections occurred in Swansea the former will outnumber the latter and so we will underestimate the 'true' R 0 , and provide a 'local' R 0 at the epicentre of the infection. Secondly it is known that the lowest level of vaccination, and the highest level of infection was amongst 10-18 year olds (Wise, 2013) ; a more accurate model would, therefore, partition the population into age groups. Age-stratified, continuous-time Markov chain SIR models are difficult to fit in general, however, and often a deterministic version of the model is used (e.g. Broadfoot and Keeling, 2015) . Finally, we treat a notification as equivalent to a removal: this is not unreasonable as once an individual has been diagnised by a GP with suspected measles they will be asked to isolate themselves.

As described in Section 2.2 we add as latent variables the number of infections at each of the reporting times, Days 30, 61, 92, ..., 212. The number of infections at times 242 and 273 must both be zero.

To understand the evolution near the peak of the epidemic and speed up inference still further, we add latent observation times during the peak of the infection, at Days 125, 130, 135, 140, 145, 155, 160, 165, 170 and 175 . This leads to 10 further latent observations of the number of infected individuals and (because of constraints) 10 further latent observations of the number removed during each reduced time period, leading to a total of 27 integer latent variables.

We use a N(log 5, 2/3) prior for log R 0 = log(β/γ), a N(log(1/15), 1) prior for log γ and, because it is very poorly identified, we set the prior for the effective population size to p(N pop = n) ∝ exp(−n/500)1 {n≥1000} .

We perform inference via a Metropolis-within-Gibbs algorithm: θ = (log β, log γ) is updated via a random walk proposal with a jump of N(0, λ 2 I 2 ), n pop via an integer-valued random walk proposal, and x latent , the latent observations via integer random walks, with physical constraints (such as the sum of all the Rs not being able to exceed n pop ) checked for automatically; see Appendix C for more details.

The basic reproduction number, R 0 , is estimated as 1.15, with a 95% credible interval of (1.01, 1.31). This fits with other information known: firstly,up until 2013, R 0 only changed gradually over time (due to year-on-year variations in infant vaccination rates) and it cannot have reached much higher than 1 in late 2012 as otherwise there would have been an outbreak in a previous year; secondly an R 0 of 1.15 if the true R 0 is 16, corresponds to a vaccination level of 93%, and R 0 = 1.3 corresponds to a 92% level, and as argued earlier, we expect to slightly underestimate R 0 . As of December 2012, the estimated coverage of one dose of MMR vaccine among 16 year-olds in Wales was 91% (Public Health Wales, 2013).

The right-hand panels of Figure 1 show the posterior median and 95% credible intervals for the number of infections at each of the monthly observation times and at the 10 additional latent times, and similar intervals for the cumulative number of infections. In any infectious disease, at any current time point, it is vital to understand the current, unknown, number of infections in order to be able to predict the future course of an epidemic.

We have shown that inference, prediction and filtering for continuous-time Markov chains with a large but finite statespace, especially those arising from reaction networks is not just conceptually straightforward when the matrix exponential is used, but it is also often practical. We have provided and demonstrated the use of robust tools for this purpose in R, which opens up the direct use of and inference for reaction-network models to a wider audience. Straightforward inference for epidemic models, such as the SIR and SEIR models is particularly apposite at the time of submission, as it might have enabled an analysis of early COVID-19 infection data by people not expert in the more complex MCMC methodology typically used.

We emphasise that we are not suggesting that the tools we provide should replace the particle MCMC, ABC and SMC methods currently employed. In our experience, inference for epidemic models coded in a fast, compiled language is often more efficient in terms of effective samples per second, for example, than the approach using matrix exponentiation. However, the matrix exponential approach is much more straightforward, and the code that uses it can be written in the simpler, interpreted language R.

As the size of the statespace increases, the efficiency of the matrix exponentiation approach decreases; however, once the statespace becomes sufficiently large, the evolution of the process is often approximated by a stochastic differential equation (e.g. Golightly and Wilkinson, 2005; Fearnhead et al., 2014) or, when the behaviour is effectively deterministic, by ordinary differential equations (e.g. Broadfoot and Keeling, 2015) .

For the scaling and squaring approach, in particular, the cost of the exponentiation of Q/K can be nearly halved by using a Padé approximant (e.g. Moler and Van Loan, 2003) , but this then requires a matrix inversion, and so, for reasons of robustness, was not pursued here. where both inequalities require < 0.04 and the latter also requires < 1−e −ρ and B > log , where

A := 2ρh m + + 1 ρ and B := − 1 2 log 4πρh m − ρ ,

and h(x) = x − 1 + x log x.

The bound (8) arises from a standard argument, whereas those in (9) and (10) 

which follows from the equivalence between at least m + 1 events of a Poisson process with a unit rate occuring by time ρ and the time until the m + 1th event being at most ρ, permit a simple but fast binary search for m (ρ).

Our binary search algorithm homes in on the required m using the upper and lower bounds of Theorem 1 together with the identity (11), the right hand side of which can be evaluated quickly and accurately using the standard C++ toolbox, boost. This is quicker than the standard implementation of the Poisson quantile function (e.g. as implemented in boost), which uses the Cornish-Fisher expansion to approximate the quantile (hence needing an expensive evaluation of Φ −1 ) and then conducts a local search.

The simple bounds for small ρ arise because e −ρ > 1 − ρ. Hence r 0 (ρ) = 1 − e −ρ < ρ and if ρ ≤ , m (ρ) = 0. Furthermore, r 1 (ρ) = 1 − e −ρ (1 + ρ) < ρ 2 , so if ρ ≤ √ then r 1 (ρ) ≤ , so m (ρ) ≤ 1.

The other bounds all use aspects of the following result.

Proof. The left hand inequality holds for x ≥ 0 and is from Boucheron et al. (2013) page 36. For the right hand inequality, set g(x) = (x − 1) 2 /2 and notice that 0 = h(1) = g(1) = h (1) = g (1), and h (x) = 1/x ≤ 1 = g (x) for x ≥ 1.

The first upper bound on m (ρ), (8), arises from a standard Chernoff argument (e.g. Boucheron et al., 2013) to the right tail of a Poisson(ρ) random variable, X. The moment generating function is M X (t) = E e Xt = exp[ρ(e t − 1)], and by Markov's inequality:

P (X ≥ m) = P e Xt ≥ e mt ≤ e −mt M X (t) = e −mt+ρ(e t −1) .

The inequality holds for all t and the right-hand side is minimised at t = log(m/ρ), giving

by Lemma 1. Setting = P (X ≥ m + 1) and y = (m +1)/ρ − 1 gives 3ρy 2 (6+ 2y) log ≥ 0, from which y ≥ − log × 1 − 18ρ/ log /(3ρ), and (8) follows on substituting for y.

The much tighter bounds in (9) and (10) use Theorem 2 of Short (2013) , which can be rewritten to state that

where m := m + 1 and where the left hand side holds provided m > ρ and the right hand side holds provided m > ρ. We first show that these conditions are satisfied. Firstly, when ρ < 1, clearly m > ρ, moreover r 0 (ρ) = 1 − e −ρ , so provided 1 − e −ρ > , we require m ≥ 1 > ρ. When ρ ≥ 1, we use the easily verified facts that r m (m) is an increasing function of m and r m (ρ) is an increasing function of ρ; thus for ρ ≥ m ≥ 1, r m (ρ) ≥ r m (m) ≥ r 1 (1) = 1 − 2e −1 > 0.04, and the tolerance condition is not satisfied. We, therefore need m > ρ (which also gives m > ρ).

Neither Φ −1 nor h −1 is tractable (functions that perform Φ −1 (p) solve Φ(x) = p iteratively), and even with the bounds on h from Lemma 1 and standard bounds on Φ in terms of φ, tractable inversion is still not possible. We use the bound (8) to create (9), and then (9) to create (10).

To prove (9), since ≤ 0.04, from the left inequality in (12),

Firstly, since m + + 1 ≥ m + 1, this ensures A > 1, so log(A − 1) is real. More importantly, it ensures that [2ρh(m /ρ)] −1/2 − [2ρh(m /ρ)] −3/2 is a decreasing function of [2ρh(m /ρ)] 1/2 and, since h (x) > 0 for x > 1, it is also a decreasing function of m . The m that we desire satisfies m ≤ m + + 1 =: m + , and hence

Since, for y > 0, Φ(−y) > (1/y − 1/y 3 )φ(y),

Combining the left inequality in (12) with the right-hand inequality in Lemma 1 gives

Equation (9) follows on rearrangement.

To show (10) we apply the right hand inequality in (12) and the bound Φ(−x) < φ(x)/x, then the fact that m ≥ m − , and finally Lemma 1 to find:

P (X > m) < , where x = m/ρ. We must, therefore, ensure that the final bound is no more than . Rearranging this gives 3ρ(x − 1) 2 − 2(B − log )(x − 1) − 6(B − log ) ≤ 0, so that when B − log > 0, x − 1 ≤ (B − log )(1 + 1 + 18ρ/(B − ))/(3ρ).

The functions in the expQ package are provided below. Each function requires a rate matrix, Q, which can be sparse or dense, and a precision, .

Unif v exp Q takes a horizontal vector, v, and calculates v exp Q via uniformisation. SS v exp Q takes a horizontal vector, v, and calculates v exp Q via scaling and squaring. v exp Q takes a horizontal vector, v, and calculates v exp Q via whichever is likely (based on empirical results on an i7-3770 CPU) to be the more efficient of uniformisation or scaling and squaring. vT exp Q takes a vertical vector, v, and calculates v exp Q via whichever is likely (based on empirical results on an i7-3770 CPU) to be the more efficient of uniformisation or scaling and squaring. SS exp Q calculates exp Q using scaling and squaring.

Our particular reduced-statespace implementation of the SIR model fit for the Swansea Measles epidemic uses 10 additional latent observation times, 5 between days 120 and 151 (at days 125, 130, 135, 140 and 145) and five between days 151 and 212 (at days 155, 160, 165, 170 and 175) . This leads to 27 latent variables: 17 unknown number of infecteds at the (true and latent) observation times and 10 (not 12 because two sums are known) unknown numbers of recovered for the time period since the previous (true or latent) observation time. We emphasise that the R latent variables are not the cumulative number of recovered individuals since the epidemic began.

When a new latent vector is proposed, we first check whether it can possibly fit with the current n pop and the known data. If it does not fit, then the proposal may be rejected without any matrix exponentiation. At the jth (true or latent) time point, denote the current number of infecteds by I j and the number removed since the previous time point by R j . Let J be the total number of (true and latent) time points. Note that S j = n pop − I 0 − I j − j i=1 R i . The following checks are performed:

1. For each j = 1, . . . , J: I j ≥ 0, R j ≥ 0 and S j ≥ 0.

2. For each j = 1, . . . , J: S j ≤ S j−1 .

3. For each j = 1, . . . , J − 1:

These constraints can hinder the mixing of the integer-valued random walk algorithm on the latent variables, so we split the latent variables into four groups, grouped by observation time. This grouping has the additional advantage that only a subset of matrix exponentiation calculations need be performed for each of the four individual proposals.

To indicate the simplicity of inference via the matrix exponential, we provide code to evaluate the log-likelihood for the Moran model. Code for the filtering distribution is very similar but there is no need to track the re-normalisation constant (in ll). 

Computing the action of a matrix exponential with an application to exponential integrators

A continuous-time hidden Markov model for cancer surveillance using serum biomarkers with application to hepatocellular carcinoma

Directly transmitted infectious diseases: control by vaccination

Stochastic Epidemic Models and Their Statistical Analysis

Particle Markov chain Monte Carlo for efficient numerical simulation

Particle Markov chain Monte Carlo methods

Concentration inequalities [electronic resource] : a nonasymptotic theory of independence

Measles epidemics in vaccinated populations. Accessed on 1st

Alive SMC 2 : Bayesian model selction for low-count time series models with intractable likelihoods

Rcpparmadillo: Accelerating r with highperformance c++ linear algebra

Inference for reaction networks using the Linear Noise Approximation

Efficient solution of parabolic equations by Krylov approximation methods

Unbiased bayesian inference for population markov jump processes via random truncations

Efficient sampling of conditioned markov jump processes

Bayesian inference for stochastic kinetic models using a diffusion approximation

Bayesian parameter inference for stochastic biochemical network models using particle Markov chain Monte Carlo

Bayesian inference for Markov jump processes with informative observations

Direct likelihood-based inference for discretely observed stochastic compartmental models of infectious disease

Back to basics: the miracle and tragedy of measles vaccine

A population Monte Carlo scheme with transformed weights and its application to stochastic kinetic models

A tutorial introduction to bayesian inference for stochastic epidemic models using approximate bayesian computation

Simulation-based Bayesian inference for epidemic models

Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later

Random processes in genetics

Markov chains

Likelihood free inference for Markov processes: a comparison

Vaccine uptake in children in wales

R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing

A stochastic model of the Eyam plague

Fast mcmc sampling for markov jump processes and extensions

Numerical transient analysis of markov models

Stochastic Processes

Analysis of some Krylov subspace approximations to the matrix exponential operator

Armadillo: a template-based C++ library for linear algebra

A user-friendly hybrid sparse matrix class in C++

Exact bayesian inference for discretely observed markov jump processes using finite rate matrices

Improved inequalities for the Poisson and binomial distribution and upper tail quantile functions

EXPOKIT: a software package for computing matrix exponentials

A numerical study of large sparse matrix exponentials arising in Markov chains

Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögel model revisited

Who position on measles vaccines

Stochastic Modelling for Systems Biology

Largest group of children affected by measles outbreak in wales is 10-18 year olds

I would like to thank Prof. Lam Ho for suggesting that the reformulation of the statespace in Ho et al. (2018) in terms of births might be applicable within the methodology presented herein. I am also grateful to Dr. Andrew Golightly for several useful discussions.

Our fast, robust and accurate method for evaluating m (ρ), as defined in Section 3.1 relies on the following new result.Theorem 1. If ρ ≤ , m (ρ) = 0, and if ρ ≤ 1/2 , 0 ≤ m (ρ) ≤ 1. More generally: m (ρ) ≤ m + , where