Sequential Monte Carlo algorithms for agent-based models of disease transmission

Agent-based models also called individual-based models are used in many fields, such as social sciences (Epstein, 2006), demographics (Hooten et al., 2020), ecology (DeAngelis and Gross, 2018) and macroeconomics (Turrell, 2016). These models describe the time evolution of a population according to a Markov process. The population comprises of $N\in\mathbb{N}$ agents who interact with one another through a probabilistic mechanism. The popularity of these models could be attributed to the ease of model building and interpretation, while allowing the representation of complex phenomena. Specialized software is available to facilitate their simulation and the visualization of their output (e.g. Tisue and Wilensky, 2004).

The primary use of these models seems to involve “forward simulations” using various parameters corresponding to hypothetical scenarios. In comparison, relatively fewer works have tackled the question of statistical inference for such models, i.e. the estimation of model parameters given available observations, also known as model calibration. Broad articles on the question of fitting agent-based models include Grazzini et al. (2017); Hooten et al. (2020); Hazelbag et al. (2020) and this is a very active research area. The question of estimation is well-posed: by viewing agent-based models as a subclass of hidden Markov models, we can define the associated likelihood function, from which maximum likelihood and Bayesian procedures can be envisioned.

Gibbs sampling or “data augmentation” Markov chain Monte Carlo (MCMC) methods, that alternate between parameter and latent agent states updates are generically applicable, see Fintzi et al. (2017); Bu et al. (2020) in the context of disease transmission. However, the mixing properties of these chains can be problematic, as we illustrate in Appendix C. Likelihood-based inference that avoid the use of data augmentation requires one to marginalize out the latent population process. This is computationally challenging as the number of possible configurations of the population of agents grows exponentially in $N$. Our contribution is to design improved particle filters to estimate the likelihood function of agent-based models, which can then be used for parameter inference procedures such as simulated maximum likelihood and particle MCMC (Andrieu et al., 2010). Particle MCMC methods, based on more standard particle filters, have been used for agent-based models before, see e.g. Kattwinkel and Reichert (2017). For indirect inference and approximate Bayesian computation applied to agent-based models we refer to (Platt, 2020; van der Vaart et al., 2016; Chen et al., 2017; Sirén et al., 2018).

We limit the scope of this article to some agent-based models of disease transmission. We next introduce a specific model which will be used to draw clear connections between agent-based models and more tractable but restrictive formulations of disease outbreak models.

The state of agent $n$ at time $t$, denoted by $X_{t}^{n}$, takes the value $0$ or $1$, corresponding to a “susceptible” or “infected” status in the context of disease transmission. We consider a closed population of size $N$, $T\in\mathbb{N}$ time steps and a stepsize of $h>0$. Initial agent states follow independent Bernoulli distributions, i.e.

$$X_{0}^{n}\sim\text{Ber}(\alpha_{0}),$$

(1)

for $n\in[1:N]=\{1,\ldots,N\}$, where we assume for the time being that each agent has the same probability $\alpha_{0}\in(0,1)$ of being infected. The state of the population $X_{t}=(X_{t}^{n})_{n\in[1:N]}\in\{0,1\}^{N}$ evolves according to a Markov process that depends on the interactions between the agents. The Markov transition specifies that, given the previous states $x_{t-h}\in\{0,1\}^{N}$, the next agent states $(X_{t}^{n})_{n\in[1:N]}$ are conditionally independent Bernoulli variables with probabilities $(\alpha^{n}(x_{t-h}))_{n\in[1:N]}$ given by

$\displaystyle\alpha^{n}(x_{t-h})=\begin{cases}h\lambda N^{-1}\sum_{m=1}^{N}x_{t-h}^{m},\quad&\mbox{if }x_{t-h}^{n}=0,\\ 1-h\gamma,\quad&\mbox{if }x_{t-h}^{n}=1.\end{cases}$

(2)

