Fast, accurate and lightweight sequential simulation-based inference using Gaussian locally linear mappings
Authors: Henrik Häggström
Pedro L. C. Rodrigues
Geoffroy Oudoumanessah
Florence Forbes
Umberto Picchini
Venue: Accepted by TMLR ()

Abstract: Bayesian inference for complex models with an intractable likelihood can be tackled using algorithms performing many calls to computer simulators. These approaches are collectively known as "simulation-based inference" (SBI). Recent SBI methods have made use of neural networks (NN) to provide approximate, yet expressive constructs for the unavailable likelihood function and the posterior distribution. However, the trade-off between accuracy and computational demand leaves much space for improvement. In this work, we propose an alternative that provides both approximations to the likelihood and the posterior distribution, using  structured mixtures of probability distributions. Our approach produces accurate posterior inference when compared to state-of-the-art NN-based SBI methods, even for multimodal posteriors, while exhibiting a much smaller computational footprint. We illustrate our results on several benchmark models from the SBI literature and on a biological model of the translation kinetics after mRNA transfection.

Published in Transactions on Machine Learning Research (06/2024)
Fast, accurate and lightweight sequential simulation-based inference using Gaussian locally linear mappings
Henrik Häggström
hehhhers@chammers.see
Dept. Mathematical Sciences
Chalmers and Universily of Gotkenburg
Pedro L. C. Rodrigues
pedro.rodrriueeehhvrii..
Unio. Grenoble Alpes, Inria, CNRS, Grenoble INP, Inserm U1216, CHU Grenoble Alpes, Grenoble Institut des Neurosciences, LJK
Geoffroy Oudoumanessah
Geoffroy Oudoumanessahh@iaaii...
Unio. Grenoble Alpes, Inria, CNRS, Grenoble INP, Inserm U1216, CHU Grenoble Alpes,
Grenoble Instifut des Neurosciences, LJK Unio. Lyon, CNRS, Inserm, INSA Lyon, UCBL, CREATIS
Florence Forbes
Horenca..orbee@inrii..
Unio. Grenoble Alpes, Inria, CNRS, Grenoble INP, Inserm U1216, CHU Grenoble Alpes, Grenoble Institut des Neurosciences, LJK
Umberto Picchini
picchini@@chabmers.es
Dept. Mathematical Sciences,
Chalmers and University of Gothenburg
Reviewed on OpenReview: //tps://opp openreview. net/ forum? id= QOnzpRcwkk
Abstract
Bayesian inference for complex models with an intractable likelihood can be tackled using algorithms performing many calls to computer simulators. These approaches are collec- tively known simallation-base inference (SBI). Recent SBI methods have made use of neural networks (NN) provide approximate, yet expressive constructs for the unavailable likelihood function and the posterior distribution. However, the trade-off between accuracy and computational demand leaves much space for improvement. In this work, we propose an alternative that provides both approximations to the likelihood and the posterior dis- tribution, using structured mixtures of probability distributions. Our approach produces accurate posterior inference when compared to state-of-the-art NN-based SBI methods, even for multimodal posteriors, while exhibiting much smaller computational footprint. We il- lustrate results several benchmark models from the SBT literature and on biological model of the translation kinetics after mRNA transfection.
1 Introduction
We propose novel methodology achieving state-of-the-art Bayesian inference for models with intractable likelihood function. We illustrate, via several benchmark models, how our methodology returns accurate in- ference, even in presence of multimodal posteriors, with a small runtime. Morcover, for numerous simulation studies and applications, machine and statistical learning solutions show an ever growing performance, but are not generally evaluated terms of resource and energy consumption. Although their environmental impact more and more documented (Schwartz et al.] 2020 Strubell et al.) 2019 Thompson al. 2022). the question of resource efficiency and energy consumption when learning model parameters not as much studied the quality of the inference itself, due to the multiple aspects to be taken into account and to the

---

Published in Transactions on Machine Learning Research (06/2024)
lack of clear metrics that limit the algorithms evaluation. When comparing with widely employed inference engines based on neural-networks, our method considerably lower the energy and memory required for the inference task, while showing competitive inference quality.
We consider simulation- based inference (SBI) for model parameters 0 in stochastic modelling. SBI refers to approaches that bypass the need to have a likelihood function p(y]0) (for data y) that is available in closed form, or one that is cheap to approximate, by instead considering a generative model f(0) or computer simulator. The simulator depends on 0 and possibly on other imputs such as covariates and stream of pseudorandom numbers. Generating dataset y by running the simulator y² f(0)) at 0 = 0" is then equivalent to simulating y" from the unavailable likelihood function. In the following this simulation is written as y² ~ p(y) Established SBI methods include approximate Bayesian computation (ABC, reviewed in [isson et al., 2018), synthetic likelihoods (Wood 2010 Price et al., 2018), particle Markov chain Monte Carlo (Andricu & Roberts 2009 Andrieu et al., 20100. and several more that we discuss in Section 2
More recently, methods exploiting neural networks (NN) have gained considerable attention. Neural- basod istimators have been used to produce approximated likelihoods or approximated posterior distributions,
typically using (discrete) normalizing flows (Rezende & Mohamed, 2015
Papamakarios et al.
2021 NN
have been used to approximate the likelihood function in implicit models (Papamakarios et al.
2019 Chen
et al.
2021). the posterior distribution (Papamakarios & Murray 2016
Greenberg et al. 2019] Durkan et al.
2020
Chen et al. 2021) or the likelihood and the posterior distribution simultaneously [iiqvist et al. 2021
Radev et al. 2023a)..
These neural-based posterior and likelihood estimation methods have produced exceptionally flexible infer- ence strategies for models with intractable likelihoods. In this work, we propose to investigate more frugal strategies (Evchenko et al. 2021] that can run with limited resources and much smaller computational footprint, while returning inference of similar quality as coural-based approaches. For SBI, we propose to investigate structured mixtures of probability distributions whose versatility has been widely recognized for a varicty of data and tasks, while not requiring excessive design effort or tuning. Their expressivity makes them good candidates to account for complex multivariate models both for the likelihood and posterior dis- tribution. Moreover, mixture models have a much smaller number of parameters & (see Table [ compared to NN, and this makes them more amenable to interpretation and efficient lcarning. More specifically, we design new simulation-base inference method named "Sequential Mixture Posterior and Likelihood Esti- mation" " (SeMPLE). SeMPLE learns approximations for both the posterior p(0ly) and the likelihood p(y) using a structured mixture model obtained via the Gaussian Locally Linear Mapping (GLLiM, Deleforge ct al. 2014) method. Interestingly, both approximations are jointly returned, simultancously, following a run of an ExxpetationnMaxxiiiaaain (EM) algorithm within GLLiM. We then use the GLLiM approximate posterior as well-informed proposal distribution in tuning-free free Settropolis-Hastinge sampler. The pro- cedure is embedded in a sequential strategy, which we compare with state-of-the-art NN-based sequential algorithms for posterior inference, namely SNL (Papamakarios et al. 2019) and SNPE-C (Greenberg et al., 2019). Other algorithms could be considered for comparison, however we stick with SNT, and SNPF-C and motivation for this is given in Section 2 Our comparisons are in terms of inference quality as well as resource and energy requirements. We did not compare against ABC samplers such as SMC-ABC except. for the example where such comparison was more suitable, as in Lueckmann et al.) 2021 it was shown that standard implementations of ABC approaches usually require a larger number of model simulations, compared to NN-based approaches, although recent research improved the SMC-ABC performance (Picchini & Tambboriioo 2024). We found that, for a given computational budget specified by the total number of model simulations, our method is cheaper to run and produces accurate inference thanks to (i) the informed proposal sampler created by ScMPLE and (ii) the fast model training. As an illustration, Figure [ shows inference for three datasets of the Two Moons model, a model with bimodal posterior that we further analyse in Section 5 The runtime, energy (consumed by the CPU and the DRAM) and memory requirements for SeMPLE are considerably lower, see Table 1 Figure 2 (right) and section 5.5 for details on how energy and memory requirements were measured. In other examples, (e.g. Figure 3 the speedup provided by SeMPLE is even more dramatic (for the model in Supplementary section
G.4 we show that SeMPLE is 74 times faster
then SNL and 26 times faster than SNPE-C). Overall SeMPLE shows low resources costs that are always of the same magnitude (<< 10 kJ and < 1 GB), while the requirements for SNL and SNPE- C vary much more
2

---

Published in Transactions on Machine Learning Research (06/2024)
with the model dimension or the inference complexity. These results highlight the advantage brought by SeMPLE, which can be run a minimal configuration while maintaining good performance.
SoMPLE
obzervation
0
diservation
observation
Figure 1: Two Moons model: inference using datasets #1 (top), #2 (middle) and #3 (bottom) from SBIBM. Posterior samples after 4 rounds of SeMPLE (column 1), 10 rounds of SNPE-C and SNL (columns 2 and 3). All methods used total budget of 10K model simulations. The true posteriors are in column 4.
SeiiLE rounds stre wih TRivesk
0.60
Cunulaaive slapsed time [mmiutes)
0.65
0.55
2500
5000
7500
10000
2500
5030
7550
10000
Number simulations
Number simulations
Figure 2: Two Moons. (left) Median C2ST (the lower the botter) and median cumulative runtime in minutes (right) for 10 runs with different data sets the number of model simulations. Shaded bands enclose the min and max values. C2ST is the two-samples classifier test, taking values in [0.5,1], see Section
Section 2 gives a brief overview of the main SBI approaches. Section describes how GLLiM can be used to construct surrogate posteriors and surrogate likelihoods, while the complete ScMPLE methodology is explained in Section 4 Results from SeMPLE, SNL and SNPE-C are then illustrated in Section 5 Test models, partly reported in the Appendices, include simulation studics comprising multimodal surfaces (two moons, multiple hyperboloid), stochastic differential equations models (Bernoulli GLM), a Markov jump process with additional measurement error (Lotka- Volterra), highly correlated posteriors with paramoters (twisted--rior), and a challenging and partially-observed biological model of the translation kinetics after mRNA transfection. Code is available at https://githuu com/ Senhagg/semple
generalised linear
3

---

Published in Transactions on Machine Learning Research (06/2024)
Table 1: Median resource requirements (across 10 independent repetitions) and number parameters & for posterior inference with budget of 10K (two moons and Bernoulli GLM) and 40K (hyperboloid and Prustein-Chlenbechh model simulations. Best (lowest) values are in bold. SNPE-C and SNL used the default SBIBM sctup described in Appendix A
Method
Energy Consumption
Peak Memory Usage
#parameters Φ
kJJ
GBB]
Two moons
SNL SNPE-C
24.5
19
28,020
15
9
60,925
SeMPLE
7.1
0.71
449
Hyperboloid
SNL SNPE-C
54
30
30,020
120
15
66,925
ScMPLE
3.2
0.88
1,479
Bernoulli GLM
SNL SNPE-C
108
75
36,100
41
10
98,025
ScMPLE
4.8
0.93
2,309
Ornstein-Uhlenbeee Uhlenbeck
SNL SNPE-C SeMPLE
66
39
41,030
134
15
102,242
6.1
0.93
30,799
2
Related work
Simulation-basee inference, also referred to as likelihood-free inference, is reviewed e.g. in Cranmer et al.) (2020) and Pesonen al. (2023). In Bayesian framework the focus on posterior estimation, which SBI directly relates to the quality of the approximations made to the unavailable likelihood, the postcrior distribution or both. This aspect raises the issue of how evaluate the quality of posterior inference when reference ground truth posterior is absent by definition. This is still an open question and is further discussed in Section 6. A second important ingredient is sampling efficiency, as simulation-based methods are typically relying on exploration, and on Monte Carlo principles that can be prohibitively expensive. Ac- way to address sample efficiency by gradually discovering and guiding where simulations should be made for a specific given observation. Our proposed method is sequential and we fo- cus comparisons with NN-based sequential methods for posterior inference, namcly SNL (Papamakarios et al. 2019 and SNPE-C (Greenberg et al. 2019 Many other approaches exist for Bayesian inference with intractable likelihoods, such as sequential neural ratio estimation (SNRE, Hermans et al. 2020), the several variants of ABC such as SMC-ABC TToni et al. 2009 Del Moral et al. 2012), synthetic likelihoods (Wood, Andricu et al., 2010). Comparisons between several the above mentioned methods are available previous works. For example, Greenberg et al. 2019) show that SNPE-C is computationally more efficient than the predecessor: SSPE-A (Papa- makarios & Murray 2016) and SNPE-B (Lueckmann et al., 2017), where efficvend means that, for given computational budget, SNPE-C approaches the true posterior using smaller number of model simulations. Lueckmann et al. (2021) show that, for nontrivial caso-studies, both SNPE-C and SNRE are more efficient than SNL and SMC-ABC. Moreover, SNPE-C and SNRE perform similarly. The above motivates why in this work we focus on the comparison with SNPE-C using the SBIBM package implementation (Lueckmann) Lotka-Volterra
live or sequential learning
2010 Price et al. 2018) and particle MCMC (Andrieu & Roberts 2009
et al. 2021 (an exception is the Lotka- Volterra Markov jump process, where we compare against SMC-ABC, as an objective comparison with NN- methods is difficult here, see section 5.3 for more details). Since
The acronyms SNPE-A and SNPE-B follow the taxonomy Papamakarios al. 2019
4

---

Published in Transactions on Machine Learning Research (06/2024)
ScMPLE includes Markov chain Monte Carlo (MCMC) step, same SNL, it is also relevant compare against SNL, using its implementation in SBIBM, which is an ameliorated one, compared to the original version of Papamakarios et al. (2019).
Note that there are crucial differences between our SeMPLE and SNPE-A, as detailed in Appendix C Briefly, in SNPF-A a component mixtare density network (MDN) used as a proposal sampler. As such, the MDN means and covariances are parametrised via neural networks. However SNPE-A does not have a built- in strategy to learn each of the K components. Instead, the GLLiM procedure embedded into SeMPLE automatically determines the means and covariance matrices for each of the K components, thus increasing flexibility to the shapes of the posteriors that can be targeted (see also Wang et al. 2018 where MDNs are employed as proposal functions within a Gibbs sampler). The fact that SeMPLE does not use MDNs, and more generally no NN architectures are employed in SeMPLE it results in an easy to train framework, with fewer parameters to learn (see Table T. Other architectures using NN within proposal samplers are e.g. the neural proposal of Stites et al. (2021) and the mixture of Gaussian copulas of Chen & Gutmann (2019)
Recent work that is relevant for SBI, but is less comparable with work, includes amortized strategies. Amortization is a property that refers to the ability of algorithms to make use of previous computational efforts, not to repeat computationally expensive steps for each new observation considered. Examples in this class are the approaches in Radev et al. (202a), Radev et al. (003b)), Dyer et al. (2022), Gao et al. (2023
Although GLLiM is by design an amortized procedure, this property is partly lost when used in a sequential approach such as SeMPLE, as amortization and sequential learning are somewhat opposite objectives. The
next section gives a brief introduction to the Gaussian local linear mapping (GLLiM) procedure, which allows to obtain amortized surrogates of the posterior and of the likelihood in closed form.
3 Amortized surrogates for likelihoods and posteriors via GLLiM
The inverse regression approach via Gaussian local linear mapping (GLLiM) was first introduced in Deleforge et al. (2014). GLLiM is parametric mixture model on the joint distribution of (y. 0), denoted by 93 and depending on a set of paramcters Φ. For interpretation or visualization, the model can also be presented as an invertible regression model as follows. GLLiM the relationship between observations y E Red and 0 E Rd is assumed to be locally-linear with K components, defined through an additional latent variable [I,.,,,K as;
y (==k) (A=0 bn + Ex),
(1)
where 1 denotes the indicator function, and A E Rdxt and by Idd are matrices and vectors, respectively, defining the affine transformations, while & E Rdd corresponds an crror term capturing both the observa- tional noise and the modelling error due to assuming an affine approximation for the data. In the following we consider Gaussian noise, ég ~ (0, Ex), and assume éç to not depend on 0, y nor Therefore, con- ditionally on component = k, an approximate generative model (while the assumed generative model is denoted p())()) is given by
43 0, k) = NNy:; A,O + by, Ex).
(2)
To complete the hierarchical model, 0 is assumed to follow a mixture of Gaussian distributions specified by
( 0 = k) = MM ēt. fx),
(3a)
U k) Te-
(3b)
The GLLiM model hierarchical construction above (eq. to 3bb defines then joint Gaussian mixture model on (y, 0):
K
23 (3. 0)
9 0 2-)0 0 =k)g
While NN could potentially be incorporated into SeMPLE for the purpose of obtaining summary statistics of the data, as done eg ABC, this is separate problem that not considered the present work.
5

---

