figurec

# Bayesian inference in Epidemics: linear noise analysis

Samuel Bronstein1 1Department of Mathematics and Applications, ENS Paris,
75005 Paris, France.  ,  Stefan Engblom2 and  Robin Marin2 2Division of
Scientific Computing
Department of Information Technology
Uppsala University
SE-751 05 Uppsala, Sweden<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Date: August 28, 2024)

###### Abstract.

This paper offers a qualitative insight into the convergence of Bayesian
parameter inference in a setup which mimics the modeling of the spread of a
disease with associated disease measurements. Specifically, we are interested
in the Bayesian model’s convergence with increasing amounts of data under
measurement limitations. Depending on how weakly informative the disease
measurements are, we offer a kind of ‘best case’ as well as a ‘worst case’
analysis where, in the former case, we assume that the prevalence is directly
accessible, while in the latter that only a binary signal corresponding to a
prevalence detection threshold is available. Both cases are studied under an
assumed so-called linear noise approximation as to the true dynamics.
Numerical experiments test the sharpness of our results when confronted with
more realistic situations for which analytical results are unavailable.

###### Key words and phrases:

Parameter estimation; Bayesian modeling; Stochastic epidemiological models;
Network model; Ornstein-Uhlenbeck process

###### 2010 Mathematics Subject Classification:

Primary: 60J70, 62F12, 62F15; Secondary: 65C30, 65C60, 92D30

Corresponding author: S. Engblom, telephone +46-18-471 27 54, fax +46-18-51 19
25, URL: http://user.it.uu.se/~stefane.

##  1\. Introduction

Computational models in epidemics are commonly relied upon to estimate the
disease spread at fairly large spatial- and temporal scales, often referred to
as _scenario generation_. With the increasing volumes and improved resolution
of data from, e.g., mobile apps and disease testing time series from hospitals
and nursing homes, predictive data-driven models formed from first principles
are within reach. The accuracy of such models is ultimately limited by the
specifics of the available disease surveillance data. In this paper we attempt
to gain a qualitative understanding of how Bayesian inference of
epidemiological parameters may be expected to perform and what the limiting
factors are.

Epidemiological models are typically formed by postulating laws for the flow
of individuals between different _compartments_ in a large population and have
been studied in this form for a long time. When connected with data in the
form of observations, the associated inference problem is also a fairly mature
field, see, e.g., [22, 29, 1]. With the increasing qualities and quantities of
data, the various incarnations of data-driven modeling allow for substantially
higher modeling resolution compared to traditional macroscopic approaches [13,
14, 4, 30, 8]. For example, individual-level contact tracing has been used to
study disease spread models at various population sizes [43, 3, 35, 33, 48].
Data-driven models have allowed epidemic and endemic conditions to be
investigated at a level of detail not previously possible [53, 27, 50].

The identification of the epidemiological parameters from data falls under the
scope of problems formally studied in _System Identification_ [41]. However, a
‘system’ viewpoint of epidemiological modeling is not yet standard, and
identification of parameters is rather more often approached through
calibration of residuals [14, 50, 16], sometimes also blending in aspects of
Bayesian arguments. Fully Bayesian approaches are rarer, albeit with some
exceptions [8, 11], typically due to the technical difficulties with
formulating suitable (pseudo-)likelihoods and the slow convergence associated
with the conditioning of the problem. Although Bayesian inference is notably
well-posed thanks to its use of prior distributions, the posterior
distribution itself is often a computationally ill-conditioned object whenever
strong parameter correlations are present, e.g., resulting from nearly
singular maps from parameters to observables. These conditions, together with
the societal importance of this modeling domain, make it relevant to reason
around the limits of fully Bayesian techniques.

Fundamental questions concerning Bayesian convergence in general settings,
including infinite-dimensional ones, have been treated [37, 46], and have also
been revisited with specific tools and applications in mind, e.g., Gaussian
processes and other machine learning algorithms [45, 26]. “Brittleness” [34,
42], or high sensitivity to small perturbations have been proposed to be a
problematic phenomenon under certain conditions. For applications, this points
to the importance of striking a balance between the granularity of the model
and the information content and level of detail of the available data.

With the specific aim of reaching qualitative conclusions for Bayesian
modeling in epidemiology, we will connect some usual modeling approaches in
infectious disease spread with a basic linear stochastic differential
equation. We analyze the pre-asymptotic Bayesian posterior in the high-quality
data regime as well as the estimate’s convergence under weakly informative
disease measurements. This is motivated by current trends in disease spread
monitoring through for example sewage water analysis and symptoms data
collection using smartphones [17, 23].

The rest of the paper is organized as follows. In §2 we suggest the Ornstein-
Uhlenbeck process as a meta-model of epidemiological models, covering various
Susceptible-Infected-Recovered (SIR)-type models locally in time, during
endemic as well as under epidemic conditions. In §3 we briefly analyze the
posterior convergence when accurate full state measurements are available. The
main results are found in §4 where we more fully develop an analysis for the
case of poorly informative data. We offer some examples of relevance in §5 and
a summarizing discussion is found in §6.

##  2\. The Ornstein-Uhlenbeck process

We connect in this section the Ornstein-Uhlenbeck process with some basic
epidemiological models via linearization and quasi steady-state arguments. We
also briefly discuss, by means of a backward analysis, how continuous and
discrete time in this setting can be connected in probability law under a
certain parameter map.

###  2.1. Disease spread modeling

Epidemiological modeling typically involves ordinary differential equations
(ODEs), e.g., the SIR-model which is often used to model the annual flu [24,
22],

(2.4)
$\displaystyle\left.\begin{array}[]{rl}S^{\prime}(t)&=-(\beta\Sigma^{-1})S(t)I(t)\\\
I^{\prime}(t)&=(\beta\Sigma^{-1})S(t)I(t)-\gamma I(t)\\\ R^{\prime}(t)&=\gamma
I(t)\end{array}\right\\}$

in terms of Susceptible, Infected, and Recovered individuals, respectively,
and for a total population size $\Sigma\equiv S+I+R$. The model parameters
$\beta$ and $\gamma$ define the transmission and recovery rates, and the basic
reproduction number is then given by $R_{0}=\beta/\gamma$, that is, the
expected number of new infections resulting from a single index case [22].

Simplifying the SIR-model by removing the $R$-compartment one obtains the SIS-
model, where the recovered state has been removed and hence effectively
identified with the susceptible state,

(2.7)
$\displaystyle\left.\begin{array}[]{rl}S^{\prime}(t)&=-(\beta\Sigma^{-1})S(t)I(t)+\gamma
I(t)\\\ I^{\prime}(t)&=(\beta\Sigma^{-1})S(t)I(t)-\gamma I(t)\\\
\end{array}\right\\}$

and with the same $R_{0}$ as the SIR-model. It is sometimes useful to add an
environmental compartment expressing the _infectious pressure_ $\varphi$,
i.e., the amount of infectious substance per unit of space. This defines the
SIS${}_{\text{E}}$-model [10], where the $E$ stands for the Environmental
compartment,

(2.11) $\displaystyle\left.\begin{array}[]{rl}S^{\prime}(t)&=-\beta
S(t)\varphi(t)+\gamma I(t)\\\ I^{\prime}(t)&=\beta S(t)\varphi(t)-\gamma
I(t)\\\
\varphi^{\prime}(t)&=\Sigma^{-1}I(t)-\rho\varphi(t)\end{array}\right\\}$

where the indicated governing equation for $\varphi$ is just a basic example.
The SIS${}_{\text{E}}$-model is convenient to adapt to spread over a network
and to detection situations involving sampling the environment [11]. Due to
the indirect transmission the basic reproduction number now scales as a square
root, with $R_{0}=\sqrt{\beta/(\gamma\rho)}$.

The use of ODEs can be justified for epidemiological models in sufficiently
large populations. In smaller populations, e.g., networks of small
communities, stochastic variants are necessary to properly capture the
underlying dynamics of the spread [7, 1]. For example, a stochastic
differential equation (SDE)-version of (2.7) was analyzed in [20] which
essentially replaces $\beta\,dt$ by a Brownian diffused version
$\beta\,dt+\eta\,dB(t)$. This has the effect of lowering the basic
reproduction number to $R_{0}=(\beta-\eta^{2}/2)/\gamma$. A first-principle
stochastic approach is to rather express the dynamics as a continuous-time
discrete state Markov chain (CTMC). The SIS-model above is then defined via
discrete transitions with exponentially distributed waiting times,

(2.14)
$\displaystyle\left.\begin{array}[]{rcl}S+I&\xrightarrow{\beta\Sigma^{-1}}&2I\\\
I&\xrightarrow{\gamma}&S\end{array}\right\\}$

meaning, e.g., that one infectious individual may infect one susceptible
individual, and such that $\beta$ and $\gamma$ are now understood as rate
parameters in the driving Poissonian processes. One can show that now
$R_{0}=(1-\Sigma^{-1}/2)\times\beta/\gamma$ such that consistency with the
previous SDE formulation would require that $\eta^{2}=\Sigma^{-1}\beta$ and
hence follows from the typical Poissonian population-dependent noise scaling
$\eta\propto\Sigma^{-1/2}$.

