Direct statistical inference for finite Markov jump processes via the matrix exponential.

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, Q , since exp ( 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 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.

Introduction

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; see Golightly and Wilkinson (2015), Golightly and Sherlock (2019), Black (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, . Whilst such inference is conceptually straightforward, it has often been avoided in practice for general MJPs, except in cases where the number of states is very small e.g. Amoros et al. (2019). The computational cost of each iteration of particle MCMC is proportional to the number of particles used and, for efficient estimation; see Doucet et al. (2015),Sherlock et al. (2015) this is approximately linear in the size of the statespace, d. In contrast, Matrix exponentiation has a computational cost of , 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, and it has been used successfully in a subclass of these situations ( e.g. Jenkinson and Goutsias (2012), see Sect. 2.2). There are three main reasons why this is possible: The sparsity of 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.

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

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.

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

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 Sect. 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 Sect. 4 benchmarks some of the algorithms and demonstrates their use for inference, filtering and prediction. The article concludes in Sect. 5 with a discussion.

Examples and motivation

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 . The off-diagonal elements of the rate matrix, , are all non-negative, and the ith diagonal element is . A chain that is currently in state i leaves this state upon the first event of a Poisson process with a rate of ; the state to which it transitions is j with a probability of . Whilst the rate matrix, , is a natural description of the process, the likelihood for typical observation regimes involves the transition matrix, , the (i, j)th element of which is exactly .

Both examples take the form of a reaction network, where from the current state , the next state change will occur according to one of the specified reactions. The state can be thought of as an integer vector, where each element of the vector indicates the numbers of a particular species that are currently present in the system. When, as here, the maximum number of each species is finite, the set of possible states can be placed in one-to-one correspondence with the natural numbers as required to define . Each reaction occurs according to a Poisson process with a rate, , and when it occurs species combine according to the reaction formula. For example, the first reaction in the SIR model, below occurs with a rate of , and when it occurs the state (S, I) changes to .

Example 1

The SIR model for epidemics. The SIR model for a disease epidemic has 3 species: those who are susceptible to the epidemic, , those both infected and infectious, , and those who have recovered from the epidemic and play no further part in the dynamics, . The non-negative counts of each species are denoted by S, I, and R. For relatively short epidemics the population, , is assumed to be fixed, and so the state of the Markov chain, represented by (S, I), is subject to the additional constraint of , with . 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, and in a population with a fixed size of . Individuals with allele reproduce at a rate of , and those with 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 switching to with a probability of u and switching to with a probability of v. Let and represent individuals with alleles and respectively and let N be the number of individuals with allele . The two reactions areSetting , the corresponding infinitesimal rates arewhere 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).

Data and likelihood calculations

Denote the statespace of the Markov chain by . Let the prior mass function across states be and define . Let the infinitesimal generator be , and suppose there are observations at times , where has a mass function of , .

Likelihood for noisy and partially observed data

For any continuous-time Markov chain with an infinitesimal generator, or rate matrix of , the element of gives the transition probability (e.g. Norris (1997)):where here and elsewhere we abuse notation by identifying the state with the corresponding index .

Defining the diagonal likelihood matrix to be and , , the likelihood for the observations is thenwhere is the d-vector of ones. Similarly, the filtering distribution after observation isConsider the required multiplication from left to right: since the likelihood vectors are diagonal, pre-multiplication by a d-vector is an operation. Pre-multiplication of the exponential of a sparse matrix by a d-vector via the uniformisation algorithm is also (see Sect. 3.1), so the entire likelihood calculation is . In the case of certain epidemic models d itself can be much smaller than might naively be assumed.

Statespace reduction for epidemic models

An alternative formulation of the statespace of the Markov chain for an SIR epidemic model (or more general models such as the SEIR), in terms of the degree of advancement (DA), was first pointed out in Jenkinson and Goutsias (2012). Instead of representing the state in terms of the number of susceptibles and the number of infecteds, given a known initial condition it is represented by the number of new infections and the number of new removals, and , neither of which can be negative and both of which are non-decreasing. Given the initial condition, the map from (S, I) to is one-to-one; however the rate matrix with the DA formulation is lower triangular, a key ingredient in the implicit Euler integration scheme used in Jenkinson and Goutsias (2012) to integrate the master equation.

When performing Bayesian inference for the SIR model using noisy, partial observations, Ho et al. (2018) points out that augmenting the state space of the MCMC Markov chain to include not just the model parameters but also the true values of S and I at each observation time can massively reduce the sizes of the statespaces that need to be considered when evolving the SIR process from one observation time to the next provided the DA formulation is used. Consider the case of exact observations and suppose, for example, that in a population of size , and for some , . Then and . The size of the statespace for evolution between time a and time , , is then reduced from the size of the full statespace, to . The exponential of the rate matrix is not used in Ho et al. (2018); instead, 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 operations, where d is the dimension of the statespace actually required. Inference is then performed for the SIR model using data from the Ebola outbreak in regions of Guinea.

We may use the DA formulation with data augmentation, provided we include an additional coffin state, , with for all . Any births that would leave the statespace (and hence contradict the observation at time ) instead go to . The aforementioned implementation, a square grid of possible states, includes “impossible” states to which the rate of entry is zero: the current number of infections can never be negative, so, throughout the time interval , . Removing these states altogether allows us to make a further reduction in the size of the statespace, by a factor of up to one half. In the example above, this reduces the statespace size still further, from 240 to 162.

Matrix exponentiation

The exponential of a square matrix, is defined via its infinite series: . As might be anticipated from the definition, for a matrix, algorithms for evaluating take operations (see Moler and Van Loan (2003), for a review of many such methods). However, for a d-vector, v, the product is the solution to the initial value problem , , and is the key component of the solution to more complex differential equations such as . 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).

