Neural Likelihood Approximation for Integer Valued Time Series Data
Authors: Luke O'Loughlin
Andrew J. Black
John Maclean
Venue: Accepted by TMLR ()

Abstract: Stochastic processes defined on integer valued state spaces are popular within the physical and biological sciences. These models are necessary for capturing the dynamics of small systems where the individual nature of the populations cannot be ignored and stochastic effects are important. The inference of the parameters of such models, from time series data, is challenging due to intractability of the likelihood. To work at all, current simulation based inference methods require the generation of realisations of the model conditional on the data, which can be both tricky to implement and computationally expensive. In this paper we instead construct a neural likelihood approximation that can be trained using unconditional simulation of the underlying model, which is much simpler. We demonstrate our method by performing inference on a number of ecological and epidemiological models, showing that we can accurately approximate the true posterior while achieving significant computational speed ups compared to current best methods.

Published in Transactions on Machine Learning Research (10/2024)
Neural Likelihood Approximation for Integer Valued Time Series Data
Luke O'Loughlin
luke. luke. edu...
Deportment of Mathemaical Sciences University of Adelaide
John Maclean
mohn. Wadelaide..ddddeee edu......ed..
Department of Mathemaical Sciences University of Adelaide
Andrew Black
andrew. black www.bbllckkaalelneo @adelaide. edu
Department of Mathemaical Sciences University of Adelaide.
Reviewed on OpenReview: https: /ttps://openreer opwnrevee.......... net/ forum? id- MMyRBe4oKF
Abstract
Stochastic processes defined on integer valued state spaces are popular within the physical and biological sciences. These models are necessary for capturing the dynamics small systems where the individual nature of the populations cannot ignored and stochastic effects are important. The inference the parameters of such models, from time series data, challenging due to intractability of the likelihood. To work all, current simulation based inference methods require the generation of realisations of the model conditional the data, which can be both tricky to implement and compulationally expensive. In this paper we instead construct a neural likelihood approximation that can be trained using unconditional simulation of the underlying model, which is much simpler. We demonstrate our method by performing inference on a number of ecological and epidemiological models, showing that we can accurately approximate the true posterior while achieving significant computational speed ups compared to current best methods.
1 Introduction
Mechanistic models where the state consists of integer values are ubiquitous across all scientific domains, but particularly in the physical and biological sciences (van Kampen 1992 Wilkinson 2018 Integer valued states are important when the population being modelled is small and so the individual nature of the population cannot be ignored. Models of these types are used extensively for modelling in epidemiology [Allen] 2017), ecology (McKane & Newman, 2005), chemical reactions [Wilkinson] 2018), queuing networks Breuer & Baum 2005) and gene regulatory networks (Shea & Ackers 1985) to name but a few examples.
Typically the main interest in learning or inferring the model parameters, rather than unobserved latent states, as the parameters dynamics. In spite of this, most methods for actually performing such inference typically still involve gencrating realisations of the full hidden process, commonly known as simulation based inference (Cranmer et al. 2020). Particle marginal Metropolis Hastings (Doucet et al.] 2015 Chopin & Papaspiliopoulos: 2020 and approximate Bayesian computation (isson et al., 2019 are two well known methods in this class.
directly related to the individual level mechanisms that drive the overall
As with most state-space models, the underlying state only partially observed, but a particular feature of many integer valued models that the subset of the state that is observed is done so accurately. The typical situation includes highly informative observations of parts of the state, while othe parts completely

---

Published in Transactions on Machine Learning Research (10/2024)
unobserved Del Moral et al. (2015). Unconditional simmlation of the model is unable to generate realisations of the process that are "close" to the data, which are effectively rare events and hence generated with low probability (rovandi & McCutchan 2016). Thus the application of simulation based inference methods for this class of models requires simulation of realisations conditional on the data to work (Bickel et al. 2008). Typically this involves the use of importance sampling or control variates to guide simulations (Lin jtt al. 2013 Black 2019). These strategies can work well, but there a fundamental limitation that better realisations (i.e. those closer to the data) are more computationally expensive to generate and hence this often offsets any gains in the overall running time. Importance sampling can also be difficult to implement correctly and has to be done on model by model basis-it is the opposite of black box algorithm (Owen 2013). In this paper we take a different approach and instead construct a neural conditional density estimator (NCDE) for modelling the full likelihood of integer valued time series data that can be trained using unconditional simulations of the model. Such simulations are both easy to implement and much faster to run.
The NDCE likelihood model is designed to be autoregressive using a convolutional neural network (CNN) architecture with an explicit casual structure to model conditional likelihoods by a mixture of discretised logistic distributions. The likelihood approximation is trained by modifying the sequential neural likelihood algorithm (Papamakarios et al., 2019) to jointly train the model and sample from an approximation of the posterior. We demonstrate our methods by performing inference on simulated data from two epidemic models, as well as a predator-prey model. The resulting posteriors are compared to those obtained using the exact sampling method particle marginal Metropolis Hastings (PMMH) (Andrieu et al. 2010 which is currently the sratt-of-thheane method (Black) 2019). We are able to demonstrate accuracy in the posteriors generated from our method, suggesting that the autoregressive model is able to learn a good approximation of the true likelihood. Furthermore, method is significantly less computationally expensive to in scenarios where PMMH struggles.
2 Problem Setup
The models we consider are continuous-ting time Markov chains (CTMC), (t), ( 0}, also known as a Markov jump process (Ross 2000 Wilkinson 2018). The state of the system, ((), is a vector of positive integers, which typically represent populations sizes or counts of events. The dynamics of such models are specified by describing the possible transitions the state and the rates at which these occur. The possible state transitions are fixed for given model, but the rates depend on small number of parameters, 0, some or all of which we wish to infer from partial time-series data, V1.., generated from the model. Our observation model is that parts of the state are observed exactly, but other parts not at all. Formally the observation likelihood can be written as a product of delta functions
p(y|=)= II 5 (()))))) JEB
(1)
where of represents the jth component of the state and B is the set of components that are assumed observed.
Although the models evolve in contimuous time, we assume that obscrvations are made at discrete time intervals and only depend on the state of the system at that time, so the inference problem fits within the a state-space model (SSM) framework DDel Moral 2004 Murphy| 2023). Denoting the latent state and the observed data at time step i as æ; and y, respectively, then the joint likelihood given by
1
p( Vmm xb) II p(a 0),
(2)
where p( 0) is the Markov transition kernel of the model and observations are assumed conditionally independent of all previous saates/observations given the current state (Del Moral 2004)
In the Bayesian setting, we wish to sample from the posterior distribution of 0 given by
(1))):.) & p(y1.n ())(())
2

---

Published in Transactions on Machine Learning Research (10/2024)
where p(0) is the prior and p(y1m 0) the likelihood. Typically this is done using Markov chain Monte Carlo (MCMC), but this is difficult for this class of models as the likelihood is intractable. This is because, although in theory the transition kernel p(æ,/2;-1,, 0) can be numerically computed as a matrix exponential, in practice this is intractable for all but the simplest or smallest systems (Black et al.) 2017 Sherlock 2021 Hence the likelihood, obtained by marginalising over #i:n in Equation is also intractable.
While the likelihood itself is intractable, simulation of CTMCs (sampling from the transition density p(a,/x,,.,, 0)) is straightforward [Gillespic 1976.. Particle filter approaches are natural and popular way of estimating the likelihood, which can then be used for inference of the parameters using particle marginal Metsopolis-Hastings [Wilkinson 2018 Chopin & Papaspiliopoulos] (2200). For the class of problems we are interested in, with perfect but incomplete observations, the difficulty implementing particle filters stems from the form of the observation likelihood Equation []] As this equation the product of delta functions, only realisations of the latent process that exactly match the observations are assigned non-zero likelihood. If the state space is naturally large, or the proposed parameters are not compatible with the observations, then producing iimulations that match the data becomes a rare event problem. Conditional simulation meth- ods (i.e. simulating realisations of the model conditional the observed data using importance sampling) can be used to guide the simulations to match the data (Black] 2019 We use these as gold standard comparison in this paper, but as we also show they can become computationally expensive in certain situa- tions. Our approach, detailed in the next section, is to instead construct a NCDE that models the likelihood p(y (0) directly. This still relies on simulations of the model for training, but these are basic unconditional simulations that are computationally cheap.
3 Methods
In this section, we construct an autoregressive model to approximate p(y1.. The simplest version of the model works with univariate time series data Vı.n and is described in Section 3.1, We break up the construction into two parts: First, in Section 3.1.1 we describe CNN architecture using causal convolutions (van den Oord et al.) 2016), which used to map the input sequence Vinn to an output sequence 0 such that o "does not look into the future" (it is only dependent on y1.1 and 0). Second, in Section 3.1.2 we describe a way to map each o onlo probability distribution supported on the set of integers Z (or a subset thereof), so that we can use o, approximate p(y++1 0). In Section 3.2] we describe how the autoregressive model can be extended to multidimensional integer valued time series y1:n. In particular, this extension allows the autoregressive model to model correlations between the features of y:- Finally, in Section 3.3 provide overview of sequential neural likelihood (SNL), which is an algorithm for simultancously training the surrogate likelihood and performing approximate inference (Papamakarios et al. 2011)). We also discussing some modifications to SNL to help increase consistency in the approximate posterior.
3.1 Autoregressive Model
Autoregressive models approximate the density of a time series yiuu by modelling a sequence of conditionals,
namely p(y) exploiting causality in the time series (Shunaay & Stoffer 2017). We will denote the autoregressive model conditionals by qe (3 where vb are the autoregressive model parameters e.g. neural network weights. Suppressing the dependence on 0 the joint density for the likelihood can be factorised as
p(y1.)) p(v1)
I IplM
then we can define the density of yıu under the autoregressive model in the same way as
qo (Vin) = qu (vn) II (olllan)))
(3)
Note that these parameters are unrelated the parameters the CTMC
3

---

