key: cord-0627532-sig2bqwk
authors: Walder, Adam; Hanks, Ephraim M.
title: A New Framework for Inference on Markov Population Models
date: 2021-01-02
journal: nan
DOI: nan
sha: 1507010a12b457950243f841564227d3873766ad
doc_id: 627532
cord_uid: sig2bqwk

In this work we construct a joint Gaussian likelihood for approximate inference on Markov population models. We demonstrate that Markov population models can be approximated by a system of linear stochastic differential equations with time-varying coefficients. We show that the system of stochastic differential equations converges to a set of ordinary differential equations. We derive our proposed joint Gaussian deterministic limiting approximation (JGDLA) model from the limiting system of ordinary differential equations. The results is a method for inference on Markov population models that relies solely on the solution to a system deterministic equations. We show that our method requires no stochastic infill and exhibits improved predictive power in comparison to the Euler-Maruyama scheme on simulated susceptible-infected-recovered data sets. We use the JGDLA to fit a stochastic susceptible-exposed-infected-recovered system to the Princess Diamond COVID-19 cruise ship data set.

Statistical inference methods for population process models are required to understand the dynamics of disease outbreaks such as the current COVID-19 global pandemic (Chen et al., 2020; He et al., 2020; Mwalili et al., 2020) . Markov population models are commonly used to model disease dynamics in small to moderate sized populations (Sun et al., 2015; Allen, 2017; Fricks and Hanks, 2018) but simulation and inference can be computationally burdensome. In this work, we construct an approximate joint Gaussian likelihood for Markov population models that helps facilitate inference and simulation. This approximation, which we call the joint Gaussian deterministic limiting approximation (JGDLA) makes inference on Markov population models easier as the joint likelihood of all timereferenced observations can be computed using only the solution to a system of ordinary differential equations (ODEs). This removes the need of for stochastic infill, as is commonly needed for inference on Markov population models.

A common framework for inference on Markov population models is to consider a Gaussian approximation of the Markov population model using the functional central limit theorem (FCLT) . This approximates the Markov population model as a system of linear stochastic differential equations (SDEs). Except in a few specific cases, systems of linear SDEs with time-varying coefficients do not have analytical solutions. In the absence of an analytic solution, numerical methods are required to solve the system of SDEs. The Euler-Maruyama scheme is the most commonly used numerical method for performing inference on SDE models (Sun et al., 2015; Allen, 2017; Eisenhauer and Hanks, 2020) . The

Euler-Maruyama scheme relies on a fixed length time-lag that requires stochastic infill estimates to perform inference or predict at unobserved locations leading to computational bottlenecks.