When is dense,can be evaluated in operations if the series is truncated at an appropriate point. However, motivated by the examples in Sect. 1.1 our interest lies in large sparse matrices, and the number of operations can then be reduced to , where r is the average number of entries in each row of .

With double-precision arithmetic, real numbers are stored to an accuracy of approximately . 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 identityapplied for a sufficiently large integer K, and evaluated via K successive evaluations of product of 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 of Sidje (1998), estimates via its projection on to the Krylov subspace of , where . 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. These and other algorithms are compared, specifically for the case of the SIR model (which has a special structure; see Sect. 2.2) in Kinyanjui et al. (2018).

Both Krylov methods and that of Al-Mohy and Higham (2011) use the fact that is sparse and that only the action of on a vector is required, but neither uses the structure of interest to us: we require where is a rate matrix for a general MJP and is a non-negative vector. Since is also a rate matrix, we henceforth set without loss of generality. Let is a Markov transition matrix, and the key observation is thatFirstly, has no negative entries so cancellation of terms with alternating signs is no longer a concern. Secondly, can be interpreted as a mixture over a random variable I, of I transitions of the discrete-time Markov chain with a transition matrix of .

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 , whereas the scaling and squaring algorithm has a cost of . 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.

The uniformisation algorithm

In many statistical applications, the most appropriate algorithm for calculating 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 , 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 , . Consider, for example, a Moran process (Example 2), where is tridiagonal. Then is also banded, with a band width of . For any given , and , set . The truncated approximation to gives a transition probability of 0 for all states above , yet, in truth there is a non-zero probability of such a transition. However, the combined probability of all transitions which have not been accounted for is, by design, at most .

From (6),Now,So the absolute relative error, or (when is a probability vector) missing probability mass, due to truncation isthe tail probability of a random variable. Of direct interest to us isthe smallest m required to achieve an error of at most , or, essentially, the quantile function for a random variable, evaluated at . Chebyshev’s inequality applied to , where gives , implying the computational cost given earlier in this section.

In many programming languages, standard functions are available to evaluate . However, for example, in R we findi.e., an inability to calculate 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 , 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 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.

Scaling and squaring

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 but is required, some further computational savings can be obtained.