Here $\lambda\in(0,1)$ and $\gamma\in(0,1)$ represent infection and recovery rate parameters respectively. An uninfected agent becomes infected with probability proportional to $\lambda$ and to the proportion of infected individuals. An infected agent recovers with probability proportional to $\gamma$. The ratio $R_{0}=\lambda/\gamma$ is known as the basic reproductive number, which characterizes the transmission potential of a disease (Lipsitch et al., 2003; Britton, 2010). The latent population process $(X_{t})$ can be related to aggregated observations $(Y_{t})$ of the entire population. For example, the number of infections reported at each time could be modelled as a Binomial distribution, i.e.

$$Y_{t}\mid X_{t}=x_{t}\sim\text{Bin}(I(x_{t}),\rho),$$

(3)

where $I(x_{t})=\sum_{n=1}^{N}x_{t}^{n}\in[0:N]$ counts the number of infected individuals and $\rho\in(0,1)$ represents a reporting rate. Equations (1), (2) and (3) define a susceptible-infected-susceptible (SIS) hidden Markov model.

As the latent state space has cardinality $2^{N}$, evaluating the likelihood function with the forward algorithm would be impractical. Statistical inference is still feasible as the latent process admits low-dimensional summaries. Since (2) assumes that all agents are interchangeable and all pairwise contact are equally likely, the dynamics depends only on the number of infected individuals $I_{t}=I(X_{t})$ and susceptible individuals $S_{t}=N-I_{t}$. More precisely, the Markov transition for each $t>0$ can be summarized as

$\displaystyle S_{t}=S_{t-h}+N^{S}_{t}-N^{I}_{t},\qquad I_{t}=I_{t-h}-N^{S}_{t}+N^{I}_{t},$

(4)

where $N^{S}_{t}\sim\text{Bin}(I_{t-h},h\gamma)$ and $N^{I}_{t}\sim\text{Bin}(S_{t-h},h\lambda N^{-1}I_{t-h})$ denote the number of new recoveries and infections, respectively. The initialization in (1) corresponds to having $I_{0}\sim\text{Bin}(N,\alpha_{0})$ and $S_{0}=N-I_{0}$. The process $(S_{t},I_{t})$ is a Markov chain on the lower-dimensional state space $\{n,m\in[0:N]:n+m=N\}$, and not in $\{0,1\}^{N}$. When combined with (3), the resulting hidden Markov model can be estimated using particle filtering strategies; see e.g. Ionides et al. (2006); Endo et al. (2019); Whiteley and Rimella (2020).

For large population sizes $N$, one can exploit the asymptotic properties of the dynamics (4) to further simplify inference procedures. As $N\rightarrow\infty$, the finite population proportions $\bar{S}_{t}^{N}=S_{t}/N$ and $\bar{I}_{t}^{N}=I_{t}/N$ admit deterministic limits $\bar{S}_{t}$ and $\bar{I}_{t}$, defined by the recursion

$\displaystyle\bar{S}_{t}=\bar{S}_{t-h}+h\gamma\bar{I}_{t}-h\lambda\bar{S}_{t}\bar{I}_{t},\qquad\bar{I}_{t}=\bar{I}_{t-h}-h\gamma\bar{I}_{t}+h\lambda\bar{S}_{t}\bar{I}_{t},$

(5)

with initial condition $(\bar{S}_{t},\bar{I}_{t})=(1-\alpha_{0},\alpha_{0})$. If fine time resolutions are desired, one can also consider the limit $h\rightarrow 0$, in which case $(\bar{S}_{t},\bar{I}_{t})$ converges to a deterministic continuous-time process $(\bar{S}(t),\bar{I}(t))$, defined by the following system of ordinary differential equations

$\displaystyle\frac{d\bar{S}(t)}{dt}=\gamma\bar{I}(t)-\lambda\bar{S}(t)\bar{I}(t),\qquad\frac{d\bar{I}(t)}{dt}=-\gamma\bar{I}(t)+\lambda\bar{S}(t)\bar{I}(t),$

(6)

with $(\bar{S}(0),\bar{I}(0))=(1-\alpha_{0},\alpha_{0})$. We refer to Allen and Burgin (2000) and Allen et al. (2008, Chapter 3) for formal links between stochastic and deterministic SIS models, in the form of limit theorems as $N\to\infty$. The dynamics in (5) or (6) combined with an observation model such as (3) leads again to a fairly tractable statistical model, which can be estimated by maximum likelihood or in the Bayesian framework (e.g. Osthus et al., 2019).