The main novel contribution of this work is the construction of the JGDLA model, which is built on the premise that the joint distribution of an Itó diffusion approximation to the Markov population model can be fully constructed from the solution to a system of ODEs. The result is an approximate inference method for Markov population models which does not require stochastic infill and offers inference familiar to ODE modeling. This is not the first presentation of such a model: Kurtz (1978 Kurtz ( , 1981 introduced the mathematical framework needed to justify converge of the approximation, and Baxendale and Greenwood (2011) used the approximation technique to investigate the behavior of stochastic population models through simulation. However, previous work only considers simulation and exploration of this approximation. In this work, we develop methods for statistical inference using the JGDLA.

The population models considered in this work have a direct application in the field of epidemiology (Allen and Allen, 2003) . We demonstrate that JGDLA outperforms a Euler-Maruyama scheme on simulated susceptible-infected-recovered (SIR) data sets of varying population sizes in terms of mean absolute prediction error at infill locations. The results of this simulation study show that JGDLA offers inference that only relies on solving an ODE system, does not require stochastic infill for predictions and inference, and provides an improved model fit in comparison to Euler-Maruyama, the most common approximation method as a framework for statistical inference.

The remainder of the manuscript is organized as follows. In Section 2 we introduce Markov population models as the sum of Poisson processes with stochastic rates. In Section 3 we use the FCLT to form a system of linear SDEs for approximating Markov population models. We discuss the Euler-Maruyama scheme for approximate inference on Markov population models in Section 4. In Section 5, we introduce the JGDLA in the statistical framework. We compare JGDLA and the Euler-Maruyama approximation on simulated SIR data sets in Section 6. In Section 7, we use the JGDLA to fit a stochastic SEIR model to the Princess Diamond Cruise COVID-19 data set. We conclude with a discussion in Section 8.

Deterministic population models are widely used in fields such as chemistry, ecology, and epidemiology to model large scale population dynamics such as disease outbreaks (Keeling and Rohani, 2011; Fricks and Hanks, 2018) . While deterministic models work well for very large populations, stochastic methods are needed to capture fine scale dynamics for populations with few individuals (Fricks and Hanks, 2018) . In this section we formulate the Markov population model in terms of stochastic reactions and rates. Our treatment follows the formulation of Kurtz (1978) and Fricks and Hanks (2018) .

Let X(t) = (X 1 (t), X 2 (t), ..., X d (t)) be a d-dimensional random vector on the nonnegative integers. X j (t) represents the number of individuals in a population of size N belonging to subpopulation j at time t (e.g. number of infected individuals in a population).

A population reaction occurs when one individual moves from one class to another. There are n possible reactions for any given model, and d subpopulations or classes. We let R i be a d-dimensional vector denoting the i th reaction. For example, if an individual can move from class 1 to 3, the reaction vector will contain −1 as the first element, 1 as the third element, and 0 for all other d − 2 elements.

Each individual reaction is assumed to occur at a stochastic rate which depends on the current state X(t) and unknown rate parameters θ. We denote the reaction rate corresponding to R i by λ i θ (X(t)). Let Y i (λ(·)) be an independent Poisson process with rate λ(·). We define the stochastic population model as the sum of Poisson processes in terms of reactions vectors and reaction rates given by

A classic example of a Markov population model is the stochastic SIR model with a closed population. The SIR model tracks the proportion of susceptible and infected individuals. We let X N (t) = (S N (t), I N (t)) , denote the scaled proportions of susceptible and infected individuals at time t. There are two possible reactions; a susceptible individual becomes infected R 1 = (−1, 1) , and an infected recovers R 2 = (0, −1) . Susceptible individuals become infected at rate of λ 1 θ (X N (t)) = βS N (t)I N (t), and infected individuals recover at a rate of λ 2 θ (X N (t)) = γI N (t). We note that θ = (β, γ), where β is the contact rate, and γ is the recovery rate. For fixed values of θ, simulation methods such as the Gillespie algorithm (Gillespie, 1977) and tau-leaping (Cao et al., 2006) can be used to generate stochastic realizations from (1). The proportions of infected and susceptible individuals from a simulated SIR data set with a population of size N = 100 are shown in 

Markov population model dynamics are controlled by the rate parameters θ which are often unknown and need to be estimated from data. The most common approaches for inference on θ are based on diffusion approximations to (1). We use the FCLT to construct a system of linear SDEs to help facilitate inference on (1). Our development follows that of Baxendale and Greenwood (2011) and Fricks and Hanks (2018) .

Let X N (t) = 1 N X(t) be the normalized population process. We scale (1) by N to obtain

We apply the FCLT for Poisson processes (see Appendix A.1) to each scaled Poisson process in (2) to obtain the Gaussian approximation

where B i (t) are independent Brownian motions with variance t. We define

and differentiate (3) to obtain

where dB(t) = (dB 1 (t), dB 2 (t), ..., dB d (t)) is a d-dimensional vector of differentiated independent Brownian motions B i (t). In some cases, the system of SDEs in (6) yields an analytic solution (Øksendal, 2003) . However, in most cases, numerical methods are needed to provide approximate solutions.

The remainder of this work focuses on approximate inference methods for Markov population models that rely on numerical solutions to the system of SDEs in (6). In the next section, we highlight the computational bottlenecks of the most commonly used numeric scheme for solving (6), Euler-Maruyama. In Section 5, we introduce the JGDLA, which is an approximate joint Gaussian likelihood for (6) constructed from the solution to a system of ODEs. We then illustrate how to perform inference on (6) using the JGDLA.

In Section 3 we showed that Markov population models can be approximated by the system of SDEs in (6). Performing inference on θ in the absence of an analytic solution consists of two steps; first (6) is approximated by a numerical method, second an approximate likelihood is constructed from the numerical solution (Doucet and Johansen, 2009; Kou et al., 2012) . In this section, we introduce the Euler-Maruyama scheme, which is the most commonly used method for approximating (6) (Sun et al., 2015; Allen, 2017) . We also highlight the most prominent drawbacks of inference methods that rely on the Euler-Maruyama scheme.

The Euler-Maruyama scheme approximates the SDE in (6) with the difference equation

where Z = (Z 1 , Z 2 , ..., Z d ) ∼ N (0, I d×d ) and t is the time lag between sequential observations of X N (t). From (7), we obtain the conditional distributions

which are used to perform inference in a likelihood framework.

Inference with the Euler-Maruyama scheme generally requires stochastic infill (Kou et al., 2012) . To see this note that the likelihoods in (8) require X(t) to be observed at all time lags t. Three common reasons Euler-Maruyama requires stochastic infill are;

1) X(t) was observed at irregularly spaced time points, 2) the time lag t is shrunk to improve inference accuracy, and 3) to predict X(t) at unobserved time points. We return to the simulated SIR data set from Section 2 to illustrate the implications of shrinking t.