Published in Transactions on Machine Learning Research (10/2024)
The parameters v are learnt from training data, using Equation to construct suitable loss function. The standard choice for the loss is the negative log-likelihood (NI..) (Shumway & Stoffer 2017], and it is what we use here. In particular, for a dataset consisting of N time series {yim))=1 the NLL loss leads minimisation of unbiased estimate of the forward KL-divereence DKL(plluu)
(3)
p
y i:n
D(pllas)=
log
N
())
q
(1:n)
N
1
log do (yim)
(3)
+ C,
where C' is constant with respect to vb. Since the KL divergence is non-negative and equal to 0 if and only if (3(m)) plyin) (Bihhop 2006), it follows that sufficiently flexible autoregressive model will be able to approximate the true density arbitrarily well in the large dataset limit. When conditioning on 0, the NLL loss is still appropriate, as if e ~ p(0) for some proposal p(0), then as pointed out by Papamakarios et al. 2201)), in the large limit
N
yi:n (0)
1 N
log
Be-p(e)
(V).m 0))]
y Vion (06)
The right hand side is non-negative and equal to zero if and only if p(y1:n 0) (3.mmm almost everywhere the support of p(0).
Since we are interested in using qu (y1m.)0)) as a surrogate within an MCMC scheme, where it is repeatedly evaluated, it is preferable to use a model architect which allows for efficient evaluation of Equation 3 For univariate time series, CNN based architectures have been proposed which allow for parallel evaluation of all of the terms in Equation 3 These CNNs, called causal CNNs, map the univariate input sequence yi:n to the sequence of log-conditionals (log qoi (y: ()) A typical example is the model WaveNet (van den) Oord et al. 2016), which uses causal CNN to synthesise audio data. Due to the efficiency of convolutions in deep learning, differentiation allows for computation of the gradient of the causal CNN surrogate (with respect to the model parameters 0), and this has the same complexity as the ivaluation (Baydin et al. 2017), so efficient MCMC samplers such as NUTS (Hoffman & Gelman 2014] can be employed to sample from the posterior rather than basic random walk Mctropolis Hastings samplers, which are much less efficient.
construct our autoregressive model using causal convolutions. Additionally, automatic
3.1.1 Causal Convolutions
A one-dimensional causal convolution can be understood as regular one-dimensional convolution with appropriate zero padding on the left of the input. For a convolution kernel with length and de output channels, say Wi.k xk then the causal convolution operation is defined by
Cs (1 Vi:n) wi:k 0,
vim)
k-1
ca (wik V1..): =
wi (ob-Jiisilscs
where we have used to denote the standard convolution operation. It is clear that Cx(Wx-jj V1:n); depends only on the input sequence only through the elements (и((); where 1(i) max max (1, (A: 1)}, hence the causal structure in the input sequence is preserved by ck. By constructing CNN layers using the operation Cs (assuming that the kernel length is l in every layer), then it is clear that composing these layers will yield a function which also preserves the causal structure of the input sequence. Moreover, denoting the output
WaveNet also dilated convolutions to increase the length receptive field of the model.
4

---

Published in Transactions on Machine Learning Research (10/2024)
sequence by Ol:ns then it is straightforward to show that O, only depends Vi.n through the elements (0)= (see Figure Th where h is the number of CNN layers and 4 max (1, i - h(k 1)}. The purpose of constructing such a CNN is that since 0,-1 only depends on elements to time step 1- then it can used to define qu (1 To do this, 0;-1 is used to compute the paramcters of a suitable parametric distribution with traciable probability mass function, which discussed Section 3.1.2) To aid with training of the model, we also use residual connections the hidden layers of the CNN, which are known help stabilise training (He et al. 2016).
01
0
03
04
05
06
07
0g
O
O
1/2
y/3
y/s
y/
V7
V8
Figure 1: An example of the dependency structure in causal convolutional neural network with 4 layers and using a convolution kernel length of 3 in every layer. Notice that paths from back to the input elements Yj only exist for j < i.
By construction, the cansal CNN can only model conditionals of the form p(u p(u)l(,,())))))) due to the finite receptive field of the neural network. particular, the receptive field grows linearly with the convolution kernel length & and the number of hidden layers in the network h. For the problems considered in this paper, there are no relevant long range correlations in the data, the finite receptive field approximation is reasonable. If long range correlations between observations were important for particular problem, then the model could be altered by dilating the convolutions as this modification produces receptive field that grows exponentially with the number of hidden layers (van den Oord et al. 2016).
To model the likelihood p(11: 0), we must augment the CNN so that 0,-1 is dependent on 0 as well as V1s-1- We perform the augmentation by using shallow neural network calculate a context vector c(0), and then adding lincar transformation of c(9) to the output of cach layer of the CNN (broadcasting across the time axis). For inference problems where some parameters are unidentifiable, or where their elfect on the dynamics should allow for greater accuracy in the autoregressive model. We chose to make every layer of the CNN dependent on 0, as this allows for greater interaction between the model parameters and observations. This is important for non-linear models where small changes in the model parameters can cause noticeable, or even substantial changes in the observation dynamics (Kuznetsov) 2004 Börgers 2017 Figure 2 shows a diagram of the details of the evaluation of the CNN.
minimal, the neural network can help to learn suitable reparameterisation of 9, which
Further details of the model architecture given in Appendix A of the supplement.

---

Published in Transactions on Machine Learning Research (10/2024)
Causal Convolutiin
0l:
(A4, Rc)
Clnnoal Convotation (du. da,)
Residual Block
X h
# wir c(0)
Causal Convolation
(ds, da)
(i)
Caasal Convolation
Skallow
c(0)
1
91:n
0
Figure 2: Evaluation of the causal CNN for one-dimensional inputs. The mumbers (a, b) in the causal channels for the convolution layer. The quantity dA
convolution blocks denote the number of denotes the number of channels used in the hidden layers, and h is the number of residual blocks used. Wo have used @ to denote addition a vector and sequence of vectors broadcasted across the time axis.
3.1.2 Discretised Mixture of Logistics
To calculate the autoregressive model conditionals from O1:, we use the components of each o to compute the parameters of a discretised mixture of logistic distributions (DMoL) (Salimans et al., 2017). The discretisod logistic (DL) is family of distributions supported Z, characterised by a shift parameter μ R and scale parameter >> A DI.(H, s) random variable Y is sampled by first sampling X (ogistic((s, and then calculating Y [X], where I-] denotes rounding down to the nearest integer. The logistic distribution has similar properties to the Gaussian distribution, however it has an analytically expressible cumulative distribution function (CDF), which makes it possible to calculate the probability mass function
Pr Pr = m) = Pr (X < m + Pr < m).
Denoting the probability mass function (PMF) of a discretised logistic with shift and scale by (DL H4 8). then a mixture of DLs has PMF
PDMoL (v)
Σ ae PDL PDL (v:
(4)
where the &; terms are weights satisfying o 2 and Ei 6: = 1. The DMoL serves as an analytically tractable family of distributions that flexible enough to approximate a wide class of distributions supported on Z, analogous to a mixture of Gaussians the continuous case (Salimans et al. 2017). It can also helpful truncate the discrctisod logistics when a problem involves data supported on subset {a < n < b Z).
This is easy, it only involves renormalising the PMF of each DL component in Equation 4 which can be done analytically, noting that the DI. CDF the same the logistic CDF restricted to Z.

---

Published in Transactions on Machine Learning Research (10/2024)
To use DMoL to represent the autoregressive model conditionals, use the outputs of the CNN, 01:: to calculate the shift, scale and mixture proportion parameters. Assuming that by (Vil 1, has e mixture components, then we must calculate 3e 1 parameters, e shifts { (1djjml....d
scales {s] and c c-1 mixture proportions {o] (0)))_,....--1 ( is redundant since Σj 0z = 1). We use a CNN with 3e output channels (i.c. O; I,, and then the parameters of the DMoL are given by
m = yi-1 + o S softplus(o) +)) +
(5a)
(5b)
(o), of) softmax 20+1 ((o, ,0,,))
(5c)
where of is the jth element of o; and is a small constant. for numerical stability which take to be 10-6 Note that Equations 5 (--c) ensures that of > E and of for all i, J, and Ej of 1 for all i, so the conditionals qo (8 (3 0) are well defined distributions.
For some models we alter the default expressions in Equation 5 to better reflect the particular structure of model. For example, in Sections 4.1 and 4.2 the epidemic models investigate produce count data which is bounded above by the total population size N, hence we scale si by (N - -V--)/N to reflect the shrinking support of the conditionals over time.
3.2 Multivariate Data
When the time series consists of multivariate observations Vi, then constructing the antoregresive model is less straightforward. The main hurdle is to replace the DMoL distribution used in the univariate case with a suitable parametric family of multivariate discrete distributions. For y, E zd, d > 1, the conditionals p(y, (1y.111111100 0) can be factorised as follows:
d II j=2
p(y 91s-1, 0) = p(u! 1 0)
plu v. 19-1 194-1110 0),
(6)
where of is the jth component. of y: and y; 1:5-1 a vector of the first j j-1 components of yi- It follows that the problem of modelling ( 1. 0) can be approached by modelling the d univariate conditionals on the [914-11
right hand side of Equation We will denote the univariate conditionals approximated by the autoregressive model by qu y, 1:j
11.-1, 0), and define the full autoregressive model conditionals by
d go 0) = qu (vt -1 0) II go (v!l 1--1 -y.--,10) 0).
j=2
The CNN architecture described in Section 3.1.1 is not immediately suitable for calculation of go(2ll 1-1 11--1, 1, 9). This is bocause to construct a suitable model for qe (o output element of the CNN o; must decompose a way such that there is a group of elements which only depend on y: through y 1:3-1 which can then be used to define (s (s/lll 1351 -1.-11,0 0). However using the CNN from Section 3.1.1 yields an output O; that is dependent on Yi in an arbitrary way, hence o, cannot be split like this. A simple solution to this issue is to notice that may 'unroll' the time series Vın into long univariate time series
y1.--1, 0), the zth
Lidn = (vi .i...., 12
Vn
M
and then the CNN architecture from Section 3.1.1 can be used on 2iddn to define
qu (aadt)//1144-111 0) = ( 1y 1-j-1 1
In practice, we do not unroll the time scries, and instead force certain weight matrices the convolution kernels of the CNN to have block lower triangular structure. The block lower triangular structure of the weight matrices is achieved with binary masking in the density model MADE (Germain et al. 2015 The
7

---

Published in Transactions on Machine Learning Research (10/2024)
block lower triangular structure ensures that o can be partitioned into blocks of elements, where the first block is not, dependent on Yi, and the jth block is dependent y, for > 1. A more detailed discussion of the masking procedure is given in Appendix A.2 For the generic case, we use 3cd output channels in the CNN, so that cach block (of size 3c) of O can be used to calculate the 3c 11 parameters of a DMoL distribution, as described in Section 3.1.2 This can be altered necessary, e.g. if one of the features is binary, then number of output channels and masks for the weight matrices can be adjusted so that only one parameter is output for the binary random variable (corresponding to the parameter of a Bernoulli distribution).
Note that by modelling the conditionals as described above, we must choose an ordering for the components of y In general, such an ordering will not be natural, and although in theory the decomposition in Equation 6 will be agnostic to the choice of ordering, in practice the neural network will be able to represent certain orderings more easily than others. Hence, the training and inference may be more stable for some orderings of the components compared to others, and problem specific knowledge may be required find a sensible choice for the ordering.
3.3 Sequential Neural Likelihood
Sequential neural likelihood (SNL.) is a SBI algorithm which uses a NCDE as a surrogate for the true likelihood within a MCMC scheme (Papamakarios et al.) 2019). The main issue that is addressed by SNL is how to choose the proposal distribution over parameters in the training dataset p(0). Idcally, the proposal would be the posterior over the model parameters so that the surrogale likelihood is well trained in regions where the MCMC sampler is likely to explore, without wasting time learning the likelihood in other regions. Obviously the posterior is not available, so instead SNL iteratively refines by sampling additional training data over a series of rounds. The initial proposal is taken to be the prior, and subsequent proposals are taken to be an approximation of the posterior, which is sampled from using MCMC with the surrogate from the previous round as an approximation of the likelihood. We refer the reader to Paapaakaariss al. (2019) for more details on SNL.
We use our autoregressive model go (11:1 I0) within SNL to perform inference for the experiments described in Section 4 We make modification to how we represent the SNL posterior, in that we run SNL over r rounds, and keep a mixture of MCMC samples over the final r' < r rounds. The resulting approximate posterior is therefore an equally weighted mixture of approximations over successive SNL rounds. We do this because somctimes the postcrior statistics exhibit variability which does not settle down over the SNL rounds, it can be difficult to assess whether the final round is truly the most accurate. Instead we use an ensemble to average out errors present in any given round, which worked well in practice. Variability across the SNL rounds occurs mostly when the length of time series is random, presumably because the random lengths inflate variability in the training dataset.
4 Experiments
We evaluate our methods on simulated data from three different models, comparing the results to the exact sampling method PMMH (Andrieu et al.] 2010).. The first two models we consider are epidemic models, namely the SIR and SEIAR model, the former being simple but ubiquitous model in epidemiology (Allen 2017), and the latter being a more complex model which admits symplomatic phases (Black, 2019). The third example the predator-prey model described in McKane & Newman (2005).
We perform two sets of experiments using the SIR model; one set considers single simulated outbreak in population of size 50, and the other considers multiple independent observations from outbreaks within houscholds (small closed populations). The experiments on the single outbrcak serve as a proof of concept for our method, as PMMH can target the posterior for this data efficiently [Black] 2019). The experiments on household data allow us to assess the scalability of our method with an increasing number of independent time series, where PMMH does not scale well due to the exponential increase in the number of particles nceded to keep the variance of the likelihood estimate low cnough (Sherlock et al. 2015
TTe ordering must. be fixed, Le. cannot depend
8

---

Published in Transactions on Machine Learning Research (10/2024)
The SEIAR model experiments serve two purposes, one being to assess the accuracy of our method when trained on data from a more complex model whose likelihood includes complex relationships between the model parameters. The other purpose is to assess the scalability with the population size, as done in Black 2019), as increasing population sizes generally causes the performance of PMMH to degrade due to combination of longer time series and a larger model stale space. The predator-prey model experiments also serve two purposes, one being to assess our method on model with multivariate observations, and the other being to assess the scalability with lowering observation errors, where PMMH again generally performs poorly due to weight degeneracy (Del Moral & Murray 2015).
We run SNL 10 times with a different random seed for every outbrcak experiment, and 5 times for the predator prey model experiments. Each run of SNL uses 10 rounds with NUTS (Holfman & Gelman 2014)) as the MCMC sampler, and using the final 5 rounds to represent the posterior. Unless otherwise stated, we used a simulation budget of 60,000 samples for SNL, beginning with 10,000 samples from the prior predictive. All neural network modcls were constructed using JAX (Bradbury et al.) 2018) and the MCMC was run using numpyro (Phan et al. 2019). We validate our method against. PMMH: for all outbreak experiments use the particle filter of Black 2019 as it is specifically designed for noise free outbreak data. For the prodator-prey model, we use a bootstrap filter (Kitagawa 1996). We tuned the paramcters of PMMH (e.g. number of particles, the Metropolis- Hanttiss proposal) based on a pilot run, and subsequently ran the tuned PMMH chain for as many MCMC iterations as necessary to generated a target effective sample size (FSS) of at least 1,000.
To assess the accuracy of our method, we compare the posterior means and standard deviations for each parameter obtained from PMMH and SNL. We compare the estimated posterior means relative to the prior standard deviation (denoted M) to account for variability induced by using wide priors. We also compare the ratio of the estimated standard deviations, shifted so that the metric is 0 when SNL is accurate (denoted S):
M
JINNT. μPMMH
S
ONNL
1.
PPrior
ØPMMH
We also report classifier 2-sample t-test (C2ST) score [Lucckmann et al. 2021 for the SEIAR model experiments, which measures the similarity of two distributions by evaluating the average score of a classifier trained to distinguish between samples from the two distributions. We do this because the 5-dimensional posteriors have complex goometries that make mcan and variance comparisons insufficient assessing the discrepancies between the two posteriors. In the case of the predator-prey experiments, achieving a large enough ESS from the PMMH posterior to train a classifier becomes far more expensive, which is why we omit the C2ST score in this case. All posterior pairs plots can be found in Appendix for a qualitative comparison between the results. We assess the runtime performance of our method by calculating the effective sample size (of the MCMC samples) per second (ESS/s) averaged over the parameters, where we take the ESS for SNL to be an average of the ESS values obtained from the MCMC chains over the final rounds. For all experiments, the chosen simulation budget ensured that we generated an ESS of approximately 1,000-3,000 per SNL round. More detailed descriptions of the models and additional details for the experiments such as hyperparameters can be found in Appendix B
4.1
SIR model
We first consider outbreak data from the SIR model (Aleen 2017). This model a CTMC whose state at time & counts the number of susceptible and infected individuals, S, and I, respectively. Equivalently, the state can be given in terms of two different counting processes: the total number of infections Z1, N S, and removals Z2, N - (S: + L:), where N is the population size. Our observed data is time series of total case numbers Vi.,, where yi Zi.i, equivalent to recording the number of new cases every day without error. Additionally, we assune that the outbreak has ended by day n so the final size of the outbreak is equal to yn- Since the outbreak never grows beyond yn. cases, we can replace the full time series by yi.r: where T is random, and defined as the smallest with y/4 = yn- The task is to infer the two parameters of the model, Ro > 0 and 1/12 > 0, using the data yi.r. To model the randomness J, we augment the autoregressive model with the binary sequence fi where f. dzy- The autoregressive model for this data is constructed by using
9

---

Published in Transactions on Machine Learning Research (10/2024)
Ro
Rg
Prior
Prior
SNL
SNL
PMMH
PMMH
Truth
Truth
1/v
1/y
2
6
8 10
2.5
7.5
10.0
2
10
2.5
7.5
10.0
Figure 3: Posterior Pairs Plots for SIR Model Experiments. Left: single outbreak experiment from a population of N = 50. Right: experiment using 500 houscholds. For the univariate posteriors, we show all 10 runs of SNL, but only of the runs for the bivariate posterior.
Vi:r as input to the CNN, and adding an additional output channel to the final layer of the CNN, which allows for calculation of (fllia 0). Using the multivariate formulation of the autoregressive model would be redundant, as by construction we always have p(y, fily1.4-1 fis-1, 0) p(y,, ft 0), and fi-1 for any i S this fact and more details concerning the autoregressive model construction in these experiments are explained in Appendix B.2
Single realisation Our proof of concept example is a single observed outbreak in population of size 50. The left of Figure 3 shows the posterior pairs plot of SNL compared to PMMH, which clearly shows that SNL provides accurate approximation of the true posterior for this experiment; see Appendix C for quantilative comparison. For this particular experiment, PMMH obtains an ESS/s of 69.38, which outperforms SNL with only 3.74. This is unsurprising, as we are not using large amount of data and the particle filter we used was specifically designed to perform well for observations of this form.
Household data Here, mimic outbreak data collected from households (Walker et al. 2017). by gen- erating h observations with household sizes uniformly sampled between 2 and 7. We use values h = 100, 200 and 500, assume that all households are independent, and that all outbreaks are governed by the same val- ues of Rg and 1/7, which have been chosen to be typical for this problem. The likelihood of an observation where y/ = 1, or where the household size is 2, can be calculated analytically, so we incorporate these terms into inference manually. We therefore train the autoregressive model to approximate p(y1:r 0 yr 1) for houschold sizes of 3 through 7. For all of these experiments, we increased the simulation budget to 100,000 as the observations are much shorter and less informative than in the previous set of experiments.
Table shows the quantitative comparison between the two posteriors, indicating a negligible bias in the SNL mcans on average. The SNL variance appears to be slightly inflated on average, mainly for the à = 500 experiments, though this is partly due to the PMMH variance shrinking with h, so the difference between the PMMH and SNL posteriors still small, as indicated by the right of Figure 3 The ESS/ of SNL. and PMMH are compared in the left panel of Figure 4 which clearly shows that SNL outperforms PMMH across each experiment, and that the drop off in SNL's performance is less significant for increasing h. This is mainly because the autoregressive model can be evaluated elficiently on the observations in parallel.
4.2
SEIAR model
This more sophisticated outbrcak model includes latent exposure period (E), symptomatic phases (I) and asymptomatic (A) phases [Black] 2019). The state of the model can be given in terms of 5 counting processes Zi,..... Zs,t: which count the total number exposures, pre-symptomatic cases, symptomatic cases, removals and asymptomatic cases respectively. In addition to Ro and 1/7, the model has three extra parameters: 1/0 2 0, K E (0, 1) and q E (0, 1). The observations are of the same form as the SIR model, where the case numbers at day given by y1 Zz,i, corresponding to counting the number of new
10