The tractability in (4) and its limits (5)-(6) hinges critically on the assumption that all agents are interchangeable and equally likely to be in contact. These assumptions are strong as it would be desirable to let the infection and recovery rates depend on individual characteristics or features, and the structure of interactions be specified by an undirected network. Adopting any one of these generalizations would result in an agent-based model that cannot be summarized by low-dimensional counts. This article is concerned with the task of inference when these assumptions are relaxed. Our approach is to return to the hidden Markov model formulation with latent space $\{0,1\}^{N}$, and seek novel sampling algorithms that have polynomial rather than exponential costs in $N$. The guiding principle of our contribution is that we can improve the performance of Monte Carlo algorithms, sometimes dramatically, by tailoring them to the models at hand, and by using available observations.

The rest of this article is organized as follows. In Section 2, we consider a “static” case with a single time step to illustrate some of the challenges and proposed solutions in a simple setting. We describe how the cost of likelihood evaluation and sampling the posterior of agent states given observations can be reduced from exponential to polynomial in $N$. In Section 3, we consider the dynamic setting with multiple time steps and design auxiliary particle filters and controlled SMC algorithms tailored to heterogeneous SIS agent-based models. We also provide theoretical support for the approximations employed to make the algorithms practical. We extend our methodology to susceptible-infected-recovered (SIR) models in Section 4 and highlight the associated difficulties. Section 5 summarizes our findings and open questions. This manuscript represents preliminary work and further numerical experiments will be added in a later version. The code is available at https://github.com/nianqiaoju/agents. The proofs and more details are in appendices.

A naive approach to compute the marginal likelihood $p_{\theta}(y)$ is to sum over all possible population configurations

$$p_{\theta}(y)=\sum_{x\in\{0,1\}^{N}}p_{\theta}(x,y)=\sum_{x\in\{0,1\}^{N}}f_{\theta}(x)g_{\theta}(y|x).$$

(8)

This requires $\mathcal{O}(2^{N})$ operations. A simple Monte Carlo approach involves sampling $X^{(p)}=(X^{(p),n})_{n\in[1:N]}\sim f_{\theta}$ independently for $p\in[1:P]$, which represents $P$ possible configurations of the population, and return the Monte Carlo estimator $P^{-1}\sum_{p=1}^{P}g_{\theta}(y|X^{(p)})$ that weights each configuration according to the observation density. The estimator is unbiased and only costs $\mathcal{O}(NP)$ to compute, but its variance might be prohibitively large for practical values of $P$, depending on the observation model $g_{\theta}(y|x)$. Another issue with this estimator is that it can collapse to zero whenever $I(X^{(p)})<y$ for all $p\in[1:P]$. Following Del Moral et al. (2015), this can be circumvented by repeatedly drawing samples $X^{(p)}\sim f_{\theta}$ independently until there are $P$ configurations that satisfy $I(X^{(p)})\geq y$, and return the estimator ${(R-1)}^{-1}\sum_{p=1}^{R-1}g_{\theta}(y|X^{(p)})$, where $R\geq P$ denotes the number of repetitions needed. The resulting estimator can be shown to be unbiased and has a random cost of $\mathcal{O}(N\,\mathbb{E}[R])$, that depends on the value of $y$. An obvious shortcoming of the above estimators is that the agents are sampled from $f_{\theta}$ without using the available observation $y$.

We can in fact reduce the cost of computing (8) without resorting to Monte Carlo approximations. The starting observation is that $I=I(X)$ under $X\sim f_{\theta}$ is the sum of $N$ independent Bernoulli random variables with non-identical success probabilities, and follows a distribution called “Poisson Binomial”. We will refer to this distribution as PoiBin($\alpha$) and write its PMF as

$$\text{PoiBin}(i;\alpha)=\sum_{x\in\{0,1\}^{N}}\mathbbm{1}\left(\sum_{n=1}^{N}x^{n}=i\right)\prod_{n=1}^{N}(\alpha^{n})^{x^{n}}(1-\alpha^{n})^{(1-x^{n})},\quad i\in[0:N].$$