We assume X(t) is observed at time points t = 0, 5, 10, 15, 20, 25, 30. We could solve the system (6) with a time step of t = 5, however, we would not be able to predict X(t) at any other time points. We consider reducing t from 5 to 1 to reduce numerical error and predict at the 24 unobserved time points t = 1, 2, 3, 4, 6, 7, 8, 9, ..., 26, 27, 28, 29 . This results in 48 latent states, 24 latent states for each of the two subpopulations tracked by the the SIR model.

There is no clear way to ignore or integrate over the latent infill states produced by Euler-Maruyama. In the next section, we construct the the JGDLA as an alternative method for approximate inference on Markov population models that does not require stochastic infill for inference or predictions. We show that the JGDLA is a joint Gaussian likelihood for the observed time points of X(t) constructed from a deterministic system.

The result is a method for approximate inference on Markov population models that relies solely on the solution to a deterministic system.

In this section we construct the joint Gaussian likelihood of the JGDLA. We follow the results of Kurtz (1978) to show that the SDE approximation of the Markov population model converges to a system of ODEs. We use the result to construct the JGDLA covariance from the solution to the system of ODEs.

We summarize the results of Kurtz (1978) which details the sufficient conditions required for the system of SDEs in (6) to converge to a deterministic system. Assume that x is an element in some compact subset K ⊂ R d and the following conditions hold:

Theorem 8.1 of Kurtz (1981) states that for each t > 0,

where

or equivalently

with initial state X † (0). That is, the diffusion approximation of the Markov population process converges to a system of ODEs as the population size N tends towards infinity.

Scientists in fields such as chemistry, ecology, and epidemiology commonly utilize deterministic population models in the form of (11) (Keeling and Rohani, 2011; Fricks and Hanks, 2018) . We defined the Markov population model of Section 2 in terms of reaction vectors R i and rates λ i θ (·). We note that the reaction vectors and rates of the deterministic system in (11) are of the same form as the diffusion approximation in (6). This allows stochastic population models to be constructed from their deterministic analogues. In the following sections, we derive the distribution of the JGDLA from the solution of (11). Kurtz (1978) constructed an approximate distribution for the system of SDEs in (6) from the limiting deterministic system in (11). We present the results of Kurtz (1978) required to obtain an approximate distribution for (1). We begin by considering the difference between X N (t) and its corresponding infinite population limit X † (t). Subtracting equations (3) and (10) gives