The basic scaling and squaring algorithm takes advantage of the identitywhere for any integer s, a square matrix is raised to the power of by squaring it s times. We set from (5). And define . First, is approximated via the uniformisation algorithm applied to a matrix, e.g. Ross (1996): . This quantity is then squared s times. The optimal value of s is chosen according to an algorithm described in Appendix B.

When evaluating via scaling and squaring with it is never most efficient to first evaluate . Let and be two integers such that . Thenwith matrix vector products. The cost of matrix squares and vector-matrix products (where the matrix is dense) is . We round the minimiser down to the nearest integer for simplicity, settingEven with this gives .

Improvements

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 there is no need to keep track of the logarithmic offset (c in Algorithm 3.1). Instead the final vector ( in Algorithm 3.1) is renormalised at the end so that its components sum to a.

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 and the early values accumulated into (Steps 6 and 10 of Algorithm 3.1) is negligible (has a relative value smaller than , say) compared with the sum of the later values. In such cases may be initialised to 0 and step 10 omitted for values of j beneath some . Proposition 1 below shows that if m is chosen such that then setting ensures that the missing probability mass is no more than . For large , , so with two-tailed truncation the cumulative cost of Step 10 dwindles compared with the other costs, which are repeated times.

Proposition 1

Given , let , and let . Then for ,

Proof

For any integer b, and ,where . Hence, if , , and so

Implementation

Our C++ implementation uses the recent basic sparse matrix functionality in the C++ Armadillo library; see Sanderson and Curtin (2016), Sanderson and Curtin (2018) to calculate , where is non-negative and is a large, sparse rate matrix. Direct function calls from the R programming language are enabled through RcppArmadillo; see 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 C.

The speed of a vector multiplication by a sparse-matrix depends on the implementation of the sparse matrix algorithm. In RR 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, , is much quicker than post multiplication, . In other languages, such as Matlab, sparse matrices are stored in row-major order and post-multiplication is the quicker operation, so should be stored and used, rather than Q.

Numerical comparisons and demonstrations

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 of 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.

When performing maximum-likelihood estimation, each iteration of the optimisation algorithm tries a new parameter value, and when performing Bayesian inference, each iteration of the algorithm proposes a new parameter value. In each case, given the parameter value, the pertinent rate matrices are created and then the matrix exponentiation function is called in turn for each of the matrices as required by (1). If the generic exponentiation function is called then the decision on whether to use Algorithm 1 or Algorithm 2 is based upon the dimension, d and the maximum absolute value on the diagonal, . Whether Algorithm 1 or 2 is called directly or via the generic exponentiation function, the first task it performs is the evaluation of . The cost of this is neglibible compared with that of the exponentiation itself, so it is essentially immaterial that is evaluated twice when the generic exponentiation function is called.

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 , and for AMH we used the double-precision option. For expAtv and for Ho et al. (2018) we use the default package setting.

Comparison with other matrix exponentiation algorithms

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, 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 Sect. 2.2; the final row is the absolute value of the largest entry of , . The rate matrix for the single jump between times 0 and 4 had , and . 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).

Table 1

Time (in units of 31 days), and numbers of susceptibles and infecteds, originally from Raggett (1982). The final rows indicates, for each pair of consecutive observations, the size of the statespace for evaluating the transition probability and the value for the associated rate matrix

Time	0	0.5	1.0	1.5	2.0	2.5	3.0	4.0
S	254	235	201	153	121	110	97	83
I	7	14	22	29	20	8	8	0
d_HCS	-	261	946	2059	1387	289	197	346
d	-	245	867	1868	1308	282	181	240
ρ	-	101.5	171.4	217.1	170.1	83.1	53.6	106.3

For the uniformisation and scaling and squaring algorithm, with , the algorithm of Ho et al. (2018) and the expAtv function from the R package expm which uses the technique of Sidje (1998) we found the CPU time for 1000 estimations of the likelihood (20 estimates for the likelihood for the jump from to ). We also recorded the error in the evaluation of the log likelihood. Since for uniformisation, using renormalisation and two-tailed 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.