(9)

Exact evaluation of (9) has been considered in Barlow and Heidtmann (1984); Chen et al. (1994); Chen and Liu (1997); Hong (2013) using different approaches. Defining $q(i,n)=\text{PoiBin}(i;\alpha^{n:N})$ for $i\in[0:N]$ and $n\in[1:N]$, we will employ the following recursion

$$q(i,n)=\alpha^{n}q(i-1,n+1)+(1-\alpha^{n})q(i,n+1),\quad i\in[1:N],n\in[1:N-1],$$

(10)

with initial conditions $q(0,n)=\prod_{m=n}^{N}(1-\alpha^{m})$ for $n\in[1:N]$, $q(1,N)=\alpha^{N}$ and $q(i,N)=0$ for $i\in[2:N]$. The desired PMF $q(i,1)=\text{PoiBin}(i;\alpha)$ for $i\in[0:N]$ can thus be computed in $\mathcal{O}(N^{2})$ operations; see Appendix A.1 for a derivation of (10).

Using the above observation, we can rewrite the marginal likelihood as

$$p_{\theta}(y)=\sum_{i=0}^{N}\text{PoiBin}(i;\alpha)\text{Bin}(y;i,\rho)\mathbbm{1}(i\geq y),$$

(11)

which costs $\mathcal{O}(N^{2})$ to compute. Using a thinning argument detailed in Appendix A.2, the above sum is in fact equal to $p_{\theta}(y)=\text{PoiBin}(y;\rho\,\alpha)$. Although the marginal likelihood will not admit such tractability in the general setup considered in Section 3, the preceding observations will guide our choice of approximations.

One can also rely on approximations of the Poisson Binomial distribution to further reduce the cost of computing (11) to $\mathcal{O}(N)$. Choices include the Poisson approximation (Hodges and Le Cam, 1960; Barbour and Hall, 1984; Wang, 1993), the Normal approximation (Volkova, 1996) and the translated Poisson approximation (Barbour and Xia, 1999; Cekanavicius and Vaitkus, 2001; Barbour and Ćekanavićius, 2002). We will focus on the translated Poisson approximation which exactly matches the mean and approximately the variance of a Poisson Binomial distribution. Let $\text{Poi}(\lambda)$ denote a Poisson distribution with rate $\lambda>0$ and write the mean and variance of PoiBin($\alpha$) as $\mu=\sum_{n=1}^{N}\alpha^{n}$ and $\sigma^{2}=\sum_{n=1}^{N}\alpha^{n}(1-\alpha^{n})$, respectively. The translated Poisson approximation of (9) is given by

$$\text{TransPoi}(i;\mu,\sigma^{2})=\begin{cases}0,&i\in[0:\lfloor\mu-\sigma^{2}\rfloor-1],\\ \text{Poi}(i-\lfloor\mu-\sigma^{2}\rfloor;\sigma^{2}+\{\mu-\sigma^{2}\}),&i\in[\lfloor\mu-\sigma^{2}\rfloor:N],\end{cases}$$

(12)

where $\lfloor\mu-\sigma^{2}\rfloor$ and $\{\mu-\sigma^{2}\}$ denote the floor and fractional part of $\mu-\sigma^{2}$, respectively. Since $\mu$ and $\sigma^{2}$ can be computed in $\mathcal{O}(N)$ operations, the translated Poisson approximation (12) and the resulting approximation of (11) only require $\mathcal{O}(N)$ operations. Hence this can be appealing in the setting of large population sizes $N$ at the expense of an approximation error that is well-studied. Indeed results in Cekanavicius and Vaitkus (2001, Theorem 2.1 & Corollary 2.1) and Barbour and Ćekanavićius (2002, Theorem 3.1) imply that the error, measured in the total variation distance, decay at the rate of $N^{-1/2}$ as $N\rightarrow\infty$.