Let $P(t)$ be some given measure of the intensity of the disease, such as the
absolute or relative disease prevalences, $I(t)$ or $\Sigma^{-1}I(t)$,
respectively. Another alternative could be the infectious pressure
$\varphi(t)$, which also informs on the current disease intensity. Gathering
data from the disease spread means collecting information about $P(t)$,
typically in the form of a time-series, $(F(P(t_{i}))_{i}$, say, for some
measurement operator $F$. Consider _endemic_ conditions first, that is, $P$ is
considered stationary. As an ansatz, suppose

(2.15) $\displaystyle P(t_{i})$
$\displaystyle\sim\mathbf{P}[P=p]=\rho_{\infty}(p)\propto\exp\left(-2/\sigma^{2}\,E(p)\right),$

for some _epidemiological potential_ $E$. Under endemic conditions we can
expect this to be a single well potential, say,

(2.16) $\displaystyle E(p)$
$\displaystyle\propto(p-\alpha)^{2}+\mbox{constant,}$

or at least can locally be approximated by this form. This is consistent with
the _Ornstein-Uhlenbeck ( $\mathcal{OU}$-) process_ [49],

(2.17) $\displaystyle X(0)$ $\displaystyle=X_{0},$ (2.18) $\displaystyle
dX(t)$ $\displaystyle=k-\mu X(t)\,dt+\sigma\,dW(t),$

where $X(t)\in\mathbf{R}$ and where the parameters of the model are
$\theta=[k,\mu,\sigma]\in\mathbf{R}_{+}^{3}$ and in (2.16), $\alpha=k/\mu$.

As a concrete example, we may approximate the continuous-time Markov chain
(2.14) by linearizing the rates around the non-trivial stationary state
$I_{\infty}=\Sigma(1-\gamma/\beta)$ of (2.7). Similarly, the noise term
$\sigma$ can be determined by inspecting the total variance of the Poissonian
rates in (2.14) around this equilibrium. This yields a _linear noise_
$\mathcal{OU}$-approximation to the state variable $I$ of the SIS-model (2.14)
with parameters

(2.22)
$\displaystyle\left.\begin{array}[]{rll}k&=\Sigma(\beta-\gamma)^{2}/\beta&=\gamma\Sigma(R_{0}-1)^{2}/R_{0}\\\
\mu&=\beta-\gamma&=\gamma(R_{0}-1)\\\
\sigma^{2}&=2\Sigma(\beta-\gamma)^{2}/\beta&=2\gamma\Sigma(R_{0}-1)/R_{0}\end{array}\right\\}$
where we use the approximation
$R_{0}=(1-\Sigma^{-1}/2)\times\beta/\gamma\approx\beta/\gamma$. Hence, _with
this specific interpretation of the $\mathcal{OU}$-process,_ (2.23)
$\displaystyle
R_{0}=\left(1-\Sigma^{-1}k/\mu\right)^{-1}=\left(1-\Sigma^{-1}\alpha\right)^{-1}.$

The quality of this approximation is exemplified in Fig. 2.1; the effect of
linearizing around the stationary state can be seen as a slightly too fast
transient compared to the Markov chain.

Figure 2.1. The approximation (2.22) exemplified for $R_{0}=[1.1,1.5,2]$
(bottom and up), with $\Sigma=1000$ and $\gamma=1$. For comparison the ODE-,
the $\mathcal{OU}$-, and the continuous-time Markov chain (CTMC)
interpretation are shown.

During _epidemic_ conditions we are rather sampling a transient of the process
and, moreover, there are also typically many kinds of feedback involved, e.g.,
from minor adjustments of individual level behavior, to major societal changes
and governmental intervention strategies. Although this means that there is
now a greater challenge in formulating a data-centric meta-model of the
situation we may still consider a window of time for which the disease spread
parameters are approximately constant and such that (2.17)–(2.18) remain a
relevant model of the situation. A difference is then that data is gathered
out of equilibrium (hence away from $\alpha$) and presumably also under more
noisy conditions with larger values of $\sigma$. More general epidemiological
potentials could be considered through the SDE in gradient formulation,

(2.24) $\displaystyle dP(t)$ $\displaystyle=-\nabla E(P_{t})\,dt+\sigma
dW(t),$

for which the stationary measure is still given by (2.15). For general SDE
models, a change of variables allows us to locally consider an SDE with
constant noise [38] and, in turn, any time-homogeneous SDE with constant
noise,

(2.25) $dX(t)=f(X_{t})dt+\sigma dW(t),$

is of course readily linearized into an $\mathcal{OU}$-process.

###  2.2. Exact sampling

The $\mathcal{OU}$-process (2.17)–(2.18) is a Gaussian continuous process with
mean and covariance given by

(2.26) $\displaystyle\operatorname{\mathbb{E}}[X_{t}]$
$\displaystyle=\frac{k}{\mu}(1-e^{-\mu t})+e^{-\mu t}X_{0},$ (2.27)
$\displaystyle\operatorname{\mathbb{C}ov}[X_{t},X_{s}]$
$\displaystyle=\frac{\sigma^{2}}{2\mu}(e^{-\mu|t-s|}-e^{-\mu(t+s)}),$

such that the stationary distribution is
$X_{\infty}\sim\mathcal{N}(k/\mu,\sigma^{2}/(2\mu))$.

Numerical simulation procedures are typically based on Euler-type
discretization methods. Assume for convenience a fixed numerical time step
$\Delta t$ and compute

(2.28) $\displaystyle X_{n+1}$ $\displaystyle=X_{n}+(k-\mu X_{n})\Delta
t+\sigma\Delta W_{n},$

with $\Delta W_{n}$ being i.i.d. normally distributed numbers of zero mean and
variance $\Delta t$. The numerical trajectory $(X_{n})$ is then an
approximation to the $\mathcal{OU}$-process $(X(t_{n}))$ at discrete times
$t_{n}=n\Delta t$ for $n=0,1,\ldots$ and is also a Gaussian vector with mean
and covariance

(2.29) $\displaystyle\operatorname{\mathbb{E}}[X_{n}]$
$\displaystyle=\frac{k}{\mu}+(1-\mu\Delta t)^{n}(X_{0}-\frac{k}{\mu}),$ (2.30)
$\displaystyle\operatorname{\mathbb{C}ov}(X_{n},X_{n+p})$
$\displaystyle=\frac{\sigma^{2}\Delta t}{1-(1-\mu\Delta t)^{2}}(1-(1-\mu\Delta
t)^{2n})(1-\mu\Delta t)^{p}.$

Actually, (2.28) forms an $AR(1)$-sequence [19]. Comparing (2.26)–(2.27) with
(2.29)–(2.30) we readily find the following useful result.

######  Proposition 2.1 (Exact $\mathcal{OU}$-samples and backward analysis).

The discrete process $(X_{n})$ given by the explicit Euler method (2.28)
follows the law of the discrete samples of an exact $\mathcal{OU}$-process
with perturbed parameters $(k_{\Delta t},\mu_{\Delta t},\sigma_{\Delta t})$,
where

(2.34) $\displaystyle\left.\begin{array}[]{rl}1-\mu\Delta t&=\exp{(-\Delta
t\mu_{\Delta t})}\\\ \frac{k}{\mu}&=\frac{k_{\Delta t}}{\mu_{\Delta t}}\\\
\frac{\sigma^{2}}{2\mu-\mu^{2}\Delta t}&=\frac{\sigma_{\Delta
t}^{2}}{2\mu_{\Delta t}}\end{array}\right\\}$

This system of equations can be solved explicitly provided $\Delta
t<\mu^{-1}$,

(2.38) $\displaystyle\left.\begin{array}[]{rl}k_{\Delta
t}&=-\frac{k\log{(1-\mu\Delta t)}}{\mu\Delta t}\\\ \mu_{\Delta
t}&=-\frac{\log{(1-\mu\Delta t)}}{\Delta t}\\\ \sigma_{\Delta
t}^{2}&=-2\sigma^{2}\frac{\log{(1-\mu\Delta t)}}{2\mu\Delta t-\mu^{2}\Delta
t^{2}}\end{array}\right\\}$

These perturbed parameters are all first order perturbations in $\Delta t$ of
the exact parameters.

The inverse of (2.38) is that

(2.42) $\displaystyle\left.\begin{array}[]{rl}k&=\mu\frac{k_{\Delta
t}}{\mu_{\Delta t}}\\\ \mu&=\frac{1-\exp(-\Delta t\mu_{\Delta t})}{\Delta
t}\\\ \sigma^{2}&=\sigma^{2}_{\Delta t}\,\frac{2-\Delta t\mu}{2\mu_{\Delta
t}}\mu\end{array}\right\\}$

It follows that, for _given_ parameters $[k_{\Delta t},\mu_{\Delta
t},\sigma_{\Delta t}]$, a forward Euler simulation using the new set of
parameters $[k,\mu,\sigma]$ defined by (2.42) will produce a sample trajectory
obeying the law of the $\mathcal{OU}$-process with the given parameters
exactly.

##  3\. Bayesian filtering with full information

We first briefly consider in this section the behavior of the Bayesian
posterior as a function of the prior and of data in the sense of full state
measurements. Hence we suppose that sampled process data
$D_{N}=(d_{i})_{i=0}^{N}$ is available at some fixed time step $h$, that is,
$d_{i}=X(t_{i})$ with $t_{i}=ih$ and $X(\cdot)$ an $\mathcal{OU}$-process with
$X_{0}$ for simplicity considered a known sure number.

It will be convenient to consider the following reparametrization of the
$\mathcal{OU}$-process

(3.4)
$\displaystyle\left.\begin{array}[]{rl}\alpha&=\frac{k}{\mu}\in\mathbf{R}\\\
\beta&=e^{-\mu h}\in(0,1)\\\
\gamma&=\frac{2\mu}{\sigma^{2}}\cdot\frac{1}{1-e^{-2\mu
h}}\in\mathbf{R}_{+}\\\ \end{array}\right\\}$

and we put $u\equiv[\alpha,\beta,\gamma]$ for brevity, also defining
$u_{0}\equiv[\alpha_{0},\beta_{0},\gamma_{0}]$, i.e., the true parameters
generating the data.

Under this reparametrization we readily find from Proposition 2.1 that the
_exact_ discrete process can be written in the $AR(1)$-sequence form

(3.5) $\displaystyle X_{i+1}$ $\displaystyle=\beta
X_{i}+\alpha(1-\beta)+\gamma^{-1/2}\xi_{i},$

where the $\xi_{i}\sim\mathcal{N}(0,1)$ are independent. Writing
$\delta:=\alpha(1-\beta)$, explicit least squares estimators for the
parameters can be found by solving

(3.6)
$\displaystyle\min_{\delta,\beta}\|A\times[\delta,\beta]^{T}-b\|_{2}^{2}$
$\displaystyle=:\min_{\delta,\beta}\left\|\begin{bmatrix}\mathbf{1}&[d_{i}]_{i=0}^{N-1}\end{bmatrix}\times[\delta,\beta]^{T}-[d_{i+1}]_{i=0}^{N-1}\right\|_{2}^{2},$
where the brackets over the data form column vectors. The residual of this
solution implies the corresponding estimator for $\gamma$, (3.7)
$\displaystyle\hat{\gamma}^{-1}$
$\displaystyle=(N-2)^{-1}\|A\times[\hat{\delta},\hat{\beta}]^{T}-b\|_{2}^{2}.$

By the Gaussian character of the $\mathcal{OU}$-process these estimators
coincide with the maximum likelihood estimators.