Using a first-order Taylor expansion of F (X N (s)) about F X † (s) , (12) becomes

The FCLT ensures that as the population size N → ∞,

in distribution to some zero-mean Gaussian process

We define Q X † (t) to be the d by n matrix with its i th column given by

We differentiate (14) to obtain the system of linear SDEs with timevarying coefficients

where dB(t) = (dB 1 (t), dB 2 (t), ..., dB d (t)) . We then approximate the finite sample process by

We note that the approximate distribution in (16), which was first constructed by Kurtz (1978) , is centered about the solution to the system of ODEs X † (t) in (11) (15), we observe that V (t) is a zero-mean Gaussian process, whose covariance depends only on the deterministic solution X † (t). All past work on this approximation has focused on simulation rather than inference (Kurtz, 1981; Baxendale and Greenwood, 2011; Fricks and Hanks, 2018) . In this work we develop methods for statistical inference on θ using the approximation (16). To this end, we construct the covariance for the JGDLA as a function of the solution to the deterministic system in (11).

To our knowledge, we are the first to construct this covariance and use it for statistical inference.

Here we define the joint Gaussian likelihood of the JGDLA. To do so, we first solve for the covariance of V (t) by solving (15). While this covariance cannot, in general, be obtained analytically, we develop a novel form for the covariance which can easily be approximated numerically. To do so, we propose a separable solution of the form

where U (t) is a d by d matrix and Y (t) is a d-dimensional Gaussian process. The solution to (17) (see Appendix A.2 for details) requires

where ∂F X † (t) and Q X † (t) are defined in (15), and are respectively a d by d matrix of partial derivatives of F X † (t) and a d by n matrix with columns q i X † (t) =

An exact solution to (18) does not exist unless ∂F X † (t) is constant or commutes. For cases in which ∂F X † (t) is not constant and does not commute, (18) is a system of time inhomogeneous linear differential equations that must be solved numerically.

We obtain an approximate solution to (18) by first solving the deterministic system in (11) for X † (t), then numerically solving for U (t) in (18). Since V (t) is a zero-mean Gaussian process, the covariance of V (t) is a linear transformation of the covariance of Y (t). We define the covariance of the zero-mean process Y (t) and obtain the covariance of V (t) by matrix multiplication. For each i = 1, 2, ..., n, we define the d-dimensional vectors a i (t) = (a i1 (t), a i2 (t), ..., a id (t)) = U −1 (t)R i . Then from

It follows from (20) that

The covariance of V (t) is then given by

It follows from (22) that for s < t,

We use the approximation in (16) and (23) to construct a joint Gaussian likelihood for observed time points t 1 , t 2 , t 3 , .., t T . Define observations of the Markov population model

, ..., X N (t T )) and the ODE solution for a given θ as X † =

The resulting JGDLA covariance is

where " * " denotes element-wise multiplication. Our novel JGDLA likelihood is then given

From (25), we see that JGDLA likelihood evaluations rely solely on solutions to a deterministic system. In summary, the JGDLA is obtained by the following procedure. For a set of parameters θ,

1. Solve the system of ODEs in (11) for X † (t). (18) for U (t).

3. Solve the integrals in (21) at all observed time points to construct Σ Y .

4. Construct Σ † θ in (24).

5. Evaluate joint likelihood π X N |θ, X † given in (25).

We note that the integrals involved in the algorithm described above are all deterministic. The integrals are solved numerically, but do not require stochastic infill. Further, predictions at unobserved time points can be obtained via conditional predictive normal distributions. Thus, the JGDLA provides a framework for inference based on the approximate joint likelihood (25) of all observed data which relies solely on the solution to a deterministic system.

In this section we formulate the stochastic SIR model directly from its popular deterministic analogue. We show that the JGDLA offerers improved predictive power in comparison to the Euler-Maruyama scheme of Section 4.