Sampling from the posterior distribution $p_{\theta}(x|y)$ by naively enumerating over all $2^{N}$ configurations is computationally impractical. A key observation is that the observation density in (7) depends on the high dimensional latent state $X\in\{0,1\}^{N}$ only through the one-dimensional summary $I(X)\in[0:N]$. This prompts the inclusion of $I=I(X)$ as an auxiliary variable. Thus the joint posterior distribution can be decomposed as

$\displaystyle p_{\theta}(x,i|y)=p_{\theta}(i|y)p_{\theta}(x|i).$

(13)

The dominant cost of sampling the posterior distribution of the summary

$\displaystyle p_{\theta}(i|y)=\frac{p_{\theta}(i)p_{\theta}(y|i)}{p_{\theta}(y)}=\frac{\text{PoiBin}(i;\alpha)\text{Bin}(y;i,\rho)\mathbbm{1}(i\geq y)}{p_{\theta}(y)},\quad i\in[0:N],$

(14)

is the evaluation of the Poisson Binomial PMF (9). Recall from Section 2.1 that this can be done exactly in $\mathcal{O}(N^{2})$, and in $\mathcal{O}(N)$ using a translated Poisson approximation. The conditional distribution of the latent state given the summary

$\displaystyle p_{\theta}(x|i)=\frac{p_{\theta}(x)p_{\theta}(i|x)}{p_{\theta}(i)}=\frac{\prod_{n=1}^{N}(\alpha^{n})^{x^{n}}(1-\alpha^{n})^{1-x^{n}}\mathbbm{1}(I(x)=i)}{\text{PoiBin}(i;\alpha)}$

(15)

is known as a conditional Bernoulli distribution, which we will write as $p_{\theta}(x|i)=\text{CondBer}(x;\alpha,i)$. Various sampling schemes have been proposed (Fan et al., 1962; Chen et al., 1994); see Chen and Liu (1997, Section 4) for an overview. We will rely on the sequential decomposition

	$\displaystyle p_{\theta}(x|i)$	$\displaystyle=p_{\theta}(x^{1}|i)\prod_{n=2}^{N}p_{\theta}(x^{n}|x^{1:n-1},i)$		(16)
		$\displaystyle=\prod_{n=1}^{N-1}\frac{(\alpha^{n})^{x^{n}}(1-\alpha^{n})^{1-x^{n}}q(i-i_{n-1}-x^{n},n+1)}{q(i-i_{n-1},n)}\mathbbm{1}(x^{N}=i-i_{N-1}),$

where $i_{0}=0$ and $i_{n}=\sum_{m=1}^{n}x^{m}$ for $n\in[1:N-1]$. A derivation of (16) can be found in Appendix A of Heng et al. (2020b). The values of $q(j,n)$ for $j\in[0:i],n\in[1:N]$ needed in (16) can be precomputed using the same recursion as (10) in $\mathcal{O}(iN)$ cost. This precomputation is not necessary if these values are stored when computing Poisson Binomial probabilities.

To reduce the cost we can employ a Markov chain Monte Carlo (MCMC) method to sample from the conditional Bernoulli distribution $\text{CondBer}(\alpha,i)$. This incurs a cost of $\mathcal{O}(1)$ per iteration and converges in $\mathcal{O}(N\log N)$ iterations, under some mild assumptions on $\alpha$, as shown in Heng et al. (2020b). On a related note, we can design MCMC algorithms to target $p_{\theta}(x|y)$. These might for example update the state of a few agents by sampling from their conditional distributions given every other variables, and alternately propose to swap zeros and ones in the vector $X\in\{0,1\}^{N}$. Each of these steps can be done with a cost independent of $N$, but the number of iterations for the algorithm to converge is expected to grow at least linearly with $N$.