The Bayesian convergence to the true parameters as $N\to\infty$ can be
characterized by either the Bernstein von Mises theorem (BvM), for $h$ fixed,
or by contrast functions convergence for a fixed time window $t\in[0,T]$. BvM
states that the Bayesian posterior converges to a normal distribution centered
at the maximum likelihood estimate with the inverse Fisher information matrix
(FIM) as covariance [5, 31, 6]. Similarly, the contrast functions approach
considers the convergence of the approximation of discretized observations,
approaching the same normal distribution for $N$ large enough [15, 18, 25]. In
the present case the FIM can be determined explicitly using the log-likelihood
$\log\mathcal{L}(X_{i+1}|X_{i};\;u)$ induced by (3.5) and the definition

(3.8) $\displaystyle\Sigma$
$\displaystyle\equiv-\operatorname{\mathbb{E}}\left[\frac{\partial^{2}}{\partial
u^{2}}\log\mathcal{L}\right](u_{0}).$ The result is that (3.9)
$\displaystyle\Sigma$
$\displaystyle=\operatorname{\text{Diag}}\left((1-\beta_{0})^{2}\gamma_{0},1/(1-\beta_{0}^{2}),1/(2\gamma_{0}^{2})\right).$

######  Proposition 3.1 (BvM Theorem for (3.5)).

Provided that the prior has a continuous and positive density in an open
neighborhood of $(\alpha,\beta,\gamma)$, as $N\to\infty$ we have that the ML-
estimators (3.6)–(3.7) converge to the true values as

(3.10)
$\displaystyle\sqrt{N}(\hat{u}-u_{0})\xrightarrow{d}\mathcal{N}(0,\Sigma^{-1}),$

(convergence in the sense of distribution) with $\Sigma$ defined by (3.9).
Alternatively, (3.10) dictates the asymptotic convergence of the mean of the
Bayesian posterior.

To summarize, the asymptotic variances of all parameters are independent of
$\alpha_{0}$. The noise term $\gamma_{0}$ mainly affects the convergence of
$\alpha_{0}$ (and $\gamma_{0}$ in an absolute sense). Finally, an increase of
$\beta_{0}$ implies a faster convergence towards $\beta_{0}$ at the cost of a
slower convergence towards $\alpha_{0}$. Under the interpretation of an SIS-
model (2.22)–(2.23), the parameter $\alpha_{0}$ is in one-to-one
correspondence with the basic reproduction number $R_{0}$, hence its central
interest here.

The asymptotic nature of both the BvM Theorem and the contrast function
convergence is a poor match in epidemiological situations where data is often
scarce and poorly informative at the time scale over which the parameters can
be considered static. This motivates our interest in also the pre-asymptotic
regime of the Bayesian posterior. We therefore consider prior densities of the
specific form

(3.11)
$\pi(u)\propto\gamma^{r}\exp{\left(-\frac{\gamma}{2}P(\alpha,\beta)\right)},$

where, for integrability, $P$ is to be a polynomial of degree $2$ in $\alpha$
and nonnegative for $\beta\in(0,1)$. This choice has the convenient property
that the posterior measure after $N$ observations is

(3.12) $\displaystyle\Pi^{(N)}(u)$
$\displaystyle\propto\gamma^{\frac{N}{2}+r}\exp{\left(-\frac{\gamma}{2}(Q_{N}+P)(\alpha,\beta)\right)},$
where (3.13) $\displaystyle Q_{N}(\alpha,\beta)$
$\displaystyle\equiv\sum_{i=0}^{N-1}(d_{i+1}-\beta
d_{i}-\alpha(1-\beta))^{2},$

and where $d_{i}=X(t_{i})$ is the observation at time $t_{i}$.

######  Remark.

To include also the case of flat priors, while avoiding technicalities for
nonintegrable densities, we note that, since the value of a constant prior has
no influence on the posterior, we can still determine a posterior after a
single initial observation, and then use this posterior as a prior for the
rest of the observations.

The following result examines the convergence of the Bayesian posterior.

######  Theorem 3.2 (Convergence of log-likelihood).

Consider the polynomial $q_{N}$:

(3.14) $q_{N}(\alpha,\beta)\equiv
N^{-1}Q_{N}(\alpha,\beta)=N^{-1}\sum_{i=0}^{N-1}(d_{i+1}-\beta
d_{i}-\alpha(1-\beta))^{2}.$

Then as $N\to\infty$ we have the almost everywhere uniform convergence on
every compact of $q_{N}$ towards the function $f$,

(3.15)
$f(\alpha,\beta)=(1-\beta)^{2}\left[\frac{1}{\gamma_{0}(1-\beta_{0}^{2})}+(\alpha_{0}-\alpha)^{2}\right]+\frac{2\beta}{\gamma_{0}(1+\beta_{0})}.$

###### Proof.

For a general point $(\alpha,\beta)$, we rely on the ergodic theory of Markov
chains [2, p. 472] to get that $N^{-1}\sum g(d_{i},d_{i+1})$ converges almost
surely towards $\operatorname{\mathbb{E}}[g(d_{i},d_{i+1})]$ for $g$
integrable against the stationary measure. Using (2.26)–(2.27) this implies
the limits:

(3.16) $\displaystyle N^{-1}\sum_{i=0}^{N-1}d_{i+1}-\beta d_{i}$
$\displaystyle\longrightarrow\alpha_{0}(1-\beta),$ (3.17) $\displaystyle
N^{-1}\sum_{i=0}^{N-1}(d_{i+1}-\beta d_{i})^{2}$
$\displaystyle\longrightarrow(\alpha_{0}(1-\beta))^{2}+\frac{1+\beta^{2}}{\gamma_{0}(1-\beta_{0}^{2})}-\frac{2\beta\beta_{0}}{\gamma_{0}(1-\beta_{0}^{2})}.$

Combined, we obtain the claimed limit in a pointwise sense. As each $(q_{N})$
is a polynomial in $(\alpha,\beta)$ of bounded degree, pointwise convergence
is equivalent to convergence of the coefficients, and hence implies the
convergence on all compact sets. ∎

An obvious extension is to consider measurements polluted by noise, say,
$d_{i}=X_{i}+\eta_{i}$, for i.i.d. $\eta_{i}\sim\mathcal{N}(0,\eta)$ and some
variance $\eta$. The posterior so obtained is readily computed via a recursive
Kalman filter but does not have a simple analytic form. To first order in
$\eta$, however, (3.11)–(3.12) still holds provided $\gamma$ is replaced with
$[1/\gamma+\eta(1+\beta^{2})]^{-1}$ and $Q_{N}$ in (3.13) is replaced with

(3.18) $\displaystyle Q^{\prime}_{N}$
$\displaystyle\equiv\sum_{i=0}^{N-1}\left[d_{i+1}-\beta
d_{i}-\alpha(1-\beta)+\eta\beta\gamma(d_{i}-\beta
d_{i-1}-\alpha(1-\beta))\right]^{2},$

where the right $\eta\beta\gamma$-term is skipped when $i=0$. Intuitively, the
first order effect of noise in data is to broaden the posterior with the
variance of this noise. We next proceed to investigate more severe truncations
of the measurements from the epidemiological process.

##  4\. Bayesian filtering of surveillance data

In the previous section we considered full process data to be available
without any extrinsic noise. This models the best possible Bayesian setup but
is also clearly unrealistic in epidemics. As a model of a more challenging
situation we consider in this section the recorded data to be some (possibly
stochastic) function of the $\mathcal{OU}$-process $X(\cdot)$, which itself is
considered a latent variable. We are initially concerned with _binary data_ of
the form $D_{N}=(d_{i})_{i=0}^{N}$ where $d_{i}=Y_{i}=\mathbf{1}_{X(t_{i})\geq
c}$, again over a uniform grid in time $t_{i}=ih$, and for a known _filter_
cut-off value $c$ (see also the related setup in [44]). That is, the
epidemiological interpretation is that the data is considered to be time-
discrete information about whether or not the prevalence $X(\cdot)$ of the
population is above or below a certain known threshold $c$. With the
prevalence process hidden one is forced to estimate it simultaneously with any
parameter estimates. Using the $AR(1)$-form (3.5) we have that the stationary
measure for the $(p+1)$ steps $(X_{i},\ldots,X_{i+p})$ is Gaussian
$\mathcal{N}(\alpha,\Sigma)$ with $\Sigma$ given by (for $p\geq 0$)

(4.1)
$\displaystyle\Sigma=\Sigma(\beta,\gamma)=\frac{\gamma^{-1}}{1-\beta^{2}}\begin{bmatrix}1&\beta&\beta^{2}&\cdots&\beta^{p}\\\
\beta&1&\beta&\cdots&\beta^{p-1}\\\ \beta^{2}&\ddots&\ddots&\ddots&\cdots\\\
\vdots&\ddots&\beta&1&\beta\\\
\beta^{p}&\cdots&\beta^{2}&\beta&1\end{bmatrix}.$

In order to be able to conduct an analysis where information is obtained only
from the observable $Y(\cdot)$, we take (4.1) as a motivation for the
following:

######  Assumption 4.1 ($p$-step Markovian stationary assumption).

We say that we work under the stationary assumption whenever we assume that
the law of the latent variable $Y(\cdot)$ can be directly inferred from the
stationary $p$-step law (4.1).

We stress that it is known that this type of clipped Gaussian processes are
not $p$-step Markov for any $p$ [40]. Assumption 4.1 rather serve as an
approximation where we approximately model $Y(\cdot)$ as if it was $p$-step
Markov with law deduced from (4.1). The pseudo-likelihood for the filtered
variable $Y$ then becomes

(4.2) $\displaystyle\mathcal{L}_{N}(u)$
$\displaystyle=\prod_{i=0}^{N-p}\varphi_{u}(Y_{i},\ldots,Y_{i+p})$

where $\varphi_{u}(e)$ denotes the probability for a Gaussian stationary
filtered process with parameters $u$ to be $e\in\\{0,1\\}^{p+1}$. We also
define the pseudo-potential

(4.3)
$\Phi_{N}^{u}\equiv\log\mathcal{L}_{N}(u)=\sum_{i=0}^{N-p}\log\varphi_{u}(Y_{i},\ldots,Y_{i+p}).$

We show below in §4.1 that the stationary assumption allows for converging
posterior estimates for the parameter $\beta$, but leaves any prior density
unchanged over a certain curve in $(\alpha,\gamma)$. In §4.2 we briefly
consider the filter cutoff value $c$ (the sensitivity) to be uncertain and
straightforwardly show that any prior density on $c$ is unaffected by data. We
next replace the sharp deterministic filter by a stochastic filter
implementing a sigmoidal response function and sharpen our results in this
more general setting in §4.3. Finally, in §4.4 we consider slightly more
informative measurements consisting of a finite discrete response and we show
that this resolves the singularity issues associated with purely binary
measurements.