---

Published in Transactions on Machine Learning Research (10/2024)
Table 1: Comparison of SNL and PMMH posterior statistics for household data experiments. The values reported are the average over the 10 SNL runs and the maximum deviation from 0.
Parameter
h
M
avg
max
avg
max
100
-0.003
-0.02
0.01
-0.08
Ro
200
-0.002
-0.02
0.02
0.13
500
-0.003
-0.02
0.04
0.11
100
0.004
-0.10
0.02
0.10
1/7
200
0.002
-0.05
0.06
0.15
500
-0.03
-0.07
0.11
0.17
SIR HH Model
SEIAR Model
PP Model
ESS/s
10a
10a.
10-1.
SNL PMMH
10-1
10-1
10-2
1000
# of Households
Population
Figure 4: Comparison of ESS/s for SNL and PMMH for all three models. The lines plotted for SNL are the averages over the runs, and the error bars correspond to the maximum and minimum. All plots are on a log-log scale. The horizontal axis of the left panel scales the observed dataset size, showing that our method has improved scaling compared to the gold standard PMMH. The horizontal axis of the middle panel scales the system size and similarly shows improved scaling. The horizontal axis of the right panel scales the observation accuracy, discussed in Section 4.3. We observe the well-known phenomenon that PMMH runtime degrades with more accurate data, but SNL improves.
symptomatic cases each day. The likelihood has a complicated relationship between the parameters Ro and which manifests as a banana shaped joint posterior. We assess our method on observations from populations of size = 350, 500, 1000 and 2000.
1.00
PMMH
SNL
0.90
0.85
0.75
0.70
Ro
Figure 5: (Ro. q) joint posterior samples for SEIAR experiment with N = 500.
Table _ shows a comparison of the posterior statistics for the parameters Ro and k. The metrics for the other parameters are deferred to Appendix C The bias in the mean value of Ro is consistently small, but the standard deviations are slightly inflated, more so with larger population sizes. The Ro posterior for these experiments quite right skewed, and SNL produces postcriors which are slightly too heavy in the tail, which explains the inflated variance estimates. For H, the posterior means do exhibit noticeable bias, but the posterior standard deviations do not average. The bias in lowest for the N 500
11

---