We set up numerical experiments as follows. We generate covariates $(w^{n})$ from $\mathcal{N}(4,1)$ for $N=1000$ individuals independently, where $\mathcal{N}(\mu,\sigma^{2})$ denotes a Normal distribution with mean $\mu\in\mathbb{R}$ and variance $\sigma^{2}>0$. Individual specific infection probabilities are computed as $\alpha^{n}=(1+\exp(-\beta w^{n}))^{-1}$ using $\beta=0.3$. We then simulate $X^{n}\sim\text{Ber}(\alpha^{n})$ independently for all $n$, and sample $Y\sim\text{Bin}(\sum_{n=1}^{N}X^{n},\rho)$, with $\rho=0.8$. We adopt a Bayesian approach and assign an independent prior of $\mathcal{N}(0,1)$ on $\beta$ and $\text{Uniform}(0,1)$ on $\rho$. We focus on the task of sampling from the posterior distribution of $\theta=(\beta,\rho)$ given $y$ and the covariates $(w^{n})$.

We consider random walk Metropolis–Hastings with exact likelihood calculation (“MH-exact”), associated with a quadratic cost in $N$. We also consider the same algorithm with a likelihood approximated by a translated Poisson (“MH-tp”), with a cost linear in $N$. Finally we consider a pseudo-marginal approach (Andrieu and Roberts, 2009; Andrieu et al., 2010) with likelihood estimated with $P=20$ particles sampled from $f_{\theta}$ (“PMMH”). Note that $P=20$ samples resulted in a variance of the log-likelihood estimates of approximately $0.3$ at the data generating parameters (DGP). These samplers employ the same random walk, based on Normal proposals, independently on $\beta$ and $\log(\rho/(1-\rho))$ with standard deviation of $0.2$. As a baseline we also consider a Gibbs sampler that alternates between the updates of $\theta$ given $X$, employing the same proposal as above, and updates of $X$ given $\theta$. These employ an equally weighted mixture of kernels, performing either $N$ random swap updates, or a systematic Gibbs scan of the $N$ components of $X$; thus the cost per iteration is linear in $N$.

We first run “MH-exact” with $100,000$ MCMC iterations (after a burn-in of $5,000$ iterations), to obtain estimates of the posterior means of $\beta$ and $\rho$. Using these posterior estimates as ground truth, we compute the mean squared error (MSE) of the posterior approximations obtained using each method with $20,000$ MCMC iterations (excluding a burn-in of $5000$) and 50 independent repetitions. Table 1 displays the MSE, as well as the relative wall-clock time to obtain each estimate. Comparing the results of MH-exact and MH-tp shows that it is possible to save considerable efforts at the expense of small differences in the parameter estimates using a translated Poisson approximation. The appeal of the PMMH approach compared to exact likelihood calculations is also clear. On the other hand, Gibbs samplers that alternate between the updates of $\theta$ and $X$ do not seem to be competitive in this example.

SMC methods (Liu and Chen, 1998; Doucet et al., 2001), also known as particle filters, provide approximations of $p_{\theta}(y_{0:T})$ and $p_{\theta}(x_{0:T}|y_{0:T})$ by simulating an interacting particle system of size $P\in\mathbb{N}$. In the following, we give a generic description of SMC to include several algorithms in a common framework.

At the initial time, one samples $P$ configurations of the population from a proposal distribution $q_{0}$ on $\{0,1\}^{N}$, i.e. $X_{0}^{(p)}=(X_{0}^{(p),n})_{n\in[1:N]}\sim q_{0}(\cdot|\theta)$ independently for $p\in[1:P]$. Each possible configuration is then assigned a weight $W_{0}^{(p)}\propto w_{0}(X_{0}^{(p)})$ that is normalized to sum to one. To focus our computation on the more likely configurations, we perform an operation known as resampling that discards some configurations and duplicates others according to their weights. Each resampling scheme involves sampling ancestor indexes $(A_{0}^{(p)})_{p\in[1:P]}\in[1:P]^{P}$ from a distribution $r(\cdot|W_{0}^{(1:P)})$ on $[1:P]^{P}$. The simplest scheme is multinomial resampling (Gordon et al., 1993), which samples $(A_{0}^{(p)})_{p\in[1:P]}$ independently from the categorical distribution on $[1:P]$ with probabilities $W_{0}^{(1:P)}$; other lower variance and adaptive resampling schemes can also be employed (Fearnhead and Clifford, 2003; Gerber et al., 2019). Subsequently, for time step $t\in[1:T]$, one propagates each resampled configuration according to a proposal transition $q_{t}$ on $\{0,1\}^{N}$, i.e. $X_{t}^{(p)}=(X_{t}^{(p),n})_{n\in[1:N]}\sim q_{t}(\cdot|X_{t-1}^{(A_{t-1}^{(p)})},\theta)$ independently for $p\in[1:P]$. As before, we weight each new configuration according to $W_{t}^{(p)}\propto w_{t}(X_{t-1}^{(A_{t-1}^{(p)})},X_{t}^{(p)})$, and for $t<T$, resample by drawing the ancestor indexes $(A_{t}^{(p)})_{p\in[1:P]}\sim r(\cdot|W_{t}^{(1:P)})$. For notational simplicity, we suppress notational dependence of the weight functions $(w_{t})_{t\in[0:T]}$ on the parameter $\theta$. To approximate the desired quantities $p_{\theta}(y_{0:T})$ and $p_{\theta}(x_{0:T}|y_{0:T})$, these weight functions have to satisfy