###  4.1. Binary measurements

We first show that under the $p$-fold stationary assumption one can estimate
the correlation term $\beta$ rather well, but increasing the gap between the
filter threshold $c$ and the mean $\alpha$ has an effect on the likelihood
which is indistinguishable in law from increasing the noise $\gamma$.

We start with a technical lemma.

######  Lemma 4.1 (Equality in law).

For $u$ in any admissible set of parameters, let $p_{u}$ be the law of
$(Y_{0},\ldots,Y_{p})$ under Assumption 4.1 corresponding to
$\mathcal{N}(\alpha_{u},\Sigma_{u})$. Then $p_{u}=p_{w}$ if and only if
$\beta_{u}=\beta_{w}$ and
$\sqrt{\gamma_{u}}(c-\alpha_{u})=\sqrt{\gamma_{w}}(c-\alpha_{w})$.

###### Proof.

$(\Longrightarrow)$ Suppose $p_{u}=p_{w}$. Consider
$(X_{0}^{u},\ldots,X_{p}^{u})$ a Gaussian vector with law
$\mathcal{N}(\alpha_{u},\Sigma_{u})$. As $\mathbf{P}(X_{0}^{u}\leq
c)=\sum_{e\in\\{0,1\\}^{p}}\varphi(0,e)$, we deduce that
$\mathbf{P}(X_{0}^{u}\leq c)=\mathbf{P}(X_{0}^{w}\leq c)$. This can be written
as an equality between standard cumulative distribution functions of
Gaussians, since $\sqrt{\gamma_{u}(1-\beta_{u}^{2})}(X_{0}^{u}-\alpha_{u})$
and $\sqrt{\gamma_{w}(1-\beta_{w}^{2})}(X_{0}^{w}-\alpha_{w})$ are standard
Gaussians. This translates into:

(4.4)
$\sqrt{\gamma_{u}(1-\beta_{u}^{2})}(c-\alpha_{u})=\sqrt{\gamma_{w}(1-\beta_{w}^{2})}(c-\alpha_{w})$

We also have $\mathbf{P}(X_{0}^{u},X_{1}^{u}\leq
c)=\mathbf{P}(X_{0}^{w},X_{1}^{w}\leq c)$, and this implies, with
$A=\sqrt{\gamma_{u}(1-\beta_{u}^{2})}(c-\alpha_{u})$,

(4.5)
$\frac{1}{\sqrt{1-\beta_{u}^{2}}}\int\limits_{(-\infty,A]^{2}}\exp\left(-\frac{x^{2}+y^{2}-2\beta_{u}xy}{\sqrt{1-\beta_{u}^{2}}}\right)\,dxdy\\\
=\frac{1}{\sqrt{1-\beta_{w}^{2}}}\int\limits_{(-\infty,A]^{2}}\exp\left(-\frac{x^{2}+y^{2}-2\beta_{w}xy}{\sqrt{1-\beta_{w}^{2}}}\right)\,dxdy.$

This last expression shows that the map
$\beta_{u}\mapsto\mathbf{P}(X_{0}^{u},X_{1}^{u}\leq c)$ is locally analytic.
As a corollary of Slepian’s lemma [21], this map is also increasing and hence
strictly increasing. It follows that the last equality implies
$\beta_{u}=\beta_{w}$.

$(\Longleftarrow)$ The law $p_{u}$ is uniquely determined by its $2^{p+1}$
values. All these values are of the type

(4.6) $I(Q;\;u)=\int_{Q(A)}(\det
M(\beta_{u}))^{-(p+1)/2}\times\exp\left(-xM(\beta_{u})x^{T}\right)\,dx,$

where $Q(A)$ is a product of $p+1$ intervals, each being either $(-\infty,A]$
or $[A,+\infty)$. And so, if $\beta_{u}=\beta_{w}$ and
$\sqrt{\gamma_{u}}(c-\alpha_{u})=\sqrt{\gamma_{w}}(c-\alpha_{w})$, the laws
$p_{u}$ and $p_{w}$ are indeed equal. ∎

######  Theorem 4.2 (Non-identifiability).

Assume $p\geq 1$. For any integrable prior $\pi$, non-zero on $\mathcal{E}$,
and for $f$ bounded, as $N\to\infty$,

(4.7) $\int f\mathcal{L}_{N}d\pi\rightarrow\int_{\mathcal{E}}fd\tilde{\pi},$

where
$\mathcal{E}\equiv\\{u;\;\sqrt{\gamma}(c-\alpha)=\sqrt{\gamma}_{0}(c-\alpha_{0}),\,\beta=\beta_{0}\\}$
and $\tilde{\pi}$ is the prior restricted to the set $\mathcal{E}$.

###### Proof.

At first, thanks to Lemma 4.1, the law $\varphi(Y_{i}\ldots,Y_{i+p})$ is not
characterized by $u_{0}=[\alpha_{0},\beta_{0},\gamma_{0}]$, but rather by
$[\sqrt{\gamma_{0}}(c-\alpha_{0}),\beta_{0}]$. As the $Y_{i}$’s are binary,
the quantity $\varphi(Y_{i},\ldots,Y_{i+p})$ can only take $2^{p+1}$ values,
which we denote by $\varphi_{e}$ for $e\in\\{0,1\\}^{p+1}$. This means that
the log-likelihood under Assumption 4.1 can be written as a finite sum:

(4.8)
$\log\mathcal{L}_{N}=\sum\limits_{e\in\\{0,1\\}^{p+1}}N(e)\log\varphi_{e}.$

where

(4.9) $N(e)=\sharp\\{i\in\\{0,N-p\\};\;(Y_{i},\ldots,Y_{i+p})=e\\}$

Now, as the discrete $\mathcal{OU}$-process is an $AR(1)$-sequence [19], we
may infer the convergence of $N(e)/N$ towards a stationary value:

(4.10)
$\frac{\log\mathcal{L}_{N}}{N}\rightarrow\sum\limits_{e\in\\{0,1\\}^{p+1}}\varphi_{e}^{u_{0}}\log\varphi_{e}^{u},$

where $\varphi_{e}^{u}$ denotes the law of a stationary filtered process with
parameter $u$. We note that

(4.11)
$\lim_{N\to\infty}\frac{\log\mathcal{L}_{N}}{N}=\sum\limits_{e\in\\{0,1\\}^{p+1}}\varphi_{e}^{u_{0}}\log\frac{\varphi_{e}^{u}}{\varphi_{e}^{u_{0}}}+\sum\limits_{e\in\\{0,1\\}^{p+1}}\varphi_{e}^{u_{0}}\log\varphi_{e}^{u_{0}}.$

Up to a constant this limit is the negative of the Kullback-Leibler divergence
between the laws $\varphi^{u}$ and $\varphi^{u_{0}}$, which vanishes if and
only if the two distributions are equal. We already know that the laws of
$\varphi^{u}$ and $\varphi^{u_{0}}$ are equal if and only if $\beta=\beta_{0}$
and $\sqrt{\gamma}(c-\alpha)=\sqrt{\gamma_{0}}(c-\alpha_{0})$, i.e., if
$u\in\mathcal{E}$. Consider the associated relation:

(4.12) $u=(\alpha,\beta,\gamma)\sim
u^{\prime}=(\alpha^{\prime},\beta^{\prime},\gamma^{\prime})\text{ iff
}\beta=\beta^{\prime}\text{ and
}\sqrt{\gamma}(c-\alpha)=\sqrt{\gamma^{\prime}}(c-\alpha^{\prime}).$

It is straightforward to show that this relation is reflexive, symmetric, and
transitive, thus forming an equivalence relation. Thanks to Lemma 4.1, we know
that the map $u\mapsto\mathcal{L}_{N}(u)$ factorizes through the equivalence
relation to a map $w\mapsto\tilde{\mathcal{L}}_{N}(w)$ where $w$ ranges over
the different equivalence classes of $u$.

Lemma 4.1 also allows us to state that the map $\tilde{\mathcal{L}_{N}}$ is
injective. Consider $w_{0}$ in the equivalence class of $u_{0}$. Then for any
neighborhood $W$ of $w_{0}$, there is a $\delta>0$ such that for $N$ large
enough, $w\not\in W$ implies

(4.13)
$\frac{\log{\tilde{\mathcal{L}}_{N}(w_{0})}}{N}\geq\frac{\log{\tilde{\mathcal{L}}_{N}(w)}}{N}+\delta,$

which means that the posterior measure of $W$ converges towards 1\. Since the
factorized map $\tilde{\mathcal{L}}_{N}\rightarrow\delta_{w_{0}}$, the theorem
follows. ∎

The same argument implies that, if we consider the case $p=0$, the posterior
will be even more degenerate. For $p=0$ one gets the convergence

(4.14) $\int
f\mathcal{L}_{N}d\pi\rightarrow\int_{\mathcal{E}^{\prime}}fd\tilde{\pi},$

where now $\tilde{\pi}$ is the prior restricted to $\mathcal{E}^{\prime}$:

(4.15)
$\mathcal{E^{\prime}}\equiv\\{u;\;\sqrt{\gamma_{u}(1-\beta_{u}^{2})}(c-\alpha_{u})=\sqrt{\gamma_{0}(1-\beta_{0}^{2})}(c-\alpha_{0})\\}.$

This case is relevant whenever we consider large gaps of time in between the
measurements such that they can be considered practically independent.

At this point, let us remind that the quantities $\mathcal{L}_{N}(u)$ and
$\Phi_{N}^{p}(u)$ are defined in equations (4.2)–(4.3).

######  Theorem 4.3 (Rate of convergence of pseudo-potential).

Under the $p$-Markovian stationary Assumption 4.1 there exists a mapping
$f_{p}$ and a constant $\sigma_{p}>0$ such that, as $N\to\infty$, we have the
convergence in law:

(4.16)
$\sqrt{N}\bigg{(}\frac{\Phi_{N}^{p}(u)}{N}-f_{p}(u)\bigg{)}\rightarrow\mathcal{N}(0,\sigma_{p}).$

###### Proof.

This is a direct consequence of the Central Limit Theorem, with $f_{p}(u)$
being the mean value of $\Phi_{N}^{p}(u)$ and $\sigma_{p}$ a nonnegative
constant. ∎