Published in Transactions on Machine Learning Research (10/2024)
experiment, where the posterior has mode near O. The parameter only has weak relationship with the data, being the most informative near 0 or 1, hence the effect of k the likelihood is presumably difficult to learn. This would explain the biases in the A posteriors, and would also explain why the N = 500 experiment had the lowest bias on average. We also report the C2ST scores in Table 3 All outputs from the classifier are very close to 0.5, which is obtained when the classifier cannot learn distinguish between the two distributions. Additionally, p-value can be calculated from the classifier scores under the hypothesis that the bwo distributions are the same (Lopez-Paz & Oquabel 2017), and we see that all p-values are larger than 0.5 except for N = 2,000 (which is still larger than 0.05), suggesting that the bulk of the posterior masses must overlap significantly. For some visual confirmation of this, Figure 5 shows that SNL manages to correctly reproduce the banana shaped posterior between Rg and in the N 500 experiment. In terms of runtime, the centre panel of Figure | shows that the performance of SNL drops off al a slower rate than PMMH with increasing N, and is noticeably faster for the = 1,000 and 2,000 experiments.
Table 2: Comparison of SNL and PMMH posterior statistics for Ro and in the SEIAR model.
Parameter
N
M
avg
max
avg
max
350
-0.003
-0.02
0.04
0.08
500
0.005
0.09
-0.02
-0.09
Ro
1000
0.05
0.08
0.13
0.20
2000
0.02
0.04
0.20
0.31
350
0.20
0.32
0.01
0.04
500
-0.08
-0.20
-0.06
-0.14
1000
-0.16
-0.29
-0.0004
0.03
2000
-0.25
-0.35
0.02
0.07
Table 3: C2ST scores and p-values for the SEIAR model experiments.
N
Classifier Score
p-value
350
0.4996
0.9514
500
0.5031
0.6670
1000
0.4975
0.7235
2000
0.5127
0.0834
4.3 Predator-Prey model
Here, consider data generated from the predator-prcy model described in McKane & Newman (2055 which exhibits resonant oscillations in the predator and prey populations. The model is CTMC with state (Pt,Q.) where Pt and Qi are the numbers of predators and prey at time t respectively. The model has 5 parameters, 6, di, P1.P2 > 0, which are the prey birth rate, predator death rate, prey death rate, and two predation rates respectively. The task is to infer these parameters from predator and prey counts collected at discrete time steps, when any individual predator/pree has a probability r (assumed to be fixed and known) of being counted. We use the values 0.5, 0.7 and 0.9 for r, and all noisy obscrvations are generated from single underlying realisation, where we simulate the predator and prey populations over 200 time units, recording the population every 2 time units; Appendix C for plot of the data with 0.9. We order the predators before the prey in the autoregressive model conditionals in Equation 0 For these experiments, we only use a simulation budget of 25, 000, because the observations are longer and more informative than in previous sets of experiments.
Table 4 shows the comparison between the posterior statistics for the experiments with r 0.9, which sug- gests negligible differences between the posterior mean and standard deviations. It also worth noting that SNT correctly reproduces highly correlated posterior between dj and p1- Furthermore, the d2 posterior
12

---

Published in Transactions on Machine Learning Research (10/2024)
is almost indistinguishable from the prior, suggesting that the autoregressive model can learn to reparam- eterise 0 into set of identifiable parameters and suppress dependence on unidentifiable parameters. The accuracy of SNL increases with increasing for this model (see Appendix C for the other metrics), which is unsurprising, as there is a higher "signal-to-noise ratio' in the training data for larger т, which lcads to training data which more informative for p(y1:100) This is direct contrast to PMMH, where more noise makes the particle filter less prone to weight degeneracy. The right of Figure 4 shows that SNL an order of magnitude more efficient than PMMH for these experiments, and the ESS/s does not drop off with increasing r.
Table 4: Comparison of SNL and PMMH posterior statistics in the predator prey model with 0.9.
Parameter
M
S
avg
max
avg
max
-0.03
-0.07
-0.04
0.11
d1
0.01
0.03
-0.002
-0.05
d2
-0.002
-0.12
0.01
0.10
Pt
0.01
0.02
-0.006
-0.05
P2
-0.05
-0.10
-0.05
-0.08
5 Discussion
In this paper we have constructed surrogate likelihood model for integer valued time series data. Such data often arises from counting processes and models of small populations and presents particular challenges for performing inference due to the low noise in the observation process. Low noise leads to highly informative obscrvations of parts of the state, which means that naive application of simulation based inference methods usually fails; more sophisticated approaches, lypically utilising bespoke importance sampling simulations, do work, but are computationally expensive and difficult to implement. The main feature of our neural model is that it can be trained using simple unconditional simulations of the dynamical model. Unconditional simulation both much casier to implement and faster to run and our experiments demonstrated significant reduction in the overall computational expense of performing inference in scenarios where exact sampling methods struggle.
Our method produces good approximations of the parameter posterior, accurately reproducing many features of the posteriors obtained from exact sampling methods. There are negligible biases in the posterior means across experiments, with the exception of the parameter λ in the SEIAR model. This parameter has a very weak relationship with the data except near к = or 1, so it unsurprising that it difficult to learn. Our method always yielded posterior standard deviations which were either very close to that of РMMH, slightly too large, suggesting that we generally avoid overconfident estimates, which can be a problem in simulation based inference (Hermans et al. 2022). Furthermore, approximate posterior that has slightly inflated variance can still be used as an importance sampling distribution in more computational expensive exact method, which can significantly improve efficiency (Papamakiros & Murray] 2015)) Our method is also relatively robust to randomness within the algorithm, not showing any major differences over independent runs of SNL.
In the right-skewed Ro posteriors in the SEIAR model experiments, we saw that the SNL posteriors were slightly heavy in the tail. Improved training of the likelihood in these tail regions requires scheme for generating SNL proposals that better cover these, e.g. by jittering the proposed samples. A particularly simple alteration would be to use the geometric mean of the SNL posterior and the prior, which leads to more sampling in the tails and better performance in sequential ABC methods (Alsing el al. (018)). Bayesian active learning has also been employed in SNL before [Lueckmann et al. 2018). and this kind of approach has the advantage that allows for assessing the uncertainty in the surrogate likelihood, so might be possible to derive an uncertainty score for the approximate posterior samples based on the surrogate uncertainty. In general, more robust. SNL training routines are an important topic of current. research (Kelly et al. 2024
13

---

Published in Transactions on Machine Learning Research (10/2024)
Neural density estimators have been the topic of research for some time now (Germain et al. 2015 Papa-
makarios et al. 2017). The work most similar to ours is described in the original SNL. paper Papamakarios et al. (2019), but a limitation of that is the use of summary statistics (as in an ABC approach) for training the model. Instead, our autoregressive approach allows us to fit directly to the raw time scries data, not just statistics derived from it. This important for our applications involving very small populations (house- holds) where it is difficult to form mcaningful statistics due to the length and sparsity of the data (Walker eet al., 2017). The convolutionial architecture of our model also means that the approximate likelihood can be evaluated in single pass, rather than the sequential evaluation of recurrent based architecture which is slower (van den Oord et al. 2016: Vaswani et al. 2017).
A neural surrogate for the likelihood has a number of additional advantages over other simulation based inference approaches such as PMMH. Our likelihood approximation is differentiable with respect to the parameters of the model, so a fast Hamiltonian Monte Carlo (HMC) kernel (Neal 2011] can be used rather than a naive random walk MCMC kernel. Our method also has specific advantages over PMMH for problems involving longer time series or tall data, i.e. large numbers of iid observations [Bardenet et al. (2017), such as seen in Section 5.1. The performance of PMMH depends strongly on the variance of the likelihood estimate (Doucct et al. 2015 Sherlock et al. 2015 with the chain becoming 'sticky' if this is too high. The variance grows linearly with the number of observations, but in practice this will lead to a quadratic growth in the computational expense of running the filter as more simulated particles are needed for the additional observations and the number of particles for all observations must also be increased to reduced the variance of each. Our neural approximation is not an estimate, so does not suffer this problem, but because it is only an approximation it may introduce a bias in the posterior samples.
In practice, we found that some care is needed in choosing the support of the priors. For example, in the predator prey experiments, without truncating the tails of the log normal priors sometimes the MCMC chain would end up stuck in non-physical regions of parameter space. For these non- non-physical parameter values, the simulated realisations are more or less trivial (cxtinction/saturaal the surrogale has difficulty in learning the behaviour of the observations in these regions of the parameter space. Many SBI methods can suffer from similar problems if the simulated data and observation are not similar enough. Truncating the priors then works well to solve this issue because it prevents the surrogate from needing to be evaluated in these non-physical regions. Prior predictive checks are standard way of checking the validity of priors and can also help diagnose these issues (Gelman 2014). Additionally, physical interpretations of the parameters (as in the predator prey example) can often be used to find reasonable bounds for the model parameters.
occurs very quickly), so prosumably
A potential drawback of our method on time series with large number of features is that the autoregressive model requires an ordering in the features. Such an ordering may not be natural and therefore could require significant amounts of tuning. Furthermore, in general the number of hidden channels in the CNN grows linearly with the input dimension, so the computational complexity of evaluation generally grows quadratically with the input dimension. Choosing a feature ordering is an issue for the model MADE as well (Germmiin et al., 2015), and the authors address this by randomly permuting the feature orderings and creating an ensemble, so a similar technique could be employed for our model. We also point out that for many physical models choosing a feature ordering is unlikely to be a problem. In epidemiological settings, multidimensional time scries of counts would be rare, and for other dynamical systems of interest there is often additional causal structure that creates an obvious choice of ordering. Furthermore, for high dimensional data, it is likely the case that continuum approximation can be employed, or the use of summary statistics becomes viable choice, and in both cases the method developed in this paper would not be required.
This work asuuess that conditioning a finite number of steps back in time is sufficient. to model the condi- tionals. If longer range dependencies were important, then dilated causal convolutions (van den Oord et al. 2016) or self attention (Vaswani et al. 2017)) could be employed to construct the autoregressive model. Self attention could also be beneficial for high dimensional time series, since it has linear complexity in the num- ber of hidden features compared to the quadratic complexity of convolution layers. Another line of enquiry would be to develop better diagnostics to assess the results of SNL. The diagnostics used in PPapamakarios lt al. (2019) did not fit well with our integer valued data, nor did it fit well with the random length time series in the epidemic models. Studying how robust our method is to model misspecification would also
14

---

Published in Transactions on Machine Learning Research (10/2024)
be important future work, any real data is unlikely to truly follow specific model (Kelly et al. 2024 Finally, one could use family of distributions other than the DMoL. to represent the surrogate likelihood conditionals. For example, the household model experiments in Section 4.1 could possibly be treated by a simple multinomial/ categorical distribution. This works because the population sizes are small enough so that the observations yi can only take a few possible values, which can be treated as distinnt calegories.
References
Linda J.S. Allen. A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis. Infectious Disease Modelling, 2:128-144 142, 2017. doi: 10.1016/j.idm.22 201mmm..01.... 201..3.000.. 7.03.001.
Justín Alsing, Benjamin D. Wandelt, and Stephen M. Feeny. Alsing, j., wandelt, b.d. and feeney, s.m., 2018. optimal proposals for approximate bayesian computation. arXiv preprint arXiv:1808 10X00:1808.0644, 2018.
Christophe Andrieu, Arnaud Doucet, and Roman Holenstein. Particle Markov chain Monte Carlo methods.
Journal of the Royal Stalislicul Sociely: Series B (Stalisticul Methodology), 72(3):269 342, 2010. doi: 10.1111/j.1467-0 11611///446--96888 2009000090...
R. Bardenet, A. Doucet, and C. Holmes. On Markov chain Monte Carlo methods for tall data. Journal of Machine Learning Research, 18:1-43, 2017.
Atilim Günes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Au- tomatic Diffcrentiation in Machine Learning: A Survey. Journal of Machine Learning Research, 18(1):
January 2017.
Peler Bickel, Bo Li, and Thomas Bengtsson. Sharp failure rates for the bootstrap particle filter in high dimensions. In Instātute of Mathematical Statistics Collections, pp. 318-329, Institute of Mathematical Statisties, Beachwood, Ohio, 2008. doi: 10. 10.1214/074921110000000033
Christopher M. Bishop. Pottern Recognition and Machine Learning. Information science and statistics. Springer, New York, 2006. ISBN 978-0-387-31073-2.
Andrew J. Black. Importance sampling for partially observed temporal epidemic models. Stabistics and Computing, 29:617-630, 2019. doi: 10.1007/s112221 018 9827-1
Andrew J. Black, Nicholas Geard, James M. McCaw, Jodie McVernon, and Joshua V. Ross. Characterising pandemic severity and transmissibiliit from data collected during first few hundred studies. Epidemics, 19:61 73, 2017. doi: 10.1016/j.epidemm.22 2017.00000.01700000000...
Christoph Börgers. Canard Explosions, pp. 105-110. Springer International Publishing, Cham, 2017. ISBN 978-3-319-51171-9 doi: 10.1007/0088-33319--1777 15.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milnn and Qiao Zhang. JAX: compos- able transformations of Python+NumPy programs, 2018. URL http:///gittuub..mmooggggg//////
L. Breuer and D. Baum. An Introduction to Queueing Theory and Matriz-Analytic Methods. Springer- Verlag, 2005. doi:
Vpringer-Verlag,
Nicolas Chopin and Omiros Papaspiliopoulos. An Introduction to Sequential Monte Carlo. Springer, 2020. doi: 100007997--3000-0000
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation- inference. Proceedings of the National Academy of Sciences, 117(48):30055 2020. doi: 10.1073/pnas. p00.3///as.1111789111 1912789117.
P. Del Moral, P. Jasra, A. Lee, C. Yau, and X. Zhang. The alive particle filter and its use in particle Markov chain Monte Carlo. Stoch. Anal. Appl., 33:943 974, 2015.
Pierre Del Moral. Feymman- Kac Formulae. Springer, New York, 2004. doi: 10.1007 1---9797---664-9333--
15

---

Published in Transactions on Machine Learning Research (10/2024)
Pierre Del Moral and Lawrence M. Murray. Sequential Monte Carlo with Highly Informative Observations. STAM/ASA Journal on Uncertainty Quantification, 3:969-999 997, 2015. doi: 10.1137/15M1012114
A. Doucet, M. K. Pitt, G. Deligianmidis, and R. Kohn. Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood cstimator. Biometrika, 102:295-311 -13,, 2015. doi: h.1099//boomc aaannt/assu75...
C. C. Drovandi and R. A. McCutchan. Alive SMC?: Bayesian model selection for low-count time series models with intractable likelihoods. Biomelries, 72:344 353, 2016. doi: 10.1111/biom.122 12449.
Andrew Gelman. Bayesian Data Analysis. Chapman & Hall/CRC texts in statistical science. CRC Press, Boca Raton, third edition, 2014. ISBN 978-1-4398-4095-5.
Mathicu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle. MADE: Masked Autoencoder for Distribution Estimation. In International Conference on Machine Leurning, volume 37, pp. 881 889, 2015.
Daniel T Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics, 22(4):403-433 -34, 1976. doi: 10.1016/0021-000 (76) 90041-3.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep Residual Learning for Image Recognition, Tn IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 770-778, 2016. doi: 20.100/CVPR.2016.900
Dan Hendrycks and Kevin Gimpel. Gaussian Error Linear Units (GELUs), arXiv preprint or orXix:1606.08// 2016.
Joeri Hermans, Arnaud Delaunoy, François Rozet, Antoine Wehenkel, Volodimir Begy, and Gilles Louppe. A Crisis In Comulation-Based Inference? Beware, Your Posterior Approximations Can Be Unfaithful. Transactions on Machine Learning Research, 2022.
Matthew D. Hoffman and Andrew Gelman. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(47):1593-1623 2014.
Samuel Karlin and Howard M. Taylor. A First Course in Stochastic Processes. Academic Press, second edition, 1975.
Ryan P. Kelly, David J. Nott, David T. Frazier, David J. Warne, and Chris Drovandi. Misspecification- robust Sequential Neural Likelihood. 2024. doi: 1.48550/arxiv...0 wttp:/i....1113300 13366.
URL
abs/2301.13365
Genshiro Kitagawa. Monte Carlo Filter and Smoother for Non-Gaussian Nonlinear State Space Models.
Journal of Computational and Graphical Statistics, 5:1-25, 1996. doi: 10.1080 11008//10116600.19.1000774444000
Y. A. Kuznetsov. Elements of applied bifureation theory. Number v. 112 in Applied mathematical sciences. Springer, New York, 3rd ed edition, 2004. ISBN 978-0-387-21906-6.
Ming Lin, Rong Chen, and Jun S. Liu. Lookahead strategies for sequential monte carlo. Statistical Science, (1):69-94, 94, 2013.
David Lopez-Paz Paz and Maxime Oquabe. Revisiting classifier two-sample tests. In 5th International Conference on Learning Representations, 2017.
llya Loshchilov and Frank Hutter. Decoupled Weight Decay Regularization. arXiv preprint 1rXiv:1711.05100 2017.
Mat-Matthis Lucekmann, Giacomo Bassetto, Theofanis Karaletsos, and Jakob H. Macke. Likelihood inference with emulator networks. In Symposium on Advances in Approzimate Bagesian Inference, 2018.
16

---

Published in Transactions on Machine Learning Research (10/2024)
Man-Matthis Lueckmann, Jan Bocits, David Greenberg, Pedro Goncalves, and Jakob Macke. Benchmarking Simulation-Based Inference. In Arindam Banerjee and Kenji Fukumizu (eds.), ) Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pp. 343-354 1511. PMLR, 2021.
A. J. McKane and T. Newman. Predator-ppl prey cycles from resonant amplification of demographic stochas- ticity. Phys. Rev. Lett., 94:218102, 2005. doi: httee//ehhseeeeee....2........
Kevin Murphy. Probabilistic Machine Learning: Advanced Topics. 2023.
Radford Neal. MCMC Using Hamiltonian Dynamics. Chapman and Hall/CRC, first edition, 2011, doi: 10.1201/b10905.
Art B. Owen. Monte Carlo theory, methods and ezamples. https://artowen.. http://oowwsssu.oo.ooo... www.ins/mm/m///mmmmoo https///artooooo.uu..ommmiiii/// 2013.
George Papamakarios, Theo Pavlakou, and lain Murray, Masked Autoregressive Flow for Density Estimation. In Aduances in Newral Information Processing Systems, volume 30, 2017.
George Papamakarios, David Sterratt, and lain Murray. Sequential Neural Likelihood: Fast Likelihood- free Inference with Autoregressive Flows. In Proceedings of Machine Learning Research, volume 89, pp. 837-848, 2019.
George Papamakiros and lain Murray, Distilling Intractable Generative Models. NeurIPS workshop on Probabilistic Integration, 2015.
Du Phan, Neeraj Pradhan, and Martin Jankowiak. Composable Effects for Flexible and Accelerated Prob- abilistic Programming in NumPyro. arXiv preprint Xrxi:: 11554, 2019.
S. M. Ross. Introduction to Probability Models. Academic Press, 7th edition, 2000.
Tim Salimans, Andrej Karpathy, Xi Chen, and Dicderik P. Kingma. PixcICNN++: A PixeICNN Implemen- tation with Discretized Logistic Mixture Likelihood and Other Modifications. In International Conference on Learning Representations, 2017.
Madeline A. Shea and Gary K. Ackers. The OR control system of bacteriophage lambda. Journal of Molecular Biology, 181:211-230, 1985. doi: (0.1016/0022 (00000002220222222 (88336685)00066444
C. Sherlock, A. H. Thiery, G. O. Roberts, and J. S. Rosenthal. On the efficiency of pscudo-marginal random walk metropolis algorithms. Ann. Stats., 43:238 275, 2015.
Chris Sherlock. Direct statistical inference for finite Markov jump processes via the matrix exponential. Computational Statistics, 36:2863 2887, 2021. doi: 10.10077/000150002--0000
Robert H. Shumway and David S. Stoffer. Time Series Analysis and Its Applications: With R Examples. Springer, 2017. doi: 39119999777--99525555-7
S. A. Sisson, Y. Fan, and M. Beaumont. Handbook of Approzimate Bayesian Computation. Chapman and HaR//CRC, 2019.
Aaron van den Oord, Sander Dieleman, Hciga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukeuoglu. WaveNet: Generative Model for Raw Audio. arXiv preprint (rXit:1609.03/000 2016.
Aaron van den Oord, Nal Kalchbrenner, and Koray Kavukeuoglu. Pixel Recurrent Neural Networks. Maria Florina Balcan and Kilian Q. Weinberger (eds.), 33rd International Conference on Machine Learn- ing, volume 48 of Proceedings of Machine Learning Research, pp. 1747 1756, June 2016.
N. G. van Kampen. Stochaslic processes in physics and chemistry- Elsevier, Amsterdam, 1992.
17

---

Published in Transactions on Machine Learning Research (10/2024)
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Milia Polosukhin. Attention is All You Need. In Advances in Nerual Information Processing Systems, pp. 6000 6010, 2017.
J. N. Walker, J. V. Ross, and A. J. Black. Inference of epidemiological parameters from household stratified data. PLoS ONE, 12:c0185910, 2017. doi: httal///ooaall.on.....
Darren J. Wilkinson. Stochastic Modelling for Systems Biology, Third Edition. Chapman and Hall/CRC third edition, 2018. doi:
18

---

Published in Transactions on Machine Learning Research (10/2024)
A Details of Autoregressive Model
A.1 One dimensional inputs
A causal convolution layer & Raxdn + Baxdowt with convolution kernel (Wi,.., Wa) (where W, E dow xdu) and bias E Rdou given by
k-1 2i- -
f(zim)) Σ 1iss Wk j=0
To make the layer O-dependent, we augment it with an additive transformation of 0. We first calculate a context vector from using a shallow neural network (with activation o)
c(0) WWo be).
(7)
We then use linear transformation of c(0) to add 0 dependence to t:
(Z c(0)); Σ
+oc(() b.
Woc(f)
(8)
The size of the hidden layer and the output in Equation are both hyperparameters. Using activation function I I,, residual block in the CNN is given by
Cres (Zi c(0)): = z, +((2) (o (00 (o(2)s),c(())) ,(()),
(9)
where ((2) and (2) are of the form in Eqqation both having the same number of input/output channels for simplicity.
From Equation it can be seen that ((Z (())),((0)), is dependent on 2 Defining f( 0) as the composition of causal convolutions L it follows that f(ym c())); dependent on Vi-N(k-1)ss Evaluation of the CNN (V 0) and calculation of the logistic parameters from the output given in Algorithm Note that we assume that the initial observation yo is known, so that the shifts of (v 0, yo) can be caleulated in the same form the other shifts.
A.2 Multidimensional inputs
Assune y; has d components, and for simplicity assume that the output of each residual block has dp output channels, where p is positive integer. For j 1,2, define the following d x block matrix
1pxd; Opxd,
Opxd, Opxd,
Lpxd, Lpxd,
M, =
Ipxd;
1pxd,
1 1xxds_
where Maxx/Imxx are X matrices of ones/zeros and di =, de = p. In other words, M, block lower triangular matrix where all entries and below the block diagonal are 1. The masked causal convolution appropriate for the hidden layers of the CNN is then defined by
? (Zz1.. 0) Woc(O) + (M, O We)
Σ T=1
WE-1Zi-i+b.
(10)
where j 1 for the first layer, and j = for the causal convolutions within residual blocks. The mask ensures that the first block of p channels in the hidden layers of the CNN only dependent on y the
4The composition include residual blocks, where one residual black for cansal convolutions
19

---

Published in Transactions on Machine Learning Research (10/2024)
Algorithm Evaluation of (yp 0)
Input: Time series with initial value yo.., first CNN layer ł, r residual blocks 4 final CNN layer 4, neural network parameters ub. Output:
(()
Sequence of shift, scale and mixture proportions for a mixture of c discretised logistic conditionals.
Calculate context c W2 bo).
2:
Apply first layer Z1:m < E(v1m. c).
3: For k= do
Apply the 6th residual block Zin +- lres (2im, c). Apply the final layer Oun +- (Z c).
6: Calculate the logistic parameters.
7: For = 2, n do
For j = c do
Calculate the shifts in 1 8:-1.
10:
Calculate the scales sy + softplus( + 10-6
Calculate the mixture proportions way x exp(o/=?0) (
11:
12:
Return M2 82m: W2n.
second block is dependent on (u).v?) etc. In other words, if 2imm is the output of any hidden layer, then the 1+(111) ..... elements z;
S , are only dependent on through the elements ....... vo. The CNN will output 3ed output channels, where each block of channels corresponds to 3c mixture logistic paramcters. To remove dependence the jth block of output channels on y we apply a d block lower triangular mask matrix similar to My to the final layer, but where the block diagonal set to zeros:
Ozexp
0gexp 03exp
0zexp
1zexp
0zcxp
M
(11)
1aexp
1sexp
Osexp_
The same procedure holds when z/ contains categorical components, though number of output channels and the block sizes in the mask Equation 11 will need to be adjusted so that the correct number of parameters are obtained to model the categorical distributions rather than DMoLs.
B Additional Experiment Details
B.1 Hyperparameters and Implementation Details
All experiments were performed using a single NVIDIA RTX 3060TT GPU, which provides significant speed up over a CPU due to the use of convolutions. All instances of PMMH were written in the Julia programming language, and the particle filters were run on 16 threads to improve performance.
We used kernel length of for all experiments, and since all observations have relatively low to no noise, this was sufficient to ensure that p(ylyi-mm-i,, 0) p(g:ly 0) for all experiments. We used 5 mixture components in all experiments. We found that more mixture components resulted in slower training without significant increase in performance, and less would slightly worsen accuracy without much increase in speed. For the SIR/SEIAR model experiments, we used 64 hidden channels and 2/3 residual blocks. For the predator-prey model experiments, we used 100 hidden channels and 3 residual blocks. All of these hyperparameters were chosen initially to be relatively small, and were increased when necessary based on an pilot run of training, adding depth to the network over incrcasing hidden channnels, since computational cost grows quadratically with the number of hidden channels. We also used the same activation function for all experiments, namely the Gaussian linear error unit (GeLU) (Hendrycks & Gimpel 2016).
20

---

Published in Transactions on Machine Learning Research (10/2024)
For all experiments, we used 90% of the data for training and the other 10% for validation. We used Adam with weight decay (Loshchilov & Hutter 2017] as the optimiser, using a weight decay parameter of 10-5 and a learning rate of 0.0003. The only other form of regularisation which we used was early stopping, using a patience parameter of 30 for all experiments. For the SIR/ LIR/SIAR// TeRR///ddatorprrgee model experiments, we used a batch size of SIR model experiments, because general each observation was quite short and therefore sparse of data. For the SEIAR model experiments, the autoregressive model was prone to getting stuck at local minimum in the first round of training with too high batch size, presumably the stochasticity in the gradients helped to break through the local minimum.
to calculate the stochastic gradients. We used large batch size for the
For the SIR model experiment with a single observation and the SEIAR model experiments, we used a simulation budget of 60,000 samples, using 10,000 from the prior predictive in the first round, hence we drew 5,000 MCMC samples per round. For the household experiments, we used a simulation budget of 100,000 due to the shorter observations, hence we drew 9,000 MCMC samples per round. For the prodator-prey model, we used simulation budget of 25,000, using 10,000 from the prior predictive, hence we drew 1,500 MCMC samples per round. We dropped the simulation budget in the predator prey experiments, because the time series were much longer on average, and MCMC more expensive for this data. Additionally, we always drew an extra 200 samples at the beginning of MCMC, which were used in the warmup phase to tune the hyperparameters of NUTS, and then were subsequently discarded.
For the classifier 2-sample I-test values, we use the same procedure employed by Lueckmann et al. (002) namely we use a two layer network with 50 hidden units (10 times the SEIAR model parameters), a ReLU activation, and single output unit, which after applying a sigmoid function gives the probability that the sample comes from the true (PMMH) posterior. To generate the training data, we used one instance of SNL and generated an ESS of approximately 10,000 using PMMH, thinning so that the number of raw PMMH and SNL samples are approximately equal. We randomly permuted the full dataset, and saved 5000 samples to calculate the classifier score, while 10% of the remaining samples were used for validation and the rest for training. After training the classifier using the Adam optimiser, we evaluate the classifier on the 5000 samples set aside and recorded their average. We performed this 10 times with different random seed to account for biases that might have been introduced through the datasct splits, reporting the average over these 10 instances for each experiment. The p-values reported were calculated by taking the classifier score and calculating z 2/5000(s 0.5), which is standard normal random variable under the hypothesis that the two distributions are the same (Lopez-Paz & Oquabe 2017), so the p-value is simply Pr(|Z] > z) where Z 1).
B.2
SIR Model
B.2.1
The Model
The SIR model is a CTMC which models an outbrcak in homogencous population of fixed size N, classing The model has state space
each individual either susceptible, infectious removed (Allen 2017
1S, e 2 881, R. S ++I+M0 N).
and transition rates
BSI
(5,1, - (5- 1,1,,11 at rate N - 1
(12)
R), (5,1 at rate 71,
(13)
where B> 0 is the transmissibilitt and > 0 is the recovery rate. The SIR model can be recast in terms of two counting processes by using the transformations
Z S. S,
(14)
Ze R.
(15)
where Z1 and Ze count the number of times Equation 12 and Equation 13 occur respectively. The model
parameters (8,7) can be reparameterised into more interpretable quantities Ro .//y and 1/7, where
21

---

Published in Transactions on Machine Learning Research (10/2024)
the former is the average number of infections caused by single infectious individual (when most of the population is susceptible), and the latter is the mean length of the infectious period (Allen 2017
B.2.2 Autoregressive Model Tweaks
We assume that there is initial infection, and that subsequently observe the number of new cases every day (without any noise), equivalent to observing Z1 at frequency of once per day. We also assume that the outbreak has reached completion, hence the final size of the outbreak Ny can be measured. Using the data Wim,np), we can calculate the day at which the final size is reached:
T min mp).
Clearly,
p (yim, NF = nF 0) - (Wirr NF = nF 0) since case numbers can never grow beyond the final size. We can take this further and define indicator variables
f1 = 8ymups
from which it is casy to see that
plyer: NF = =p = fo:
where we need to include fo in case the outbreak never progresses from the initial infection. In the context of the autoregressive model, we decompose the conditionals p(ys, filyr. fi.s, 0) as follows
ply,, filii foss; 0) p(filyia, fo-i 0)p(3 11 foi-1,,,
p(filuia fi-1, 9 1p
since fi-1 = 0 implies that fa = for all k i-- 1. Furthermore, since fi follows that fi-1 = 0 for all S т. It would therefore be redundant to make (v fi- explicitly dependent on S1-1, since f4-1 to the CNN. We do require tweak to the output layer of the CNN, that it outputs 16 parameters per conditional, 15 for the DMoL parameters, and one for (filla, fi-1, 0). Since fi should be conditioned on Vi, we use mask on the final causal convolution layer which has a similar structure to Equation I1 in Appendix A Explicitly, the mask matrix has 15 rows of zeros, followed by a single row of ones, ensuring that only the parameter corresponding to qu (falb fi-1, is dependent on ys. and the other 15 are not.
and only if = then it
always by construction. As result, we only need to use V1rr
input
In practice, we do not need include fo in our calculation with the autoregressive model, p( can be caleulated analytically using the rates Equation 12 Equuaion 13 Karlin & Taylor 1975):
p(fo fo + Ro
(16)
Ro +
We therefore train the autoregressive model on samples conditioned on NF > 1, and define the surrogate likelihood to be
ge (V ny (0) = p (it
This prevents us from wasting simulations on time series with 0 length, which are likely when Rg small. In principle, we could also calculate p(NF t)0) for some threshold size t, and only train on simulations which cross this threshold, which would further reduce the number of short time series in the training data. This is possible largely due to the simplicity of the SIR model, and such procedure generally becomes more difficult. for CTMCs with a higher dimensional state space.
The final augmentation that we make to the autoregressive model concerns masking out certain terms of the form log qu (f Vy1 fi-1, 0), which forces (f fi-1, 9) = p(f.lyiss ft = 1. There are three scenarios where this is appropriate:
If y = Vt-1 and fa 0. as we know that NF * ys-1, and the outbreak does not process over the interval [ i, fa with probability 1.
22

---

Published in Transactions on Machine Learning Research (10/2024)
2. If 1/1 1, since we always have NB by construction.
3. If y: N, since this corresponds to the whole population becoming infected, therefore we have Np N and fi with probability 1.
We calculate the context vector c((8,)) 7)) using shallow neural network with 64 hidden units and 12 output units, using GeLU non-linearity. For the household experiments, we additionally concatenate one- encoding of the household size to (3, 7) before calculating c, so that the autoregressive model can distinguish between houschold sizes. To obtain the logistic paramcters from the output of the ncural network 01 we calculate
1 = oi_ N (1 1111----1100 o'-))u--:
5f = (N yi- yo-rssofpllus( +5)
This encourages the shifts to lie in the support of the conditional (vi-1..... NJ (since yi is increasing), and allows the scales to be learnt in proportion to the width of the support of each conditional.
B.2.3 Training and Inference Details
For inference, we used a Uni1orm 10) prior on Ro and a Gamma (10, 2) prior on 1/7, where we truncate the prior of 1/y so that it is supported on [1, oc). For the single observation experiment, we simulate the data using the Gillespie (1976) algorithm, running for a maximum length of 100 days, which ensures that approximately 99.9% of samples (estimated based on simulation) from the prior predictive reach their final size. We construct binary mask for the simulations where T 100 so that values following 7 do not contribute to the loss during training.
For the houschold experiments, we did not need to generate training data from houscholds of size 2, as the likelihood of such observations can be calculated analytically. The case of Nr = 1 is handled by Equation 16 and if NF = 2 it can be shown that
ply or = 2,0, > = e (1 (1-e - -3--)))))
The training data generated from five realisations for a single value of 0, using a unique household size (between 3 and 7) for each of the realisations. We found that it was useful to increase the variance of the proposals generated from MCMC, as the proposals had small variance and sometimes this seemed to reinforce biases over the SNL rounds. We inflated the variance the proposals by fitting a diagonal covariance Gaussian to the proposed samples (centred to have mean 0), and added noise from the Gaussian to the proposals, effectively doubling the variance of the proposals over the individual paramcters.
For all experiments on the SIR model we, used the particle filter given in Black (2019) within PMMH. For the single observation experiment, which was sufficient for this experiment since the particle filter works particularly well on the SIR model. For the household experiments, we used 20, 40 and 60 particles for 100, 200 and 500 households respectively. We ran the MCMC chain for long enough as necessary to obtain an ESS of approximately 1000 or above for both parameters. Mixing of the MCMC chain slows down with incrcasing houschold numbers due to increasing variance in the likeliool estimate that must be offset by also increasing the number of particles. The covariance matrix for the Eatropolis-Hasttnng proposal was tuned based pilot run for all experiments.
used 20 particles, and the MCMC chain was run for 30,000 steps,
B.3 SEIAR Model
This model is another outbreak model, and is described in Black 2019). It is a CTMC with state space
((4 72 Z3, Zz) 7 Z
Zg >
Zz
2 0, 7 > 2 2 0, 7, + + Z5 =N},
5 Note that we input y0.r inta the neural network we can caleulate %s (v 10).
23

---

Published in Transactions on Machine Learning Research (10/2024)
and transition rates
rate (N Z)p Zz3) B. Za))
(2,, (2,,,,2,
43)
(17)
+ 1,,,,,,,,,,,,,85
at rate ga - -2- -Zs).
(18)
(Z1, Z2, Z3, 73, 75)
(2,,2,23+1,24,23
rate 7 7z);
(19)
73, 75)
74,
rate 7(3g 74);
(20)
(21, L1, [s)
+1 at rate (1 () Zz -)
(21)
The event Equation 17 represents a new exposure, Equation [8| represents an exposure becoming infectious (but pre-symptomatie) Equation 19 represents a pre-symptomtic individual becoming symptomatie, Equa- tion 20 represents recovery, and Equation 21 represents an exposure who becomes asymptomatio (and not infectious). The parameters of the model are the transmissibilitie attributed to the pre-symptomatic and
symptomatic population Sp, 3. > 0, the rale of leaving the exposed class > 0, the rate of leaving the (pre)symptomatic infection class y > O, and the probability of becoming infectious q E (0, 1). We fol- low Black (2019) and reparameterise the model in terms of the parameters Ru = (3p + 3.)/g7, 1/a, 1/7, pp B.) and q- We use the same parameters to generate the observations and priors that were used
in Black 2019 Rg 2.2, 1/y 1, 1/or 1, 0.7 and 0.9, and priors
Ro Uniform (0.1, 8),
1
Gamma 10, 10),
~ Gamma 10, 10),
K ~ Uniform (0, 1),
q ~ Uniform 1).
truncating the priors for 1/0 and 1/4y to have support [0.1, and [0.5, oo) respectively.
The observations are given by the values of Z3 at discrete time steps, corresponding to the number of new symptomatic individuals per day, as well as the final size Np. We handle the final size data in the same way discussed for the SIR model, but we include fu for this model, p(Nr 1|0) is more difficult to calculate due to the asymptomatic individuals. The initial state is taken to be (Zi, Zz, Za, Z4, Zs) (1, 1, 0, 0), which we simulate forward using the Gillespie algorithm until Z3 = 1, after which we begin recording the values of y at integer time steps. For population sizes of 350, 500, 1000 and 2000, we generate time series of maximum length 70, 80, 100 and 110 respectively. This ensured that at least 99% of samples from the prior predictive reached their final size (estimated by simulation). The autoregressive model follows the exact same construction as given for the SIR model experiments, except that we rescale the time series by the population size before inputting to the CNN. We also divide Vir by N before inputting to the CNN, which helps to stabilise training.
For PMMH, we used the particle filter given in Black 2019)) using 70, 80, 100 and 120 particles for population sizes of 350, 500, 1000 and 2000 respectively. The MCMC chain was run until all parameters had an ESS of approximately 1000 or above, and the proposal covariance matrix was estimated based on a pilot run.
B.4 Predator- Prey Model
The predator-prey model that use is described in McKane & Newman 2005 For large populations in predator-prcy models, it is often sufficient to consider deterministic approximation via an ODE. For this model, the corresponding ODE has a stable fixed point, which makes it incapable of exhibiting the cyclic behaviour seen in e.g. the Voltera model (Wilkinson) 2018). In the stochastic model however, there is a resonance effect which causes oscillatory behaviour. Unlike predator-prey models with limit cycle behaviour, in this model the dominant frequencies in the oscillations have a random contribution to the trajectory, which causes greater variability among the realisations for fixed set of parameters.
24

---

Published in Transactions on Machine Learning Research (10/2024)
The model is CTMC which represents the population dynamics of predator population of size P and prey population of size Q. The state space given by
((P.Q) E z² P.OK, P+ KK]
where 2 0 is a carrying capacity parameter.
The transition rates are given by
(P₁Q) (P rate di P.
(22)
Q
(P,Q) + (600)((P)) (P.Q+1) at rate 26
(K-
P Q),
(23)
K
PQ
(P,Q) - (P.Q at rate 2p2
d2Q,
(24)
K
PQ
(P,Q) (P 1, 1, al rate 2pt
(25)
K
where b, dj, 11.22.11 P1.PP P2 > are the parameters of the model. To simulate the data, we use parameters = 0.26, d = 0.1, d2 = 0.01, P1 = 0.13, P2 = 0.05. We use informative priors on b, d1 and d2, assuming that reasonable estimates can be obtained based on observing the two individual species. In particular we use
log N (log 0.25, 0.25),
log d1 N (log 0.1, 0.5),
log d2 (log 0.01, 0.5),
truncating each prior to have support (-o0, 0] in the log domain. The effect of this truncation minimal: the widest prior, for log(b), has upper bound -1.2 for the 99% confidence interval. For the two predation parameters, we use priors
P1 (x0.nential(..1),
P2 Exponential(0.05).
where we truncate the prior of pi that it is supported on [0.01, xo). This truncation, which removes only 5.4% of the mass from P1, rules out population explosion in the prey due to unphysically low predation. Furthermore, when p1 is close to 0, 25 becomes negligible and the resulting model does not behave like a true predator-prey model, so truncating the prior rules out this possibility. To simulate the realisation for the observations, we use a carrying capacity of K 800 and initial condition (250, 250). We then run the Gillespie algorithm over 200 time steps, taking y (P2., Q2i), and apply observation error to generate the obscrvations, namely
D 7)
1
V Binomial(y'.::7)
for 0.5, 0.7, 0.9.
To generate the training data, we firstly infer the initial condition from the noisy initial observation yo- In particular, we assume p(P)) x 1 and p(Qo) x 1, from which it is easy to derive the posterior p((Po, Qo)33)):
(Po vo) yo ~ Nogativeeinoomillgu. + 1,r),
(Qo vd) yo NegativeBinnomia(pp. + 1,r).
We then generate the training data in the same way we produced the observations. No tweaks are made the autoregressive model, though we do divide 91:100 by 60 before inputting to the CNN, and multiply the shifts (before adding B:-1) and scales by 60. We do this because the standard deviation across the time axis for samples from the prior predictive was estimated to be approximately 60. We also tried rescaling the imput/outpnn by larger values, but this generally resulted in less stable training and slower MCMC iterations. We do not centre the inputs, bocause if y! 0 for multiple time steps in a row, then this implies that extinction has likely occurred, hence zeroing the weights associated with yf in the first layer
25

---

Published in Transactions on Machine Learning Research (10/2024)
should correspond to the model learning the dynamics of the individual species. If centred the data, then extinctions would not correspond to zeroing the weights. We also attempted to order the prey before the predators in the autoregressive model, but found that there was not large difference between the two orderings.
For PMMH, we used a bootstrap particle filter (Kitagawa 1996). using 100, 200 and 250 particles for 0.5, 0.7 and 0.9 respectively. We used 180,000 MCMC iterations for each experiment, which produced ESS of greater than 500 for all parameters. Generating large ESS for this data is difficult due to the extremely high posterior correlation between d1 and P1, and the computational time per MCMC iteration is quite high for these experiments. We estimated the covariance matrix for the MCMC kernel from prior SNL run, and tuned further based on a pilot run. Without the SNL run to inform the initial covariance matrix, tuning the MCMC proposal would likely take huge amounts of computational time, since a proposal with diagonal covariance would mix extremely slowly with the highly correlated parameters in the posterior.
C
Additional Results
C.1 Posterior Metrics
Table 5: Comparison of SNL and PMMH posterior statistics for the SIR model.
Parameter
M
S
avg
max
avg
max
Ro
0.00007
-0.0034
0.036
0.068
1/4
-0.029
-0.087
0.019
0.073
Table 6: Comparison of SNL and PMMH posterior statistics for 1/0, 1/7 and in the SEIAR model.
Parameter
N
M
avg
max
avg
max
350
0.04
0.10
0.03
0.07
500
0.01
0.08
-0.01
-0.05
1/a
1000
0.02
0.08
0.02
0.05
2000
-0.007
-0.05
0.04
0.07
350
0.07
0.20
0.03
0.06
500
0.08
0.30
0.03
0.08
1/7
1000
0.05
0.17
0.03
0.06
2000
0.04
0.10
0.03
0.10
350
-0.03
-0.06
0.02
0.06
500
-0.05
-0.12
-0.02
-0.12
1000
-0.07
-0.12
-0.02
-0.06
2000
-0.02
-0.07
0.06
0.08
26

---

Published in Transactions on Machine Learning Research (10/2024)
Table 7: Comparison SNL and PMMH posterior statistics predator prey model experiments.
Parameter
M
S
avg
max
avg
max
0.5
-0.09
-0.22
-0.02
-0.10
b
0.7
-0.05
-0.11
-0.009
-0.06
0.5
-0.08
-0.18
-0.01
-0.08
d1
0.7
-0.003
-0.02
0.03
0.08
0.5
-0.00001
-0.14
-0.08
-0.16
d2
0.7
0.04
0.12
0.05
0.17
0.5
-0.06
-0.13
-0.01
-0.09
0.7
-0.0003
-0.01
0.03
0.08
0.5
0.02
0.13
0.02
0.07
P2
0.7
-0.06
-0.08
-0.11
-0.17
C.2
Runtimes
Table 8: ESS/s values from PMMH versus SNL for household data experiments. Recall that we report the ESS/s value for SNL by averaging the ESS over the final rounds.
h
PMMH ESS/s
SNL ESS/s
avg
min
max
100
2.33
3.67
3.37
4.00
200
0.70
3.28
3.03
3.44
500
0.25
2.38
2.25
2.57
Table 9: ESS/S values from PMMH versus SNL for SEIAR experiments. Recall that report the ESS/s value for SNL by averaging the ESS over the final 5 rounds.
N
PMMH ESS/s
SNL ESS/ /s
avg
min
max
350
1.94
2.55
2.24
2.96
500
1.65
1.78
1.60
1.90
1000
0.50
1.50
1.29
1.62
2000
0.20
1.23
1.08
1.39
Table 10: ESS/s values from PMMH versus SNL for predator- prey experiments. Recall that we report the ESS/s value for SNL by averaging the ESS over the final rounds.
r
PMMH ESS/s
SNL ESS /
avg
min
max
0.5
0.034
0.33
0.29
0.37
0.7
0.027
0.51
0.48
0.56
0.9
0.015
0.45
0.34
0.50
27

---

Published in Transactions on Machine Learning Research (10/2024)
D Additional Plots
340
Predator
Prey
320
300
Counts
280
OK
260
240
220
0
25
50
75
100 Time
125
150
175
200
Figure 6: Predator and prey counts for = 0.9. Both populations oscillate around an equilibrium in somewhat irregular fashion.
Prior
SNL PMMH
Ro
Truth
1/y
10
2
8
Figure 7: Posterior pairs plot for single observation from SIR model experiment.
28

---

Published in Transactions on Machine Learning Research (10/2024)
Prior
SNL
Ro
PMMH
Truth
1/y
10
Figure 8: Posterior pairs plot for observed data from the SIR model across 100 households.
Prior
SNL
Ro
PMMH
Truth
1/y
10
2
8
Figure 9: Posterior pairs plot for observed data from the SIR model across 200 households.
29

---

Published in Transactions on Machine Learning Research (10/2024)
Prior
SNL
Ro
PMMH
Truth
1/y
10
Figure 10: Posterior pairs plot for observed data from the SIR model across 500 households.
Rp
PMMH
Prior
Truth
1/o
1/v
0.5
1.0
15
0.00
0.25
0.50
0.75
1.00
0.6
0.8
Figure 11: Posterior pairs plot for observed outbreak from SE1AR model with population size N 350.
30

---

Published in Transactions on Machine Learning Research (10/2024)
SNL
PMMH
Prior
Truth
1/a
1/v
K
0.5
1.0
1.5
2.0
0.00
0.25
0.50
0.75
1.00
0.6
0.8
1.0
Figure 12: Posterior pairs plot for observed outbreak from SEIAR model with population size N = 500.
Rp
PMMH
Prior
Truth
1/o
1/v
K
0.5
1.0
15
0.00
0.25
0.50
0.75
1.00
0.6
0.8
Figure 13: Posterior pairs plot for observed outbreak from SEIAR model with population size N 1000.
31

---

Published in Transactions on Machine Learning Research (10/2024)
SNL
PMMH
Prior
Truth
1/a
1/v
K
0.5
1.0
1.5
2.0
0.00
0.25
0.50
0.75
1.00
0.6
0.8
1.0
Figure 14: Posterior pairs plot for observed outbreak from SETAR model with population size = 2000.
SNL
PMMH
Prior
Truth
d
d2
0.0
0.1
0.2
0.3
0.00
0.02
0.04
0.2
0.0
0.1
Figure 15: Posterior pairs plot for predator prey experiments with r 0.5.
32

---

Published in Transactions on Machine Learning Research (10/2024)
SNL
b
PMMH
Truth
d1
d2
Pj
0.0
0.2
0.4
0.0
0.1
0.2
0.3
0.00
0.02
0.04
0.2
0.4
0.0
0.1
0.2
Figure 16: Posterior pairs plot for predator prey experiments with r = 0.7.
SNL
b
Prior
Truth
d1
d2
0.0
0.1
0.3
0.00
0.02
0.04
0.2
0.0
0.1
Figure 17: Posterior pairs plot for predator prey experiments with 0.9.
33


Reviews:
[{"review_id": "HPP1SzSxI8", "reviewer": "Reviewer_KPNs", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strength**\n\n- The SBI literature on time-series data with discrete observation spaces is still limited while providing important problems. Therefore,\nI believe that providing a neural architecture for this specific set of tasks is a good addition to sequential neural likelihood estimation to alleviate the need to represent inter-valued time-series by summaries. \n\n- The tasks are well suited and show that SNL with the proposed architecture approximates the PMMH posterior well. In particular, SNL does not drop performance in terms of the ESS/s when dealing with increasing number of observations and longer time-series.\n\n**Weaknesses**\n\n- Given that one of the paper's major contributions is the introduction of an autoregressive model for modeling the CTMC, the experiments lack an ablation study to demonstrate the advantages of this autoregressive approach compared to a standard feed-forward network.\n\n- The extension to multivariate data seems not to be straight forward and and introduces potential biases related to the ordering of observations.\n\n- The evaluation metrics chosen are somewhat limited. Given that the posterior results of PMMH are treated as reference data, it would be beneficial to include additional metrics [1] such as the classifier 2-sample test (C2ST) [2], which is widely used. Moreover, considering that the pair plots for the SEIAR task reveal skewed marginal densities with high pair-wise correlations, relying solely on posterior means and standard deviations may not provide the most accurate assessment.\n\n[1] Lueckmann, J. M., Boelts, J., Greenberg, D., Goncalves, P., & Macke, J. (2021, March). Benchmarking simulation-based inference.\n[2] Lopez-Paz, D., & Oquab, M. (2016). Revisiting classifier two-sample tests.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents an enhancement to simulation-based inference methods for discrete time-series data through the introduction of an autoregressive neural likelihood model. This model leverages causal convolutions to maintain the autoregressive nature of a continuous-time Markov chain and utilizes a discrete mixture of logistics to model the intractable likelihood.\nThe authors evaluate their approach on two epidemiological tasks and a predator-prey model involving discrete time-series observations. When compared to particle marginal Metropolis Hastings (PMMH) in a conditional inference setting, the proposed method demonstrates comparable posterior expressivity while offering superior performance in terms of the ESS/s metric.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- I would appreciate experiments that highlight the benefits of the introduced autoregressive model/causal convoluations over a standard MLP.\n- Eq(4): To my understanding, the DMoL operates with a probability mass function rather than a density function, as does the DL. In general, it should be made more clear when random variables are discrete or continuous. \n- I understand the inclusion of ESS/s as a performance metric, but I believe that ESS is a particularly significant aspect of particle-based inference methods and should be reported explicitly. Additionally, it is difficult to estimate inference time directly from the ESS/s plots. Therefore, I would appreciate it if both of these quantities could be presented.\n-  > The work most similar to ours is described in the original SNL paper Papamakarios et al. (2019), but a limitation of that is the use of summary statistics (as in an ABC approach) for training the model.\n\n   While I generally agree with your statement, it’s worth noting that results of SNL on discrete time-series observations for the SIR task were reported in [1]. Despite the network architecture not being specifically tailored for discrete data, the posteriors closely matched the reference posterior data. I would be interested if you deem a comparison to the vanilla SNL be worth and if so, I would be curious how the proposed architecture compares against it.\n- Please define $\\mathbb{Z}$.\n- It would be great if the authors could report on the simulation budget of the PMMH.\n\n**I further list my open questions here:**\n1. As someone with limited experience in CTMCs, I'm curious about your decision to use the DMoL parameterization instead of the more commonly used tokenization and embedding methods typical in transformer-based architectures or in TCN [3] in which the output would be a softmax over the state-space. Could you elaborate on this choice, especially considering the complexities involved in handling multivariate data?\n2. > adding depth to the network over width since computational cost grows quadratically with width\n\n   Shouldn’t width be only increasing in space complexity, while depth should increase the time complexity?\n3. Could you clarify how the ensemble of SNL posterior is utilized? Specifically, is the ensemble treated as a mixture of posterior distributions from which samples from one of the posteriors are drawn with equal probability?\n4. How can the ESS/s of PMMH and SNLE be compared if the codebases are different (Julia vs Jax)?\n6. > Studying how robust our method is to model misspecification would also be important future work, as any real data is unlikely to truly follow a specific model.\n\n   Since vanilla SNL has been demonstrated to struggle with model misspecification [2], I’m curious to know if you anticipate greater robustness with the added structure in your approach.\n\n[1] Lueckmann, J. M., Boelts, J., Greenberg, D., Goncalves, P., & Macke, J. (2021, March). Benchmarking simulation-based inference.\n[2] Kelly, R., Nott, D. J., Frazier, D. T., Warne, D., & Drovandi, C. (2024). Misspecification-robust Sequential Neural Likelihood for Simulation-based inference.\n[3] Bai, S., Kolter, J. Z., & Koltun, V. (2018). An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.", "broader_impact_concerns": "I don't have any broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "bEAd8vGIaG", "reviewer": "Reviewer_vfTp", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The paper is well-written. The method is explained clearly and appears sound.\n- The method is intuitive and outperforms previous state-of-the-art in terms of scalability, and thus could be a useful addition for such applications.\n\nWeaknesses:\n\n- The behavior of the method depends on an ordering of observation dimensions.\n- While this is not necessarily a weakness given the TMLR policies, I am somewhat uncertain about the relevance of the submission in ML venues.  The method appears obvious given the constraints of the problem, and while the model studied can certainly be important in some application fields, to me it is less clear if they are of major importance to a general ML audience, so it is not clear how they could benefit from seeing this paper in this venue.  I would note that I do not work in this specific area and am open to the opinions of other reviewers.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work studies the inference of model parameters in integer-valued state space models. The authors proposed a method based on the combination of neural conditional density estimation (NCDE) and sequential neural likelihood (SNL) inference, and designed a causal convolution model for NCDE.  Simulation experiments shows that the method outperforms particle marginal Metropolis-Hastings (PMMH), a previous state-of-the-art for such models, in terms of computational efficiency.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "As the behavior of the method depends on the ordering of output dimensions, it may be helpful to study the sensitivity empirically.\n\nMinor issues:\n\n- There are several places where \\cite or \\citet should be \\citep.  E.g. the Bickel et al (2008) on p.1.\n- Appendix D: title should be capitalized consistently.\n- p.5, \" rounding to down to the nearest integer\" -- should be \"rounding down to\"?", "broader_impact_concerns": "There is no such concern.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "2AkYgqJ1Sr", "reviewer": "Reviewer_7xDi", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\nS1: the overall approach seems both clever and relatively novel approach to likelihood approximation in small population models, and the differentiability aspect makes it amenable to using plenty of exciting MCMC methods.\n\nS2: The method demonstrates good scalability with increasing data size and system complexity, cleanly outperforming the state-of-the-art (at least wrt. PMMH and ESS, which are good baselines).\n\nS3: The experimental section showcases a good degree of versatility, as SIR, SEIAR, and predator-prey are fairly different settings and overall the proposed method seems to be straightforwardly performant.\n\nS4: I strongly appreciate that the focus is on raw time series data. Summary statistics are often extremely hard to get right in real-world data, and are a significant drawback of more commonly used methodologies in the literature.\n\nS5: The neural likelihood approximation is differentiable with respect to model parameters, enabling the use of efficient MCMC methods like Hamiltonian Monte Carlo.\n\n## Weaknesses\n\nW1: The choice of priors seems to have a significant impact on the method's performance, and the paper doesn't provide clear guidelines on how to select appropriate priors. The authors mention that:\n\n>In practice, we found the choice of priors to be important. For example, in the predator prey experiments, sometimes the MCMC chain would wander into regions of very low prior density, resulting in a degenerate proposal. This was easily fixed by truncating the priors to only allow parameters in a physical range.\n\nI would have appreciated some kind of somewhat extensive sensitivity analysis wrt. this aspect. I agree that in principle it's usually not too hard to pick the right prior models, however I wonder whether there are clear failure cases that can be reported and (at least superficially) studied.\n\nW2: There's a similar argument (as W1) to be made for the ordering of the components of $y_{i}$. Are there simple heuristics that could help alleviating the somewhat arbitrary hyperparameter optimization process here?\n\nW3: SIR, SEIAR, and predator-prey don't really feel high-dimensional (though admittedly they seem to be extremely relevant within their own application domains, at least for SIR and SEIAR). Even though the manuscript makes some statements about scalability wrt. model parameters, it's not clear to me whether there are any traps that would arise from living in a highly-multivariate model.\n\nW4: While the method shows good performance in reproducing posterior means, there are instances of inflated posterior variances (which the authors acknowledge), particularly for some parameters in the SEIAR model. Discussing how to more accurately capture the tail behavior of skewed posteriors would nicely enhance the manuscript. This also seems to be a problem that shows up primarily in SEIAR, so I wonder whether we'd see more of this for more complex models.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This manuscript proposes a neural likelihood approximation method for performing Bayesian inference on integer-valued time series data, particularly focusing on stochastic processes in physical and biological sciences (with relevant constraints regarding dynamics of model population). The authors construct an autoregressive model using CNNs to approximate the likelihood of the data. This approach allows for efficient inference using unconditional simulations of the underlying model, which is simpler and computationally less expensive than current simulation-based inference methods that require conditional simulations.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "C1: I'd love to see a section on guidance on prior selection wrt. SNL, and its impact on the method's performance, possibly including a sensitivity analysis.\n\nC2: I think it'd be great to probe the method for performance (and issues!) on significantly more complex models. I don't know what the literature in epidemiology looks like, but I presume there exist some. if not, might be good to look at earth-models from the econometrics / ecology / climate literature -- e.g. the AYS AYS model from Kittel et al., 2017, though that still somewhat simple)\n\nC3: I think an analysis on the dynamics of feature learning wrt. feature ordering wrt. the presented problems in the manuscript would be neat to see.\n\n### Typos\n\nT1 (p2): becasue -> because\n\nT2 (p7): autoregresssive -> autoregressive", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Based on the positive recommendation of all three reviewers, i suggest acceptance.
Decision: Accept (poster)