$\displaystyle w_{0}(x_{0})\prod_{t=1}^{T}w_{t}(x_{t-1},x_{t})=\frac{\mu_{\theta}(x_{0})\prod_{t=1}^{T}f_{\theta}(x_{t}|x_{t-1})\prod_{t=0}^{T}g_{\theta}(y_{t}|x_{t})}{q_{0}(x_{0}|\theta)\prod_{t=1}^{T}q_{t}(x_{t}|x_{t-1},\theta)}.$

(21)

Given the output of the above simulation, an unbiased estimator of the marginal likelihood $p_{\theta}(y_{0:T})$ is

$\displaystyle\hat{p}_{\theta}(y_{0:T})=\left\{\frac{1}{P}\sum_{p=1}^{P}w_{0}(X_{0}^{(p)})\right\}\left\{\prod_{t=1}^{T}\frac{1}{P}\sum_{p=1}^{P}w_{t}(X_{t-1}^{(A_{t-1}^{(p)})},X_{t}^{(p)})\right\},$

(22)

and a particle approximation of the smoothing distribution is given by

$\displaystyle\hat{p}_{\theta}(x_{0:T}|y_{0:T})=\sum_{p=1}^{P}W_{T}^{(p)}\delta_{X_{0:T}^{(p)}}(x_{0:T}).$

(23)

In the latter approximation, each trajectory $X_{0:T}^{(p)}$ is formed by tracing the ancestral lineage of $X_{T}^{(p)}$, i.e. $X_{0:T}^{(p)}=(X_{t}^{(B_{t}^{(p)})})_{t\in[0:T]}$ with $B_{T}^{(p)}=p$ and $B_{t}^{(p)}=A_{t}^{(B_{t+1}^{(p)})}$ for $t\in[0:T-1]$. Convergence properties of these approximations as the size of the particle system $P\rightarrow\infty$ are well-studied; see for example Del Moral (2004). However, the quality of these approximations depends crucially on the choice of proposals $(q_{t})_{t\in[0:T]}$ and the corresponding weight functions $(w_{t})_{t\in[0:T]}$ that satisfy (21).

The BPF of Gordon et al. (1993) can be recovered by employing the proposals $q_{0}(x_{0}|\theta)=\mu_{\theta}(x_{0}),q_{t}(x_{t}|x_{t-1},\theta)=f_{\theta}(x_{t}|x_{t-1})$ for $t\in[1:T]$ and the weight functions $w_{t}(x_{t})=g_{\theta}(y_{t}|x_{t})$ for $t\in[0:T]$. Although the BPF only costs $\mathcal{O}(NTP)$ to implement and has convergence guarantees as $P\rightarrow\infty$, the variance of its marginal likelihood estimator can be too large to deploy within particle MCMC schemes for practical values of $P$ (see Section 3.3). Another issue with its marginal likelihood estimator, for this particular choice of observation equation, is that it can collapse to zero if all proposed configurations have less infections than the observed value, i.e. there exists $t\in[0:T]$ such that $I(X_{t}^{(p)})<y_{t}$ for all $p\in[1:P]$. With increased cost, this issue can be circumvented using the alive particle filter of Del Moral et al. (2015), by repeatedly drawing samples at each time step until there are $P$ configurations with infections that are larger than or equal to the observed value.