###  4.2. Propagation of filter uncertainty

A relevant variation of the theme is to consider the parameter $c$ (the “test
sensitivity”) an uncertain parameter. Mathematically, this means considering
the likelihood $\mathcal{L}_{N}$ a function of $(\alpha,\beta,\gamma,c)$, and
priors and posteriors depending also on $c$. However, the following result
shows that this setup will not produce any more information about the
parameters.

######  Theorem 4.4 (Translation of filter uncertainty).

Consider a prior of the form $\mu(\alpha)\pi(c-\alpha,\beta,\gamma)$ and the
pseudo-likelihood in §4.1. Then the resulting posterior measure will be of the
form $\mu(\alpha)\Pi(c-\alpha,\beta,\gamma)$, where, for fixed $c$,
$\Pi(c-\cdot,\cdot,\cdot)$ is the pseudo-posterior from the prior
$\pi(c-\cdot,\cdot,\cdot)$.

###### Proof.

This result is immediate once one realizes that the pseudo-likelihood has the
property that, for any $t\in\mathbf{R}$,

(4.17)
$\mathcal{L}_{N}(\alpha,\beta,\gamma,c)=\mathcal{L}_{N}(\alpha+t,\beta,\gamma,c+t).$

∎

In other words, the only information we can infer on the parameters
$(c,\alpha)$ is the gap $c-\alpha$. Any uncertainty of $c$ can be understood
as an uncertainty on $\alpha$, and vice-versa.

###  4.3. Non-perfect binary measurements

Rather than a sharp cut-off value $c$, most environmental sampling methods
obey some kind of _sensitivity response_ , e.g., of sigmoid character, with a
quick rise in the detection probability as one progresses through some
threshold region. Examples here could include sampling and subsequent analysis
of sewage water or animal droppings, but this would also be a relevant model
in the case of statistical regression estimates using data obtained via self-
reporting smartphones applications.

A general ansatz to capturing this situation is to consider

(4.18) $Y_{i}\sim\mathcal{B}(s(X_{i})),$

where $s$ is a map from $\mathbf{R}$ to $[0,1]$ and $\mathcal{B}$ denotes the
Bernoulli law. Typically, the map $s$ is sigmoidal with a gradient around the
threshold $c$ which depends on the sensitivity and specificity of the test. We
naturally ask that the efficiency of the filter does not depend on its
previous use, i.e., that $(Y_{i}|X_{i})_{i}$ is an independent family. To
construct such an object, one may consider an i.i.d. family of uniformly
distributed variables $(\xi_{i})$ on $(0,1)$, independent from $(X_{i})$, and
set

(4.19) $Y_{i}=\mathbf{1}_{\xi_{i}\leq s(X_{i})}.$

This object clearly fulfills all properties mentioned and is general enough to
capture also quite specific situations.

We will now establish our result in this more general setting. Consider a
bounded and continuous map $f:\mathbf{R}^{p+1}\rightarrow\mathbf{R}$. We
observe the $\mathcal{OU}$-process $X(t)$ through the map $f$, so the
observations are $Y_{i}=f(X_{i},\ldots,X_{i+p})$ where as before
$X_{i}=X(t_{i})$. To define a pseudo-likelihood, we reason as if the observed
data can be regarded as stationary. We thus denote by $f_{u}$ the law of a
_stationary_ $\mathcal{OU}$-process with parameters $u$ filtered via the
measurement map $f$ and we remind ourselves that the latent process $X(t)$ is
not necessarily stationary.

We consider a likelihood of the kind

(4.20) $\mathcal{L}_{N}\propto\exp\left(\sum_{i=0}^{N-p}\log
f_{u}(Y_{i})\right),$

normalized to mass one, and we denote the potential by

(4.21) $g_{N}(u):=N^{-1}\sum\limits_{i=0}^{N-p}\log f_{u}(Y_{i}).$

This setup simply corresponds to observing the data portion
$(X_{i},\ldots,X_{i+p})$ via the filter $f$ and, assuming stationarity,
building a pseudo-likelihood. Unfortunately this does not directly include the
intended case of a sigmoid filter, as this would rather involve measurements
$h(\xi_{i},X_{i},\ldots,\xi_{i+p},X_{i+p})$ with $(\xi_{i})$ an i.i.d.
sequence independent from $(X_{i})$ and uniform on $(0,1)$. However, once the
proposed case is examined, one can get to the latter case by exchanging the
order of integration using the Fubini property, i.e., studying the behavior of
$h(t_{0},X_{i},\ldots,t_{p},X_{i+p})$ and then integrate over
$(t_{0},\ldots,t_{p})$. For example, consider the case $p=1$. Then the
intended map $h$ would be

(4.22) $h(t_{0},x_{0},t_{1},x_{1})=(\mathbf{1}_{t_{0}\leq
s(x_{0})},\mathbf{1}_{t_{1}\leq s(x_{1})})$

For a fixed $u$, $f_{u}$ can take at most $4$ values, e.g.,

(4.23) $f_{u}(1,1)=\operatorname{\mathbb{E}}[s(Z_{0})s(Z_{1})],$

where $(Z_{0},Z_{1})$ is a stationary $\mathcal{OU}$-process of parameter $u$,
and similarly for $f_{u}(0,1)$, $f_{u}(1,0)$, and $f_{u}(0,0)$. As this
technique would render the proof lengthy, we decided to put it aside.

In order to state our result, we need to specify some minimal set of
regularity conditions. Convergence of the potential is required as well as
definiteness in the Kullback-Leibler divergence ($\operatorname{D_{KL}}$).
Additionally, we shall also require a separation condition in the large data
limit.

######  Assumption 4.2 (Regularity).

We assume the following specifics:

  1. (1)

The potential $g_{N}$ converges uniformly on the compacts of $u$.

  2. (2)

$\operatorname{D_{KL}}(f_{u_{0}},f_{u})=0\Longrightarrow u=u_{0}$.

  3. (3)

There is a $\delta>0$ and a compact neighborhood $K$ of $u_{0}$ such that, for
$N$ large enough and for $u\in K^{\complement}$,
$g_{N}(u_{0})>g_{N}(u)+\delta$.

######  Theorem 4.5 (Weak convergence).

Under Assumption 4.2, we have the weak convergence

(4.24) $\mathcal{L}_{N}\rightarrow\delta_{u_{0}}.$

This theorem can be adapted to cover other situations:

  * •

If $\operatorname{D_{KL}}(f_{u_{0}},f_{u})=0$ does not imply $u=u_{0}$, one
can try to factorize the map through an equivalence relation (as in the proof
of Theorem 4.2) to get the convergence towards the indicative function of the
set $\\{u;\;\operatorname{D_{KL}}(f_{u_{0}},f_{u})=0\\}$.

  * •

If we only want to consider parameters within a subset $U$ of the set of
parameters, one can always consider the assumptions restricted to the set $U$,
and then one would have the convergence only for priors with support in $U$.

###### Proof.

We divide the proof in two steps as follows.

#### Step 1

Let us call $g$ the limit of $(g_{N})$, which by assumption is not random. The
same argument as in the proof of Theorem 3.2, the ergodicity of the process
$(X_{t})$, allows us to state that the limit function must be
$\operatorname{\mathbb{E}}[\log f_{u}]$ (up to a constant, the limit is
actually $-\operatorname{D_{KL}}(f_{u_{0}},f_{u})$). As adding a constant is
equivalent to multiplying all posteriors by a nonnegative constant, we may
assume the limit function to be equal to
$u\mapsto-\operatorname{D_{KL}}(f_{u_{0}},f_{u})$.

#### Step 2

As $\mathcal{L}_{N}$ is proportional to $\exp(g_{N})$, the assumptions allow
us to deduce that the mass of any neighborhood of $u_{0}$ converges towards
$1$. This implies the weak convergence of the posteriors towards the Dirac
$\delta_{u_{0}}$. ∎

###  4.4. Non-binary measurements

As a final variation on the theme we now show how Theorem 4.5 can be applied
in such a way as to overcome the issues with the non-identifiability of
Theorem 4.2. The idea is that the test with one filter $c$ evidently at best
gives us a curve containing the parameter $u_{0}$, namely the one satisfying
$\beta=\beta_{0}$, and
$\sqrt{\gamma}(c-\alpha)=\sqrt{\gamma_{0}}(c-\alpha_{0})$. Suppose instead
that we have two kinds of tests, one with cut-off $c_{1}$ and one with cut-off
$c_{2}$. Since the intersection of the two curves implied by the respective
filters $c_{1}$ and $c_{2}$ is exactly the point $u_{0}$, it is natural to
assume that two filters are enough to get the convergence of the posterior
towards $\delta_{u_{0}}$. We show that this is indeed the case.

######  Theorem 4.6 (Trinary filter).

Let $c_{1}<c_{2}$ and consider the filtered values $Y_{i}$ to be

(4.25) $Y_{i}=\mathbf{1}_{X_{i}>c_{1}}+\mathbf{1}_{X_{i}>c_{2}},$

and the associated pseudo-likelihood

(4.26) $\mathcal{L}_{N}=\prod_{i=0}^{N-p}\varphi(Z_{i},\ldots,Z_{i+p}),$

where $\varphi$ is the multivariate cumulative distribution function of a
Gaussian $\mathcal{N}(\alpha,\Sigma)$ (4.1). Then the sequence of pseudo-
likelihoods converges weakly:

(4.27) $\mathcal{L}_{N}\rightarrow\delta_{u_{0}}.$

###### Proof.

It is sufficient to check that the sequence $(Y_{i})$ verifies Assumption 4.2.
Considering $g_{N}:=N^{-1}\log{\mathcal{L}_{N}}$, one has

(4.28) $g_{N}=\sum\limits_{e\in\\{0,1,2\\}^{p}}\frac{N_{e}}{N}\log{p_{u}(e)},$

where $N_{e}$ is the number of times $(Y_{i},\ldots,Y_{i+p})$ matches $e$, and
$p_{u}$ is the cumulative distribution function of the underlying Gaussian
process $\mathcal{N}(\alpha,\Sigma)$. From this we get directly conditions (1)
and (3) of Assumption 4.2. For condition (2), if the Kullback-Leibler
divergence is zero, then we must have $p_{u}(e)=p_{u_{0}}(e)$ for all $e$.
This implies the equalities