Table 2

Timings for estimating the full log-likelihood (1000 repeats) and the log-likelihood for the jump from the initial to the final observation (20 repeats) for the Eyam data set, number of sparse vector-matrix multiplications for one repeat, and the accuracies of the estimates. Results are given for the method of Ho et al. (2018) (HCS), the expAtv function in the expm package, which uses the Krylov subspace techniques of Sidje (1998), the method of Al-Mohy and Higham (2011) (AMH), the uniformisation algorithm (Unif) and the scaling and squaring algorithm (SS). The timing for SS on the jump likelihood was estimated from a single repeat

Algorithm	Full likelihood			Jump likelihood
	Time (secs)	Mult	Accuracy	Time (secs)	Mult	Accuracy
HCS	45.3	–	5.7×10^-8	9.7	–	4.3×10^-9
expAtv	558.5	–	1.6×10^-10	323.2	–	8.2×10^-11
AMH	–	3701	<1×10^-15	-	14300	<4×10^-14
Unif	18.72	1596	<1×10^-15	15.2	3921	<6×10^-14
SS	1678	–	1.1×10^-13	8940^1	-	<6×10^-14

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 but , and the scaling and squaring algorithm outperforms uniformisation by orders of magnitude.

Since 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 accuracy. There is a gain in efficiency by using two-tailed-truncation because is larger ( and ), 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.

Maximum likelihood inference, filtering and prediction

We now consider the Moran model, which has four unknown parameters: and . Setting , we simulate a path of the process for time units. We then sample 51 observations at times , by taking the value of the process at each of these times and adding independent noise with a distribution of .

We then perform inference on by maximising the likelihood based on all the data and, separately, based on the data up to . In each of these two data scenarios we find the filtering distribution, , at time T via (2); finally we predict forward from T in steps of 200 for a further time of by repeatedly multiplying the current distribution vector by . The true values, observations and filtering and prediction distributions are shown in Fig. 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 , only 10 lines of code are required to calculate the log-likelihood, and fewer than this to produce the filtering distribution (see Appendix E).

Fig. 1

Moran model (left): true values (), observations (), filtering/prediction mean (solid lines) and quantiles (dashed and dotted lines) for a further time of 5000 from data up to and data up to . Swansea measles SIR model (right): posterior median (solid line, with to show the positions) and posterior quantiles for the number of infected people at each real or latent observation time (top) and the cumulative number recovered by that time (bottom)

Bayesian inference for the Swansea measles epidemic of 2013

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, , 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, . For measles, 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 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 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’ , and provide a ‘local’ 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, see 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 diagnosed by a GP with suspected measles they will be asked to isolate themselves.

Table 3

Number of measles notifications in the Swansea Local Authority area by month (from http://www.wales.nhs.uk/sitesplus/888/page/66389, February 10th 2020)

Month	2012			2013
	Oct	Nov	Dec	Jan	Feb	Mar	Apr	May	Jun	Jul
Day number	0	30	61	92	120	151	181	212	242	273
Notifications	0	10	27	34	59	183	278	56	17	0

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 prior for , a prior for and, because it is very poorly identified, we set the prior for the effective population size to .

We perform inference via a Metropolis-within-Gibbs algorithm: is updated via a random walk proposal with a jump of , via an integer-valued random walk proposal, and , the latent observations via integer random walks, with physical constraints (such as the sum of all the Rs not being able to exceed ) checked for automatically; see Appendix D for more details.

The basic reproduction number, , is estimated as 1.15, with a credible interval of (1.01, 1.31). This fits with other information known: firstly,up until 2013, 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 of 1.15 if the true is 16, corresponds to a vaccination level of , and corresponds to a level, and as argued earlier, we expect to slightly underestimate . As of December 2012, the estimated coverage of one dose of MMR vaccine among 16 year-olds in Wales was ; see Public Health Wales (2013).

The right-hand panels of Fig. 1 show the posterior median and 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.

Discussion

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 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.