Deterministic SIR population models are widely used in epidemiology to model disease outbreaks (Keeling and Rohani, 2011 ). The SIR model tracks the deterministic proportions of susceptible S † (t), infected I † (t), and recovered R † (t) individuals in a closed population of size N . In large populations, systems of ODEs are commonly used to model disease dynamics (Keeling and Rohani, 2011) . The differential equations governing the deterministic system are given by

where 1 = I † (t) + S † (t) + R † (t) due to the assumption of a constant population size. Note that since R † (t) = 1 − I † (t) − S † (t), the system can be reduced to equations (26) and (27).

The two unknown parameters θ = (β, γ), are the direct transmission rate β ∈ (0, ∞) and the recovery rate γ ∈ (0, ∞).

We construct the stochastic SIR model from the system of ODEs by defining the reaction vectors and their corresponding deterministic rates. Recall the two reactions that may occur previously defined in Section 2, a susceptible individual becomes infected R 1 = (−1, 1) , or an infected recovers R 2 = (0, −1) . We denote the deterministic class proportions X † (t) = S † (t), I † (t) , and define the reaction rates

We let X N (t) = (S N (t), I N (t)) denote the random proportions of susceptible and infected individuals at time t. We obtain the stochastic SIR model by swapping X N (t) and X † (t) in (29-30).

We use the Gillespie algorithm (Gillespie, 1977) 

We first fit the simulated data to the fully deterministic system as a point of comparison for the stochastic models. The fully deterministic model is fit with likelihood given by

Maximum likelihood estimates from the deterministic SIR model likelihood in (31), with mean X † (t) being the numerical optimization solution of (11) and a function of θ = (β, γ), which are estimated. We use ode from the deSolve package in R with a time step of 0.1 to solve for X † (t) and obtain maximum likelihood estimates for (σ, β, γ) using the built-in optimizer optim in R. π X obs N |X † , θ , where X obs N = X obs N (5), X obs N (10), ..., X obs N (30) and X † = X † (5), X † (10), ..., X † (30) . We use ode from the deSolve package in R with a time step of 0.1 to solve all the integrals required to form the JGDLA likelihood. We again obtain maximum likelihood estimates for the JGDLA model optim in R. Details on the construction of the JGDLA are contained in Appendix A.3.

For comparison, we also make inference based on an Euler-Maruyama approximation.

We implement the Euler-Maruyama scheme (EM) for the stochastic SIR model with a time-lag of t = 1 by defining (8) of Section 4. We also consider the independent Euler-Maruyama (EM Ind) model, which is obtained by setting the off-diagonal elements in G θ (X N (t)) of (8) to be zero. Taking G θ (X N (t)) to be diagonal allows for fast inversions and likelihood evaluations, and is common in some settings (Wikle and Hooten, 2010) .

Both Euler-Maruyama schemes require 48 latent states to be estimated at T pred . We fit both Euler-Maruyama schemes via MCMC. We specify iid half normal scale 1 priors for π(β) and π(γ), giving π(θ) = π(β)π(γ). Let X pred denote the collection of X N (t) at infill time points T pred . We estimate the latent infill states by drawing all-at-once Metropolis-Hasting samples from the conditional distribution

where π(X N (t i ) |X N (t i−1 ) , θ) is the conditional normal likelihood obtained from the Euler-Maruyama scheme in (8). We draw 300,000 burn-in states from (32) to allow for the Markov chain to reach a stationary distribution and store an additional 100,000 post burn-in states to perform the analysis.

These three approaches to inference are all compared in terms of mean absolute prediction error (MAPE) on the predicted infected. The MAPE at the 24 infill time points for the infected is computed as M AP E = t∈T pred |I N (t) − I truth N (t)|/24. To compute the MAPE for the JGDLA, we take the expectation of the conditional predictive multivariate normal distribution for X pred conditioned on X obs N and the maximum likelihood estimates for θ. For the Euler-Maruyama schemes, we compute the MAPE as We summarize the MAPE for each approximation method in Table 1 and note that all 95% credible/confidence intervals contained the true value of θ. We see that the JGDLA offers improved predictive power in comparison to both Euler-Maruyama approximations and the deterministic fits. We note that the Euler-Maruyama scheme can be improved by Table 1 : MAPE for the predicted infected proportions for the four models; Euler-Maruyama (EM), Euler-Maruyama with off-diagonal covariance terms set to 0 (EM Ind), JGDLA, and the ODE system fitted to the four simulated data sets with N = 100, 300, 500, 1000.