(4.32) $\displaystyle\left.\begin{array}[]{rl}\beta&=\beta_{0}\\\
\sqrt{\gamma}(c_{1}-\alpha)&=\sqrt{\gamma_{0}}(c_{1}-\alpha_{0})\\\
\sqrt{\gamma}(c_{2}-\alpha)&=\sqrt{\gamma_{0}}(c_{2}-\alpha_{0})\end{array}\right\\}$

As $c_{1}\neq c_{2}$ we have $u=u_{0}$, and so the assumptions are fulfilled.
∎

##  5\. Illustrations

We devote this section to some illustrations of selected results from §§3 and
4. We shall do this in the intended epidemiological setting and thus no longer
assume the $\mathcal{OU}$-process, but rather the SIS- and
SIS${}_{\text{E}}$-models from §2 in the form of continuous-time Markov chains
over a discrete state-space. In §5.1 we investigate the precision of the
predicted posterior uncertainty under full state measurements, while in §5.2
we offer a demonstration of the singular behavior under filtered measurements.
Finally, in §5.3 we highlight the use of synthetic data when approaching more
realistic problems defined over a network.

The software for the numerical experiments is available for download via the
corresponding author’s web-page111Refer to the BISDE-code at
https://user.it.uu.se/$\sim$stefane/freeware.html.

###  5.1. Asymptotic uncertainty

We first consider the Bayesian uncertainty under accurate measurements and
take the continuous-time Markov chain version of the SIS-model (2.14) as an
example. Using the $\text{SIS}\leftrightarrow$ $\mathcal{OU}$ approximate
interpretation (2.22) we have from (3.4) the relations

(5.7)
$\displaystyle\left.\begin{array}[]{rl}\alpha_{\text{$\mathcal{OU}$}}&=\Sigma(1-R_{0}^{-1})\\\
\beta_{\text{$\mathcal{OU}$}}&=\exp(-\gamma h(R_{0}-1))\\\
\gamma_{\text{$\mathcal{OU}$}}&=\Sigma^{-1}R_{0}\cdot[1-\beta_{\text{$\mathcal{OU}$}}^{2}]^{-1}\end{array}\right\\}\quad\left.\begin{array}[]{rl}R_{0}&=1+\alpha_{\text{$\mathcal{OU}$}}(1-\beta_{\text{$\mathcal{OU}$}}^{2})\gamma_{\text{$\mathcal{OU}$}}\\\
\Sigma&=\alpha_{\text{$\mathcal{OU}$}}+[(1-\beta_{\text{$\mathcal{OU}$}}^{2})\gamma_{\text{$\mathcal{OU}$}}]^{-1}\\\
\gamma&=-\log\beta_{\text{$\mathcal{OU}$}}\cdot[h(R_{0}-1)]^{-1}\end{array}\right\\}$

where $\\{R_{0},\Sigma,\gamma\\}$ are the SIS-model parameters.

We generate synthetic data from the Markov chain as illustrated in Fig. 2.1
for $\Sigma=1000$ and with ranges of values $\gamma\in[0.1,1]$ and
$R_{0}\in(1,3.5]$. We let $I(0)=0.01\times\Sigma$ and sample exact values of
$I(t_{i})$ for $t_{i}=ih$, $i=1,\ldots,N$, and $h\equiv 1$. This corresponds
to 100 perfect samples in a closed population at a rate equivalent to between
one to one tenth the disease period unit ($=1/\gamma$).

We evaluate the posterior over a grid of values in the $(R_{0},\gamma)$-plane
by simply normalizing the likelihood for the Markov chain given the synthetic
data discussed previously. The likelihood of the Markov chain is formally
obtained by solving the associated master equation which, however, is
inconvenient to do except for small populations $\Sigma$. A more general
approach is via a local linear Gaussian approximation and a Kalman filter. Put
$I_{0}=I(0)$ and define

$\displaystyle I_{k+1}$ $\displaystyle=\left(1+\Delta
t\beta\Sigma^{-1}(\Sigma-I_{k})-\Delta t\gamma\right)I_{k}+w_{k},$ (5.8)
$\displaystyle w_{k}$ $\displaystyle\sim\mathcal{N}(0,[\Delta
t\beta\Sigma^{-1}(\Sigma-I_{k})+\Delta t\gamma]I_{k}),$

that is, this is the forward Euler discretization of the Langevin equations
approximating the Markov chain. The Kalman filter associated with (5.8)
computes a likelihood for each data point, albeit for a perturbed model. The
relative error in the Langevin approximation generally scales with the inverse
of the population size $\Sigma$, and can be expected to be rather small in the
present context (see [12, Ch. 11.3]). Further, Proposition 2.1 suggests
analyzing the Euler discretization via backward analysis as a parameter
perturbation, but unfortunately this is not generalizable to non-additive
noise [36]. For a resolved discretization, however, $\gamma\Delta t\ll 1$,
weak first order convergence can be expected under broad conditions. We use
the constant Kalman resolution $\Delta t=h/4$ and next focus on the estimation
error.

We have already evaluated the asymptotic covariance matrix under accurate data
in (3.10). Using the linear uncertainty transformation $Q^{\prime}\approx
JQJ^{T}$, where $J$ is the Jacobian of the parameter map (5.7), and where $Q$
is the (diagonal) covariance matrix in (3.10), we can estimate the posterior
variance

(5.9) $\displaystyle\text{Var}(R_{0})$
$\displaystyle\approx\frac{R_{0}^{3}}{\Sigma
N}\times\frac{1+\beta_{\text{$\mathcal{OU}$}}}{1-\beta_{\text{$\mathcal{OU}$}}}.$

A similar formula can be worked out for the variance of $\gamma$ as well,
although a bit more involved. For small enough $h$, the denominator
$1-\beta_{\text{$\mathcal{OU}$}}\sim\gamma h(R_{0}-1)$ in (5.9), and so the
_relative_ uncertainty in any consistent estimator of $R_{0}$ can be expected
to depend weakly on $R_{0}$ itself. This effect is seen for $R_{0}>1$ in Fig.
5.1 (top), where it can also be seen that the approximation (5.9) derived from
the $\mathcal{OU}$-approximation is somewhat optimistic. Similarly, we find
that the relative uncertainty of $R_{0}$ goes down with increasing values of
$\gamma$, or, which by (5.7) is the same thing, with decreasing correlation
$\beta_{\text{$\mathcal{OU}$}}$ (cf. Fig. 5.1, bottom).

Figure 5.1. Marginal posterior uncertainty ($\pm 2$ SD) for the SIS-model and
a range of parameters. Top: with $\gamma\equiv 1$ fixed, bottom: with
$R_{0}\equiv 1.5$ fixed. Crosses: MMSE-estimators $\hat{R}_{0}$ (i.e.,
posterior means), dotted: estimated uncertainty according to (5.9), red:
posterior width ($\pm 2$ SD).

###  5.2. Limits of convergence

The SIS-model investigated previously was dependent on two parameters only and
hence the singularity detected in §4 is not likely to be limiting any
convergence. While various model modifications naturally lead to additional
independent parameters, e.g., an extra transition $S\to I$ modeling external
infectious events, the most immediate modification is to simply consider the
population size $\Sigma$ uncertain. For instance, this is a possible setup for
inference relying on sewage water analysis, where the data is binary according
to whether the infectious substance is above or below some known threshold
value $c$, but the sewage uptake area is populated with an unknown number of
individuals.

We remark that this is a considerably challenging task, and although we are
able to demonstrate the sharpness of our negative results from §4 in this
setting, the fact that this problem can at all be approached is quite
remarkable.

As ground truth we use the same parameters as in the previous section, but
with $\gamma=1/10$ (time-1) and $R_{0}=1.5$ fixed, and we need to sample more,
$N=1000$ points equispaced with $h=1$. The data is then taken to be the
filtered sequence $Y_{i}=\mathbf{1}_{X_{i}\geq c}$ with $c=0.9\times
I_{\infty}=0.9\times\Sigma(1-R_{0}^{-1})$.

The measurement map is strongly nonlinear and so the Kalman filter needs to be
extended in some way. We took an immediate approach by simply discretizing the
state variable $I$ into $M=100$ Gaussian particles, distributed according to
the percentiles of the stationary measure for a given proposed set of
parameters. Each particle is evolved in time $[0,h]$ using steps of size
$\Delta t$ according to the Kalman filter (5.8), after which the prior
distribution is formed by aggregating the probability mass in the vicinity of
each particle. This yields the likelihood for a single data point after which
the posterior distribution for the state is obtained by setting selected
particles’ mass to zero (according to the data point) and rescaling
appropriately.

From Theorem 4.2 we have the singular curve
$\sqrt{\gamma_{\text{$\mathcal{OU}$}}}(c-\alpha_{\text{$\mathcal{OU}$}})=\mbox{constant}$,
which gets transferred via the map (5.7) into a surface in
$(R_{0},\Sigma,\gamma)$-space. After arbitrarily fixing $\gamma$ we thus
obtain a curve in the $(R_{0},\Sigma)$-plane. Since the SIS-to-$\mathcal{OU}$
map (2.22) is an approximation and, moreover, the likelihood is approximated
via a Kalman filter, this analytical curve can be expected to be a
perturbation of the observed numerical posterior level curves. As shown in
Fig. 5.2 the match is quite remarkable.

Figure 5.2. Left: (log-)posterior for the SIS-model under filtered data and
conditioned on the true value $\gamma=1/10$ (time-1). Also indicated is the
singular curve as predicted by theory passing through the true parameter
generating the data (circle). Right: marginal distribution for $R_{0}$
together with a normal fit (dashed).

###  5.3. Network epidemics

Epidemic models on networks can give rise to phenomena not observed in single-
node systems, e.g., the _rescue effect_ [22], where the infection is “rescued”
from extinction through the network structure. Here we consider both the SIS-
and the SIS${}_{\text{E}}$-model, respectively, cf. (2.7)–(2.11). We assume
these models at each node in a network implicitly defined by pre-recorded
movements of individuals between the nodes. As a concrete example this would
be an appropriate model for estimating disease parameters in a monitoring
program for bovine animals using cheap, but low-informative, tests collected
on a weekly basis for a subsample of _sentinel_ nodes.