Published in Transactions on Machine Learning Research (06/2024)
where the fulll vector of mixture model parameters is
の TT CR2 Ts, Ar, bs, Ex)
(5)
From this joint model, we can deduce straightforwardly both its conditional distributions qs (5(v))) ()) and ( 10/ y) in closed-form. First, it comes that,
( 0)
ii (0) + bu. En), Σ k=
(6)
with
TNi(0;; Cz I'n) NgN((0; 2, r,
(0)
(7)
and similarly,
K
( = Σ
Th (B) Auy + bk. Ex).
(8)
with
TNNa Fs) Ty
(1)
(9)
where Φ (Tk: Ck. Ig A bg. Σ deduced from & (Deleforge et al.
is just a convenient equivalent reparameterization that can easily be
2014 as
Anc + bi
(10a)
T₂ Σn A A An
(10b)
Σ, (Tx) As DS Av)
(10c)
A, Σ
A Σx
(10d)
b. Σn (
Ag 5-
(10e)
The modelling power of GLLiM comes essentially from the fact that both conditional distributions above are mixtures of Gaussian affine experts (MoE), which are known for their modelling and approximation capabilities (Nguyen et al.; 2019 2021 2023). As such, they are good candidates to provide surrogate likelihoods and posteriors. All the more so as, given a training dataset made of many pairs of (0... ya h==1t surrogate models specified via GLLiM can be estimated via an EM algorithm, see Deleforge et al. 22014 4) and Forbes el al. (2022) for details. In fact surrogale likelihood and posterior can oblained with ( (O., yn. ) J prior-predictive realizations simulated from p(0)p(3) 0), but more generally, as exploited in SeMPLE, GLLiM be used in sequential way, by changing the training set to include posterior draws.
Once estimate for the parameter vector & (or equivalently Φ) is available, we have an amortized procedure providing a surrogate likelihood function 93 (y (), which can be evaluated for any data y. In particular, if
we specify y to be given observed dalasel yo, and plug yo into (6), we obtain an approximate (surrogate) likelihood function 4% (yo ()) for observation yo- Similarly, and using the same training data, an amortized surrogate posterior 26 0 y), for a generic y, can also be derived as mixture of Gaussian experts distribu- tions from If we plug yo into (8), then (0) y.) is an approximate posterior based on observed data. In summary, once ò is estimated via EM, we obtain surrogates of both the likelihood and the posterior in the same GLLiM run, simultaneously and using the same training data. The way this is performed and integrated with the rest of the ScMPLE method is explained in Section 4 Note that it is possible to assume a different distribution for éµ than the Gaussian one. For example [1] can be a mixture of generalized Student distributions, resulting the so-called SLLiM method (Perthame et al.) 2018 Both GLLiM and SLLiM are implemented in the R package xLLiM (Perthame et al. (022) which we employ our experiments, see Section 5 for details. Notice that GLLiM has already been used SBI, see sce Forbes al. (2022). namely they produced amortized sampler whose posterior draws then "refined" via single ABC step. However,

---

Published in Transactions on Machine Learning Research (06/2024)
in Forbes et al. (2022) only surrogate posterions were investigated and, in addition, they did not sequentially refine the surrogate posterior and surrogate likelihood. In contrast, our SeMPL.E approach makes use of both posterior and likelihood surrogates, to then provide an automatically "tuned" proposal function for MCMC (as described in Scction 4 but also to provide a surrogate likelihood, which finds excellent use in the MCMC we construct in Section Then, SeMPLE sequentially refines, multiple rounds, both the surrogate likelihood and the surrogate posterior (and hence the proposal function), without requiring an ABC step, and hence without the critical tuning/ /stuup decisions that ABC inference entails (eg the choice of distances, threshold values, the determination of suitable proposal functions), and without the large number of simulations that ABC requires.
Another approach that, similarly to GLLiM, could be considered for closed-form approximations of densities without involving neural networks, is the "elliptical process" Bånkestad et al. 20200. in particular when using their piecewise constant mixing distribution.
Algorithm SeMPLE
1: (0) t =-(0), set starting K and number of rounds R.
2: for r == 0: R do
3:
r = or r 1 then
sample [0,]mm) iid 0n ~ (0) (without MCMC), I,... ........ N.
else
sample [0 via 0, (0) using Algorithm 2
end if
8:
Conditionally on every On, obtain fyn J by sampling (implicitly) p(v 0,,) (n I...,N). N).
9:
if = then
10:
Collect Dp = (in, Yn I Train 9ı. (y ()) on Do. Update K. Obtain 9dh (0 y).
11:
12:
13:
Set 51 (0) 400 (0 y yo).
14:
else
15:
Collect D. = Dr-1 U (On, Yn J.M Doo
16:
Train 9z. (5.(y11 ()) on Dr. Update K.
17:
Obtain 9s. (0) y).
18: Set Pr+1(() - Lz. (yol ()p(O)
19:
((optional) pr+ (0)
p(0) P. to. 00 y.))
20: end if
21:
end for
Algorithm Independence MH with GLLiM surrogale posterior proposal distribution
For the initial value 01, pick the last accepted value from the currently available Markov chain.
2: for 2: N do
3: Proposc 0" ~ 9 (0 yoo,
min (1, (0 1] G0 0 ye)
(0")
90.-, , 1 11o)
Sample 00
6: if < a then
7:
0, 0-
8: elsc
0, 0j-1-
10: end if
11: end for

---

Published in Transactions on Machine Learning Research (06/2024)
4 SeMPLE: sequential learning of likelihoods and posteriors
We have shown how GLLiM can provide amortized estimations of both the likelihood and the postcrior distribution at additional cost, given training data. This GLLiM feature is exploited in our proposed Sequential Mixture Posterior and Likelihood Estimation (SeMPLE), to obtain a sequence of increasingly improved samplers, to more accurately and more efficiently sample from the targeted posterior of observed data Yo- Similarly to SNL (Paaaaakaarios et al. 2019), given a likelihood approximation (V. 0). SeMPLE targets an approximate posterior p(B) x V3 ()p(() with posterior sampling being carried out using a MCMC procedure. Typically, the MCMC scheme used in SNL is slice sampling (Neal, 2003), which has the advantage to be Luming-fres free but is likely to have difficulties in exploring high dimensional posteriors with complex multimodalities. As matter of fact, Papamakarios et al. (2019) recommend SNL "at most for tens of parameters that do not have pathologically strong correlations" " In ScMPLE, since GLLiM is able to learn a closed-form approximation of the posterior distribution (in addition to a likelihood function), such posterior is used as a proposal distribution within a Metropolis- Hastingg (MH) algorithm. SeMPLE, by using an efficient EM procedure to fit a flexible proposal distribution, is expected to benefit from a MH step that better guides simulations towards more likely paramcter values. To ease the comparison with SNL, the latter is recalled in Appendix D while SeMPLE is summarized in Algorithms 1
Algorithm shows that the targeted approximate posterior for 0 is a function denoted with p~ x 93. (yo 0 () which is soquentially obtained, where r = O, R denotes the ScMPLE rounds for the current p) @pproximation In the following we denote by ((r) and / a parameter and an observation simulated in round r. respectively. At first (r = 0) (o) is initialised at the prior p(0). An initial batch of N pairs yn IN produced from the prior-predictive distribution (step 4) constitutes the starting train- ing data for GLLiM, and the parameters de and do are obtained via EM within GLLiM to produce the initial surrogate likelihood and posterior Lz. (y/ 0) and qou (0 yy) respectively (steps 11-12). The initial surrogate posterior becomes the current target (step 13) p1 (0) = Ps0 (0 yo). In the next round (r = 1), the latter can be sampled from without MH (step 4) since (0) is defined as Gaussian mixture, from which it is trivial to sample N parameters (0(² j These N draws are then plugged into the model simulator to produce y p(y 2 and define Dj (0.)) yy J==I ) as the new training data (we discard Do for future rounds as it is decmed too uninformative for training). GLLiM is then fitted to Dj and the corresponding surrogate likelihood and surrogate posterior are obtained. From this step onwards, the surrogate posterior is employed as a proposal distribution within a tuning-free free MH algorithm (Algorithm 2 since, for an arbitrary prior p(0), when r 2 2 we do not explicitly know which distribution is proportional to the product p(0)ga, (y= ()) (for r = we know by construction that the targeted posterior is a Gaussian mixture), except for the special case of Gaussian and Uniform priors (see Papamakarios & Murray 2016 In Algorithm 2 the acceptance
3
Do = (0 ( (0)
probability of a proposed 0 is (for r > 1)
= mine 1, (0') - (0,-1 yo) (r 9 (f- yoo
Samples from the posterior p, are obtained by first proposing parameters 0 from Gaussian mixture model. From the current proposal function go. ,0 yo), which is the mixture
K
L (0) yo)
=
(r-1) (r-t) (v.)N(( m
- yo + b (r-1) 26
1=1
a component A* is randomly sampled from the set {1,..., K} with associated probabilities (r-1) (gw)). Then 0* is sampled from the k -th component Ni(0; A y + 6 (r-1) . - (r-1))
(m) (r-1) (vo).... Ux
by (optionally) slightly inflating its covariance matrix with a factor 2 1, which we set to = 1 (no in- flation) or 7 E [1.1, in our experiments. The inflation factor is to ensure that the proposal is sampled from a broud-enough distribution covering the support of the target (no inflation factor is used when r = 1).
We the term round instead of the most common ieeraton, as the latter could be confused with the MCMC iterations that are run within given SeMPLE round.
8

---

Published in Transactions on Machine Learning Research (06/2024)
Also, note that the N retained draws via MH are N post-burnin draws, where the burnin consists of 100 iterations, and the last accepted draw from round r of SeMPLE is used to initialize the chain at round r+ 1.
Importantly, our MCMC procedure is tuning-free. We use an independence sampler (Robert & Casella 2004], relying on GLLiM to provide an informative proposal function. This proposal function is trained on draws accepted in previous rounds which are, reasonably, more spread than the currently targeted posterior. In addition, being mixture model, it is suitable to deal with multimodalities in the target. This greatly simplifies the construction of ScMPLE no complex tuning is necessary for the MCMC sampler. While K could in principle be set arbitrarily large, in practice this may slow the EM algorithm convergence. In ScMPLE, the value of K is allowed to decrcase between rounds, as expressed with "update K " in lines 11 and 16 of Algorithm | In fact, mixture components having probabilities Гіk smaller than a threshold are discarded (ie these gets removed from the mixture model). In most cases we find that setting a starting value of K 20 or 30 produces good results, with a threshold for ii typically set to 0.005 (multiple hyperboloid, Lotka-Volterra Volterra, biological SDE modcl), 0.03 (Bernoulli GLM model) or simply 0 for no deletion (Two Moons model). This is to prevent having too many irrelevant components, that could be problematic when fitting the mixture to training data produced via MCMC, as the latter may not manage to visit tiny modes. It is also possible to choose the starting value of K in principled way. In [Deletorge et al.) (2014h, the Bayesian Information Criterion (BIC) is used to select K, as we do in some of our examples, see Appendix B An appropriate value for K can be learned by repeatedly computing the BIC on realizations of the prior-predictiv distribution, before SeMPLE is started. This use of the BIC is an amortized procedure, as the selected values can be reused regardless the specific observed data you The impact on the overall computational budget is negligible. We have also experimented with version of SeMPLE that learns corrected posterior (sce the optional part in step 19 of Algorithm D using the correction formula from Papanakarios & Murray (2016). However, found that learning surrogale likelihood was much more beneficial, providing better inference coupled with a higher acceptance rate. Notice that none of our reported results use the optional corrected--osterrio approach.
In summary, a major difference with the several variants of sequential neural posterior estimation (where the most used one is SNPE-C, GGeenberg et al.] 2019), and with SNL, is that SeMPLE uses Gaussian mixtures learned via EM to approximate both the likelihood function and the postcrior distribution of 0, as opposed to neural density estimation. Thus NN architecture has to be designed, nor special stochastic gradient descent. (SGD) setup has to be considered to fit. NNs. Only the number of components K and the structure of the covariance matrices Σn (as the Σ3 are obtained as an algebraic by-product), according to the options offered by the xLLiM package") have to be set. As we show in our examples, inference based on procedures using SGD can be much more computationally demanding and time consuming.
5 Numerical illustrations
We illustrate ScMPLE on several simulation studies, see also the Appendices for extra results and additional models. SeMPLE is written in R and uses the GLLiM implementation from the xLLiM package (Perthame
ct al. )022). Comparisons with neural-based likelihood and posterior estimation methods (SNL and SNPE- C) use the Python implementations in the SBIBM package for SBI benchmarking (Lueckmann et al. 2021 and architecture specifications for SBIBM are in Appendix A To compare approximate posteriors with a suitable reference posterior, we consider the classifier two-samples test (C2ST), as implemented in SBIBM, and Wasscrstein distances via the POT Python package (Flamary et al. 2021). C2ST varies in the interval [0.5,1), the lower the better. It equals 0.5 when samples from two distributions appear isatistiallly indistin- guishable and equals 1 when the classifier is able to perfectly separate samples belonging to two different distributions. We also report runtimes in figures and resource requirements in Table D For ScMPLE, the time required to select K via BIC (when the BIC procedure is performed, which is not the case for all studies) is reported separately and not taken into account in the figures below. In fact, for some case studies, such as Lotka-Vo Volteraa, we merely picked a starting value of K without running the BIC procedure. When we used the BIC, adding this overhead time does not change the overall conclusion that SeMPLE is much
4The Ey's 's can he isotropic, diagonal or full matrices, and set all equal or varying with
9

---

Published in Transactions on Machine Learning Research (06/2024)
more time-efficient than the other methods. Moreover, the selection of K is amortized procedure, hence it does not need to be repeated for different observations of the same model.
In the next two sections 5.1 and 5.2 we illustrate results from two models having multimodal posteriors, in section 5.3 we summarise the Lotka- Volter- experiment, and in section [.6] and Table 2 we summarise results from other experiments detailed the Appendices. For comparison, all methods are given the same budget of model simulations. For SNL and SNPE-C the SBIBM default of distributing these simulations uniformly across R 10 rounds is used, but to case comparison we also run SNL and SNPE-C with the same number of rounds used for SeMPLE, which is R = 4, except for the Bernoulli GLM model where R = 2 suffices. In fact, we observe that ScMPLE does not need as many sequential steps as the other methods. It may be argued that SNL could perform better if a larger number of MCMC iterations was used, at the cost of a smaller mumber of rounds R. Results in the Appendices show that both SNL and SNPE-C perform worse under this alternative configuration. Regarding SeMPLE, when the computational budget limits the number of model simulations, we found it beneficial to keep the number of rounds R between and 3, where R=2 2 implies that we have exactly one round where MCMC is used, and this is typically already enough for satisfying inference, though usually additional round with MCMC (hence R = 3) produces further refined inference. In our work we show results up to R = 4, but in general, four rounds are not recommended. The algorithm performance relies first on the GLLiM ability to learn reasonably accurate likelihood and posterior approximations and then on the MCMC sampler to explore the parameter space. Both steps require a large enough number of model simulations, which depends mainly on the model dimension.
5.1
Two Moons
The Two Moons model is a standard benchmark example in SBI (Greenberg et al., 2019 Lueckmann et al. 2021 Picchini & Tamborrino (024)) to illustrate how algorithms deal with multimodality and the local crescent shapes of the posterior distribution. The model itself is defined in Appendix D Our experiments are run on the same 10 datasets provided in SBIBM, each of them being vector yo of length two. We use K = 30, value we determine as described in Appendix B (this determination only required 90 seconds). We run ScMPLE, SNL and SNPE-C algorithms for total budget of 10a model simulations, uniformly distributed across inference rounds, and we used both R = 4 and R= 10 rounds. The obtained samples are then compared with the corresponding posterior reference samples (obtained via MCMC, since the likelihood is tractable) provided in SBIBM. The performance metric results are averaged over the 10 different data sets. Figure 2 shows C2ST and runtime values. ScMPLE provides the best approximation while SNL performs the worst. The comparison with SNL is particularly interestingg as while both SNL and SeMPLE incorporate tuning-fres free MCMC step, SeMPLE is able to visit both modes with acceptance rate of 60-70%, while SNL often struggles in efficiently exploring a multimodal posterior as demonstrated by the higher C2ST values. Moreover, in this simple example the runtime for SeMPLE is about half the SNPE-C runtime, and is about three times lower than the SNL runtime. In other examples, (e.g. Figure 3), the speedup provided by SeMPLE is even more dramatic. All these benefits from SeMPLE come with a very small impact on the memory usage, see Section 5.3 for details.
5.2 Multiple hyperboloid model
The multiple hyperboloid model (Forbes et al. 2022 is constructed from a sound source localization problem in audio processing. This is a challenging inference task with a complex symmetrical and multimodal posterior distribution. The model is provided in Appendix [F] The observation is vector of length 10 and 0 has dimension 2. SeMPLE is run with K 40, value that is selected via BIC in 4.5 minutes. The total number of model simulations is set to 4 x 10 C2ST and runtime results are reported in Figure 3 An exemplary posterior is shown in Figure 4 SeMPLE shows the best performance according to C2ST, and its runtime is about half that of SNL, and is about 10 times more efficient than SNPE-C. Adding the BIC overhead does not change the conclusion that ScMPLE is much more efficient. The runtime comparison with SNL not particularly relevant, since the latter does not generally manage to explore the multimodal posterior (see Figure 4 and the posteriors in Appendix
10

---

Published in Transactions on Machine Learning Research (06/2024)
sutie with Creents SVL rounds
Cumulative elapoed time [minates)
0.
120
an
0.7
60-
20
0.5
43030
simulations
Figure 3: Hyperboloid. Median C2ST (left) and median cumulative runtime in minutes (right) for 10 runs with the same data vs the number of model simulations. Shaded bands enclose the min and max values.
p(6, Oalyo)
p(l
0.3
0.2
-2
91
p(8/y9))
0.4
Tuu posterier SeMPLE with munds
SNPE-C with noanss
10
0.2
0.1
Figure 4: Hyperboloid. An example of marginal posteriors and samples from the last round of cach algorithm. The SNL marginal posteriors are not reported their inference is completely off.
5.3 Lotka-Volterra
The VolkarVolterra predator prey model is a classical example of population dynamics model. We consider the stochastic Markov jump process version of the model, which is defined in Appendix ! This model, although being relatively simple and characterized by a set of three interactions on the two species, each with a corresponding rate 02 (j = 3),, encompasses many of the difficulties associated with larger, more complex systems. The observations consider the addition non-negligible measurement error to the simulated dynamics, and the standard deviation of this noise term inferred as an additional parameter, resulting in total of four model parameters to infer (0,, 03, For this model, are not primarily focusing the performance comparison to other neural-based methods, this version of the model, as Markov jump process, is not included in SBIBM. Morcover, running times for a single model call are random even for repeated calls using the same parameter value. We use the Gillespie algorithm,, as explained in Appendix ! and the number of "reactions" at any given stochastic, Due to the reasons above, the primary focus is to show that high quality inference can be obtained with ScMPLE using a comparably small computational effort, and compare our results against two versions of sequential Monte Carlo ABC (SMC-ABC). "reference posterior" is obtained with an efficient version of SMC-ABC called blockedopt in Picchini & Tamborrino
11

---

Published in Transactions on Machine Learning Research (06/2024)
2024), which is run for very many model simulations, around 3.3 million, see Appendix for details. We also compare it with more standard off-the-shelf SMC-ABC named standard. SeMPLE inference is obtained with only 30,000 model simulations. In Figure 5 we also include the corresponding inference obtained with around 30,000 model simulations for both blockedopt and standard SMC-ABC. With only 30,000 model simulattins, SeMPLE manages lo produce inference for parameters that is similar to the inference returned after 3.3 million model simulations from blockedopt SMC-ABC, a 100-fold computational saving against the most efficient version of SMC-ABC. We observe that, unlike SMC-ABC, SeMPLE produces accurate posterior inference for all model parameters, including the noise parameter d, but with much smaller number of model simulations compared to SMC-ABC. The latter fact, namely the SMC-ABC inefficiency in terms of model simulations, was not unexpected and known from Greenberg et al. (2019); Lueckmann et al. (2021). However it is reassuring to notice the excellent inference quality brought by SeMPLE in this challenging case study.
SPCAAC Staaadic
SHC-MCC Joclecptt SAAPLESEE
4-
3.
1-
(0g)
log(0z)
CACC shstedor XIME UPC. MECiLa
2.5
OMGASC Exciedee ME
Lac. dasted IBC
2.0
1.5
0.5
0.0
-6
-4
-2
2
-2
0
log(0g)
log(a)
Figure 5: Lotka- Volterra. Marginal posteriors from SMC-ABC and SeMPLE. SeMPLE is run for a total of 30,000 model simulations. We report posterior samples from blockedopt SMC-ABC after a total of 32,865 and 3,337,648 model simulations. We also report the posterior samples from standard SMC-ABC after a total of 30,436 model simulations. True parameter values are indicated with black vertical lines.
5.4 Biological model of the translation kinetics after mRNA transfection
We consider two-dimensional stochastic differential equation (SDE) model of the translation kinetics after mRNA transfection, as studied in Pieschner et al. 2022). Only one of the two coordinates of the process is observed (with measurement noise). Details about the model are in Appendix K Pieschner et al. (2022) perform several experiments, where data are either exact observations the system or are perturbed with measurement error, and inferences for both an ODE and an SDE model are compared. However, here we only study the SDE model with data perturbed with measurement error. The parameter to infer is (we consider pa- rameters log-scale, same for their priors) 0 (log 8, log 7, log k log mo, log scale, log fo, log offset, log o). We ran inference for five datasets, each consisting of time-scries of length 60, cach generated with different parameter values. Details about datta-generating paramelers, and other set I lings, are in Appendix of the five datasets, in Figure we show the marginal posteriors from the third round of SeMPLE, in compar- ison to the parameter priors, while Appendix [ we also show comparisons with SNL and SNPE-C. For all datasets, SeMPLE accurately infers (6.7, 7, fo, offset, σ), and we confirm the findings in Pieschner et al. (2022) (which they obtain using exact. Bayesian inference via the Hamiltonian Monte Carlo method implemented
K For one
12

---

Published in Transactions on Machine Learning Research (06/2024)
in Stan), namely that k, my and scale suffer from identifiability issues (see the supplementary section A.4.1 in Pieschner et al. 2022 where their parameter O2 is our k). In Appendix K we show that only SeMPLE can identify to accurately, while SNL and SNPE-C provide biased estimation of to. Moreover, SNPE-C is often unable to infer the measurement error variability a, while this is not a problem for ScMPLE.
-5
0
5
logó
logy
logk
logma
-1
2
3
-2
0
1
2
3
-4
2
logscale
logto
logaffset et
loga
Figure 6: Translation kinetics model: in green are the marginal posteriors from SeMPLE (at round r = 3), priors (solid black lines) and ground truth parameters (vertical red lines).
5.5 Resource requirement metrics
The metrics reported in Table I correspond to four models and have been measured in kiloJoules (kJJ) using the PyJoules" module found in the PowerAPI package (Bourdon et al.) 2013) Python, which makes use of the Running Average Power Limit (RAPL) technology available on Intel CPUs. It is important note that these energy measurements are based on hardware readings, and not on estimations depending eg running time. The experiments 3.60GHz and 16GB of DRAM. Memory cost and peak memory usage results have been mcasured by the tracemalloc Python library. These results highlight the advantage brought by SeMPLE, which can be run a minimal configuration while maintaining good performance. For instance, in the Two Moons case, the memory requirement (Peak Memory Usage) for ScMPLE 13 times lower than SNPE-C and 27 times lower than SNL. The gain is even larger for the other models in Table I Similarly, for the Ornstein- Uhlenbeck model the energy consumption (CPU+DRAM) more than 10 times smaller with SeMPLE.
machine with CPU Intel Core 17-4790 CPU â
5.6 Results summary
In Table 2 give summary of the main results from all experiments where objective comparison with NN-based methods was possible, including those in the Appendices. Instead, for the Lotka-Voltere Volterra model we did not run NN-based methods and we reported comparisons with SMC-ABC both in section 5.3 and in Appendix ! Generally, ScMPLE returns accurate inference, either by producing the smallest metric value (the smaller the better), as for Two Moons and Hyperboloid, or being second best after SNL, as for Bernoulli GLM and Ornstein-Uhlenbeck, also considering the larger variation of the SNPF-C results. However, even when SNL produces smaller Wasserstein distance, has to be balanced with respect the running time. For Bernoulli GLM, SNL is between 18 and 75 times slower than SeMPLE, depending the number of rounds, SNL has to perform expensive neural network training at each round. For the same reason, when looking at Bernoulli GLM, SNPE-C is between 8 and 26 times slower than SeMPLE,
http:://ppjollee.rattthoooo.
13

---

Published in Transactions on Machine Learning Research (06/2024)
Table Median metric results, across 10 independent runs (brackets include minimum and maximum values), and median runtime for posterior inference. Best (lowest) values are in bold. Notice, "Metric (final)" reports the metric at the final inference round, instead "#Model sims" gives the total number of model simulations across all rounds for single run, and similarly for "Runtime". SNPE-C and SNL used the default SBIBM setup except when the number of algorithm rounds is reduced from 10. The C2ST metric is reported for multimodal posterior distributions and the Wasserstein distance for unimodal posterior distributions according to the discussion in Section 5.6
Method
Metric (final)
# #oodel sims
Runtime (minutes)
# Rounds
Two moons (multinodal)
SNL
C2ST=0.62 [0.56, 0.92]
10K
19
[17, 31
4
SNL
C2ST=0.60
[0.56, 0.68
10K
22
211 23
10
SNPE-C
C.59
[0.54, 0.66
10K
8.1
[7.1, 12
SNPE-C
0257=0.57 00.52, 0.60]
10K
15
[13, 17]
10
SeMPLE
-.57-0.54 (0.50, 0.58)
10K
8.3
[7.9, 8.8
4
Hyperboloid (multimodal)
SNL
C25T=0.90 [0.84, 1.0]
40K
35 [30, 39
4
SNL
C2ST=0.90 [0.86, 0.95
40K
50 45, 57
10
SNPE-C
C2ST=0.76 (0.74, 0.81
40K
100
[76, 150
4
SNPE-C
C.ST=0.74 (0.70, 0.75
40K
124 [100, 151]
10
SeMPLE
C25T=0.58 [0.56, 0.60]
40K
14
[13, 17]
4
Bernoulli GLM (unimodal)
SNL
Wasserstein=0.58 [0.44, 0.71
10K
22
22, 24
2
SNL
W.sserstein==069 [0.46, 0.89)
10K
91,
92
SNPE-C
Wasserstein=0.94 [0.65, 1.6]
10K
9.5
[4.9, 15
2
SNPE-C
Wasserstein=0.71 00.49, 1.78
10K
31
[25, 40]
10
SeMPLE
0.57 [0.43, 1.0]
10K
1.2
00.98, 1.2]
2
Orastein-Uhlenbeck (unimodal)
SNL
Wasserstein=00067 [0.057, 0.23]
40K
56
[54, 60]
4
SNL
Wasserstein=0.12 [0.079, 0.20]
40K
62 [58, 71]
10
SNPF-C
Wasserstein=0.76 [0.27, 1.0)
40K
100
[83, 150)
4
SNPE-C
Wasserstein==0344 [0.19, 0.48]
40K
130
[100, 150]
10
ScMPLE
Wasscrstein=0015 j0.11, 0.28]
40K
5.9 [3.7, 14]
again depending on how many times it needs be retrained. When considering the Ornstein-Uhlentee model, the Wasserstein distances have good degree of overlap between SNL and SeMPLE, but SNL is between 9.5 and 10.5 times slower than SeMPLE, while SNPF-C is between 17 and 22 times slower than SeMPLE. In Table 2 we highlight whether experiment has a unimodal or multimodal target, as this determines the type of metric we report. Namely we use the C2ST metric for multimodal targets and the Wasserstein distance for unimodal targets. The reasons for this are now discussed, however notice that in the Appendices we report results for both metrics in each experiment. In our experiments, we found the C2ST metric (advocated in Lueckmann et al., (021)) useful evaluate the inference quality with multimodal targets, when draws from reference posterior produced contradicting results with unimodal targets. An example is the Bernoulli GLM model considered in Appendix G There, it is obvious that the marginal posteriors returned by ScMPLE are accurate and that inference improved between the first two rounds, however C2ST showed that inference was apparently deteriorating for SeMPLE from round r = to r = 2, which was not the case in reality. We also computed Wasserstein distances, which could confirm that the inference quality via ScMPLE did improve from round to round 2. However, for multimodal targets, it
available. However, we also found that times C2ST
that Wasserstein distances may not be a suitable tool, e.g. Balaji et. al. (2019). We therefore the reader from blindly relying single performance
14

---

Published in Transactions on Machine Learning Research (06/2024)
metric. It is also important to examine marginal and joint posterior densities. For example, in Bischoff et al. (2024) it is mentioned that interpreting results for C2ST requires knowledge of the classifier used and its appropriateness for the data at hand, and that its effectiveness is critically dependent on the selection and training of suitable classifier.
6 Conclusion and perspectives
We have addressed the challenging question of Bayesian inference in the absence of explicit likelihoods. We have investigated mixture models through the GLLiM method, and shown that their expressivity could account for multivariate distributions, while their much simpler structure made them more amenable to efficient learning than NN solutions. The resulting new inference strategy, named SeMPLE, yielded inference that was on par with inference from the state-of-art neural-network. (NN))based methods that we examined (SNL and SNPE-C), while being more energy and memory efficient. For the considered examples, and the considered sinulation selup, SeMPLE produced inference of superior similar quality than NN-based methods, and when it did occasionally under-perform, it was by negligible amount that can largely be compensated by the much smaller runtime. In fact, in all tested cases, SeMPLE displayed much lower running times and memory requirements. The energy gains appear to increase with the complexity of the model simulator, and we expect. to obtain even larger gains with real-life simulators. For example, for the Chlecbin model, ScMPLE achieved a 17-fold runtime spcodup compared to SNPE-C. Although very preliminary, we hope with this attempt meet the constraints of "Green Al" (Schwartz al., 2020), and open the way to the development of more methods that can achieve accurate posterior inference, balancing the environmental impact of growing energy costs with the obtained results quality.
Of course, SNL and SNPE-C could certainly be improved for each specific case study, but we did not examine or fine-tune SNL and SNPE-C on the many other possible configurations, as this can be influenced by many factors. While examining several case studies implemented in SBIBM, we restricted to simulation settings provided in the SBIBM package, including reusing the same NN architectures and training setups, see Appendix A for details.
Some open questions, that can be addressed in future research, follow. SeMPLE can benefit from the summarization of data when the dimension of y is large, just like other SBI methods. In addition to reducing the dimension, using summary statistics (as we did with the Lotke-Volterra model) can also improve the GLLiM performance. Indeed the performance of SeMPLE depends on the availability of good first approximations of the posterior and likelihood, which themselves depend on the ability of GLLiM to find probabilistic mappings between parameters and observations, in regression-like manner using a Gaussian mixture structure. The construction of informative summary statistics is important, but independent, research question that affects not only standard ABC methods, but also SBI methods based of neural networks, e.g. Chen et al. [2021); Radev et al. (2023b). Finally, following the discussion in Section 5.6 we wish to emphasize that measuring the quality of posterior inference in SBI is certainly an important question deserving further rescarch, and that multiple metrics should be considered when evaluating experiments. For a recent study on metrics to evaluate the quality of SBI, see Bischolf et al. (2024)
Acknowledgements
We thank Emcline Perthame (Institut Pasteur) for assistance with the xLLiM package. HH and UP acknowl- edge funding from the Swedish National Research Council funding from the Chalmers AI Research Centre.
Caancillleeenneeeae (Vetenskapsrádet 2019-03924). UP acknowledges
References
Christophe Andricu and Gareth O Roberts. The pscudo-marginal approach for efficient Monte Carlo com- putations. The Annals of Statistics, 2009.
Christophe Andrieu, Arnaud Doucet, and Roman Holenstein. Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society Series B: Statistical Methodology, 72(3):269 342, 2010.
15

---

Published in Transactions on Machine Learning Research (06/2024)
Yogesh Balaji, Rama Chellappa, and Soheil Feizi. Normalized Wasserstein distance for mixture distributions with applications in adversarial learning and domain adaptation. arXiv preprint arXiv:1902.0014 2019.
Maria Bänkestad, Jens Sjöhnd, Jalil Taghia, and Thomas Schön. The elliptical processes: a family of fat-tailed stochastic processes. arXiv preprint arXiv:2003.02001 2020.
Mark A. Beaumont, Jean-Marie Cornuct, Jean-Michel Marin, and Christian P. Robert. Adaptive approxi- mate Bayesian computation. Biometrika, 96(4):983 990, 10 2009.
Scbastian Bischoff, Alana Darcher, Michael Deistler, Richard Gao, Franziska Gerken, Manuel Gloeckler, Lisa Haxel, Jaivardhan Kapoor, Janne K Lappalainen, and Jakob H Macke. practical guide to statistical distances for evaluating generative models in science. arXiv preprint arXiv:2403.1220 12636, 2024.
Christopher M. Bishop. Mixture density networks. Technical report, Aston University, 1994.
Aurélien Bourdon, Adel Noureddine, Romain Rouvoy, and Lionel Seinturier. PowerAPI: A software library to monitor the energy consumed at the process-lovel. ERCIM News, 2013.
Yanzhi Chen and Michael U Gutmann. Adaptive gaussian copula ABC. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 1584-1592. PMLR, 2019.
Yanzhi Chen, Dinghuai Zhang, Michael U. Gutmann, Aaron Courville, and Zhanxing Zhu. Neural approx- imate sufficient statistics for implicit models. In International Conference on Learning Representations, 2021.
Kylc Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation- inference. Proceedings of the National Academy of Sciences, 117:201912789, 05 2020.
Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. An adaptive sequential Monte Carlo method for approx- imate Bayesian computation. Statistics and Computing, 22:1009- 2012.
Antoine Deleforge, Florence Forbes, and Radu Horaud. High-dimensional regression with Gaussian mixtures and partially-latent response variables. Statistics and Computing, 25(5):893-99) 911, mar 2014.
Conor Durkan, lain Murray, and George Papamakarios. On contrastive learning for likelihood inference. In International conference on machine learning, pp. 2771-2781. PMLR, 2020.
Ritabrata Dutta, Marcel Schoengens, Lorenzo Pacchiardi, Avinash Ummadisingu, Nicole Widmer, Jukka- Pekka Onnela, and Antonietta Mira. ABCpy: A high-performance computing perspective to approximate Bayesian computation. Journal of Statistical Software, 100:1-34 2017.
Joel Dyer, Patrick W Cannon, and Sebastian M Schmon. Amortised likelihood-ffre inference for expen- sive time-scries simulators with signatured ratio estimation. In International Conference on Artificial Intelligence and Statistios, pp. 11131 11144. PMLR, 2022.
Dirk Eddelbuettel, Romain Francois, JJ Allaire, Kevin Ushey, Qiang Kou, Nathan Russell, Inaki Ucar, Douglas Bates, and John Chambers. Repp: Seamless R and C++ Integration, 2024. URL https:
//CRAN.R-par R-projectt.....mm www.ooogggg/aaee.oo Rackage=Rcpp R package version 1.0.12.
Mikhail Evchenko, Joaquin Vanschoren, Holger H. Hoos, Marc Schoenauer, and Michèle Sebag. Frugal machine learning. arXiv preprint axXiv:211..03711. 2021.
Sarah Filippi, Chris P Barnes, Julien Cornebise, and Michael PH Stumpf. On optimality of kernels for approximate Bayesian computation using sequential Monte Carlo. Statistical applications in genetics and molecular biology, 12(1))87 107, 2013.
Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Alaya, Aurélic Boisbunon, Stanislas Cham- bon, Lactitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. POT: Python optimal transport. The Journal of Machine Learning Research, 22(1):3571-333 -3778, 2021.
16

---

Published in Transactions on Machine Learning Research (06/2024)
Florence Forbes, Hien Duy Nguyen, TrungTin Nguyen, and Julyan Arbel. Summary Statistics and Dis- crepancy Measures for Approxinate Bayesian Computation via Surrogate Posteriors. Statistics and Computing, 32(5), oct 2022.
Richard Gao, Michael Deistler, and Jakob H Macke. Generalized Bayesian inference for scientific simulators via amortized cost estimation. arXiv preprint arXiv:2305.1220 2023.
Gillespie. Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry 81 (25), 2340-2361, 1977.
Andrew Golightly and Darren J Wilkinson. Bayesian paramoter inference for stochastic biochemical network models using particle markov chain monte carlo. Interface Focus, 1(6):807 820, 2011.
David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood-free inference. In International Conference on Machine Learning, pp. 2404 2414. PMLR, 2019.
Heikki Haario, Eero Saksman, and Johanna Tamminen. Adaptive proposal distribution for random walk Metropolis algorithm. Computational Statistics, 14(3):375-38) 395, 1999.
Jocri Hermans, Volodimir Begy, and Gilles Louppe. Likelihood-free MCMC with amortized approximate ralio estimators. In International conference on machine learning, pp. 4239 4248. PMLR, 2020.
Jingjing Li, David J Nott, Yanan Fan, and Scott A Sisson. Extending approximate Bayesian computation methods to high dimensions via Gaussian copula model. Computational Statistics & Data Analysis, 106: 77-89, 2017.
lan-Matthis Lueckmann, Pedro J Goncalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher, and Jakob H Macke. Flexible statistical inference for mechanistic models of neural dynamics. Advances in neural information processing systems, 30, 2017.
Jan-Matthis Lueckmann, Jan Boclts, David Greenberg, Pedro Goncalves, and Jakob Macke. Benchmarking simulation-based inference. In International conference artificial intelligence and statistics, pp. 343 351. PMLR, 2021.
Radford M Neal. Slice sampling. The Annals of Statistics, 31(3):705-777, 2003.
Hien D Nguyen, Faicel Chamroukhi, and Florence Forbes. Approximation results regarding the multiple- output Gaussian gated mixture of linear experts model. Neurocomputing. 366:208-21 -244 2019.
Hien Duy Nguyen, TrungTin Nguyen, Faicel Chamroukhi, and Geoffrey John McLachlan. Approximations of conditional probability density functions in Lebesgue spaces via mixture of experts models. Journal of Statistical Distributions and Applications, 8:1-15, 2021.
Huy Nguyen, TrungTin Nguyen, and Nhat Ho. Demystifying Softmax Gating Function in Gaussian Mixture Experts. Thirty-sevnnth Conference on Neural Information Processing Systems, 2023.
J. J. Nott, V. M.-H. Ong, Y. Fan, and S.A. Sisson. Handbook of approximate Bayesian computation, chapter High-dimensional ABC. Chapman and Hall/CRC, 2018
Jamie Owen, Darren J. Wilkinson, and Colin S. Gillespie. Likelihood free inference for markov processes:
comparison. Statistical Applications in Genetics and Molecular Biology, 14(2):180-209, 2015a.
Jamie Owen, Darren J Wilkinson, and Colin S Gillespie. Scalable inference for Markov processes with
intractable likelihoods. Statistics and Computing, 25:145 156, 2015b.
George Papanakarios and lain Murray. Fast e-free inference of simulation models with Bayesian conditional density estimation. Advances in Neural Information Processing Systems, 29, 2016.
George Papamakarios, David Sterratt, and Iain Murray. Sequential Neural Likelihood: Fast Likelihood-free Inference with Autoregressive Flows. In Proceedings of the Twenty- Second International Conference on Artificial Intelligence and Statistics, volume 89, pp. 837 848. PMLR, 16-16 18 Apr 2019.
17

---

Published in Transactions on Machine Learning Research (06/2024)
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshmi- narayanan. Normalizing flows for probabilistic modeling and inference. The Journal of Machine Learning Research, 22(1):2617 2680, 2021.
E. Perthame, F. Forbes, and A. Deleforge. Inverse regression approach to robust nonlinear high-to-low dimensional mapping. Journal of Multivariate Analysis, 163(C): 2018.
Emeline Perthame, Florence Forbes, Antoine Deleforge, Emilie Devijver, and Melina Gallopin. xLLiM: High Dimensional Locally-Lincar Mapping, 2022. R package version 2.2.1.
Henri Pesonen, Umberto Simola, Alvaro Köhn-Lauqe, Henri Vuollckoski, Xiaoran Lai, Arnoldo Frigessi, Samuel Kaski, David T Frazier, Worapree Maneesoonthorn, Gael M Martin, and Jukka Corander. ABC of the future. International Statistical Review, 91(2):243 268, 2023.
Umberto Picchini and Massimiliano Tamborrino. Guided sequential ABC schemes for intractable Bayesian models. 2024. Forthcoming in Bayesian Analysis 161:10..214/24-BA. BA1451. Also available as arXiv preprint, arXiv:2206.1223:
Susanne Pieschner, Jan Hasenauer, and Christiane Fuchs. Identifiability analysis for models of the translation kinetics after mRNA transfection. Journal of Mathematical Biology, 84(7):56, 2022.
Leah F Price, Christopher Drovandi, Anthony Lee, and David J Nott. Bayesian synthetic likelihood. Journal of Computational and Graphical Statistics, 27(1):1 11, 2018.
Stefan T. Radev, Marvin Schmitt, Valentin Pratz, Umberto Picchini, Ullrich Koethe, and Paul Buerkner. JANA: Jointly amortized neural approximation of complex Bayesian models. In Proccodings of the Thirty-Ninth Conference on Uncertainly Artificial Intelligence, pp. 1695 1706. PMLR, 2023a.
Stefan T Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ull- rich Köthe, and Paul-Christian Bürkner. BayesFlow: Amortized Bayesian workflows with neural networks. arXiv preprint 16111 2023b.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International conference on machine learning, pp. 1530-1538. PMLR, 2015.
Christian P Robert and George Casella. Monte Carlo Statistical Methods. Springer, 2004.
Schwartz, J. Dodge, N.A. Smith, and O. Etzioni. Green AI. Commun. ACM, 63(12):54466) 63, 2020.
Yannik Schälte, Emmanuel Klinger, Emad Alamoudi, and Jan Hasenauer. AyABC: Efficient and robust easy-Lo-use approximate bayesian computation. Journal of Open Source Sofwware, 7:4304, 2022.
Scott Sisson, Yanan Fan, and Mark Beaumont (eds.). Handbook of approximate Bayesian computation. CRC Press, 2018.
Sam Stites, Heiko Zimmermann, Hao Wu, Eli Sennesh, and Jan-Willem van de Meent. Learning proposals for probabilistic programs with inference combinators. In Uncertainty in Artificial Intelligence, pp. 1056-1066. PMLR, 2021.
F. Strubell, A. Ganesh, and A. McCallum. Energy and Policy Considerations for Deep Learning in NLP. Tn 57th Meeting of Assoc. Computational Linguistics, 2019.
Neil C. Thompson, Kristjan Greenewald, Keeheon Lee, and Gabriel F. Manso. The computational limits of deep learning. ArXiv, abs/2007.00533 005558..555 2022.
Tina Toni, David Welch, Natalja Strelkowa, Andreas Ipsen, and Michael PH Stumpf. Approxinate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface, 6(31):187-20) 2009.
G. E. Uhlenbeck and L. S. Ornstein. On the Theory of the Brownian Motion. Phys. Rev., 36:823 841, Sep 1930.
18

---

Published in Transactions on Machine Learning Research (06/2024)
Tongzhou Wang, Yi Wu, Dave Moore, and Stuart J Russell. Meta-learning MCMC proposals. Advances in neural information processing systems, 31, 2018.
Darren Wilkinson, smfss: Stochastic Modelling for Systems Biology, 2024. URL org/package=snnfsb package version 1.5.
Darren J. Wilkinson. Summary stats for ABC, 2013. https://darrenjii wordpress. com/2013/09//11 /
summary- fraats-for-bbe
Samuel Wiqvist, Jes Frellsen, and Umberto Picchini. Sequential neural posterior and likelihood approxima- tion. arXiv preprint WrXiv2202..06222 2021.
Simon N Wood. Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466(7310): 1102-1104, 2010.
Appendix
SBIBM
BIC computation for selecting K
C SNPE-A
SNL
E Two Moons model
E.1 Two Moons model definition
E.2 Two moons: targeting the posterior directly vs targeting the likelihood.
E.3 Two moons: selection of K via BIC
E.4 Two moons:
Results with R
F
Multiple hyperboloid model
40
F.1
Multiple hyperboloid model definition
F.2 Multiple hyperboloid model: selection of via BIC
F.3 Multiple hyperboloid model: results with R = 4.
40
Bernoulli GLM model
G.1 Bernoulli GLM model definition.
G.2 Inference setup
G.3 Results when using
for all algorithms
G.4 Results when using R 10 for SNL/SNPE-C and R 2 for SeMPLE
Orastein-Uhlenbeck model...
H.1 Inference setup
H.2 Results......
50
H.2.1 Results when using R for algorithms...
H.2.2 Results when using R = 10 for SNT/SNPP-C and R = 4 for SeMPLE.
I Lotka-Volteras model (Markov jump process with measurement error)
56
I.1 Model definition with three reactions
56
T.2 Inference setup
56
T.3 Results
19

---

Published in Transactions on Machine Learning Research (06/2024)
J. Twisted prior model with 20 parameters
J.1 Model
J.2 Results.....
Biological model of the translation kinctics after mRNA transfection
K.1 Model definition..
K.2 Results
64
A SBIBM
Our runs of SNL (Papamakarios et al. 2019 and SNPE-C given in the Python package SBIBM (Lueckmann et al. 2021 which we ran with default parameters. However, we produced a fork of the SBIBM GitHub repository which can be found at htttt:://g// /ttp:///////..//g/////// http://iituuu....oeeaaaaaaeaeeee com/henhagg/ sbibm, and conveniently modified this output posterior samples at each round of the inference, instead of only from the last round as given in the default implementation. This made it easier display the performance of SNL and SNPE-C at intermodiate rounds. Additionally, the SNL and SNPE-C algorithms in SBIBM for comparison of runtimes. Note that the multiple hyperboloid model and the Ornstein- Uhlenbeck process had to be implemented as new tasks in SBIEM to run these models with SNL and SNPE-C. The implementation of these tasks can also be found in the SBIBM fork mentioned above.
Greenberg et al., 2019 used the implementation
added time mcasurements to
The default setting for SNPE-C in SBIBM to utilize Neural Spline Flows (NSF) as density estimator with 50 hidden features and 10 atoms. Correspondingly, SNL by default uses Masked Autoregressive Flow (MAF) for density estimation with 50 hidden features. The MCMC step of SNL utilizes vectorized slice sampling with 100 Markov chains and thinning interval of 10 by default.
B BIC computation for selecting K
In GLLiM, the number K of components in the Gaussian mixture model has to be specified prior to per- forming the EM procedure. A too large value of K may result in overfitting and unnccessary computational effort, while a too small value of K may limit the ability to represent the relationship between 0 and The Bayesian Information Criterion (BIC), used in Deleforge et al. (2014) to guide the selection of K, is given in equation [1
BIC = -2C( + D log
(11)
where ò the GLLiM model parameter, L(Φ) is the maximised value the GLLiM log-likelihood function at the MLE d, D(Ø) is the total number of parameters in the model and N is the number of observations in the training dataset. We wish to select a returning a small BIC when GLLiM is fitted with K components to a training dataset (0,, yn m==1 obtained by sampling parameters 0 ~ p(0) from the prior and simulating the corresponding Ya ~ pyy from the generative model. The GLLiM logeliieldoood is specified as (see
eq. (4) in the main body),
((Φ)
log 93 (ym, On):
(12)
with
94 (nn Q.)) = Σ
Ny,: A.o + b.. Σx) N(O.; Cx2 F&)) m
(13)
k=
The number of parameters that GLLiM needs to estimate
D = (K- - K(DL + D nbpar,; + nbparr),
(14)
where ubpary and nbparr are the number of parameters in the covariance matrices E, and F&, respectively. The covariance structure of the matrices Egg and F can be constrained to reduce the number of parameters
20

---

Published in Transactions on Machine Learning Research (06/2024)
that GLLiM needs to estimate. One possible constraint is to set the covariance matrices to be isotropic, which means that they are set to be proportional to the identity matrix, which is the default in the xLLiM package (Perthame et al. 2022).
The GLLiM implementation in the xLLiM package removes a mixture component & if the mixture probability becomes zero. Empirically it was found that when mixture probabilities becomes very small, GLLiM could occasionally crash with errors about failed internal matrix inversions. To prevent this, threshold was set within ScMPLE to remove mixture components with probabilities îi (0) below this threshold, which we set to be 0.005 for the multiple hyperboloid, Ornstein- Uhlenbeck, Lotka- Volterra and the biological SDE model, zero for Two Moons (hence no component was deleted) or 0.03 for Bernoulli GLM and SLCP models. This check was repeated in every round of SeMPLE. This shows the need for setting reasonable initial number of mixture components K, as setting a too large value inflates the problem of small mixture component probabilities. However, several cases our inference tasks worked out-of-the-bon by setting the initial value to K 20 or 30.
C SNPE-A
In Papamakarios & Murray (0166) a parametric approxination to the exact posterior is learned via sequential updates, and the method is denoted SNPF-A (that is, Sequential Neural Posterior Approximation A This approximation is acquired by training mixture density network (MDN, Bishop 1994), which is a Gaussian mixture where the means and covariances are parameterised via neural networks. However, no method of estimating the number of mixture components is provided in SNPF-A. SNPF-A is given Algorithm 3 and The MDN with parameters & is denoted by qs y) and the proposal distribution for parameters O is denoted ((0) and called "proposal prior". Therefore, parameters On 5(0) are drawn and corresponding data Yn are simulated from the generative model p(y 0n), and the MDN (0 | is trained on the data set (0n, Vn J. to obtain posterior approximation.
Papamakarios & Murray (2016) show that using preliminary posterior fits to guide future simulations dras- tically reduces the number of simulations required to learn an accurate posterior approximation. In terms of notation, this means that parameters are proposed from p(0), which can be different from the prior p(0). However, not sampling from the prior distribution requires a correction factor to target the true posterior distribution. The motivation behind this is given in Proposition 1 in Papamakarios & Murray (0166), and ultimately they suggest to estimate the posterior by
(()) ( I ( yo).
p yo) x
As mentioned in Papamakarios & Murray 2016), in A Algorithm g as long as qa (a(y)) has only one Gaussian component (K = 1) then p(0) remains a single Gaussian throughout, This Gaussian approximation can be the true posterior, or it can serve as good proposal prior in
used as a rough but cheap approximation Algorithm g for fine-tuning non-Gaussian multi-component posterior. When the second strategy is used, i.e. Algorithm 4 is employed, they reuse the single-componee neural density estimator learnt in Algorithm B to initialize ds in Algorithm 4 and in this case the weights in the final layer of the MDN are replicated K times, with small random perturbations to break symmetry. Notice, in SNPE-A the proposal prior p(0) trained in Algorithm is restricted to be a Gaussian distribution, and the prior p(0) is assumed to be a uniform or Gaussian distribution to allow analytical calculation of ((0 y = yo).
In SeMPLE we do not have such restrictions, meaning that we can consider any prior p(0), however this comes at the price of introducing MCMC sampling (as in SNI.). However, unlike for SNPE-A, in SeMPLE are not constrained to reuse the estimated one-component Gaussian and replicate this K times to obtain K-components mixture, instead the GLLiM algorithm automatically determines the means and covariance matrices for each of the K components, which as such can be different between components, thus adding flcxibility to the shapes of the postcrior distributions that can be targeted.
6Thii is the name given in Papamakarios Murray 2016 but follows the taxonomy established in Durkan et al. 2020
21

---

Published in Transactions on Machine Learning Research (06/2024)
Algorithm SNPE-A: Training of proposal prior
Algorithm 4 SNNE-A: A: Training of posterior
Initialize (0< y) with one component (K 1). p(0) < p(0).
Initialize qs ( x) with K components.
(if 46 available by Algorithm initialize by replicating its one component K times}
repeat
for N do
sample 0n (0)
for n narn=1- 1 : N do
sample yn p 0n)
sample 0n ~ p(0)
end for
sample y/ ~ ply (n.) end for
retrain (0 y) on (0.: Un ) n=11 ((0) p(0)
train (ø y) (fn, Yn |
p()] 4s (0 yo)
until p(0) has converged
p(0 y = yo) + P(O) 1sb (0 yo)
p(0)
D SNL
Sequential Neural Likelihood (SNL) Papamakarios et al. 2019) sequentially trains a conditional neural density estimator for the likelihood function, instead of training an estimator for the posterior as in SNPE- This eliminates the need for corrections stemming from the proposal distribution not being the prior, but introduces an additional MCMC step (in the form of slice sampling) to sample from the posterior distribution, by employing the estimation (3o of the likelihood. SNL is detailed in Algorithm 5 which is the same Algorithm 1 in Papamakarios ()))0)) and then plugs the latter into MCMC for posterior sampling, SeMPLE trains q3 (z(1)0) and also obtain (as an immediate algebraic by-product of the qs (y)0) training) proposal distribution qa( (0ly)) that is used in the MCMC step.
al. (2019) The main difference with SeMPLE is that, while SNL trains
Algorithm 5 Sequential Neural Likelihood (SNL)
Input: observed data yo: estimator 4 ()) number of rounds R, simulations per round N. Output: approximate posterior p(0) yo).
Set )
y.) =0(0) and D = 0
for r= = 1: R do
for N do
sample 0, ~ 1 (0 yo) with MCMC
simulate yn p(y O..)
add (On, 3/h) into D end for (re-)train ( ( 0) on D and set (0) y0) x (0 v())
end for return (0 yo)
22

---

Published in Transactions on Machine Learning Research (06/2024)
E
Two moons model
E.1 Two moons model definition
The Two Moons model, described in Greenberg et al. (009 generates y E TR22 for a given parameter E TR²2 according to
a ~ U( π π
(15)
2' 2
N....,,00000 0.013)
(16)
p (r (rcos(a) + 0.25, r rsin(a))
(17)
181 181 02
- 01 + 02
y' p
(18)
v2
v2
Same as in Greenberg et al. (2019), set uniform priors 01 ~ u(-1,,,111 02 ~ ~((11,)) By construction,
this toy example has a posterior distribution with two crescent shapes for given observation.
E.2 Two moons: targeting the posterior directly targeting the likelihood
In the main paper we discussed the two possibilities of targeting the posterior (+) x 9s. (y. 10)p(0) by first obtaining an estimate 9g. (yo O)) of the likelihood and then sampling via MCMC from (0), as opposed to the option of using the GLLiM surrogate posterior qs. (0) yo) to update the posterior 1
D. q0. (0) yo), and then sampling from the latter via MCMC. We hinted at the fact that it is preferable to follow the first route, which we used obtain all our results, and here we provide results in support of our recommendation. Figure shows the C2ST metric and the Wasserstein distance as function of the number of model simulations for different SeMPLE settings, including different choices for the covariance matrices Σ₁ ("isotropic" or "full unconstrained") of the mixture models. We first considered the "isotropic" covariances case (specified in xLLiM via the option "ntr&Sigma="""" which means that all the K covariance matrices are proportional to the diagonal matrix, however these matrices are all different between the K components. With this option, we ran SeMPLE by targeting (0) x Qz. (3o ()p(O) ("likelihood; isotropic" in the legend of Figure and also ran ScMPLE by targeting ( p(0)
p(0) 9..
(0) lo. (0 ()) ("posterior; isotropic" the legend of Fi Figure 8 Targeting (0) x Iz. (Bo ()) is clearly the preferred route according to C2ST. This is less clear by looking at the Wasserstein distance, however the computation of this particular metric may not be accurate with multimodal targets (while we recommend to use Wasserstein distances for unimodal targets), and we refer the reader to sect 5.6 the main paper for a discussion about metrics to evaluate posterior inference. Additionally, we explored the results of assuming full matrices Σn without constraints (specified in xLLiM via the option cstr$Sigma= method of sampling from (0) 9z. (yo ())( ("likelihood; unconstrained" in the legend of Figure g This ScMPL.E configuration produces the best results. This may not be surprising, as for two-moons Yo has dimension 2, therefore the covariances Σn have dimensions 2 x 2, thus not inducing a large number of parameters & to estimate via EM. Notice that the C2ST value, when using our recommended method but with isotropic covariance matrices, still lower than C2ST as obtained with SNPE-C (see the main paper). Regarding the performance of the MCMC algorithm, Figure 7 shows that the acceptance rate noticcably higher with the "likelihood" option rather than the "posterior" one, which strengthens our choice. While higher acceptance rate does not denote, per-se, any special ability for algorithm to efficiently explore the posterior, however this fact coupled with the already displayed high quality inference provided by ScMPLE for this example, means that Figure 7 further shows that we can obtain a larger number of representative posterior samples when targeting ( x 9z. (y. (0).
This setting is run together with (recommended)
Based on these results, was concluded that using the surrogate likelihood to target q,, | ()) in the MCMC step, with the surrogate posterior qo (0 yo) as proposal distribution in the Metstpolis>Hastinns Hastings algorithm, superior method overall, which we used in all the other simulation studies.
23

---

Published in Transactions on Machine Learning Research (06/2024)
target
Surrogate likelhood
Surrogate posterior
Acceplance rate
0.2
0.0
Iteration
Figure Two Moons: Median acceptance rate of the Hasstnglis-aastings different data sets. Error bars show min/man values.
Hastings step of 10 ScMPLE runs with
0.80
Settings
Settings
Posterior; Isotropic
1.0
Posterior, Isotropic
Likellhood; Isotropic
Wasserstein distance
Likelihoad; lsotropic
Likellhood: Unconstrained
Likelihood; Unconstrained
0.8
0.70
0.6
0.65
0.55
0.50
10000
2500
5000
10000
Number of simulations
Number of simulations
(a)
(b)
Figure 8: Two Moons: Median C2ST and Wasserstein distance on 10 SeMPLE runs with different data sets vs number model simulations. The three different SeMPLE settings refers to whether the surrogate likelihood (reported as "likelihood; isotropic" or "likelihood; unconstrained") or the "corrected" surrogate posterior (reported as "Posterior; isotropic") is usod in the MCMC target distribution, and constraints used for the GLLiM covariance matrices ("isotropic" or "unconstrained"), Error bars show min/mas values.
E.3 Two moons: selection of K via BIC
We illustrate how to select starting value for K using BIC. However, we recall that it is not strietly necessary to choose K in very precise way, as its value is let decrease during the SeMPLE run. All BIC values are computed on the same prior-predictive data set (6n, ) with N = 2 500, obtained prior to running SeMPLEF. Figure g shows BIC as function of K. The minimum BIC of the evaluated values is obtained at K = 50 but decreases sharply from K = 10 to K = 30. By the Occam's razor principle, 30 was used as input to SeMPLE. Empirical experiments showed that a larger value of did not improve results. The runtime compute all BIC values in Figure g was seconds (on a desktop computer "amortised", that
with 6-core (12 threads) AMD Ryzen 5 2600 CPU). Note that the BIC computations is are independent of the observed data set you and only have to be performed once before running SeMPLE, and can thus be recycled whenever a different observed data set is considered.
E.4 Two moons. Results with R =
In the main paper, we have considered the inference results based on total 10* model simulations, uniformly distributed across R 10 rounds for SNPE-C and SNT. (as per SBIBM defaults), and R 4 rounds
24

---

Published in Transactions on Machine Learning Research (06/2024)
-6000
-0000
BIC
-8000
-0000
10
20
30
40
50
60
70
K
Figure 9: Two Moons. BIC with respect to the number of mixture components K.
SeMPLE
SNPE-C
SNL
Reference
1.0
1.0
Number of points per pixel
9
7
155
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
100
-0.5
-0.5
-0.5
-0.5-
-10
-1.0
-1.0
-1.0
0
Figure 10: Two Moons using R rounds for every method. Density scatter plot of posterior samples from the last round (r=4) of cach algorithm. Observed data set number in SBIBM.
for SeMPLE. Results showed that SeMPLE performed better and in partieular that SNL (which is also MCMC-based procedure, similarly to SeMPLE) was returning the worst performance. Here we show the performance of both SNPE-C and SNL with R = 4, again using total of 10 model simulations for cach method. That is, for each method we show the inference obtained after 2500, 5000, 7500 and 10000 model simulations. Figures 10 19 shows scatter plots (colored by density value) of the posterior samples from the last round (r=4) of each algorithm and for each the 10 observed data sets in SBIBM. SeMPLE consistently matches the reference posterior, unlike SNL and SNPE-C.
Moreover, it
interesting note that SeMPLE performs already well after single round, see Figures 20 29 showing posterior samples from every round of each algorithm and for each of the 10 observed data sets. While SNL improves across rounds, this is not the case for SNPE-C.
Figure 30 shows C2ST and runtime values in the case when R = used for all algorithms. Comparod to the results in the main paper with It = 10 rounds, SNPE-C and SNL perform slightly worse with rounds in terms of the C2ST metric. The SNPE-C and SNL runtimes are however reduced, with SNPE-C having essentially identical runtime SeMPLE, although with significantly higher variance. SNL is approximately twice slow the other algorithms.
25

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
Reference
Number of points per pivel
150
0.5
0.5
0.5
0.5
6
-100
0.5
0.5
0.5
-0.5-
-1.0
-1
0
Figure 11: Two Moons using R = rounds for every method. Density scatter plot of postcrior samples from the last round (r=4) of each algorithm. Observed data set number in SBIBM.
SeMPLE
SNPE-C
SNL
Reference
1.0
1.0
Number of points per pixe)
0.5
0.5
0.5
0.5
8
-0.5
-0.5
-0.5
-0.5
-1.0
-1.0
-1.0
-1
-1
61
Figure 12: Two Moons using rounds for every method. Density scatter plot of posterior samples from the last round (r=4) of each algorithm. Observed data set number in SBIBM.
SeMPLE
SNPE-C
SNL
Reference
1.0
Number of points per pixe!
0.5
0.5
0.5
0.5
0.0
-100
-0.5
-0.5
-0.5
-0.5
50
-1
Figure 13: Two Moons using R rounds for every method. Density scatter plot of posterior samples from the last round (r=4) of each algorithm. Observed data set number in SBIBM.
SeMPLE
SNPE-C
SNL
Reference
1.0
1.0
Number of points per pixe)
0.5
0.5
0.5
0.5
300
-0.5
-0.5
-0.5
-0.5
100
-1.0
-1
Figure 14: Two Moons using R rounds for every method. Density scatter plot of postcrior samples from the last round (r=4) of each algorithm. Observed data set number in SBIBM.
26

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
Reference
200
Number of points per pivel
0.5
0.5
0.5
0.5
-150
100
0.5
0.5
0.5
-0.5
50
-1.0
-1
Figure 15: Two Moons using R = rounds for every method. Density scatter plot of postcrior samples from the last round (r=4) of each algorithm. Observed data set number in SBIBM.
SeMPLE
SNPE-C
SNL
Reference
1.0
1.0
200
Number of points per pive)
0.5
0.5
0.5
0.5
150
-100
-0.5
-0.5
-0.5
-0.5
-1.0
-1.0
-1.0
-1
-1
Figure 16: Two Moons using rounds for every method. Density scatter plot of posterior samples from the last round (r=4) of each algorithm. Observed data set number in SBIBM.
SeMPLE
SNPE-C
SNL
Reference
1.0
Number of points per pixe!
0.5
0.5
0.5
0.5
-15
100
-.5
-0.5
-0.5
-0.5
-1
Figure 17: Two Moons using R rounds for every method. Density scatter plot of posterior samples from the last round (r=4) of each algorithm. Observed data. set number in SBIBM.
SeMPLE
SNPE-C
SNL
Reference
1.0
1.0
Number of points per pixe
0.5
0.5
0.5
0.5
110
-0.5
-0.5
005
-0.5
50
-1
Figure 18: Two Moons using R rounds for every method. Density scatter plot of postcrior samples from the last round (r=4) of each algorithm. Observed data set number in SBIBM.
27

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
Reference
1.0
1.0
1.0
1.0
Number of points per pixel
0.5
150
0.0
&
0.0
0.0
0.0
100
50
-1.0
-1.0
-1.0
-1.0
-1
0
1
-1
91
6
01
Figure 19: Two Moons using R = rounds for every method. Density scatter plot of posterior samples from the last round (r=4) of each algorithm. Observed data set number 10 in SBIBM.
28

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
25.00
,
5000 simulations
C
2
7500 simulations
10000 simulations
Figure 20: Two Moons using R = rounds for every method: inference using observed dataset number 1 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
29

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
25.00
5000 simulations
,
c
7500 simulations
10000 simulations
Figure 21: Two Moons using R = rounds for every method: inference using observed dataset mumber 2 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
30

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
?
25.00
e
5000 simulations
"
e
7500 simulations
10000 simulations
Figure 22: Two Moons using R = rounds for every method: inference using observed dataset number 3 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
31

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
25.00
5000 simulations
C
7500
10000 simulations
Figure 23: Two Moons using R = rounds for every method: inference using observed dataset number 4 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
32

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
25.00
5000 simulations
7500 simulations
10000 simulations
Figure 24: Two Moons using R = rounds for every method: inference using observed dataset number 5 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
33

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
?
25.00
5000 simulations
6
7500
10000 simulations
e
Figure 25: Two Moons using R = rounds for every method: inference using observed dataset number 6 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
34

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
3
25.00
3
5000
3
7500
3
10000 simulations
Figure 26: Two Moons using R = rounds for every method: inference using observed dataset number 7 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
35

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
e ?
25.00
?
5000
?
7500
10000 simulations
Figure 27: Two Moons using R = rounds for every method: inference using observed dataset number 8 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
36

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
25.00
5000 simulations
c
7500 simulations
10000 simulations
Figure 28: Two Moons using R = rounds for every method: inference using observed dataset number 9 from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
37

---

Published in Transactions on Machine Learning Research (06/2024)
SeMPLE
SNPE-C
SNL
True
25.00
5000 simulations
??
7500 simulations
10000 simulations
??
Figure 29: Two Moons using R 4 rounds for every method: inference using observed dataset number
from SBIBM. Posterior samples from each algorithm round, starting from the first round (upper row) to the fourth (bottom row).
38

---

Published in Transactions on Machine Learning Research (06/2024)
ScHPLE CWtt
30
0.80
24, with iounds
Curmuaatvve elapsed time [minutes]
0.75
0.70
15
0.65
10
2500
5000
75:00
10000
2500
5,000
7500
10000
Number simulations
Number of simulations
Figure 30: Two Moons using R rounds for every method. Median C2ST (left) and median eumulative runtime in minutes (right) for 10 runs with different data sets the number of model simulations. Shaded bands enclose the min and max values.
39

---

Published in Transactions on Machine Learning Research (06/2024)
F Multiple hyperboloid model
F.1 Multiple hyperboloid model definition
This model was introduced by Forbes et al. 2022). With reference to applications in audio processing, the unknown parameter 0 E TR.2 can be interpreted as the coordinates of a sound source location on the plane, while the two pairs of 2-dimensional parameters m (m,, m!)) and m (mi. m²) can be interpreted the location of two pairs of microphones aimed at detecting the sound source. An observation y E R has the likelihood function
1 S(y; Fm o') + (200))2((00 (m)) o² Ia,
pyy ())
(19)
2
where
Fn (110 millm 110 m2)l2), with m m2).
(20)
The likelihood equation 19 mixture of two Student's t-distributiins with a d-dimensional location parameter having all dimensions equal to Fm)(0), (0), Fm= (0) respectively, a diagonal scale matrix Lla and v degrees of freedom. In example, the observed data vector yo has length d 10, and 3, = 0.01, m = (-0.5,0) m = (0.5,0) m = (0, 0.5)* and m² = (0,0.5)² The prior is set to be a uniform 2() for both components of 0, and 0 is the only parameter to infer. Data generating parameter is 0" 0 01=51,,,1 1) as in Forbes et al. 2022) which represents the coordinates of the location of a sound source.
Since the multiple hyperboloid model is not implemented in SBIBM, there exists no observed data sets or reference posterior samples in SBIBM. To produce an observed data sot yo: the model simulator p(y (*) was run with the true parameter value 0- (1.5, 1). Since the posterior shape varies widely depending the true parameter value, and to keep the connection to the experiment in Forbes et al. (2022), only a single observed data set was produced. Instead, the metric results were averaged over 10 independent runs with this same observed data set. To produce a reference posterior sample corresponding to the produced observed data set, the Metropolis-Hastinng algorithm was used. The sample size of the reference posterior sample was set to 10a.
F.2 Multiple hyperboloid model: selection of K via BIC
Figure 31 shows the BIC for different values of K computed by fitting GLLiM, prior to rumning SeMPLE, using the same prior-predictive data set {θn. yn J.Y with N = 10000. Isotropic covariance matrices were used in GLLiM as it was found that using fully-specified covariance matrices did not improve the results. The decrease in the BIC slows down significantly at K=30, which motivates a starting value of K=40 in SeMPLE. This is close to the value of K-38 used in Forbes et al. (2022) (which was also selected for the same multiple hyperboloid model via BIC). The runtime to compute all the BIC values in Figure 31 was 270 seconds 4.5 minutes. Note once again that the BIC computation is independent of the observed data set and has to be performed only once before running ScMPLE.
F.3 Multiple hyperboloid model: results with R = 4
In the main paper we have considered the inference results based on a total of 4 x 104 model simulations, uniformly distributed across
rounds for SNPE-C and SNT (as per SBIBM defaults), and R= 4 rounds for SeMPLE. Results showed that SeMPLE performed better and in particular that SNL (which is also an MCMC-based procedure, similarly to SeMPLE) was returning the worst performance. Here we show the performance of both SNPE-C and SNL with R = 4, again using total of 4 x 10 model simulations for each method. Figures 34howw posterior samples from the last round (r=4) of cach algorithm for each of the 10 independent algorithm runs. We can clearly notice that SNL fails completely and therefore, for ease of reading, the SNL marginal posteriors are not reported. However, the differences between SNPE-C and ScMPLE are more subtle. It appears that occasionally ScMPLE samples from central arcas, where there should not be any posterior mass, while SNPE-C seems to undersample some of the branches. Therefore, There, it appears that the
we also plot the corresponding kernel-smoothed marginals in Figures 32 34
40

---

Published in Transactions on Machine Learning Research (06/2024)
-6e-04
BCC
-80-04
-1e105
K
Figure 31: Multiple hyperboloid. BIC with respect the number of Gaussian mixture components K.
marginals from SeMPLE overall closer the reference marginals, which in agreement with the better C2ST attained by SeMPLEE (see the main text).
Figure 35 shows C2ST and runtime values in the case when R = 4 rounds for all algorithms. The C2ST results of SNPE-C and SNL are very similar with R = 4, meaning that SNL fails completely and SeMPLE still performs much better than SNPE-C. In terms of runtime, SNPE-C is still about 7 times slower than SeMPLE, disregarding SNL as the inference from this method not useful.
41

---

Published in Transactions on Machine Learning Research (06/2024)
p(o, Oslvo)
Ozlyo)
0.5
0.5
0.3
0.2
0.2
0.1
0.0
81
61
01
0.5
(o)
61
pl0 lyo)
0.4
0.4
The postertor
Tree posterior
SAMPLE with rounds
0.3
SaHf with minds
SKIE Citt 4 rounds
SMPE C with rounds
with rounds
SML with rounds
0.1
0.1
0.0
-2
-2
82
62
(a) Repetition 1
(b) Repetition 2
p(8. 02lyo)
plo,. 0zlyu)
p(0 (ya)
p (o)
0.4
0.4
0.1
0.1
0.0
2
8,
-2
6
0.5
p(0:lya)
61
0.5
p(0:
0.3
0.3
SNL with rounds
SNL mounds
0.2
0.2
0.0
0.0
6a
82
(c) Repetition
(d) Repetition
Figure 32: Multiple hyperboloid model. Pair plots posterior samples from the last round (r = 4) of cach algorithm. The SNL marginal posteriors are not reported for of reading, since the SNL inference fails.
42

---

Published in Transactions on Machine Learning Research (06/2024)
p(o, OslVo)
Ozlyo)
0.5
((9119))
0.5
((0))o)
0.3
0.2
0.2
02
0.1
0.0
81
81
01
0.5
p(02/9ra)
61
p(0:
0.4
0.4
The postertor
Tree posterior
SWMPLE with rounds
0.3
SaHf with minds
SKIE Citt 4 rounds
SMPE C with rounds
with rounds
SML with rounds
0.1
0.1
0.0
-2
-2
82
62
(a) Repetition 5
(b) Repetition
p(8. 02lyo)
plo,, 0zlyo)
p(0 [ba)
p()/1o))
0.4
0.4
0.1
0.1
0.0
2
8,
-2
6:
0.5
((8:lyr)
61
0.5
p(0.lvo)
0.3
0.3
SNL with rounds
SNL mounds
0.2
0.2
0.0
0.0
-2
0a
82
(c) Repetition
(d) Repetition
Figure 33: Multiple hyperboloid. Pair plots of posterior samples from the last round (r 4) of cach algorithm. The SNL marginal posteriors are not reported for ease of reading, since the SNL inference fails.
43

---

Published in Transactions on Machine Learning Research (06/2024)
p(o,. OalMo)
plo.. 9zlvo)
0.5
((0//Vo)
0.5
0.4
0.4
0.2
0.2
0.1
0.14
61
81
0.5
p(0)/tra)
61
p(0:br))
0.4
0.3
0.3
SNPE rounds
SHPE with rounds
SKL with rounds
SNL with rounds
0.2
0.1
0.0
0.0
82
(a) Repetition 9
(b) Repetition 10
Figure 34: Multiple hyperboloid. Pair plots of posterior samples from the last round (r = 4) of each algorithm. The SNL marginal posteriors are reported for ease of reading, since the SNL, inference fails.
SENPLE rourtb
Curulaative elspred time [minutess]
3.0
1.00
0.8
0.5
20000
Number of simulations
Number of simulations
Figure 35: Hyperboloid. Modian C2ST (left) and median cumulative runtime in minutes (right) for 10 runs with the same data vs the number of model sinulations. Shaded bands enclose the min and max values.
44

---

For all inference algorithms, we executed runs with the 10 different observed data sets provided in SBIBM. The total simulation budget always set to 10, 000 model simulations, but two different setups were used for the number of algorithm rounds. In the first sctup all algorithms used R = 2 algorithm rounds and in the second setup SNL and SNPE-C used R 10 rounds while SeMPLE still used R = 2 rounds. The reason for running R = rounds with SeMPLE was that for this model the output from r = 2 produced accurate inference but had some stickiness in the Markov chain that caused problems in the GLLiM learning The GLLiM Published in Transactions on Machine Learning Research (06/2024)
G Bernoulli GLM model
G.1 Bernoulli GLM model definition
For the model definition, please see section T.5 in the appendix of Lueckmann et al. (2021 The prior distribution and model parameters can also be found there.
G.2 Inference setup
in round r 3. Hence, the best inference results with ScMPLE were obtained with R covariance struct
set to be full unconstrained matrix as our experiments showed that this improved the inference results compared to constrained covariance matrix structure. Figure 36 shows BIC for several values of K, computed by fitting GLLiM using the same prior- predictive data set (On, Vn J with N = 2500 prior to running SeMPLE. The minimum value is obtained at K = 2 but because of this small value, A was instead set to 10 as this gives more flexibility to the model without resulting in significantly longer runtimes. The runtime to compute all BIC values Figure 36 was 16 seconds. Note that the BIC computation is independent the observed data set and only has to performed once before running SeMPLE.
5X00
2000
2000
1500
1900
5000
Figure 36: Bernoulli GLM model. BIC with respect to the number of Gaussian mixture components K.
G.3 Results when using R R== for all algorithms
Here we compare results via SNL, SNPE-C and SeMPLE with the true posterior using = 2 rounds for all algorithms. Since we consider a total budget of 10,000 model simulations, Figure 38 shows the C2ST and Wasserstein performance evaluated at 5,000 an 10,000 model simulations. We found that C2ST can be problematic with unimodal targets, and we this also with the has unimodal posterior), where the performance of the several algorithms seems to be better evaluated using Wasserstein distances. In fact, as from Figure 37 it is clear that the results from ScMPLE at = 2 are noticeably improved compared to r = 1, and this is captured by the Wasserstein distances in Figure bS,,,, whereas C2ST produces a contradicting result in Figure 38(a). We therefore alert the reader that with unimodal targets C2ST at times not always consistent with other measures, such as posterior plots (sce also the Constin-UUlenbee section), while we found C2ST useful with multimodal targets. Moreover, the inference from SeMPLE is not only accurate, but the computational gain provided by SeMPLE is striking,
Uhlenbeck results (which also
45

---

Published in Transactions on Machine Learning Research (06/2024)
Figure 37: Bernoulli GLM model. Posterior plots obtained with ScMPLE and obtained round after 5,000 model simulations, and at after 10,000 model simulations.
as illustrated in Figure 39 where ScMPLE 8 times faster than SNPE-C and 19 times faster than SNL. These figures increase dramatically for SNL and SNPE-C next section, where these methods are run with the default SBIBM setup.
JWCCC
34
distrnee
14-
2.5-
66
39W3 04-
106
1999
199
simulatiors
(b)
Figure 38: Bernoulli GLM model with = for all algorithms. (a) Median C2ST of 10 runs with different data sets vs number of model simulations. (b) Median Wasserstein distance of 10 runs with different data sets vs number of model simulations. Error bars show min/ /ax values. For unimodal posteriors (as in this case) Wasserstein distances should be preferred to C2ST.
46

---

Published in Transactions on Machine Learning Research (06/2024)
Alpesttm SAMPLI SUTE-C
20
SNL
1020
Figure 39: Bernoulli GLM model with R 2. Median running times of 10 runs using several mothods. Error bars show min/mas values.
68-
Wasserraeinn
03-
Number of simulations
Number simulations
(a)
(b)
Figure 40: Bernoulli GLM model with R= 10 (SNL/SNPE-C) and = (SeMPLE). (a) Median C2ST of runs with different data sets vs number of model simulations. (b) Median Wasserstein distance runs with different data sets vs number of model simulations. Error bars show min max values. For unimodal posteriors (as in this case) Wasserstein distances should be preferred to C2ST.
G.4 Results when using R = 10 for SNL/SNPE-C and R = for SeMPLE
We executed the experiment anew, but this time by allocating the 10 model simulations uniformly across 2 rounds for ScMPLE and 10 rounds for SNL and SNPE-C (as SBIBM defaults). Results are in Figure 40 and Figure 41 Once more, we recommend the reader to consider the Wasserstein distances and not C2ST for this example, since that posterior is unimodal, as motivated in the previous section. The distances show that on average SeMPLE is as accurate as SNL but not before 10a model simulations are executed, while SNPE-C shows large variability in the results. Once more, SeMPLE comes with the benefit of much faster inference (Figure 41 in this case SeMPLE is 26 times faster than SNPE-C and 74 times faster than SNL.
47

---

Published in Transactions on Machine Learning Research (06/2024)
Cumuuame fime
3000
4002
5000
8003
10200
Numbher crstmnintionss
Figure 41: Bernoulli GLM model with R= 10 (SNL/SNPE-C) and R= 2 (SeMPLE). Median running times of 10 runs using several methods. Error bars show min/ /max values.
48

---

Published in Transactions on Machine Learning Research (06/2024)
H Ornstein-Uhlenbeck
The Ornstein- Ubieeinn-hheebeck process has original applications in physics to model the velocity of a Brownian par- ticle (Uhlenbeck & Ornstcin 1930), and has found application in many contexts such as neuronal modelling and finance, to mention a lew. In our work, it serves as an application where the data is a time series, while the likelihood function is tractable and therefore it is possible to obtain samples from the true posterior, to use as reference. The (SDE)
process {Xr) defined by the stochastic differential equation
dx.- = - aadt + adW,
(21)
where / > 0, a R.
0, Xo zo and W, denotes the Wiener process (Brownian motion). The solution
to equation 21 is
X, (mo ale a)e-
8t +
e-Bttts)dW,
(22)
The process has known Gaussian transition densities: for a generic time t and an arbitrary state xo at time to = 0, the transition density is given as
(x, Xo 20) N
N((2 a (No aae
o² 23 (1-
- -88))
(23)
The transition density 23 allows exact simulation of a solution path to the -eattinnUhhlenbeck SDE. Starting at lo 0 one can simulate 1 further points of the solution at instants (,,..., ln-1, letting A = t - t,-1, by using the transition recursively as following
(X+ Xt, =x, I4, 1)_~~~~ N a (No (Na(NNNN(
a)e- SA,
2 23 (1-e e -8811) 1= 1.
0²
(24)
This makes the Ornstein-Uhlenbeck the exact likelihood function. The likelihood function for the parameters = (ax, 3, a) is written as (due to the Markovianity of the solution of an SDE)
Uhlenbeck model one of the very few examples of SDEs where it is trivial to evaluate
C(zo.... In 0)) = II p(æ x)) c-
(25)
where p(z; 0) is the transition density in equation 24 (we assumee zo fixed constant throughout,
otherwise a multiplicative term p(zol0) would have to appear in (25) The interpretation of the parameter a is the asymptotic mean. This can be seen by letting L + co in equation (23 which shows that the process'
stationary distribution is X, Naa, )²) The paramcter can be interpreted as the variation the size of the noise, while 3 can be interpreted as the growth-rate or how strongly the system reacts to perturbations.
H.1 Inference setup
The fixed initial value was set to 2g = O. We set uniform priors a U(0, 10), B ~ (0, and o ~ U(0,2) 2). In this study, 50 points were simulated from the process to obtain a time series of length n = 51 (including the fixed starting value ιо) using 0 = (a,3,a) = (3, 1, 0.5) as ground-truth. The time frame of the process was set to T] = [0, 10], and the discrete time points f1, fn-1 set to be equally spaced in this interval. The resulting observed data is shown in Figure 42 Similarly to the multiple hyperboloid model, a single observed data set was produced by running the simulation model once with the true parameter values 0 = (3, 1, 0.5).
The covariance structure in GLLiM was set to be full unconstrained matrix, which is reasonable given the time-dependent structure, and in your experiments it was found that this improved the inference results substantially. Figure 43 shows the BIC for different values of K, computed by fitting surrogate likelihoods using the same prior-predictive data set (0n, yn with N = 10 000 prior to rumning SeMPLE. The minimum value is obtained at K but because of this relatively small value, K was instead set to 20 to have some buffer remove unneccssary mixture components if needod. The runtime to compute all BIC values in Figure 43 was 973 seconds 16 minutes. Note that the BIC computation is independent of the observed data set and only has to be performed once before running SeMPLE.
49

---

Published in Transactions on Machine Learning Research (06/2024)
8-
8-
Figure 42: Ornstein-Uhlenbed Observed time series data of the process.
-1000000
BIC
-1150000
-1300000
2
12
K
Figure 43: Ornstein-Uhlenbeck:: BIC as a function number of Gaussian mixture components K.
H.2 Results
The posterior unimodal and therefore, following the results for the Bernoulli GLM example, we find it useful report both the C2ST metric and Wasserstein distances. For all inference algorithms, we executed 10 independent runs using the same observed data set. We experiments using two configurations by always allocating the same total number of 4 x 10 model simulations for all methods: the first configuration uses R 4 rounds for all methods (section H.3 The second configuration uses R = 4 rounds for SeMPLE and R = 10 for both SNL and SNPE-C (section H.4), as R 10 is the default implementation in SBIBM for SNL and SNPE-C. is instructive to look first at the case where R = 4 is used for all algorithms, as there we can make useful remarks about the behaviour C2ST and the Wasserstein distance.
H.3 Results when using = for all algorithms
Figure shows the C2ST metric and Wasscrstein distance a function of the number of model simula- tions for each algorithm. At round r-2 (20000 000 model simulations) SNL consistently shows an unexpected behaviour where the posterior sample is clearly worse than in the first round (r=1). This happened in some scenarios with SNL where the number of simulations in the first round low enough. Apart from this, SNL scores best in both metrics after the full 40 model simulations (R=4 algorithm rounds). However, notice that the difference between SeMPLE and SNL is much smaller in terms of the Wasserstein distance.
50

---

Published in Transactions on Machine Learning Research (06/2024)
In fact, by looking at the marginal posteriors (Figures 46-45 49 approximate the true posterior, SeMPLE is strong second best alternative, and provides clearly better posteriors compared to SNPE-C. However the latter result is not apparent when looking at the C2ST met- ric, where ScMPLE and SNPE-C do not seem strikingly different. Thercfore, for this experiment with a unimodal larget, similarly to the Bernoulli GLM example, we find that the Wasserstein distance provides better summary of the quality of the posterior inference. The Prnstein-Uhlenbee process could be a simulation model where it is easier estimate the likelihood compared the posterior, which could explain the outcome of metric scores since both SNL and SeMPLE uses surrogate likelihood to obtain posterior samples.
it is clear that while SNL posteriors better
1.0
Algerithm
Alqoritim
+ SuMPLE
SNPE-C
-9 SNIE-C
-
SAMPLE
SNL
0.9
Wasserstein distance
1.50
1.25
0.8
C25
1.00
0.75
0.50
0.00
20000
30000
Number of simulations
Number of simulations
(a)
(b)
Figure 44: Oristein-Uhlenbekk, results when R = for all algorithms: median C2ST and Wasserstein distance across 10 independent runs (with the same observed data set) as a function of the number of model simulations. Error bars show min/me /max values. For unimodal posteriors (as in this case) Wasserstein distances should be preferred to C2ST.
In terms of runtime, Figure shows that SeMPLE has a much much smaller runtime compared to both SNT, and SNPE-C, even if were to include the 16 minutes of BIC computations into account. Especially SNPE-C has much longer runtimes with this model, combined with the worst performance metric results out of all three algorithms. Plots of the posterior marginals and posterior draws from the last round (r=4) of each algorithm are in Figures 46 47 for each of the 10 independent runs. Acceptance rates of the Mctropolis- Hastings step of SeMPLE are in Figure 48 showing that, particularly for the third and fourth round, the acceptance rate is fairly high (compatibly with being an MCMC algorithm).
H.4 Results when using R = 10 for SNL/SNPE-C and R = 4 for SeMPLE
We executed the experiment anew, but this time by allocating the 4 x 10 model simulations uniformly across 4 rounds for SeMPLE and 10 rounds for SNL and SNPE-C (as per SBIEM defaults). Results are in Figure more trustworthy to consider
49 and Figure 50 As explained in previous sections, for this experiment the Wasserstein distances, hence Figure 19 in this case. Once again we show that, even by comparing with SNL and SNPE-C trained over many inference rounds (as per SBIBM defaults), SeMPLE is competitive, Moreover, in terms of runtime, if we include the 16 minutes required to run the BIC procodure to select K (recall this is amortized it can be reused with other observations). the median runtime for SeMPLE is around 21 minutes, while SNPE-C requires around 136 minutes to complete the 10 rounds on average. Hence SeMPLE implies more than 6-fold acceleration, which becomes around 27-fold if we exclude the BIC procedure.
51

---

25002 Published in Transactions on Machine Learning Research (06/2024)
150
Algeeinm
Cumukkove depuoll Line (ni-itol
80
60-
60-
20
12000
15000
20203
23020
10000
40200
Kunbher instinintionss
Figure 45: Ornstein-Uhlenbeck, results when R for all algorithms: median cumulative runtime of 10 independent runs as a function of number of model simulations. Error bars show min/ max values.
52

---

Published in Transactions on Machine Learning Research (06/2024)
SHAPLZ
2.5
(a) Repetition
(b) Repetition
SHIPLE
1.0
(c) Repetition
(d) Repetition
235
200
1.0
(e) Repetition
(f) Repetition
Figure 46: Ornstein-Uhlenbeek. results when R 4 for all algorithms, runs 1 6.; marginal posteriors and posterior draws from the last round (R=4) of each algorithm.
53

---

Published in Transactions on Machine Learning Research (06/2024)
94 SANPLE
(a) Repetition
(b) Repetition 8.
SAA 20
33
00
(c) Repetition 9.
(d) Repetition 10.
Figure 47: -blenben results when R 4 for all algorithms, runs 7-10: marginal posteriors and postcrior draws from the last round (R=4) of each algorithm.
0.8
Acceptance rate
0.5
0.4
0.2
Iteration
Figure 48: Orastein-Uhlenbek: Median acceplance rate of the Metropolis-Hastings step of 10 independent SeMPLE runs with the same data set. Error bars show min/max values.
54

---

Published in Transactions on Machine Learning Research (06/2024)
2.0
0.1
009
srice
19000
25000
1at 35000
300K
1550
250en
3950 19000
1003
simulations
simulations
(a)
(b)
Figure 49: Orustein-Uhlenbeed results with R 4 (SeMPLE) and R = 10 (SNL/SNPE-C): median C2ST and Wasserstein distance across 10 independent runs (with the same observed data set) as function of the number of model simulations. Error bars show min/max values. For unimodal posteriors (as in this case) Wasscrstein distances should be preferred to C2ST.
140
Algoritmm
Cumurtivee dlased Sane (nsriiaa)
20-
20
80
15000 20200 Numberr cfstmilations
Figure 50: Orastein-UUlenbbeck, results with R = 4 (SeMPLE) and R = 10 (SNI//NPFF-C: median cumu- lative runtime of 10 independent runs as a function of number of simulations. Error bars show min/max values.
55

---

Published in Transactions on Machine Learning Research (06/2024)
(Markov jump process with measurement error)
The Lotka- Volterra predator prey model is population dynamics model with two specios. Its dcfinition varies, depending whether it is defined via ordinary differential equation (as in Lueckmann et al. 2021 continuous time Markov Jump Process (MJP, Owen et al. 2015a bj or a continuous time approximation to the MJP given by a stochastic differential equation (as in Golightly & Wilkinson 2011). Moreover, it
could be seen as a model involving three (as in Owen et al. 2015a] or four (as in Papamakarios & Murray
2016) reactions. We employ a stochastic MJP version of this model as described in Owen et al. (2015a) characterized by a of 3 reactions on the two species, with additional mcasurement crror.
1.1 Model definition with three reactions
Denoting the two species, preys in number X1 and predators in number X2, evolution of the system is governed by the following three reactions:
R
X1
+ 2X1
R₂ X1
X22 7 2X₂
(26)
R3: X2
01
These reactions can be interpreted respectively as a prey birth, a predator prey interaction resulting in the death of prey and predator birth, and predator death respectively. Let X, = (2,,) denote the number of both species at time t. Each reaction R, is assumed to have an associated constant 0 and hazard function hx,X,, 0:)) which gives the propensity for reaction event of type i at time t to occur. Reaction events are dependent on the current state of the system as well as the reaction rale parameters. Hence the evolution of species counts corresponds to a Markov process on a discrete state space. This reaction network is summarised by its stoichiometry matrix, S and hazard function h(X,0).
-1 0 S =
0 1 )
h(X,0) 0) = (0, 0,,
(27)
1.2 Inference setup
To simulate the exact dynamics, we use the smfsb package (Wilkinson 2024). implementing the "Gillespie algorithm" " (Gillespic 1977). Whilst the Gillespic algorithm allows exact simulation of the time and type of each reaction event that occurs, observed data Y (yo- 81. 31...., .......)) (r) are typically noisy observations at discrete time intervals, where each y/ is 2-dimensional with of x(X,,,, 0) ii.dd for j = 1,2. We corrupt each X5.4 with itidd Gaussian error with mean and variance o², n T Xj,,,²) ~ NNXX,,,, The true parameter values generating the observed data is set to (1, 0.005, 0.6) and the starting values Xo- = (50, 100) are used throughout. When generating the observed data let 30 and later infer along with the model parameters 0, resulting in a total of 4 parameters.
When computing model realisation with the Gillespie algorithm, we let T = 30 and observe the Markov jump process at uniformly spaced intervals with step length 0.2, resulting in a time series of length 150 for each species. After corrupting the time series data with the previously described Gaussian error we compute a total of 9 summary statistics that is used when performing the inference. These summary statistics, first suggested in Wilkinson (2013), consist of the mean, (log) variance, and the lag 1 and lag two auto- correlations of the time series of cach species. Additionally the cross-correlation between the two time series is included in the summary stalistic.
Following the setup in Owen et al. (2015a), we place uniform prior information on log(!),
log(0,) (((-6,2)) i= 1,
(28)
and log(a)
log(a) ~ (((oo((0.5) log(50))
(29)
56

---

Published in Transactions on Machine Learning Research (06/2024)
Prior to running the inference, for each of the nine summary statistics we computed trimmed means, where 1. 1.25% of the values were trimmed from each end and corresponding trimmed standard deviations, obtained from 10,000 prior predictive simulations, and we used these to standardise both the observed summaries and the summarics simulated during the inference. For each SeMPLE run we set a budget of 80000 model sirmulation, uniformly distributed across algorithm rounds. The number of mixture components is set to K = 15 (a number which was picked without any specific tuning) and we set the covariance matrices in GLLiM to be full unconstrained matrices, all different across the K components. Moreover, we do not inflate the covariance in the Meesspgs step (ie we set Y 1). The mixture probability threshold to remove negligible mixture components during the SeMPLE run was set to 0.005.
For a MJP Lotka- Volterra model, we did not find a ready implementation in SBIBM (there, an ODE-based Lotka- Voltera model is provided, not a MJP version), and this is likely because with such a model it is impossible to perfectly control the computational effort, as for a given value of the time spent to complete a trajectory is random. For this reason, we found more practical to compare ScMPLE with approximate Bayesian computation (ABC) experiments using the code from Picchini & Tamborrino (2024). For SMC- ABC, we used the same observed data and the same standardised summary statistics as with SeMPLE. We consider sequential Monte Carlo ABC (SMC-ABC), and also use the posterior from SMC-ABC as "reference posterior". SMC-ABC is the state-of-art approach to ABC inference and is the algorithm of choice in recent inference platforms such as ABCpy (Dutta et al. 2017 and PYАBC (Schälte et al., 2022 We consider two versions of SMC-ABC that differ in terms of proposal function: the first one is the "classic" SMC-ABC, in that it uses as proposal distribution a Gaussian random walk of the type l NN 22) (Beaumont et all 2009 Filippi et al. 2013), where 0* is a parameter "particle" randomly sampled from the most. recent set of "paramcter particles" " (these providing an ABC approximation to the posterior), and E is the empirical weighted covariance matrix of these "parameter particles". This first sampler is named standard in [Picchini & Tamborrino (2024), as it is the most common SMC-ABC proposal sampler, as found in most literature and in available software packages. A second SMC-ABC sampler is the one denoted blockedopt in PPichhini & Tamborrino (2024), which is also a Gaussian proposal sampler but not of random-walk type, and is much more resource-effective as it more rapidly "guide" the parameters towards the bulk of the posterior region (see Picchini & Tamborrino 2024 for details). We ran both SMC-ABC versions with 1,000 particles.
1.3 Results
For this case study, we consider the posterior returned by the blockedopt SMC-ABC as "reference posterior", as it is obtained after more than three million model simulations (specifically 3,337,648 ,64 simulations), that is compulational budget more than 100 times larger than the budget used for SeMPLE. In Figure 51 we show that, while the "guided" blockedopt SMC-ABC can infer (01, 02 03) for a smaller budget than millions of simulations, however inferring also a does requires many additional simulations. When we use the standard SMC-ABC, results after 3,110,349 model sirmulation are in Figure 52. and it is clear that the variability in the inference is much larger, and in particular standard struggles in inferring a.
It is of interest to compare results for ScMPLE and SMC-ABC for a similar number of model simulations. We anticipate to observe something that is well known already from Greenberg et al. (2019): Lueckmann et al.) (2021), , namely the SMC-ABC inefficiency in terms of model simulations. What we mostly wish to verify is the ability for SeMPLE to return accurate inference. In Figure 53 we compare results when using around 30,000 model simulations, however notice that it is impossible to compare the mothods for the exact same budget, as SMC-ABC employs while-loops that keep repeating model simulations until a prefixed number of parameter proposals accepted. It clear that when the number of model simulations limited to around 30,000, neither of the SMC-ABC versions can produce accurate inference. The SMC-ABC methods struggle especially with estimating the noise parameter a when the number of model simulations is limited. The ScMPLE posterior after 30,000 model simulations is more similar to the blockedopt SMC-ABC posterior obtained after approximately 3.3 million model sinulations, showing the computational efficiency of SeMPLE with a much smaller simulation budget.
To give an idca of the variation the ScMPLE posterior sample output, Figure (54 we report 5 independent SeMPLE runs with the same observed data as in previous analyses. We see that ScMPLE consistently provides accurate inference of the model parameters 0 in multiple repetitions. There is larger variation
57

---

Published in Transactions on Machine Learning Research (06/2024)
loge,
logê
2
0
-2
-4
-6
10
20
30
39
10
20
30
39
iteration
iteration
log0g
loga
2
3
0
-2
10
20
30
39
10
20
30
39
iteration
iteration
Figure 51: Lotka- Volterra: evolution of ABC posteriors using the "guided" blockedopt SMC-ABC across
39 rounds and 3,337,640 model simulations. Ground-ruthh parameters are marked with horizontal lines.
loge
loge,
2
0
0
-2
-6
1
10
20
30
40
50
64
10
20
30
40
50
64
iteration
iteration
log0
logo
2
3
0
-2
0
1
10
20
30
40
50
64
10
20
30
40
50
64
iteration
iteration
Figure 52: Lotka-Volterra: evolution of ABC posteriors using standard SMC-ABC across 64 rounds and 3,110,349 model simulations. Ground-truth parameters marked with horizontal lines.
58

---

Published in Transactions on Machine Learning Research (06/2024)
SeC ABC Rnciestee 17ME D4C ABC soaeurd
SMCARC enoceedoa 1:44
DMGARC cotandrrd
log10.)
lag(62)
04: Anectetie LIMC
2.5
OPC-3EC Instereet
(40 ABC condur I.M
2.0
0.5
log(
log
Figure 53: Lotka- Volterro marginal postcrior plots from SMC-ABC and ScMPLE. The ScMPLE postcrior sample based 10,000 posterior draws, and returned after lotal of 30,000 model simulations. The blockedopt SMC-ABC posterior samples from round 10 and 39 used total 32,865 and 3,337 640 model simulations respectively. The posterior samples from standard SMC-ABC used a total of 30,436 model simulations.
in terms of the noise parameter a, suggesting that this is more difficult to infer with the limited budget of model simulations. However, when comparing with the SMC-ABC posterior distributions for a, see Figure
the SeMPLE results much more informative for the given model simulation budget of 30,000.
59

---

Published in Transactions on Machine Learning Research (06/2024)
SAPPLE Pepue
2-
0
0-
-6
-4
-6
-2
log(8.)
log(82)
2.5
SenF-5 mpeet
S6n-1 mpaet
1.5
1.0
0
0.0
2
2
lag(02)
lag(a)
Figure 54: Lotka-Volt Volterra: marginal posterior plots from 5 repeated runs of SeMPLE using the same observed data. The posterior sample based on 10,000 posterior draws, and was returned after a total of 30,000 model simulations. The x-axes are adjusted to the support the uniform prior of each parameter.
60

---

Published in Transactions on Machine Learning Research (06/2024)
J Twisted prior model with 20 parameters
This case study is merely illustrative of the ability of ScMPLE to rapidly return sensible inference for a posterior with highly correlated components. Unlike with other examples, here we did not perform a full comparison with other methods, as we found it difficult to specify the non-standard prior used here to work within the framework provided by SBIBM.
J.1 Model definition
We now consider a model with 20 parameters and a challenging posterior charactorised by a strong correlation between some of the parameters. This case study is particularly interesting, as the likelihood only provides location information about the unknown parameter while the dependence structure in the posterior comes mostly from the prior. For this reason, the posterior dependence changes direction depending on whether the likelihood locates the posterior in the left right tail of the prior (see Nott et al. 2018 for graphical illustration). This case study was analysed in ABC context in Li et al. (2017) and Picchini & Tamborrino 2022h. The model assumes (0,,.... Os,) and diagonal covariance matrix V = diag .i.,()) The prior is the "twisted-normal prior
(U1,..., Dd,) drawn from dg-dimensional Gaussian y N(0, V) with
of Haario et al. (1999), with density function proportional to p(0) x exp - -09 - (02 - 608 + 1006)2/2
1 (do22)) E5=3 / where Lg denotes the indicator function of the set B. This prior essentially a product
of independent Gaussian distributions with the exception that the component for (61, 02) is modifiod produce a "banana shape", with the strength of the bivariate dependence determined by the parameter b. Simulation from p(0) is achieved by first drawing 0 from dg-dimensional multivariate Gaussian O~.N(0, A), where A diag(100, I, 1), and then placing the value 02 + 60} 1006 in the slot for 0z. We set b = 0.1, a value inducing a strong correlation in the prior between the first two components of 6. We consider dg = 20, i.e. both y and 0 have length 20, and (same as Li et al. 2017 Picchini & Tamborrino, 2024) we set 0g = as fixed and known, with observations given by the 20-dimensional vector yo = (10, 0, .... where only the first entry is non-zcro.
J.2 Results
We now show how a single round of ScMPLE (r = 1), hence without any MCMC, can in this case return excellent inference (notice, generally we recommend to run at least 2 rounds of SeMPLE, as MCMC is typically needed ensure that samples from r do not "leak" outside the prior's support). Figure 55 compares the banana-shaped prior for the first two components of 0 with draws from the true posterior obtained via MCMC (panel (a)), and then with draws from single round (r 1) of SeMPLE initialised with K = 10 components (panel (b)). Here SeMPL.E has been trained on 10,000 prior predictive simulations. Despite the true posterior occupies tiny region within the prior's support, SeMPLE immediately identifies the bulk of the posterior one round, that is Figure 55(b) shows samples via the learned GLLiM posterior model. Then, in Figure 56 (a) we zoom on the posterior region to better show the quality of the SeMPLE approximation for both latter being an illustrative example of the non-correlated components (inference behaves the same for all the remaining components, hence further marginals are not reported). It appears that, for this model, the correction brought by MCMC at = 2 is not necessary.
and r 2: see Figure 56 a) for (01,02) and Figure 56 for (0s, 04), the
61

---

Published in Transactions on Machine Learning Research (06/2024)
prior draws and true posterior
prior draws and SeMPLE posterior (1 iteration)
150
150
100
100
thela2
50
440
-20
20
-40
-20
0
20
theta1
(a) Draws from prior p(0,, 02) (black) and exact posterior (b) Draws from prior p(0,, 02) (black) and from ScMPLE (green)
(brown)
Figure 55: Twisted prior example: prior, true posterior and draws from SeMPLE at round
True posterior SeMPLE r=1 ScMPLE r=2
True posterior SeMPLE ScMPLE r=2
(a) Posteriors for (01, 02)
(b) Posteriors for (0s, 04)
Figure 56: Twisted prior example: (a) Contour plots of the exact posterior for (01, 02) (green) and SeMPLE posterior at = 1 (brown) and r = 2 (blue). (b) Same as (a) but for (0s. 02).
62

---

Published in Transactions on Machine Learning Research (06/2024)
K
Biological model of the translation kinetics after mRNA transfection
This model serves an example of a multidimensional SDE model with higher number of model parameters, compared to the Ornstein-Uhlent (only one coordinate observed), and observations incorporate measurement measurement error. Reference posterior samples cannot be obtained exactly, as the transition densities are not known for this model.
model. Moreover here the model is challenging as it is partially observed
K.1 Model definition
We consider the SDE model from Pieschner et al. (2022) in a simulation study. In Picschner et al. (2022) the model was used to describe translation kinetics after mRNA transfection. Specifically, in a real data scenario, time lapse microscopy images of the cells are taken for several hours, and for the first hour, the cells would be incubated with mRNA lipoplexes. Released mRNA gets translated into green fluorescent protein (GFP) which causes the cells to fluoresce. The production of GFPs is the translation process. The SDE below is used to study the translation kinetics of one cell, based on one observed fluorescence trajectory, where the two-dimensional stochastic process (m((),p(t)) has m(t) to represent the amount of mRNA molecules at time t, and p(t) is the amount of GFP molccules at time t. The SDE model is given by
m (m))
- -d-m(() (1) (.tm(t) - r p(t) ) di + (Vo. 0
dBc
(30)
Vk m(t) +Y p(t)
where dB, is a two-dimensional standard Brownian motion. It is assumed that all mRNA molecules (within one cell) are released at once from the lipoplexes and denote this initial time point by lu- Before lus there are neither mRNA nor GFP molecules, and at lo; an amount mp of mRNA molecules released, i.e
mll)) p(t)
for t << t0
while m(ta) = mo and p(0) = 0.
We take as observable mapping
v(44) = log(scale: p((4) + offset) e(4))
i 1,..,,
where e(tx) is iid Gaussian measurement error =(t,) N(0, o²), and ti,..,, tn are observational time instants (see further below). Note that the observations depend only indirectly on process (m(()), and only {p(t)) is observed. The parameter vector we wish to infer is 0 = (6, 7, k, mo, scale, la offiset, o).
It has been shown by Pieschner et al. (2022) that SDE modelling improves identifiability of parameters, compared to using a corresponding ODE model. Nevertheless, some parameters are still unidentifiable, see in particular supplementary section A.4.1 in Pieschner et al. (2022).
Solutions to the SDE 30) are simulated using Euler- Maruyama scheme implemented in Rcpp Eddel- bucttel et al. 2024) from t = to to = 30, with step size 0.01. The observed time series is interpolated from the Enler-Maruyama approximation at = 60 equidistant time points from t1 = 0.5 to t, 30. Note that this implies that observations can be made prior to ( = to where p(4:) = 0, l < lo-
Inference is performed on the parameters log-scale, and the prior distributions are all Gaussian as follows
log 8 (V-0.694, 0.5),
log ~ W3(3, 0.5)
log A (0.0277, 0.5),
log ~ ((5.704
log scale ~ (0.7551 0.5),
logto~>V (o(_..66001 0.5)) loga~~ N(-2,
log offset ~ (/ 0.5),
ScMPLE, SNL and SNPE-C were run for five different observed data sets, each generated with different set of parameters given in Table 3. For each method, the total simulation budget was set 30k model simulations that were uniformly distributed across the number of algorithm rounds. For SeMPLE, the
63

---

Published in Transactions on Machine Learning Research (06/2024)
number of algorithm rounds set to R 3, the number of starting mixture components set to K 10, and the covariance matrices in GLLiM were set to be full unconstrained matrices. We did not
inflate the covariance in the Hitropolis-Hastinng step (7 =). 1). The mixture probability threshold to remove
negligible mixture components was set to 0.005. For SNPE-C and SNL the number of algorithm rounds was set to 10.
log
log7
log k
log mg
log scale
log ln
log offset
log -2
Data set 1
-0.694
-3
0.027
5.704
0.751
-0.164
2.079
Data set 2
-1.827
-3.138
-0.426
6.869
0.697
0.465
2.787
-1.993
Data set 3
-1.153
-3.501
-0.181
4.601
0.386
-2.121
2.185
-2.578
Data set 4
-1,297
-2.134
1.118
5.577
0.938
0.12
1.122
-2.854
Data set
-0.243
-4.596
0.168
5.573
-0.086
1.234
1.392
-1.655
Table
Translation kinetics model: paramoter values used generate observations for the five data sets.
K.2 Results
compare the marginal posteriors obtained by SeMPLE, SNPE-C and SNT,, where each figure corresponds to a different dataset. All results are obtained after 30k model simulations but the number of algorithm rounds 10 for both SNPE-C and SNL, while we ran 3 rounds for SeMPLE. Inference very similar across the three methods, however it is clear that only SeMPLE and SNPE-C provide biased estimation of to (an exception being the inference for dataset 3 in Figure 59 where we deliberately challenged the inference by setting the ground truth value for to far in the left tail of the prior). Moreover, SNPE-C is often unable to infer the measurement error variability (datasets 2, 4, 5). As previously found in Pieschner et al. (2022), the lack of identifiability of parameters (k, mo, scale) is confirmed our analyses. However the other parameters are correctly inferred by SeMPLE. Finally, we focus on displaying the performance of SeMPLE for each round.
Figures 57 61
identify to accurately, while SNT.
Figures 62 66
5
3
log
log
log Mo
SNPE-C
SeMPLE
3
-2
-1
0
3
4
-3
22
-1
log scale
log to
log offset
log D
Figure 57: Translation kinetics model, observed data set. #1 out of 5: marginal posteriors from SeMPLE, SNPE-C and SNL (after 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines). For SNPE-C and SNL show results after 10 algorithm rounds, and after 3 rounds for SeMPLE.
64

---

Published in Transactions on Machine Learning Research (06/2024)
-2
-1
5
4
-3
-2
-1
11
2
3
4 5
log
log Y
log
log mo
SNPE-C
SNL
SeMPLE
-2 -1
log
log to
offset
log
Figure 58: Translation kinetics model, observed data set #2 out of 5: marginal posteriors from SeMPLE, SNPE-C and SNL (after 30k model simulations), priors (black lines) and ground truth paramcters (vertical
red lines). For SNPE-C and SNL we show results after 10 algorithun rounds, and after rounds for SeMPLE.
-2
-1
5
4
-3
-2
-1
2
2
3
4
5
log
log Y
log
log mg
SNPE-C SNL
SeMPLE
-2
-1
log scale
log to
offset
log
Figure 59: Translation kinetics model, observed data set #3 out of 5: marginal posteriors from SeMPLE, SNPE-C and SNL (after 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines). For SNPE-C and SNL we show results after 10 algorithm rounds, and after rounds for SeMPLE.
65

---

Published in Transactions on Machine Learning Research (06/2024)
-2
11
5
4
-3
-2
-1
2
-1
2
3
4
5
log
log Y
log
log mo
SNPE-C
SNL
SeMPLE
log -2 -1
44
log to
offset
log
Figure 60: Translation kinetics model, observed data set #4 out of 5: marginal posteriors from SeMPLE, SNPE-C and SNL (after 30k model simulations), priors (black lines) and ground truth paramcters (vertical red lines). For SNPE-C and SNL we show results after 10 algorithun rounds, and after rounds for SeMPLE.
-2
-1
5
4
-3
-2
-1
2
2
3
4
5
log
log Y
log
log mg
SNPE-C SNL
SeMPLE
-2
-1
log scale
log to
offset
log
Figure 61: Translation kinetics model, observed data set #5 out of 5: marginal posteriors from SeMPLE, SNPE-C and SNL (after 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines). For SNPE-C and SNL we show results after 10 algorithm rounds, and after rounds for SeMPLE.
66

---

Published in Transactions on Machine Learning Research (06/2024)
-2
-1
-
4
-3
-2
-1
2
2
3
4
5
log
log Y
log
log mo
SeMPLE round
SeMPLE rounds
SeMPLE round
log -2
-1
4
log to
offset
log
Figure 62: Translation kinetics model, observed data set #1 out of 5: marginal posteriors from SeMPLE for every algorithm round (corrcsponding to 10k, 20k and 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines).
2
-1
-5
4
-3
-2
-1
2
2
3
4
5
6
8
log
log Y
log
log mg
SeMPLE round
SeMPLE rounds
SeMPLE round
-2
-1
log scale
log to
offset
log
Figure 63: Translation kinetics model, observed data set #2 out of 5: marginal posteriors from SeMPLE for every algorithm round (corresponding to 10k, 20k and 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines).
67

---

Published in Transactions on Machine Learning Research (06/2024)
-2
-1
0
-5
4
-3
-2
-1
2
3
4
5
6
7
log
log Y
log
log mo
SeMPLE round
SeMPLE rounds
SeMPLE round
-2
-1
4
log
log to
offset
log
Figure 64: Translation kinetics model, observed data set #3 out of 5: marginal posteriors from SeMPLE for every algorithm round (corrcsponding to 10k, 20k and 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines).
-2
-1
-5
4
-3
-2
-1
2
2
3
4
5
8
log
log Y
log
log mg
SeMPLE round
SeMPLE rounds
SeMPLE round
-2
-1
log scale
log to
offset
log
Figure 65: Translation kinetics model, observed data set #4 out of 5: marginal posteriors from SeMPLE
for every algorithm round (corresponding to 10k, 20k and 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines).
68

---

Published in Transactions on Machine Learning Research (06/2024)
-2
-1
-
4
-3
-2
-1
3
5
log
log Y
log
log mo
SeMPLE round
SeMPLE rounds
SeMPLE round
2 -1
-3
log
log to
offset
log
Figure 66: Translation kinetics model, observed data set #5 out of 5: marginal posteriors from SeMPLE for every algorithm round (corresponding to 10k, 20k and 30k model simulations), priors (black lines) and ground truth parameters (vertical red lines).
69


Reviews:
[{"review_id": "INcflII2As", "reviewer": "Reviewer_1VPa", "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": "This is a good paper, with a single claim and sufficient evidence to support it. Reducing computational cost of amortized inference is definitely of interest to practitioners in the TMLR audience. Both the method description and the experimental sections are clear and easy to read. To the extent I'm familiar with the relevant literature, the paper described it accurately.\n\nI don't see any major weaknesses. The paper could be further improved by tuning the baselines more carefully and benchmarking against alternative classes of methods, but I don't really think any of that is necessary.", "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 proposes to use a linear Gaussian mixture model instead of distributions parameterized by neural networks for amortizing simulation based inference. The effort is motivated by reducing the associated computational costs, and as a result the energy expenditure. The paper provides a complete description of the training algorithm and a set of experiments on popular benchmarks proving that the proposed method achieves good performance, while being significantly cheaper than state of the art neural network-based alternatives.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I'm happy with the paper as is.", "broader_impact_concerns": "No concerns. The only foreseeable broader impact is reducing energy consumption, which seems unequivocally good.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cmrpS9URwu", "reviewer": "Reviewer_iTZL", "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- The presentation is clear, both in terms of language and structure.\n\n**Weaknesses:**\n- I struggle to identify the main contributions of the work given that the mixture model has been described in detail before, and has even been adapted for ABC previously. Further, if the main outcome is to have a good proposal for the MCMC method, that could be used directly as a learning objective, see e.g. Wang et al. (2018), _\"Meta-Learning MCMC Proposals\"_ and Stites et al. (2021), _\"Learning Proposals for Probabilistic Programs with Inference Combinators\"_.\n- GLLiM assumes a fixed, and relatively small, number of components and I'm missing a discussion of the pros and cons of alternative but related approaches such as (i) using a continuous mixture distribution, e.g. Deep Kernel Processes (Aitchison et al., 2021) or Elliptical Processes (Bånkestad et al., 2023), or (ii) an amortized model that uses e.g. neural networks to directly output the parameters (akin to a VAE).\n- I question the premise that the cost of inference is a major concern in SBI, I would expect that in most cases this cost is negligible compared to running the actual simulation. Please provide real-world counter-examples if you want me to reevaluate this view.\n- Experiments are only performed on toyish datasets. Moreover, I'm sceptical that the baselines, which are both from 2019, are still accurate reflections of the state of the art. \n\n\n**Comments and suggestions:**\n\nHow do you measure energy consumption? I imagine it must be dependent on the hardware.\n\nA graphical illustration of GLLiM would be helpful.\n\n**Nitpicking:**\n- Fig 1: a very close look is needed to discern the differences. Consider overlaying the true posterior or highlighting the differences in some other way.\n- Fig 2: C2ST is not defined until much later. Consider showing quartiles instead of min-max.\n- Table 1: CPU+DRAM heading should be energy.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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 authors propose a method for simulation-based (a.k.a likelihood-free) inference that uses a parametric mixture model as an approximate likelihood/posterior, that they then use as the proposal in a Metropolis-Hastings sampler. The mixture model, Gaussian locally linear mapping (GLLiM), is the same as in Deleforge et al. (2014) and was adapted for Approximate Bayesian Computation (ABC) in Forbes et al. (2022). My interpretation is that the authors make two principal claims:\n1. _\"our methodology returns accurate inference, even in presence of multimodal posteriors, with a small runtime.\"_\n2. _\"When comparing with widely employed inference engines based on neural-networks, our method considerably lower the energy and memory required for the\ninference task.\"_\n\nElaborating on 1, they state that _\"for a given computational budget specified by the total number of model simulations, our method is cheaper to run and produces accurate inference thanks to (i) the informed proposal sampler created by GLLiM and (ii) the computationally faster model training.\"_", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "See above.", "broader_impact_concerns": "No.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "8ZW3Y9PyOJ", "reviewer": "Reviewer_YJPN", "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:  \nS1. The paper is well-written.  \nS2. I am not an expert in SBI, but it seems to me that the evaluation methodology is relatively comprehensive (follow an established/published SBI benchmarking methodology - uses appropriate problem settings and metrics), well-motivated, and fairly narrated.  \n\nWeaknesses:  \nW1. The novelty is limited but in my opinion appropriately narrated.", "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 authors propose a method for likelihood-free/simulation-based inference (SBI), similar to Sequential Neural Likelihood (SNL) [Papamakarios et al., 2019], but utilizing a (far more) tractable approximate probabilistic model $\\hat{p}(\\theta, y)$ that allows for analytic determination of both the approximate likelihood $\\hat{p}(y | \\theta)$ as well as the approximate posterior $\\hat{p}(\\theta | y)$ under the model $\\hat{p}(\\theta, y)$.  \nThe advantage of this change is that the resulting method can be considerably computationally cheaper & faster than the neural-density-estimation based alternatives (on the included experiments, between 1-20x), while achieving comparable quality of inference across a range of benchmark tasks in these problem settings.  \nFor me it's particularly interesting that such a simple probabilistic model allows to achieve a competitive performance on those benchmark tasks.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "R1. In the abstract, \"they do not generally achieve an optimal trade-off between accuracy and computational demand\" - I think it's not clear what \"optimal\" means here. Could you clarify or rephrase, please?  \nR2. I'd love to see how does the quality of the approximation/inference improve with the number of rounds R for both SeMPLE and other methods. (although not crucial for acceptance). \n\nNitpick:  \nR3. maybe just “parameters” instead of “hyperparameters”?  \nR4. Algorithm 1: lines 3 - use bold \\theta for consistency and maybe rewrite this \"if\" statement with indentation such that it's easier to read?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper meets TMLR's evaluation criteria and is supported by all reviewers.
Decision: Accept (poster)