On 5 February 2020 a cruise ship hosting 3711 people docked for a 2-week quarantine in Yokohama, Japan after a passenger tested positive for the coronavirus disease (COVID-19) (Mizumoto et al., 2020) . Random testing of passengers and crew members began on February 5 th and continued through February 20 th . We note that on February 11 th and 14 th no testing occurred and denote the fourteen observed times T obs = {1, 2, 3, .., 6, 8, 9, 11, 12, . .., 16}.

The number of tests administered per day n t , positive tests per day y t , and total number of individuals remaining on the ship are shown in Table 2 . Once exposed to COVID, susceptible individuals experience an incubation prior to becoming infectious (Chen et al., 2020) . To account for the incubation period, we fit a stochastic susceptible-exposed-infected-removed (SEIR) model to the Princess Diamond cruise ship COVID-19 outbreak data set. We use the JGDLA to form a joint distribution for the latent states of the SEIR model. We specify the system of ODEs governing the SEIR model

where 1 = I † (t) + S † (t) + E † (t) + R † (t). We note that the system contains three com-

The contact rate β, incubation rate α, and the recovery rate γ all have support on the positive reals. Let X † = X † (1), ..., X † (6), X † (8), X † (9), X † (11), ...X † (16) denote the ODE solution at times T obs . We build a joint distribution for the latent proportions X N (t) = (S N (t), E N (t), I N (t)) using JGDLA at the fourteen observed time points X N = (X N (1), ..., X N (6), X N (8), centered at X † .

We estimate θ = (β, α, γ) and the initial conditions X N (0) = (S N (0), E N (0), I N (0)) ,

where t = 0 denotes February 4 th . We account for the disembarkment of susceptible passengers with the inclusion of a fixed time-varying rate α S (t) estimated from passenger records. We model the number of observed seropositive individuals on day t as y t ∼ Binom n t , P (t) = I N (t)

We obtain the stochastic model from the deterministic SEIR by defining the four reaction vectors; a susceptible becomes exposed R 1 = (−1, 1, 0) , a susceptible disembarks from the ship R 2 = (−1, 0, 0) , an exposed individual becomes infected R 3 = (0, −1, 1) , or an infected individual recovers R 4 = (0, 0, −1) . We define the corresponding reac-

We elect to take a Bayesian approach to inference and fit the model via MCMC. We

, center c, and scale parameter d. We draw 100,000 all-at-once Metropolis

Hastings samples of θ and X N (0). We discard 10,000 samples as burn-in states and use the remaining 90,000 samples to perform the analysis.

The likelihood in (33) must be approximated due to its dependence on the latent states X N . We perform a Monte Carlo approximation of (33) by drawing 1,000 samples from the JGDLA joint distribution X (l) N ∼ π X N |θ, X † , assigning P (l) (t) = I 

We use the approximate likelihood in (34) to draw all-at-once Metropolis-Hasting samples from the posterior conditional distribution π(θ, X N (0)|y) ∝π(y|θ)π(X N (0))π(θ), where y is the vector of y t at T obs . Further details on model fitting are included in Appendix A.4.

From Figure 3 , we see that the posterior mean estimate for the probability of being infected captures the mean of the observed proportions well. All parameter estimates are summarized in Table 3 . We note that α is the rate at which individuals move from exposed to infected. Studies have found that symptoms take roughly five to six days to appear (Chen et al., 2020) . Our estimate is most likely lower due to the fact that several tests were administered to a small population (i.e. many were tested positive before symptoms began). We also note that the recovery/removal rate is roughly 1.14 days. This estimate is lower than the known time to recover from SARS-COV2, as individuals on the ship were not tested once found positive, and likely confined to their rooms (i.e. quarantined), and are thus functionally removed from the population, even though they are still infectious.

Our stochastic SEIR model, assisted by the use of the JGDLA, allows us to understand the mean behavior of the small population of the Princess Diamond cruise ship. Table 2 are plotted in black. The deterministic curve for P (t) fitted from the posterior mean estimate of θ is shown in blue, with 95% credible intervals in red.

In this work we proposed the JGDLA as a method for approximate inference on Markov population models. We showed that the JGDLA produces a joint Gaussian distribution that relies solely on the solution to a system of ODEs. We showed that, unlike the Euler-Maruyama scheme, the JGDLA does not require stochastic infill to perform inference or predict at unobserved times. We also illustrated the connection between SDE and ODE modeling by demonstrating how to perform inference on Markov population models directly from a deterministic system in our simulation study and data analysis.

There are other frameworks for approximating the joint likelihood of the data, such as particle filtering (Doucet and Johansen, 2009 ) and iterated filtering (Ionides et al., 2015) .

These methods are commonly used when the Markov population process is assumed to be unobserved. Filtering methods rely on sequentially drawing samples from the filtering and predictive distributions of the latent states. This task requires recursively defining conditional densities from a Gaussian distribution. We fit a hidden Markov model in the COVID-19 data analysis of Section 7 in which the population process was unobserved.

The JGDLA does not require filtering methods to approximate the likelihood. We instead performed a Monte Carlo approximation of the likelihood using samples from the JGLDA's full joint distribution for the latent states.

Inference for JGDLA models is not restricted to Bayesian approaches. In Section 6, we considered a directly observed population process. This allowed for maximum likelihood estimates to be obtained by optimizing the JGDLA joint Gaussian likelihood. In Section 7, the JGDLA produces a joint distribution for the latent states at the time points at which seropositive data counts were observed. Maximum likelihood estimates could have been obtained by optimization, however we elected to use MCMC and found the results to be robust to initial conditions.

In summary, we have constructed the JGDLA as a new method for performing approximate inference on Markov population models. We showed that likelihood evaluations of the JGDLA depend solely on the solution of a system of ODEs. In turn, the JGDLA does not require stochastic infill to predict or perform inference. We saw the JGDLA approximation is built directly from a system of ODEs, connecting ODE and SDE modeling. We also observed that the JGDLA outperformed the Euler-Maruyama approximation in the SIR simulation study considered in Section 6. In conclusion, we suggest the use of the JGDLA as a new framework for performing inference on Markov models.

with rate λt, denoted Y (λt), as n → ∞ we have

where B(t) is a standard Brownian motion, and " ⇒ " denotes convergence in distribution (Kurtz, 1978; Van Kampen, 1992) . We use (35) to perform a Gaussian approximation of a Poisson process

for sufficiently large n.

We seek solutions to (15) of the form

and

where K(t) is d by n matrix that must be solved for. Using (37) and (38), we have

We ensure (39) agrees with (15) by taking

The solution to the SDE in (15) is obtained by sequentially solving (37) and

We denote the observed time points X obs N = X obs N (5), X obs N (10), ..., X obs N (30) . The deterministic SIR model was given in Section 6

We first solve (41-42) for X † (t) = S † (t), I † (t) using initial conditions X N (0) = (0.95, 0.05) .

We used the built-in numerical solver ode in the deSolve package of R with a step size of 0.1, solver setting lsoda, and initial guess θ 0 = (β 0 , γ 0 ). We note that the results were robust to the choice of θ 0 .

We solve dU (t) = ∂F X † (t) U (t), U (0) = I 2×2 ,

where

We solve (43) numerically again using ode with a time step of 0.1 and solver setting vode.

We interpolate a i (s) = U −1 (s)R i across the time domain (0,30) and construct Σ Y by numerically solving V ar (Y 1 (t)) = t 0 βa 2 11 (s)I † (s)S † (s) + γa 2 21 (s)I † (s) ds,

V ar (Y 2 (t)) = t 0 βa 2 12 (s)I † (s)S † (s) + γa 2 22 (s)I † (s) ds,

Cov (Y 1 (t), Y 2 (t)) = t 0 βa 2 11 (s)a 12 (s)I † (s)S † (s) + γa 21 (s)a 22 (s)I † (s) ds, (46) at observed time points t = 5, 10, 15, 20, 25, 30. Next, we construct Σ † θ in (24) and evaluate the JGDLA likelihood

where X † = X † (5), X † (10), ..., X † (30) . We use optim in R to find maximum likelihood estimates by performing the process above iteratively for differing values of θ.

We denote the fourteen times at which COVID tests n t were administrated yielding y t positive tests T obs = {1, 2, 3, .., 6, 8, 9, 11, 12, ..., 16}. We denote the solution to the system of ODEs in Section 7 at time points T obs , X † = X † (1), ..., X † (6), X † (8), X † (9), X † (11), ...X † (16) .

We detail the construction of the JGDLA for the latent proportions X N = (X N (1), ..., X N (6), X N (8), X N ( We first solve for X † conditioned on X N (0) = (S N (0), E N (0), I N (0)) and θ = (β, α, γ)

using the built-in numerical solver ode in the deSolve package of R with a step size of 0.1 and solver setting vode. Next we define

and solve for U (t) in (18) using ode with a step size of 0.1 and solver setting vode. We interpolate a i (s) = U −1 R i on (0,16] and use the result to numerically solve each in integral of (21) using ode with a time step of 0.1 and numerical solver setting vode for time points t ∈ T obs . We then form Σ † θ in (24) to obtain our JGDLA density π X N |X † , θ . We fit the stochastic SEIR model of Section 7 via MCMC. We note that for each iteration of MCMC, we propose θ * and X * N (0) all-at-once. Note that we propose S * (0) > 0 and I * (0) > 0 under the constraints S * N (0) + I * N (0) < 1, R N (0) = 0, and E * N (0) = 1 − S * N (0) − I * N (0). For each proposed value of θ * and X * N (0), the JGDLA must be constructed to obtain π X N | X † * , θ * , where X † * is the ODE solution resulting from initial conditions X * N (0) and θ * . We then perform a Monte Carlo approximation forπ(y|θ * ) by drawing 1,000 samples from π X N | X † * , θ * and following the details of Section 7. A Metropolis-Hasting accept/reject step is then performed with posterior conditional distribution π(θ, X N (0)|y) ∝π(y|θ)π(X N (0))π(θ).

A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis

A comparison of three different stochastic population models with regard to persistence time

Sustained oscillations for density dependent markov processes

Efficient step size selection for the tau-leaping simulation method

A mathematical model for simulating the phase-based transmissibility of a novel coronavirus

A tutorial on particle filtering and smoothing: Fifteen years later

A lattice and random intermediate point sampling design for animal movement

Stochastic population models

Exact stochastic simulation of coupled chemical reactions

Seir modeling of the covid-19 and its dynamics

Inference for dynamic and latent variable models via iterated, perturbed bayes maps

Modeling infectious diseases in humans and animals

A multiresolution method for parameter estimation of diffusion processes

Strong approximation theorems for density dependent markov chains

Approximation of population processes

Estimating the asymptomatic proportion of coronavirus disease 2019 (covid-19) cases on board the diamond princess cruise ship

Seir model for covid-19 dynamics incorporating the environment and social distancing

Stochastic differential equations

Parameter inference and model selection in deterministic and stochastic dynamical models via approximate bayesian computation: modeling a wildlife epidemic

Stochastic processes in physics and chemistry

A general science-based framework for dynamical spatio-temporal models

We state the FCLT for Poisson processes used in Section 3 to perform a Gaussian approximation of the Markov population model in (1). The FCLT states that for a Poisson process