The nodal model is replicated across 1,600 nodes, populated with 196,168
individuals, and the nodes are connected using 466,692 prescribed movements of
individuals over four years, see Figure 5.3. The system is not well stirred on
the aggregate level, but events occur frequently; the average # of events per
sample node and day $=0.20\,[0.19,0.21]$ with 50% credible interval (CrI).
This particular network was constructed by anonymizing a set of recorded
cattle movements and can be accessed through the publicly available R-package
SimInf [51, 11]. We extract model measurements from the same 100 randomly pre-
selected sentinel nodes every 7th day for a total of 4 years. Each measurement
is the outcome of a binary test: if the prevalence ($P=\Sigma^{-1}I$) in the
node is above a threshold value ($c=$ 30% or 4% for the SIS- or the
SIS${}_{\text{E}}$-model, respectively), where all nodes are seeded at 10% or
2% initially. In Figure 5.4, we illustrate the sampling output in time for the
SIS${}_{\text{E}}$-model.

Figure 5.3. Illustration of the transport network; the red points are the
sentinel nodes, and the grey points are the latent ones. Red/black lines are
transport events into a sentinel- or latent node, respectively. Figure 5.4.
The population-weighted average prevalence (red) is unobservable, but the
pseudo prevalence (blue) is obtained from weighting together multiple binary
measurements. The least squares $\mathcal{OU}$-fit for the pseudo prevalence
is used as summary statistics (a few samples in grey are shown).

A challenging aspect of many data-driven inference problems, e.g., including
network dynamics is that the likelihood function is intractable and must be
estimated through repeated simulations. Bayesian inference in this setting is
termed Likelihood-free inference (LFI), or Approximate Bayesian Computations
(ABC); see [28, 39] for reviews. In this example, we consider the Sequential
Monte Carlo (SMC) adaptation of ABC (SMC-ABC) implemented in SimInf and
described in [47].

Our SMC-ABC implementation determines proposal rejections per generation $n$
using the normalized Euclidean kernel
$K_{\varepsilon}(x,y)=\sqrt{\sum_{i}((x_{i}-y_{i})/x_{i})^{2}}<\varepsilon_{n}$
for statistics of the simulation proposal $y$ and observation $x$ and with a
series of decreasing tolerances $\varepsilon_{n}$. The statistics are computed
as follows. Each simulation generates a time series of _pseudo prevalences_ ,
i.e., a population-weighted sum of positive samples. We interpret the time
series as an $\mathcal{OU}$-process and select the least square estimates
(3.6)–(3.7) as _indirect_ summary statistics [9]. A word in favor of this
particular choice of statistics for other ABC implementations, e.g., synthetic
likelihoods [52], is that least square estimates are asymptotically normally
distributed under broad assumptions [32]. Notably, we get away with using
statistics with one more or equal dimension as the parameter set, suggesting
that this characterization is indeed very fitting.

For the inference we use a single initial simulation with the true parameters
$(\beta,\gamma,R_{0})_{\text{SIS}}=(0.16,0.1,=1.6)$ and
$(\beta,\gamma,\rho,R_{0})_{\text{SIS}_{\text{E}}}=(0.054,0.1,0.44,=1.108)$,
respectively, and we infer all parameters simultaneously. For priors, since we
have no likelihood and thus cannot easily produce a strictly non-informative
prior, we take uniform distributions over quite large intervals in parameter
space: $\beta$ and $\gamma\sim\mathcal{U}(0,1)$ in both cases and
$\rho\sim\mathcal{U}(0.4,0.5)$. We use the decreasing ABC tolerances
$\varepsilon_{n}=100\exp(-0.25(n-1))$, for $n=1\ldots 15$ and 1,000 SMC
particles.

We found that $\rho$ requires a tighter prior than the others for the
computations to complete in a reasonable time: the SIS${}_{\text{E}}$-model is
considerably more challenging but is also more realistic, particularly so in
our setting on a network where rather large prevalences are required for the
SIS-model to not simply die out.

The results for the $R_{0}$-marginals are displayed in Fig. 5.5, where we also
investigate the concentration effect of data through the relative change in
quartile coefficient of dispersion (QCD); a small concentration factor
indicates an accurately identifiable parameter. Although the
SIS${}_{\text{E}}$-model is clearly more challenging, $R_{0}$ is still well
reconstructed for both models.

Binary data implies identifiability when given multiple observations at the
same time, e.g., over a small collection of nodes in a network rather than on
a single node, as indicated in a qualitative sense by Theorem 4.6.
Additionally, $R_{0}$ is identifiable with quite high accuracy even when the
dependent parameters covaries. The example is prototypical of using synthetic
data to evaluate the feasibility of an intended setup. Since the posteriors
are robust with respect to capturing the synthetic truth, we have good reasons
to also have some faith in the design when approaching real data.

| $\beta$ | $\gamma$ | $\rho$ | $R_{0}$
---|---|---|---|---
SIS${}_{\text{E}}$ | 0.61 | 0.59 | 0.70 | 0.0086
SIS | 0.089 | 0.093 | - | 0.016

Figure 5.5. Top: posterior and prior distributions for $R_{0}$ with the true
values indicated. Bottom: recorded QCD concentration factors for all the
parameters (see text for details).

##  6\. Conclusions

Throughout this work we have employed the Ornstein-Uhlenbeck process as a
meta-model of more involved epidemiological models. We indicated in §3 a
convergence analysis of the Bayesian posterior under direct process
observations. Since this is an unrealistic setup in most epidemiological
applications, one can think of these results as _best possible_.

We next took the opposite standpoint and considered data to be severely
filtered such that, literally, each data point contributed only a single bit
of information. For instance, this could be a model of pooled data obtained
through environmental sampling and subsequent analysis. To obtain a closed
framework we added a fairly general stationary $p$-step Markov assumption and
worked out conditions on the data to obtain a non-singular inverse problem.

With increasing compute power and improving possibilities for gathering data,
fully Bayesian first-principle epidemiological models can be realized. As a
minimum standard, we propose, such methods should be preceded by a proof of
self-consistency: data generated from the model itself and chosen “nearby” the
actual data should allow for accurate parameter identifiability. Under this
basic standard, Bayesian epidemiological modeling with both short-term
prediction and generation of forecasting scenarios can be included as an
integrated part of the public health’s methodological arsenal.

## Acknowledgment

This work was financially supported by the ENS Paris (S. Bronstein) and in
part by the Swedish Research Council Formas (S. Engblom) and the Swedish
Innovation Agency Vinnova (S. Engblom, R. Marin). Several helpful comments and
kind suggestions on earlier versions of the manuscript were offered by friends
and colleagues. These helped to improve the presentation considerably.