Alternatively, one can construct better proposals with supports that respect these observational constraints. One such option is the fully adapted APF (Pitt and Shephard, 1999; Carpenter et al., 1999) that corresponds to having the proposals $q_{0}(x_{0}|\theta)=p_{\theta}(x_{0}|y_{0}),q_{t}(x_{t}|x_{t-1},\theta)=p_{\theta}(x_{t}|x_{t-1},y_{t})$ for $t\in[1:T]$ and the weight functions $w_{0}(x_{0})=p_{\theta}(y_{0})$ and $w_{t}(x_{t-1})=p_{\theta}(y_{t}|x_{t-1})$ for $t\in[1:T]$. At the initial time, computing the marginal likelihood $p_{\theta}(y_{0})$ and sampling from the posterior of agent states $p_{\theta}(x_{0}|y_{0})$ can be done exactly (or approximately) as described in Sections 2.1 and 2.2 respectively for the static model. More precisely, we compute

$\displaystyle p_{\theta}(y_{0})=\sum_{i_{0}=0}^{N}\text{PoiBin}(i_{0};\alpha_{0})\text{Bin}(y_{0};i_{0},\rho)\mathbbm{1}(i_{0}\geq y_{0})$

(24)

and sample from

$\displaystyle p_{\theta}(x_{0},i_{0}|y_{0})=p_{\theta}(i_{0}|y_{0})p_{\theta}(x_{0}|i_{0})=\frac{\text{PoiBin}(i_{0};\alpha_{0})\text{Bin}(y_{0};i_{0},\rho)\mathbbm{1}(i_{0}\geq y_{0})}{p_{\theta}(y_{0})}\text{CondBer}(x_{0};\alpha_{0},i_{0}),$

(25)

which admits $p_{\theta}(x_{0}|y_{0})$ as its marginal distribution. For time step $t\in[1:T]$, by conditioning on the previous configuration $x_{t-1}\in\{0,1\}^{N}$, the same ideas can be used to compute the predictive likelihood $p_{\theta}(y_{t}|x_{t-1})$ and sample from the transition $p_{\theta}(x_{t}|x_{t-1},y_{t})$, i.e. we compute

$\displaystyle p_{\theta}(y_{t}|x_{t-1})=\sum_{i_{t}=0}^{N}\text{PoiBin}(i_{t};\alpha(x_{t-1}))\text{Bin}(y_{t};i_{t},\rho)\mathbbm{1}(i_{t}\geq y_{t})$

(26)

and sample from

	$\displaystyle p_{\theta}(x_{t},i_{t}|x_{t-1},y_{t})$	$\displaystyle=p_{\theta}(i_{t}|x_{t-1},y_{t})p_{\theta}(x_{t}|x_{t-1},i_{t})$		(27)
		$\displaystyle=\frac{\text{PoiBin}(i_{t};\alpha(x_{t-1}))\text{Bin}(y_{t};i_{t},\rho)\mathbbm{1}(i_{t}\geq y_{t})}{p_{\theta}(y_{t}|x_{t-1})}\text{CondBer}(x_{t};\alpha(x_{t-1}),i_{t})$

which admits $p_{\theta}(x_{t},|x_{t-1},y_{t})$ as its marginal transition.

An algorithmic description of the resulting APF is detailed in Algorithm 1, where the notation $\text{Cat}([0:N],V^{(0:N)})$ refers to the categorical distribution on $[0:N]$ with probabilities $V^{(0:N)}=(V^{(i)})_{i\in[0:N]}$. As the weights in the fully adapted APF at time $t\in[1:T]$ only depend on the configuration at time $t-1$, note that we have interchanged the order of sampling and resampling to promote sample diversity. The cost of running APF exactly is $\mathcal{O}(N^{2}TP)$. To reduce the computational cost to $\mathcal{O}(N\log(N)TP)$ one can approximate the above Poisson binomial PMFs with the translated Poisson approximation (12), and employ MCMC to sample from the above conditioned Bernoulli distributions.