## References

  * Andersson and Britton [2012] H. Andersson and T. Britton.  _Stochastic epidemic models and their statistical analysis_ , volume 151.  Springer Science & Business Media, 2012.  doi:10.1007/978-1-4612-1158-7.
  * Athreya and Lahiri [2006] K. B. Athreya and S. N. Lahiri.  _Measure theory and probability theory_.  Springer Science & Business Media, 2006.  doi:10.1007/978-0-387-35434-7.
  * Bajardi et al. [2012] P. Bajardi, A. Barrat, L. Savini, and V. Colizza.  Optimizing surveillance for livestock disease spreading through animal movements.  _J. R. Soc. Interface_ , 9(76):2814–2825, 2012.  doi:10.1098/rsif.2012.0289.
  * Balcan et al. [2009] D. Balcan, V. Colizza, B. Gonçalves, H. Hu, J. J. Ramasco, and A. Vespignani.  Multiscale mobility networks and the spatial spreading of infectious diseases.  _Proc. Natl. Acad. Sci. USA_ , 106(51):21484–21489, 2009.  doi:10.1073/pnas.0906910106.
  * Basawa and Rao [1980] I. V. Basawa and B. P. Rao.  Chapter 10 - Bayesian Inference for Stochastic Processes.  In I. V. Basawa and B. P. Rao, editors, _Statistical Inference for Stochastic Processes_ , pages 255–293. Academic Press, London, 1980\.  doi:10.1016/B978-0-12-080250-0.50017-8.
  * Bishwal [2000] J. Bishwal.  Rates of convergence of the posterior distributions and the Bayes estimations in the Ornstein-Uhlenbeck process.  _Random Oper. Stoch. Equ._ , 8(1):51–70, 2000\.  doi:10.1515/rose.2000.8.1.51.
  * Britton [2010] T. Britton.  Stochastic epidemic models: a survey.  _Math. Biosci._ , 225(1):24–35, 2010.  doi:10.1016/j.mbs.2010.01.006.
  * Brooks-Pollock et al. [2014] E. Brooks-Pollock, G. O. Roberts, and M. J. Keeling.  A dynamic model of bovine tuberculosis spread and control in Great Britain.  _Nature_ , 511:228, 2014.  doi:10.1038/nature13529.
  * Drovandi et al. [2011] C. C. Drovandi, A. N. Pettitt, and M. J. Faddy.  Approximate Bayesian computation using indirect inference.  _J. Roy. Stat. Soc.: Ser. C (Appl. Stat.)_ , 60(3):317–337, 2011.  doi:10.1111/j.1467-9876.2010.00747.x.
  * Engblom and Widgren [2017] S. Engblom and S. Widgren.  Data-driven computational disease spread modeling: from measurement to parametrization and control.  In C. R. Rao, A. S. Rao, and S. Payne, editors, _Disease Modeling and Public Health: Part A_ , volume 36 of _Handbook of Statistics_ , chapter 11, pages 305–328. Elsevier, Amsterdam, 2017.  doi:10.1016/bs.host.2017.05.005.
  * Engblom et al. [2020] S. Engblom, R. Eriksson, and S. Widgren.  Bayesian epidemiological modeling over high-resolution network data.  _Epidemics_ , 32:100399, 2020.  doi:10.1016/j.epidem.2020.100399.
  * Ethier and Kurtz [1986] S. N. Ethier and T. G. Kurtz.  _Markov Processes: Characterization and Convergence_.  Wiley series in Probability and Mathematical Statistics. John Wiley & Sons, New York, 1986.  doi:10.1002/9780470316658.
  * Eubank et al. [2004] S. Eubank, H. Guclu, V. A. Kumar, M. V. Marathe, A. Srinivasan, Z. Toroczkai, and N. Wang.  Modelling disease outbreaks in realistic urban social networks.  _Nature_ , 429(6988):180, 2004.  doi:10.1038/nature02541.
  * Ferguson et al. [2005] N. M. Ferguson, D. A. Cummings, S. Cauchemez, C. Fraser, S. Riley, A. Meeyai, S. Iamsirithaworn, and D. S. Burke.  Strategies for containing an emerging influenza pandemic in southeast Asia.  _Nature_ , 437(7056):209, 2005.  doi:10.1038/nature04017.
  * Florens-Zmirou [1989] D. Florens-Zmirou.  Approximate discrete-time schemes for statistics of diffusion processes.  _Statistics: J. Theor. Appl. Stat._ , 20(4):547–557, 1989.  doi:10.1080/02331888908802205.
  * Fournié et al. [2018] G. Fournié, A. Waret-Szkuta, A. Camacho, L. M. Yigezu, D. U. Pfeiffer, and F. Roger.  A dynamic model of transmission and elimination of peste des petits ruminants in Ethiopia.  _Proc. Natl. Acad. Sci. USA_ , 115(33):8454–8459, 2018.  doi:10.1073/pnas.1711646115.
  * Galani et al. [2022] A. Galani, R. Aalizadeh, M. Kostakis, A. Markou, N. Alygizakis, et al.  SARS-CoV-2 wastewater surveillance data can predict hospitalizations and ICU admissions.  _Sci. Total Environ._ , 804:150151, 2022.  doi:10.1016/j.scitotenv.2021.150151.
  * Genon-Catalot [1990] V. Genon-Catalot.  Maximum contrast estimation for diffusion processes from discrete observations.  _Statistics_ , 21(1):99–116, 1990.  doi:10.1080/02331889008802231.
  * Ghosh et al. [2016] A. P. Ghosh, W. Qin, and A. Roitershtein.  Discrete-time Ornstein-Uhlenbeck process in a stationary dynamic environment.  _J. Interdiscip. Math._ , 19(1):1–35, 2016.  doi:10.1080/09720502.2013.857921.
  * Gray et al. [2011] A. Gray, D. Greenhalgh, L. Hu, X. Mao, and J. Pan.  A stochastic differential equation SIS epidemic model.  _SIAM J. Appl. Math._ , 71(3):876–902, 2011.  doi:10.1137/10081856X.
  * Joag-Dev et al. [1983] K. Joag-Dev, M. D. Perlman, and L. D. Pitt.  Association of normal random variables and Slepian’s inequality.  _Ann. Probab._ , pages 451–455, 1983.  doi:10.1214/aop/1176993610.
  * Keeling and Rohani [2011] M. J. Keeling and P. Rohani.  _Modeling infectious diseases in humans and animals_.  Princeton University Press, 2011.  doi:10.1086/591197.
  * Kennedy et al. [2022] B. Kennedy, H. Fitipaldi, U. Hammar, et al.  App-based COVID-19 syndromic surveillance and prediction of hospital admissions in COVID Symptom Study Sweden.  _Nature Commun._ , 13(2110), 2022.  doi:10.1038/s41467-022-29608-7.
  * Kermack and McKendrick [1927] W. O. Kermack and A. G. McKendrick.  A contribution to the mathematical theory of epidemics.  _Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci._ , 115(772):700–721, 1927.  doi:10.1098/rspa.1927.0118.
  * Kessler [1997] M. Kessler.  Estimation of an ergodic diffusion from discrete observations.  _Scand. J. Stat._ , 24(2):211–229, 1997.  doi:10.1111/1467-9469.00059.
  * Latz [2020] J. Latz.  On the well-posedness of Bayesian inverse problems.  _SIAM/ASA J. Uncert. Quant._ , 8(1):451–482, 2020\.  doi:10.1137/19M1247176.
  * Liu et al. [2018] Q.-H. Liu, M. Ajelli, A. Aleta, S. Merler, Y. Moreno, and A. Vespignani.  Measurability of the epidemic reproduction number in data-driven contact networks.  _Proc. Natl. Acad. Sci. USA_ , 115(50):12680–12685, 2018.  doi:10.1073/pnas.1811115115.
  * Marin et al. [2012] J.-M. Marin, P. Pudlo, C. P. Robert, and R. J. Ryder.  Approximate Bayesian computational methods.  _Stat. Comput._ , 22(6):1167–1180, 2012.  doi:10.1007/s11222-011-9288-2.
  * McKinley et al. [2009] T. McKinley, A. R. Cook, and R. Deardon.  Inference in epidemic models without likelihoods.  _Int. J. Biostat._ , 5(1), 2009.  doi:10.2202/1557-4679.1171.
  * Merler et al. [2011] S. Merler, M. Ajelli, A. Pugliese, and N. M. Ferguson.  Determinants of the spatiotemporal dynamics of the 2009 H1N1 pandemic in Europe: implications for real-time modelling.  _PLoS Comput. Biol._ , 7(9):e1002205, 2011\.  doi:10.1371/journal.pcbi.1002205.
  * Mishra and Prakash Rao [1987] M. Mishra and B. Prakash Rao.  Rate of convergence in the Bernstein-Von Mises theorem for a class of diffusion processes.  _Stock.: Inter. J. Probab. Stoch. Proc._ , 22(1):59–75, 1987.  doi:10.1080/17442508708833467.
  * Newey and McFadden [1994] W. K. Newey and D. McFadden.  Chapter 36 Large sample estimation and hypothesis testing.  In _Handbook of Econometrics_ , volume 4, pages 2111–2245. Elsevier, 1994.  doi:10.1016/S1573-4412(05)80005-4.
  * Obadia et al. [2015] T. Obadia, R. Silhol, L. Opatowski, L. Temime, J. Legrand, et al.  Detailed contact data and the dissemination of Staphylococcus aureus in hospitals.  _PLoS Comput. Biol._ , 11(3):e1004170, 2015.  doi:10.1371/journal.pcbi.1004170.
  * Owhadi et al. [2015] H. Owhadi, C. Scovel, and T. Sullivan.  On the brittleness of Bayesian inference.  _SIAM Review_ , 57(4):566–582, 2015.  doi:10.1137/130938633.
  * Salathé et al. [2010] M. Salathé, M. Kazandjieva, J. W. Lee, P. Levis, M. W. Feldman, and J. H. Jones.  A high-resolution human contact network for infectious disease transmission.  _Proc. Natl. Acad. Sci. USA_ , 107(51):22020–22025, 2010.  doi:10.1073/pnas.1009094108.
  * Shardlow [2006] T. Shardlow.  Modified equations for stochastic differential equations.  _BIT Numer. Math._ , 46(1):111–125, 2006.  doi:10.1007/s10543-005-0041-0.
  * Shen and Wasserman [2001] X. Shen and L. Wasserman.  Rates of convergence of posterior distributions.  _Ann. Stat._ , 29(3):687–714, 2001.  doi:10.1214/aos/1009210686.
  * Shoji [1998] I. Shoji.  Approximation of continuous time stochastic processes by a local linearization method.  _Math. Comput._ , 67(221):287–298, 1998.  doi:10.1090/S0025-5718-98-00888-6.
  * Sisson et al. [2018] S. A. Sisson, Y. Fan, and M. Beaumont.  _Handbook of approximate Bayesian computation_.  CRC Press, 2018.  doi:10.1201/9781315117195.
  * Slud [1989] E. Slud.  Clipped Gaussian processes are never M-step Markov.  _J. Multivariate Anal._ , 29(1):1–14, 1989.  doi:10.1016/0047-259X(89)90072-9.
  * Söderström and Stoica [1989] T. Söderström and P. Stoica.  _System identification_.  Prentice-Hall International, 1989.
  * Sprungk [2020] B. Sprungk.  On the local Lipschitz stability of Bayesian inverse problems.  _Inverse Probl._ , 36(5):055015, 2020.  doi:10.1088/1361-6420/ab6f43.
  * Stehlé et al. [2011] J. Stehlé, N. Voirin, A. Barrat, C. Cattuto, V. Colizza, et al.  Simulation of an SEIR infectious disease model on the dynamic contact network of conference attendees.  _BMC Medicine_ , 9(1):87, 2011.  doi:10.1186/1741-7015-9-87.
  * Stoltenberg and Hjort [2021] E. A. Stoltenberg and N. L. Hjort.  Models and inference for on-off data via clipped Ornstein-Uhlenbeck processes.  _Scand. J. Statist._ , 48(3):908–929, 2021.  doi:/10.1111/sjos.12472.
  * Stuart and Teckentrup [2018] A. Stuart and A. Teckentrup.  Posterior consistency for Gaussian process approximations of Bayesian posterior distributions.  _Math. Comput._ , 87(310):721–753, 2018.  doi:10.1090/mcom/3244.
  * Stuart [2010] A. M. Stuart.  Inverse problems: a Bayesian perspective.  _Acta Numerica_ , 19:451–559, 2010.  doi:10.1017/S0962492910000061.
  * Toni et al. [2009] T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. P. Stumpf.  Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems.  _J. Roy. Soc. Interface_ , 6(31):187–202, 2009\.  doi:10.1098/rsif.2008.0172.
  * Toth et al. [2015] D. J. Toth, M. Leecaster, W. B. Pettey, A. V. Gundlapalli, H. Gao, J. J. Rainey, A. Uzicanin, and M. H. Samore.  The role of heterogeneity in contact timing and duration in network models of influenza spread in schools.  _J. R. Soc. Interface_ , 12(108):20150279, 2015.  doi:10.1098/rsif.2015.0279.
  * Uhlenbeck and Ornstein [1930] G. E. Uhlenbeck and L. S. Ornstein.  On the theory of the Brownian motion.  _Phys. Rev._ , 36(5):823, 1930.  doi:10.1103/PhysRev.36.823.
  * Widgren et al. [2018] S. Widgren, S. Engblom, U. Emanuelson, and A. Lindberg.  Spatio-temporal modelling of verotoxigenic Escherichia coli O157 in cattle in Sweden: Exploring options for control.  _Veterinary Res._ , 49(78), 2018.  doi:10.1186/s13567-018-0574-2.
  * Widgren et al. [2019] S. Widgren, P. Bauer, R. Eriksson, and S. Engblom.  SimInf: An R package for data-driven stochastic disease spread simulations.  _J. Stat. Software_ , 91(1):1–42, 2019.  doi:10.18637/jss.v091.i12.
  * Wood [2010] S. N. Wood.  Statistical inference for noisy nonlinear ecological dynamic systems.  _Nature_ , 466(7310):1102, 2010.  doi:10.1038/nature09319.
  * Zhang et al. [2017] Q. Zhang, K. Sun, M. Chinazzi, A. P. y Piontti, N. E. Dean, et al.  Spread of Zika virus in the Americas.  _Proc. Natl. Acad. Sci. USA_ , 114(22):E4334–E4343, 2017.  doi:10.1073/pnas.1620161114.