Preconditioned Neural Posterior Estimation for Likelihood-free Inference
Authors: Xiaoyu Wang
Ryan P. Kelly
David J Warne
Christopher Drovandi
Venue: Accepted by TMLR ()

Abstract: Simulation-based inference (SBI) methods enable the estimation of posterior distributions when the likelihood function is intractable, but where model simulation is feasible.  Popular neural approaches to SBI are neural posterior estimation (NPE) and its sequential version (SNPE).  These methods can outperform statistical SBI approaches such as approximate Bayesian computation (ABC), particularly for relatively small numbers of model simulations.  However, we show in this paper that the NPE methods are not guaranteed to be highly accurate, even on problems with low dimension.  In such settings the posterior cannot be accurately trained over the prior predictive space, and even the sequential extension remains sub-optimal. To overcome this, we propose preconditioned NPE (PNPE) and its sequential version (PSNPE), which uses a short run of ABC to effectively eliminate regions of parameter space that produce large discrepancies between simulations and data and allow the posterior emulator to be more accurately trained.  We present comprehensive empirical evidence that this melding of neural and statistical SBI methods improves performance over a range of examples including a motivating example involving a complex agent-based model applied to real tumour growth data.

Published in Transactions on Machine Learning Research (09/2024)
Preconditioned Neural Posterior Estimation for Likelihood- free Inference
Xiaoyu Wang
x811. wal.gaan@@hr aut.eddu.ou
School of Mathematical Sciences Centre for Data Science
Queensland University of Technology
Ryan P. Kelly
r21. kelly@lhhhrr qut. edu
School of Mathematical Sciences
Centre for Data Science
Queensland Universily of Technology
David J. Warne
david. warddaarneeeeuuee edu..n
School of Mathematical Sciences
Centre for Data Science
Queensland University of Technology
Christopher Drovandi
c.drovandi@qut edu
School of Mathematical Sciences Centre for Data Science
Queensland University of Technology
Reviewed on OpenReview: https: // openreview. net/ net/ forum? id- 0gIAOKKAA
Abstract
Simulation-based inference (SBI) methods enable the estimation of posterior distributions when the likelihood function intractable, but where model simulation is feasible. Popular neural approaches SBI are neural postcrior estimation (NPE) and sequential version (SNPE). These methods can outperform statistical SBI approaches such as approximate Bayesian computation (ABC), particularly for relatively small numbers of model simulations. However, show in this paper that the NPE methods accurate, even problems with low dimension. In such settings the posterior cannot be accurately trained over the prior predictive space, and even the sequential extension remains sub-optimal. To overcome this, we propose preconditioned NPE (PNPE) and its sequential version (PSNPE), which uses a short run ABC to effectively eliminate regions of parameter space that produce large discrepancies between simulations and data and allow the posterior emulator to be accurately trained. We present comprehensive empirical evidence that this melding of neural and statistical SBI methods improves performance over range of examples including motivating example involving complex agent-based model applied to real tumour growth data.
not guaranteed to be highly
1 Introduction
Computational models, frequently termed as simulators, are typically governed by stochastic processes. When provided with set of parameter values, these simulators output synthetic data that inherently capture the stochastic nature the simulated phenomena. However, a substantial challenge arises in performing posterior inference for the parameters of these simulators, as the corresponding likelihood function often intractable. An example is agent- modelling tumour growth (c.g. Jenner al. (2020): Aylett-
Bullock et al. 2021 Warne et al. (2022 where cell proliferation, movement and invasion governed by

---

Published in Transactions on Machine Learning Research (09/2024)
probabilistic rules, which depend on key biological parameters. Consequently, standard statistical inference methods that rely on closed-form expression for the likelihood function are inapplicable in this scenario.
To address this issue, simulation-based inference (SBI) methods have been developed, which approximate the posterior based only on simulations from the model. The most thoroughly examined SBI method in
the statistical literature is approximate Bayesian computation (ABC) (Sisson et al.] 2018). Advancements
in deep neural networks have led to the emergence of neural SBI methods (Cranmer et al. 2020). Their
widespread application spans various fields, including biology (Wehenkel et al. |002), neuroscience |Fenller
et al. 2021 West et al. 2021 and astronomy (Mishra- -haarra----aaaaa
2022 Dax et al. 2021
Statistical SBI methods, such as ABC methods, are well-developed and boast strong theoretical guarantees
of convergence to the true posterior (Beaumont et. al.) 2009; Blum 2010 Biau et al., 2015 Lintusaari el al. 2017 Sisson et al., 2018 Beaumont 2019). ABC approaches compare observed to simulated data using a discrepancy function and prefer parameter values that generate discrepancics below a pre-defined threshold, 6. However, often a small value of is required to obtain accurate posterior approximations, which can significantly increase the number of model simulations required, and hence the computational cost (Csilléry
et al.
2010 Turner& & Van Zandt, 2012). Advanced ABC samplers such as adaptive sequential Monte Carlo (SM ABC, c.g. [Sisson et al. (2007); Drovandi & Pettitt (2011)) have been developed to mitigate this However, the efficiency of ABC methods may decrease as e becomes smaller, requiring a significantly larger number of simulated datasets. Euhancing the efficiency of ABC algorithms remains an active area of research (Lintusaari et al. 2017 Sisson et al. 2018
ABC
issue.
With the rapid advances in machine learning methods, more efficient approaches based on neural networks have been developed (Papamakarios ct al.) 2017 2021 A popular approach is neural posterior estimation (NPE) and its sequential version, the sequential NPE (SNPE) (Papamakarios & Murray 2016 Lueckmann
et al. 2017 Greenberg et al.) (2019). These methods use a set of training pairs of parameter values and simulated datasets to fit a neural conditional density estimator (NCDE), such as a conditional normalizing
flow (Rezende& & Mohamed| 2015 Winkler et al., 2019 Durkan et al. 2019 Dolatabadi et al. 2020), to approximate the posterior. NPE uses parameter values drawn from the prior, whilst SNPE uses parameter values drawn from previous NPE approximations for a given number of rounds. The idea of the sequential approach is that more accurate emulator of the posterior can be achieved when more (parameter value, simulated data) training pairs are generated in higher density regions of the posterior. Lueckmann et al. (202I) have shown that NPE can outperform ABC, particularly for a relatively small number of model simulations.
However, in real-world problems, there may be little known about the parameters a priori, so that a vague prior may be employed. For example, a uniform distribution with wide constraint range may be used as prior distribution. For some problems in this setting we find that it is difficult to construct an accurate NCDE across wide parameter space, which leads to NPE producing inaccurate posterior approximation. We find that even SNPE may not be able to recover from such an initially deficient approximation, even with a relatively large number of rounds, and hence model simulations. One way to mitigate this kind of unstable NCDE training is to clip the extreme simulated datasets (Shih el al. 2023 de Santi et al. 2023). However, this approach is ad-hoc and it is not clear how much clipping is required for a given problem and may require extensive experimentation, and cach level of clipping requires refitting of the NCDE.
This paper contains three key contributions. Firstly, we explore several examples where NPE methods fail to produce highly accurate posterior approximations, even in relatively low dimensional problems. Our second contribution is the development of preconditioned NPE (PNPE), and its sequential extension (PSNPE), which combines the strengths of statistical and ML approaches to SBI. The preconditioning step involves applying an ABC algorithm for efficiently discarding parts of the parameter space that lead to large discrep- ancies, which then subsequently permits NPE methods to perform well. In a sense, our preconditioning step acts as a principled clipping method. Our third contribution shows via an extensive empirical study that our preconditioned NPF approaches outperform NPF approaches when the latter performs sub-optimally, and is competitive when it performs well. Our motivating example involves a complex agent-basod model applied to real tumour growth data.
2

---

Published in Transactions on Machine Learning Research (09/2024)
2 Simulation-based Inference
Consider a simulator that takes parameters 0 E Rd where d is the number of parameters and generates a simulated dataset 2 E R where is the dimension of the data, but its density p(z(0) intractable. The objective of SBI to accurately estimate the posterior density of conditional on the observed dataset Zo EDD based only on simulating data from the model and not requiring evaluation of the intractable likelihood, p(r.)0). Two popular SBI methods are ABC and NPE, which are summarized bclow.
2.1 Approximate Bayesian computation
Statistical SBI (Sisson al. 2018); such the ABC rejection algorithm, is based Monte Carlo rejection sampling. That is, it keeps only the parameter values simulated from the prior that generale simulated data such that ((tr,tr) < 6, where To) is a user-defined discrepancy function between the simulated and observed data, and is a user-defined threshold often referred to as the ABC tolerance.
SMC ABC algorithms (e.g. Sisson et al. (0077: Drovandi & Pettitt (2011) aim to be more efficient by sampling sequence of ABC posteriors with decreasing c's, updating the importance distribution at each
iteration. More specifically, SMC ABC algorithms define a sequence of non-increasing ABC thresholds c1 62 2..22 ET, such that
Pe, (0 x p(0) [ T(p(Fo; x) < (c)p(x)()dr,
for t aalt=,.... I....,T.
(1)
Here, ET represents the target ABC posterior.
In many real-world applications, o E R.D are considered high-dimensional data, necessitating a mapping to a lower-dimensional space for computational efficiency. This is typically done using summary statistics S(.). If summary statistics are required, we use S(x) and S(ra) instead of the full datasets x and dg-- The choice of appropriate summary statistics is a subject of ongoing research and is discussed extensively in the literature (sce Sisson al. (018)).
However, even sophisticated ABC algorithms can require a significant number of model simulations to achieve
suitably small value of e [Biuu et al. 2015 Csillery et al. 2010 Beaumont et al. 2009 Blum 2010
2.2 Neural posterior estimation
NPE uses N training pairs of simulator parameter values and simulated datasets, {0,, x, )
to estimate the
posterior distribution p(0 (x) (Papamakarios & Murray 2016 Lueckmann et al. 2017 Greenberg ot al. 2019)
Once the NPE is trained on the sinulated datasels, the posterior distribution 2()))))) can be computed by inputting the observed dataset Fe-
A conditional neural density estimator (F(z.6) (0), utilizing a neural network F and its adjustable network weights v, is often used NPE. In order to train (0), the following loss is minimized:
d arg min
Σ (og
(2)
over network weights vh. For sufficiently expressive dF, (F(..)) (0) p(0)z)) as N - xc.
SNPE aims to improve the accuracy of the approximate posterior for particular observed dataset xo iteratively by sampling parameter values from a previous NPE approximation for given number of rounds. The current NPE approximation is treated as a proposal distribution p(0) for the next round. However, training fr using parameter values drawn from p(0) will not converge to the true posterior distribution, but rather to
p(0) p(0x..) & p(0 (..)) p(0)
Many approaches have been developed overcome this limitation, such as Papamakarios Murray (2016) Lueckmann et al. (2017) Greenberg al. (2019). Among all these approaches, the automatic
3

---

Published in Transactions on Machine Learning Research (09/2024)
posterior transformation (APT, also known as SNPE-C), as proposed by Greenberg et al. (2019) which has been reported to significantly outperform the others (Tueckmann et al. 2021). For simplicity, we refer to our specific implementation as SNPE hereafter, but we note that other implementations of SNPE can be used with our method.
The general leakage problem can occur when the proposal distribution assigns non-zero probability density outside the support of the prior distribution (Deistler al.) 2022a). In other words, the loss function of SNPE-C (Greenberg et al.) 2019) does not force the NCDE to place density within the prior support, causing mismatch between the posterior and prior support. This can result in extreme values for the proposals, such as generating proposals outside prior bounds, or generating a relatively large number proposals near prior bounds when & suitable transformation is applied. Another possible 'leakage' might occur because the neural networks ignore some extreme or invalid data to stabilize the training, leading to unexplored areas of parameter space. This general leakage problem can result in poor performance of SNPE. This is an espocially scrious problem in situations where leakage is more likely to occur, such as when the training dataset itself contains extreme values of the summary statisties. The truncated SNPE (TSNPE) method (Deistler et al. 2022a) aims to overcome the leakage problem by using truncated proposal distribution and rejecting proposals if they fall outside of certain quantiles of the prior.
2.2.1 Illustrative Example
We consider a sparse vector autoregrossive (SVAR) model that has been considered previously in the SBI literalure (Thomas el al. 2020 Drovandi el al. 2023 The SVAR model is given by:
yt yt-- + Ec:
(4)
where y/ E R& represents the k-dimensional observation of the time series at time 4, X E Baxk is the transition matrix, and & ~ (N(00=1) is a k-dimensional noise vector with o² being the noise parameter. The model considers a sparse transition matrix X where the only off-diagonal entries that are non-zero must satisfy the following conditions: if variable i is coupled with variable j. then X1 # and X,, # 0 (note that Xij is not necessarily equal to Xz.)) and each variable is coupled to only one other variable. Under this condition, each column will only have one element that is non-zero, and this will also be an off dia-onal clement of the matrix. To ensure the stability of the SVAR, the diagonal elements of X are set to -0.1. The parameter space of SVAR can easily scale to higher dimensions by increasing k. In this study, the model parameters f E R are the non-zero off-diagonal entries of X and its variance and we consider = 6, which leads to 7 parameters. This choice is based on the assumption that if SNPE does not produce highly accurate approximations in this low-dimensional case, it is unlikely to be accurate in a higher-dimensional parameter space. We generate an observed dataset of length T = 1000 using the summary statistics to
true parameter value 0 (0.579, -..143, 0.836, 0.745, -06660. 00.254, 0.1). We reduce the dimension of the data. Following Thomas et al. (2020); Drovandi el al. (2023). we use the lag autocovariance + vivi as the summary statistics, where vi is the tth observation of the ith time series. We use the sample standard deviation of the k time series to inform σ. Thus there is single summary statistic that is intended to be informative about each parameter.
We employ a uniform distribution as the prior, constrained between -1 and 1 for the k parameters and between 0 and 1 for a. We denote the generated dataset as x E (R+(+) and its corresponding summary as 5(x) E (R)+)), We find that extreme values of the summary statistics can be produced by parameter values away from the true parameter value. To stabilize the training, we clip simulated datasets with summary statistic outliers (any simulated values greater than 10, around 3% of training datasets). This leads to some regions of the parameter space being unexplored and results in leakage after several rounds of SNPE training. Note that some experimentation was required to obtain clipping value that led to reasonable results for SNPE.
For illustrative purposes, we run three rounds of SNPE to avoid the leakage problem that occur several rounds after and compare the results with exact posterior samples as the likelihood is casily computable in this example since it is Gaussian. Thus in this case the reference distribution estimated from exact posterior samples. Ideally, we would expect performance to improve when increasing the number of SNPF
4

---

Published in Transactions on Machine Learning Research (09/2024)
rounds. However, as shown in Figures ! and even with datasets clipped for every round, SNPE does not improve the accuracy of the estimates as the number of rounds increases.
15
15
15
15
10
10
10
10
-1
0
-1
0
-1
0/
0g
0
15
15
400
Reference
10
10
NPE
200
SNPE
True
-1
-1
0
0.1
0.2
Figure 1: Comparison of marginal posterior distributions between reference distribution (orange), NPE (dashed pink) and SNPE (red), with grey solid lines representing the true values. The SNPE results are based on three rounds.
600
600
100
100
400
400
50
50
200
200
-0.02
0.02
-0.02
0.02
0.05
0.1
0.05
0.1
X1
X
X
X
600
600
80
Reference
60
400
400
NPE
SNPE
True
200
200
20
0
0
-0.02
0.02
-0.02
0.02
0.1
0.2
X5
X
Xy
Figure Comparison of posterior predictive distribution of the summary statistics of observation datascts between reference distribution (orange), NPE (dashed pink) and SNPE (red), with grey solid lines repre- senting the true values. The SNPE results are based three rounds.

---

Published in Transactions on Machine Learning Research (09/2024)
3 Method
We find that NPE can perform sub-optimally when the prior predictive distribution of the data is complex and has significant variability. When this occurs, the NCDE may not be sufficiently accurate, especially in regions of high posterior support. The sequential version of NPE was originally designed to overcome this issue, however if the initial NPE approximation is not substantially better than the prior, then subsequent rounds of SNPE may suffer from the same issue.
Other approaches to overcome this issue may be increase the training sample size or to try different configurations of the neural network, but both of these may increase the computational cost substantially and may not address the issue. Instead, we propose the preconditioned NPE (PNPE) method, and its sequential extension below, in order make NPE methods more reliable.
3.1 Preconditioned NPE
For a vague prior distribution p(0), parameter values drawn from it might be very far from the true posterior. We suggest using a short run of ABC to refine those parameter values so that they are closer to the true posterior p(0 ro), i..,
pe(0)a)) x p(0) / (( < e)de,
(5)
where I[-) is the indicator function and e can be chosen considerably larger than what might typically be used in an ABC algorithm. Then we fit a density estimator to those parameter samples. The key idea is to use an efficient ABC algorithm to quickly discard poor regions of the parameter space that generate unusual datasets relative to the observed data, which provides belter quality training datasets for NPE. We note that any ABC algorithm could be employed here, but we use the SMC ABC algorithm of Drovandi & Pettitt (2011) in this paper (see Appendix for a full description of this method). The SMC ABC algorithm generates samples from sequence of ABC posteriors based on decreasing ABC thresholds, e1 where 6T is the target ABC threshold. The sequence of tolerances is determined adaptively, by, at each iteration of SMC, discarding a proportion of the samples, a n, with the highest discrepancy, where a is a tuning parameter. Then, the population of samples is rejuvenated through a resampling and move step. During the move step, a Markov chain Monte Carlo (MCMC) ABC kernel is employed to maintain the distribution of particles based on the current value of the tolerance. The number of MCMC steps R to apply to each particle is determined adaptively based on the overall MCMC acceptance rate, that is R, log(e) 1og(1-2t00 where peee is the estimated MCMC acceptance probability at the SMC itcration t and is a tuning parameter of the algorithm that can be interpreted as the probability that a particle not moved in the R2 MCMC itcrations. A natural stopping rule for the algorithm when the MCMC acceptance rate becomes intolerably small.
cr:
Based on parameter samples from the ABC posterior, fit unconditional normalizing flow 4G (note that other density estimators could be used). Then we can use qa the initial importance distribution for the NS) process. We call this melding of ABC and (S)NPE as the preconditioned (S)NPE method. The method is summarized in Algorithm
Algorithm Precondition SNPE
1: Choose preconditioning ABC algorithm and SNPE implementation
2: Obtain [0:)_1 from the preconditioning ABC algorithm
Sct o* <- arg min Σ -logga(0))
Perform SNPE using initial importance distribution (0)
If we obtain a well-trained unconditional normalizing flow, this unconditional normalizing flow can act as the initial importance distribution, i.e., p(0) in Equation 3 and to draw samples for training the NPE. Following Papamakarios & Murray (2016), it is noted that. given expressive enough conditional normalizing flow,
6

---

Published in Transactions on Machine Learning Research (09/2024)
NPE converges to the true posterior p(0)zo) as N - co, with an appropriate importance re-weight if p(0) # p(0).
Choosing suitable value of e for our method requires some thought. A smaller value of will focus in on more promising regions of the parameter space, but will increase the computational time of the preconditioning step. A larger value of will lead to a fast preconditioning step, but may not. eliminate enough of the poor parts of the space to improve the training of the NCDE. In this paper we use an MCMC acceptance rate of 10% (unless otherwise specified) as the stopping criteria for the SMC ABC algorithm in the preconditioning step. For our examples find that this choice is effective at balancing the aforementioned objectives. We note that other choices are possible.
Furthermore, once these poor sinulations have been removed, we find NPE to be more effective than ABC, since ABC requires an exponentially increasing number of simulations to drive e to O. To avoid the scaling problem of ABC, the preconditioning step only takes a short run of ABC, and thus we are not interested in driving to O.
3.2
Computational cost
We now consider computational cost for P(S)NPE and compare it with SNPE. The preconditioning step can be considered as the initial round of NPE where the total number of simulated datasets generated during SMC ABC is denoted as NABC- Hence, it is worth noting that P(S)NPE, like SNPE, is not amortized since it requires running an ABC algorithm for each observation datasets zo.
Furthermore, for complex real-world problems, the simulation time may depend the parameter values, and parameter values with very low posterior support can produce substantially longer simulation times. For such problems, is important from a computational perspective to quickly eliminate such regions from the parameter space, as is the motivation the preconditioning ABC algorithm. Thus, for problems where SNPE does not perform well, we find PSNPE to substantially more computationally efficient in terms of compute time.
We perform analysis of the trade-off between computational cost and estimation accuracy by using differ- eut choices for the % of stopping rule for the preconditioning step in the SVAR example. For reprodueibility, we use 20 different random seed values for each choice of the stopping rule. We compute and record the average values we choose for the % and its corresponding total number of simulations in Table
MCMC ACCEPTANCE SIMULATIONS IN ABC
15%
12%
10%
8%
5%
15603
33585
62402
100398
244667
Table 1: Computational cost for preconditioning step: The first row of the table refers to the % stopping rule selected to obtain samples from ABC, while the second row indicates the average of total number of simulations ABC used for the corresponding %.
To investigate the performance of this trade-off, we use maximum mean discrepancy (MMD) as metric to compare the approximate distribution from ABC with the reference distribution. We also run NPE with the same number of simulations as ABC and compute its MMD. In Figure 3 we plot the total number of simulations versus MMD for ABC, NPE and PNPE. The solid lines represent the average MMD values basod on 20 MMD values from 20 different scod values. The number of simulations for PNPE includes an additional 10k simulations generated after the corresponding preconditioning step. It is surprising that with increasing number of simulations, the accuracy of NPE can vary significantly. For some random seed values, the accuracy of NPE even decreases when using more than approximately 70,000 model simulations. Compared with ABC and PNPE, which have narrow uncertainty intervals, the accuracy of NPE seems to highly depend on the random seed value.
We record that 3- 5% of training datascts are elipped before training NPE, which mcans around 2, 2,100-3,500 simulated datasets were ignored, and the corresponding parameter space is not explored well. This might indicate the reason why the accuracy of NPE might reduce under a large number of simulations for some
7

---

Published in Transactions on Machine Learning Research (09/2023)
0.6
NPE
ABC
PNPE
0.5
MMD
0.4
0.3
0.2
50000
100000
150000
200000
250000
Total Number of Simulation
Figure Computational cost versus MMD for ABC, NPE and PNPE for SVAR example: Comparison of MMD values for ABC (blue), NPE (red) and PNPE (green) with corresponding color bands indicating the 95% uncertainty interval. The x-axis shows increasing number of simulations required for each algorithm.
random seed values. It is evident that an increasing number of simulations can cause a decrease in MMD values for both ABC and PNPF. However, when a large number of simulations are used in the preconditioning step, the following NPE step does not provide a large improvement in accuracy. The small improvement in accuracy may not be worth the cost incurred by the large number of model simulations. Having a percentage too large will stop the ABC too early, and we will lose the benefit of the preconditioned step. Having percentage too small will result in ABC running too long and using too many model simulations.
3.3 Illustrative Example Revisited
We apply PNPE to the illustrative example shown in Section 2.2.1 For the preconditioning step, we use the adaptive SMC ABC algorithm proposed by Drovandi & Pettitt 2011). with tuning parameters n = 1k, 0.5, and e 0.01. We employ an unconditional normalizing flow as the unconditional density estimator. For this, use the state-of-the-art neural spline flow implemented in the Pyro package.
In order to make fair comparison between P(S)NPE and SNPE, as in the SMC ABC algorithm, denoted NABC, to train the initial round of NPE. We run the SMC ABC algorithm ten times to obtain the average number of simulations it requires, which is NABC = 54k. For illustrative purposes, we only run two rounds of SNPE and compare it with PNPE.
the same mumber of simulations
The estimated marginal posterior plots are displayed Figure 4 where the black solid lines represent the true parameter values. It is evident that our PNPE method (in single round) produces a substantially sharper approximation of the posterior compared to that of SNPE.
8

---

Published in Transactions on Machine Learning Research (09/2024)
15
15
15
15
10
10
10
10
5
5
5
0
0
0
-1
-1
0
-1
0
0
0
0
0
2
4
15
15
400
Reference
300
SNPE
10
10
-reecondiiooingg
200
PNPE
5
5
True
100
0
0
0
-1
1
-1
0.05
0.1
0.15
00
Figure 4: Performance on SVAR model. Comparison of marginal posterior distributions between refer- ence distribution (orange), SNPE (red), proconditioning step (blue dash) and PNPE (green solid), with grey solid lines representing the true values.
It is evident that short of the ABC algorithm gives reasonable posterior approximation to train unconditional normalizing flow as the unconditional density estimator, which then generates samples for NPE training. Figures g and show that the improved parameter posteriors leads to more accurate postcrior predictive distributions the sunmaries compared to SNPE. This indicates that with good starting point. NPE can further improve accuracy.
In Appendix C.1 we show the results for different datasets, each with different random soeds for reproducibility purposes. We use the stopping rule for the preconditioning step across all datasets. It is evident that the performance between ABC and NPE close under the same number of simulations. In some cases, NPE performs better than ABC it has sufficient data to learn while not too much parameter space is ignored. However, PNPE performs better than SNPE as the leakage problem occurs. We also compare NPE and PNPE with the same number of simulations and find that the performance of PNPE is still better than NPE.
4 Further Experiments
We present five benchmarking tasks (two from Lueckmann et al. (2021) and three from the SBI literature) and two additional examples, including our motivating example, where SNPE, perhaps surprisingly, does not produce highly accurate posterior distributions. To fairly compare our method with vanilla SNPE and potentially TSNPE (if the leakage issue is encountered), we run the ABC algorithm 10 times and compute the average total number of simulations that the ABC algorithm requires. We then this number of

---

Published in Transactions on Machine Learning Research (09/2024)
800
800
80
80
600
600
60
400
400
40
40
200
200
20
20
0
0
0
-0.01
0
0.01
-0.01
0
0.01
-0.05
0
0.05
-0.05
0
0.05
X
X
X
X
800
800
80
Reference
600
600
60
SNPE
PPreconditioning
400
400
40
PNPE
True
200
200
20
0
0
-0.01
0.01
-0.01
0.01
0.1
0.15 0.2
X
X
X
Figure 5: Performance SVAR model. Comparison of posterior predictive distributions of the summary statistics of observation datasets between reference distribution (orange), SNPE (red), preconditioning step (blue dash) and PNPE (green solid), with grey solid lines representing the true values.
simulations for the initial NPE in both SNPE and TSNPE. For all experiments, utilize the adaptive SMC ABC algorithm proposed by Drovandi & Pettitt (001 for the preconditioning part with tuning parameters = 1k, a = 0.3, 0.01, and use an unconditional normalizing flows as the unconditional density estimator. For SNPE, we use 10k samples for each round of training.
We find that even with single round PNPE, there can be a significant improvement performance. The code reproduce the results has been included supplementary material.
4.1 Benchmarking example
We compare PSNPE and SNPE across five popular benchmarking simulators (see Appendix B.1 for a detailed description). We use the sbibm package [Lucckmann al. (021] for the two-moon and SL.CP simulators, and the ELFI package (Lintusaari et al., 2018) for the other benchmark simulators. Since the ground-truth posteriors are available, sample test (C2ST), where score of 0.5 indicates that the approximate posterior is indistinguishable from the true posterior, and a score of significs that the classifier can completely separate the approximate posterior from the true posterior. Additionally, we use the maximum mean discrepancy (MMD) between ground truth posterior and approximate posterior to quantify performance. We refrain from using the negative log probability of true parameters as performance metric because our method is not amortized. It evident from Table 2 that SNPE and PSNPE achieve similar results, except for the g-and-k example, the results of which describe in more detail next.
use specific performance metrics for comparison. We use classifter two-
10

---

Published in Transactions on Machine Learning Research (09/2024)
MODEL
C2ST
MMD
Two MOON
0.527 (0.528
0.00067,0.000
SLCP
0.664 0.664)0.691
0.00065.00000
MA(2)
0.887 \0.856
0.104 /0.0966
G-AND-K
0.975 (0.727
0.553 00.311
TOAD EXAMPLE
0.899 00.887
0.120 10.111
Table 2: Performance on popular benchmarking simulators. Classification accuracy (CS2T) and maximum mean discrepancy (MMD) computed for SNPE) with same number simulation for ABC and NPE followed by 2 rounds SNPE, each round with N 10k simulations. The bold indicates better performance metric value.
For the benchmark examples in Table 2 where the performance of SNPE and PSNPE similar, the initial training dataset is not extreme and the NCDE fits reasonably well. In more realistic examples, the prior predictive distribution of the summaries may be much more complex and contain extreme values. In simple examples, we can create more complex prior predictive distributions by widening the prior. Here we reconsider the g-and-k benchmark example but widen the prior distribution.
We describe the details of the g-and-k model and its prior distribution Appendix We use full datasets instead summary statistics for both PNPE and SNPE. For visualization purposes, plot the marginal approximate posterior distribution for cach method of the estimated results. It is evident that ABC performs better than NPE under the same number of sinulations (Figure 6a), where the ABC posterior is able to concentrate on the true value of B and k. As shown in Figure 6b PNPE performs better than SNPE.
(a)
A
g
k
0.15
0.6
0.10
1.0
NPE
Density
0.10
Density
0.4
Density
Density
ABC
True
0.05
0.5
0.05
0.2
0.00
0.0
0.00
0.0
0
10
10
0
10
10
(b)
A
g
k
0.6
0.15
1.5
SNPE
Density
0.4
Density
0.4
Density
0.10
Density
1.0
True
PNPE
0.2
0.2
0.05
0.5
0.0
0.0
0.00
0.0
10
0
10
10
10
Figure 6: Performance on G-and-K example: (a) Comparison of marginal posterior distributions be- tween the ABC preconditioning step (red) and NPE (blue), with black dashed lines representing the true values; (b) Comparison of marginal posterior distributions between PNPE (green) and SNPE (red), with black dashed lines representing the true values. The result of SNPE, uses rounds with same number of simulation as the ABC proconditioning step.
11

---

Published in Transactions on Machine Learning Research (09/2024)
4.2 High-dimensional SVAR model
To investigate how our method scales to higher dimensional problems, take the illustrative SVAR example from before and consider A: 20, which leads to 21 parameters. We detail the experimental settings in Section B.2 To ensure a fair comparison, we run the SMC ABC algorithm ten times using a 10% acceptance rate as the stopping critcrion and calculate the average number simulations it takes. We then use the same number of simulations, approximately RABC 45k, to train the initial NPE. Hence, the total number of simulations for both PNPE and SNPE the same (55k total). To stabilize the training, we apply the same clipping technique used in the previous low-dimensional case, which results in approximately 11% of the training samples being removed in the initial round of NPE and about 1% to 2% in the second rounds. Starting from the third round, SNPE is unable to sample any parameter values from the neural networks due to severe Icakage issue (Deistler et al. 2022a).
The estimated marginal posterior plots are displayed in Figure 7 where the black solid lines represent the true parameter values. is evident that as the number of parameter dimensions increases, training the unconditional normalizing flows becomes more challenging. With well-trained unconditional normalizing flows, PNPE outperforms SNPE in high-dimensional cases (in this example, except for parameter 03.).
10
10
5
5
5
0
0
005
0.5
-0.4-0.2
0.4
0
2
9
0
05
0
10
10
5
5
0
0
0.5
0.5
-0.5
0
-0.4
0
0.4
0
0
0.5
1
07
08
0g
0
1
11
0
12
10
10
5
5
5
-1
-0.5
-1
-0.5
-1
-0.5
0
-1
-0.5
-1
-0.5
0
-0.4
0,4
015
0,6
0,7
0
13
10
10
500
Reference
SNPE
5
0
PNPE
True
0
0
0.5
1
0.09 0.1 0.11
0
19
20
0
Figure Performance on SVAR model with 21 parameters. Comparison of marginal postcrior distributions between reference distribution (orange), SNPE (red) and PNPE (green), with grey solid lines representing the true values. The result of SNPE uses 2 rounds.
12

---

Published in Transactions on Machine Learning Research (09/2024)
400
400
400
400
500
500
200
200
200
0
0
0
-0.02
-0.02
0
0
0.02
-0.02
-0.02
0
0.02
X
X
X2
X4
X
X
400
400
400
500
500
500
200
200
200
0
0.02
0.02
-0.02
0.02
-0.02
0
0
0.02
x
X
Xg
X10
111
X
500
500
200
200
100
100
100
100
0
0
0
-0. 02
-0. 02
-0,02
0
-0.02
0
-0.02
-0.02
X
X 14
X15
X16
X
X 8
Reference SNPE
50
50
100
50
PNPE
True
0.02
0.04
0.02
0.04
0.12
0.14
X19
X20
X21
Figure 8: Performance SVAR model with 21 paramcters. Comparison of posterior predictive distributions of the summary statistics of observation dalasets between reference distribution (orange), SNPE (red) and PNPE (green), with grey solid lines representing the true values. The result of SNPE uses 2 rounds.
4.3 Biphasic Voronoi cell-based model
Finally, we consider challenging real-world problem in cancer biology: calibrating the biphasic Voronoi cell-based model (BVCB) (Wang al. 2024) that models tumor growth. The model uses a parameter to divide the tumor into two growth phases. Here, the term 'growth phase' refers to the different growth patterns of the tumor. There are four parameters that govern tumor growth during each phase, namely Ppsc: dmax ), where and Pase are the probability of cell proliferation and invasion, respectively, dmax is the maximum distance between cell and nutrient, and Jage is the time taken for cell to be able to divide. Thus there are nine parameters total, four parameters each of two phases, and the parameter I which the growth phase changes. In this paper, we calibrate to two real-world panercatic cancer datasets Wade (2019), which describe tumor growth time series data. The datasets span 26 and 32 days, respectively, with measurements taken each day. While the ground truth posteriors unknown for those datasets, we compute posterior predictive distributions to assess if SNPE and PNPE can effectively calibrate the model to the data.
We employ vague prior distributions for all parameters. Specifically, we use a Uniform distribution con- strained between and 24 hours the number of days for Sage during both growth phases. Additionally, we use Uniform distribution constrained between and the number of days minus for The prior distributions for the remaining parameters detailed in Appendix B.3
13

---

Published in Transactions on Machine Learning Research (09/2024)
As reported by Waan al. (0024) CPU times for model simulation range from 1.76 to 137.27 seconds per simulation when using samples from the prior distribution. This implies that 10k simulations for the first round of SNPE take approximately 2 hours. Consequently, the initial stages of SNPE are computationally expensive. In contrast, the ABC preconditioning step takes around 10-15 minutes. This is duc to the fact that the longer simulation times lend to also lead to large discrepancies with the observed data, and such samples are quickly rejected by ABC. For the ABC part, around 18k and 16k simulations are used for the 26-day and 32-day pancreatic cancer datasets, respectively.
We observed a leakage problem occurring in rounds 8 and 5 for the 26-day and 32-day datasets, respectively. Hence, we perform 10 rounds of TSNPE for all three datasets. At round 10, the acceptance rate of SNPE for the 26-day measurement dataset is above 50%, leading us to utilize rejection sampling for sample generation. While a 50% acceptance rate for rejection sampling is acceptable, it is less efficient compared to direct simulation from the trained conditional normalizing flows in the previous round. Hence, we employ TSNPE for this dataset. For the 32-day measurement datasct, only 0.000% of samples are accepted at round 8, making it computationally expensive. Consequently, we use TSNPE for this dataset.
To estimate the posterior predictive distributions, we sample 1k parameter values from (T)SNPE and PNPE, using them to simulate datasets. We then plot these data in the form of credible intervals. As a bascline, we show the prior predictive distributions in Appendix C.2 The top and bottom rows of Figure g displays the posterior predictive distribution for the 26-day and 32-day datasets, respectively (the same plots but on the log scale are shown in Appendix C.2 It is evident that both SNPE and TSNPE provide biased estimations, as the posterior predictive does not capture the observed data well. The posterior predictive distribution for the preconditioning step (middle column) shows that the ABC step can capture the data reasonably well since the observed data lie within the 90% posterior predictive interval. Our method (third column) provides a better fit, as the variance of the posterior predictive distribution is tighter than that of the preconditioning step and still captures the observed data. This demonstrates that even one round of PNPE can perform more accurate cstimations based on our results.
5 Discussion
We present a neural SBI method that is both simple and easy to deploy, designed to enhance the accuracy of SNPE methods. Our method, termed preconditioned neural posterior estimation (PNPE) and its sequential version, PSNPE, employs an ABC algorithm for the initial step. This algorithm is used to efficiently filter out poor regions of the parameter space. Additionally, unconditional densily estimator qa, enabling 4G to serve as the initial proposal distribution for SNPE. The core concept is that an improved starting point can significantly enhance the accuracy of SNPE estimations. Indeed, we obtained very good results with PNPE. Lueckmann et al.] [017]]; Deistler et al. (2022b) have proposed similar ideas to improve accuracy of SNPE. They train a classilier to predict samples that fall inlo rejection criterion based on a certain distance metric. However, such a classifier can falsely reject samples. Our method can avoid this problem and hence is more principled compared to their method.
use the ABC posterior samples to train
We showcase several examples where either SNPE failed to perform inference effectively, such in the SVAR case, or produced biased results, as observed in the BVCBM. For the SVAR example, SNPE mot struggle due to the impact of low-quality samples from certain parameter space regions, adversely affecting the training process. The ABC method can efficiently climinate these bad samples, thereby enhancing the training. For cases where SNPE results in biased estimations, our methods were effective at accurately fitting observed data (real data for BVCBM example). This is substantiated by our empirical results for the posterior predictive distribution discussed in the experimental section.
Although method demonstrates the capability to enhance estimation accuracy, it does have some limi- tations. Firstly, our method requires model sinulations in the ABC preconditioning step, which may lead to greater computational demands in situations where SNPE methods perform well. However, by perform- ing the preconditioning step, significantly fewer model simulations may be required in the SNPE part to achieve high accuracy. In this paper we usod SMC ABC for the proconditioning step, but we note that other ABC algorithms or SBI methods could be used. We do not recycle the simulations performed in the ABC preconditioning step for the SNPF phase, but it could be possible to modify our method to exploit
14

---

Published in Transactions on Machine Learning Research (09/2024)
5000 (a)
5000 (b)
4500
(c)
4000
Tumour size (mm3)
4000
Tumour size (mm
4000
Tumour size (mm3)
3500
3000
3000
3000
2500
2000
2000
2000
1500
1000
1000
1000
500
0
0
5
Time (days)
10
15
20
25
5
Time (days)
10
15
20
25
0
5
Time (days)
10
15
20
25
4500 (d)
4500
(e)
4500
(f)
4000
4000
4000
Tumour size (mm²)
3500
Tumour size (mm>)
3500
Tumour size (mm3
3500
3000
3000
3000
2500
2500
2500
2000
2000
2000
1500
1500
1500
1000
1000
1000
500
500
500
0
20 25 Time (days)
10
0
5
15 20 25
0
5 15
25 30 Time (days)
Time (days)
Figure 9: Posterior predictive distributions for two pancreatic cancer datasets. 90% Posterior predictive interval plots for SNPE (a), preconditioning step (b) and PNPE (c) for the 26-day dataset; 90% Posterior predictive interval plots for SNPE (d), proconditioning step (c) and PNPE (f) for the 32-day dataset.
these model simulations. Here we used the acceptance rale as the stopping criterion for the preconditioning step, but another option could be to check after each iteration of SMC ABC and stop the preconditioning if a suitable NCDE is found. Secondly, the choice of an unconditional density estimator necessitates careful consideration. In scenarios involving low-dimensional parameter spaces, a kernel density estimator might be preferable option compared to unconditional normalizing flows.
It is important to note that the preconditioning step of our method requires manual selection of the summary statistics or discrepancy function. However, this choice is not as critical as in typical ABC application, since we aim to remove poor parts of the parameter space rather than achieve highly accurate posterior approximation. If it is difficult to manually select summary statistics, the ABC preconditioning step could use discrepancy functions that do not rely on summary statistics such as Wasscrstein ABC (Bernton et al. 2019) and K2-ABC (Park al., 2016), followed by an NPE method that incorporates a summary network (Jiang et al. 2017). This approach would have the benefit that the training of the summary network will not be adversely affected by extreme simulated data. In this paper we used the same summary statistics in the preconditioning step for the subsequent rounds of NPE. However, note it would be possible to use different summaries after the preconditioning step, and possibly use automated summary statistic selection methods [Fearnhead & Prangle 2012 Chen et al. 2023
In this paper we considered the well-specifted scenario, where the model is either known to be correct or can provide a good fit to the data with suitable choice of paramcter values. Standard neural SBI mothods are known potentially perform poorly under model misspecification Bon et al. (2023); Cannon et al. (2022): Schmitt et al. (0233.. Our preconditioning method may be useful in the misspecified scenario, since ABC
15

---

Published in Transactions on Machine Learning Research (09/2024)
are known to perform reasonably well under model misspecification. That is, ABC still converges onto the pseudo-true parameter value Frazier et al. (2020) The preconditioning step could be followed by a robust
neural SBI method such as Kelly et al. (0024); Huang et al. (2024); Gloeckler et al. 2023); Ward et al. 2022 We plan to explore this in future research.
Overall, PNPE employs a preconditioning step to focus on important parts of the parameter space, thereby creating a good starting point for training SNPE and enhancing estimation accuracy. We have empirically demonstrated that PNPE is capable of producing more accurate estimations in complex real-world problems.
Acknowledgement
We thank the computational resources provided by QUT's High Performance Computing and Research Support Group (HPC). Xiaoyu Wang, Ryan P. Kelly and Christopher Drovandi were supported by an Australian Research Council Future Fellowship (FT210100260).
References
Joseph Ayleu-Bullock, Carolina Cuesta-Lazaro, Arnau Quera-Bofarull, Miguel Icaza-Lizaola, Aidan Sedgewick, Henry Truong, Aoife Curran, Edward Elliott, Tristan Caulfield, Kevin Fong, et al. June: open-souree individual-based epidemiology simulation. Royal Society Open Science, 8(7):210506, 2021.
Mark A Beaumont. Approximate Bayesian computation. Annual Review of Statistics and Its Application, 6:379 403, 2019.
Mark A Beaumont, Jean-Marie Cornuet, Jean- Michel Marin, and Christian P Robert. Adaptive approximate Bayesian computation. Biometrika, 96(4))983399)) 2009.
Espen Bernton, Pierre E Jacob, Mathieu Gerber, and Christian P Robert. Approximate Bayesian computa- tion with the wasserstein distance. Journal of the Royal Slatistionl Sociely Series B: Statislical Methodology, 81()):235-26) 2019.
Gérard Biau, Frédéric Cérou, and Arnaud Guyader. New insights into approximate Bayesian computation. In Annales de L'IHP Probabilités et Statistiques, volume 51, pp. 376 403, 2015.
Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D Goodman. Pyro: Deep universal probabilistic programming. The Journal of Machine Learning Research, 20(1):973 978, 2019.
Michael GB Blum. Approximate Bayesian computation: a nonparametric perspective. Journal of the
American Statistical Association, 105(491):1178-11))) 1187, 2010.
Joshua Bon, Adam Bretherton, Katie Buchhorn, Susanna Cramb, Christopher Drovandi, Conor Hassan, Adrianne L Jenner, Helen Mayfield, James M McGree, Kerrie Mengersen, al. Being Bayesian in the 2020s: opportunities and challenges in the practice of modern applied Bayesian statisties. Philosophical Transactions of the Royal Society A, 381 2023.
Patrick Cannon, Daniel Ward, and Sebastian M Schmon.
Investigating the impact of model misspecification in neural simulation-base based inference. arXiv preprint arXi::2209.0188 01845, 2022.
Yanzhi Chen, Michael U Gutmann, and Adrian Weller. Is learning summary statistics necessary for likelihood-free inference? In International Conference on Machine Learning, pp. 4529-4544 PMLR, 2023.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117(45):30055-30022 2020.
Katalin Csiller,, Michael GB Blum, Oscar E Gaggiotti, and Olivier François. Approximate Bayesian com- putation (ABC) in practice. Trends in Ecology & Fvolution, 25(7):410-118, 2010.
16

---

Published in Transactions on Machine Learning Research (09/2024)
Maximilian Dax, Stephen R Green, Jonathan Gair, Jakob H Macke, Alessandra Buonanno, and Bernhard Schölkopf. Real-time gravitational wave science with neural posterior estimation. Physical Review Lellers, 127((4):241103, 2021.
Natalí SM de Santi, Francisco Villacscusa-Navarro, Raul Abramo, Helen Shao, Lucia A Percz, Tiago Castro, Yueying Ni, Christopher C Lovell, Elena Marttnne,,,,,, Federico Marinacci, et al. Field- level simulation-based inference with Galaxy catalogs: the impact of systematic effects. arXiv preprint wrXiv:231011525) 2023.
Michael Deistler, Pedro J Goncalves, and Jakob H Macke. Truncated proposals for scalable and hassle-free simulation-based inference. Advances in Neural Informalion Processing Syslems, 35:23135 23149, 2022a.
Michael Deistler, Jakob H Macke, and Pedro J Gonçalves. Energy-efficient network activity from disparate circuit parameters. Proceedings of the National Academy of Sciences, 119(4)):c2206632111 2022b.
Hadi Mohaghegh Dolatabadi, Sarah Erfani, and Christopher Leckie. Invertible generative modeling using lincar rational splines. In International Conference on Artificial Intelligence and Statistics, pp. 4236-4246.
PMLR, 2020.
Christopher Drovandi and David T Frazier. A comparison of likelihood-fre methods with and without summary statistics. Statistics and Computing, 32(3):42, 2022.
Christopher Drovandi, David Nott, and David T Frazier. Improving the accuracy of marginal approxima- tions in likelihood-free inference via localisation. Journal of Computational and Graphical Statistics, pp. 1-19, 2023.
Christopher C Drovandi and Authony N Pettitt. Estimation of parameters for macroparasite population evolution using approximate Bayesian computation. Biometrics, 67(1):225-233 2011.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows Advances in Newral Information Processing Systems, 32, 2019.
Paul Fearnhead and Demnis Prangle. Constructing summary statistics for approximate Bayesian computa- tion: semi-automatic approximate Bayesian computation. Journal of the Royal Statistical Society Series B: Statistical Methodology, 74(3):419 474, 2012.
Alexander Fengler, Lakshmi N Govindarajan, Tony Chen, and Michael J Frank. Likelihood approximation networks (L.Nss) for fast. inference of simulation models in cognitive neuroscience. Elife, 10:e65074, 2021.
David T Frazier, Christian P Robert, and Judith Rousseau. Model misspecification in approximate Bayesian computation: consequences and diagnostics. Journal of the Royal Statistical Society Series B: Statistical Methodology, 82(2))421 444, 2020.
David T Frazier, Christopher Drovandi, and David J Nott. Better Together: pooling information in
likelihood-free inference. arXiv preprint 2xXx:2212.00658, 2022a.
David T Frazier, David J Nott, Christopher Drovandi, and Robert Kohn. Bayesian inference using synthetic likelihood: asymptotics and adjustments. Journal of the American Statistical Association, pp. 1-12, 2022b.
Manuel Gloeckler, Michael Deistler, and Jakob H Macke. Adversarial robustness of amortized Bayesian inference. pp. 11493-11552 2023.
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.
Daolang Huang, Ayush Bharti, Amauri Souza, Luigi Acerbi, and Samuel Kaski. Learning robust statistics for iinilattio-baseee inference under model misspecification. Advances Neural Information Processing Systems, 36, 2024.
17

---

Published in Transactions on Machine Learning Research (09/2024)
Adrianne L Jenner, Federico Frascoli, Adelle CF Coster, and Peter S Kim. Enhancing oncolytic virotherapy: Observations from Voronoi Cell-Based model. Journal of Theorelical Biology, 485: 485:110052, 2020.
Bai Jiang, Tung-yu Wu, Charles Zheng, and Wing H Wong. Learning summary statistic for approximate Bayesian computation via deep neural network. Statistica Sinica, pp. 1595-1618, 2017.
Ryan Kelly, David J Nott, David Tyler Frazier, David Warne, and Christopher Drovandi. Misspecifeattinn-obass sequential neural likelihood for simulation based inference. Transactions on Ma- chine Learning Research, 2024.
Pyung- Hwan Kim, Joo-Hyuk Sohn, Joung- Woo Choi, Yukyung Jung, Sung Wan Kim, Seungjoo Haam, and Chac-Ok Yun. Active targeting and safety profilc of PEG-modified adenovirus conjugated with herceptin. Biomalerials, 32(9))2314 2326, 2011.
Jarno Lintusaari, Michael U Gutmann, Ritabrata Dutta, Samuel Kaski, and Jukka Corander. Fundamentals and recent developments in approximate Bayesian computation. Systematic Biology, 66(1):e66 e82, 2017.
Jarno Lintusaari, Henri Vuollekoski, Antti Kangasrääsiö, Kusti Skytén, Marko Järvenpää, Pekka Marttinen, Michael U Gutmann, Aki Vehtari, Jukka Corander, and Samuel Kaski. ELFI: Engine for likelihood-free inference. Journal of Machine Learning Research, 19(16):1-7, 2018.
Jan-Matthis Lueekmann, Pedro J Goncalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher, and Jakob H Macke. Flexible statistical inference for mechanistic models of neural dynamics. Advances in Newral Information Processing Systems, 30, 2017.
Jan-Matthis Lueckmann, Jan Boelts; David Greenberg, Pedro Goncalves, and Jakob Macke. Benchmarking simulation- based inference. In International Conference on Artificial Intelligence and Statistics, pp. 343- 351. PMLR, 2021.
Philippe Marchand, Morgan Boenke, and David M Green. A stochastic movement model reproduces patterns of site fidelity and long-distance dispersal in population of fowler's toads (anaxyrus fowleri). Ecological Modelling, 360:63-60 69, 2017.
Frank A Meineke, Christopher S Potten, and Markus Loeffler. Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Proliferation, 34(4):253 266, 2001.
Siddharth dishra-Sharman Inferring dark matter substructure with astrometric lensing beyond the power spectrum. Machine Learning: Science and Technology, 3(1):01LTO03, 2022.
George Papamakarios and lain Murray. Fast s-free inference of simulation models with Bayesian conditional density estimation. Advances in Neural Information Processing Systems, 29, 2016.
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. Advances in Newral Information Processing Systems, 30, 2017.
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): 2021.
Mijung Park, Wittawat Jitkrittum, and Dino Sejdinovic. K2-ABC: Approximate Bayesian computation with kernel embeddings. Tn Artificial intelligence and slatistics, pp. 398 008-407. PMLR, 2016.
Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathicu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Seikit-learn: Machine learning in Python. the Journal of Machine Learning Research, 12:2825 2830, 2011.
Leah F Price, Christopher Drovandi, Anthony Lee, and David Nott. Bayesian synthetic likelihood. Journal of Compulational and Graphical Statistics, 27(1):1-11, 2018.
18

---

Published in Transactions on Machine Learning Research (09/2024)
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International Con- ference on Machine Learning, pp. 1530-1538. PMT.R, 2015.
Marvin Schmitt, Paul-Christian Bürkmer, Ulirich Köthe, and Stefan T Radev. Detecting model misspecifi- cation in amortized Bayesian inference with neural networks. pp. 541-557, 2023.
David Shih, Marat Freytsis, Stephen R Taylor, Jeff A Dror, and Nolan Smyth. Fast parameter inference on pulsar timing arrays with normalizing flows. arXiv preprint MXXi:23100112000 2023.
Scott A Sisson, Yanan Fan, and Mark M Tanaka. Sequential Monte Carlo without likelihoods. Proccedings of the National Academy of Sciences, 104(6):1760 1765, 2007.
Scott A Sisson, Yanan Fan, and Mark Beaumont. Hlandbook of approzimate Bayesion computation. CRC Press, 2018.
Alvaro Tejero- Tejero-Cantero, Cantero, Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pedro Gonçalves, David S Greenberg, and Jakob Macke. SBI: toolkit for simulation-based inference. Journal of Open Source Software, 5(52):2505, 2020.
Owen Thomas, Henri Pesonen, Raquel Sá-Leão, Hermínia de Lencastre, Samuel Kaski, and Jukka Corander. Split-BOLFF for for misspecification-robne likelihood free inference in high dimensions. arxiv preprint arxiv:2002. 09377, 2020.
Brandon M Turner and Trisha Van Zandt. A tutorial on approximate Bayesian computation. Journal of Mathemalical Psychology, 56(2):69 85, 2012.
Samantha Jane Wade. Fabrication and preclinical assessment of drug cluting wet spun fibres for pancreatic
cancer treatment. 2019.
Xiaoyu Wang, Adrianne L Jenner, Robert Salomone, David J Warne, and Christopher Drovandi. Cali- bration of agent based models for monophasic and biphasic tumour growth using approximate Bayesian computation. Journal of Mathematical Biology, 88(3):28, 2024.
Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian Schmon. Robust neural
posterior estimation and statistical model criticism. Advances Neural Information Processing Systems, 35:33845-338555 33859, 2022.
David J Warne, Ruth E Baker, and Matthew J Simpson. Rapid Bayesian inference for expensive stochastic models. Journal of Compulational and Graphical Statistics, (1(2):112-524 2022.
Antoine Wehenkel, Jens Behrmann, Andrew C Miller, Guillermo Sapiro, Ozan Sener, Marco Cuturi, and Jörn-Henrik Jacobsen. Simulation-based inference for cardiovascular models. arXiv preprint urXiv:2307.19116 2023.
Timothy O West, Lue Berthouze, Simon F Farmer, Hayriye Cagnan, and Vladimir Litvak. Inference of brain networks with approximate Bayesian computation-assessing face validity with an example application in Parkinsonism. Newroimage, 236:118020 2021.
Christina Winkler, Daniel Worrall, Emiel Hoogeboom, and Max Welling. Learning likelihoods with condi- tional normalizing flows. arXio preprint 191.::11220000 2019.
19

---

Published in Transactions on Machine Learning Research (09/2024)
A Further background for SMC ABC
We provide a detailed description of the adaptive SMC ABC algorithm we used in this paper and provide pseudocode in Algorithm 2 for reference. This algorithm starts by drawing N independent samples from the prior distribution p(0), represented as {0, [0-) For each sample 0g (known as a particle), the algorithm simulates a dataset z from the stochastic model and calculates the corresponding discrepancy pi = p(z,, to); resulting in the pair set {0s, pr) These pairs are then arranged in order of increasing discrepancy such that P P2 <............................<.. < PN- The first tolerance threshold, 61, is set as the largest discrepancy, PN- To move through the target distributions, the algorithm adjusts the tolerance dynamically. The next tolerance, e, is set as PN-N., where N, [NaJ, and is a tuming parameter. Essentially, in each step, the algorithm discards the top x 100% of particles with the highest discrepancies. After discarding these particles, only N- No particles remain. To replenish the set back to N particles, the algorithm resamples Na times from the 'alive' particles, copying both the parameter and discrepancy values. This process, however, leads to duplicates in the particle set. To add variely to the set, the algorithm applies an MCMC ABC kernel each resampled particle. The parameters for the MCMC proposal distribution ge(-l-) are derived from the current particle set. For instance, if using a multivariate normal random walk proposal, its covariance Eg is based the particle set's sample covariance. The acceptance of proposed parameter (assuming symmetric proposal) and simulated data is determined by the equation:
Pt m--min min 11 p()) p()) (p(Ea, (p(cor))<))) <)<())))))))))) <ei))
(6)
where Ø ~ q(-10) and ă ~ p(-10) are proposed parameter values and dataset, respectively. However, proposals may be rejected, lcaving some particles unchanged. To address this, the algorithm performs RG iterations of the MCMC kernel on each particle, where Rg log(c) log(1-pt-5) where e tuning parameter of the algorithm that can be interpreted as the probability that particle is not moved in the R, iterations. The acceptance probability pt is estimated from trial MCMC ABC iterations and used to compute RG for the next set of MCMC ABC iterations. For this adaptive SMC ABC algorithm, two stopping rules can be used. The first stopping rule halts the ABC algorithm when the maximum discrepancy is below a set tolerance, 6T- The second stopping rule terminates the algorithm when the MCMC acceptance probability pt falls below a predefined threshold Pace- Here we choose the second rule, and since we only require short run of ABC, set Pace to be higher than what is typically used in an ABC analysis.
20

---

Published in Transactions on Machine Learning Research (09/2024)
Algorithm 2 Adaptive SMC ABC
The observed data To: the stochastic model p(z]0), distance function p(:,), prior distribution p(0), number of particles N, tuning parameters and e. for adaptive selection of discrepancy thresholds and selecting the number of MCMC iterations in the move steps, target tolerance initial number of trial MCMC itcrations Sinit, minimum acceptable MCMC acceptance rate Pmin for == ,....,N do
Input:
Simulate x2 ~ p(x((6)) where 0, ~ p(0)
Compute Pi - (i))
end for Sort {0.)
N by {pili=1 such that #1 P2 S ... S PN Set N₂ = [NN], =
= Na-Na: 61 PN, S. = Smitt p 1
while 1 > ET or pi > Pmin do
Compute E as the sample covariance matrix of [1 N.
Generate {0; I. Po.Ji-NNAA N-N++1 by resampling from {0:) for i= N N₂ + +1....,N do
N-N. with replacement
for j = : S. do
Simulate 전 ~ p(x(0) based on proposal Ø ~~ No E) Compute p( 3)
Compute v: min 11, v(6) [
[( < ti)
p(0)
With probability pt set 0, and Ps = pú otherwise, retain the current values of 0g and Pt
end for end for
pe E Li=N-N+H 2==1 11=1 (S, - -Na)
R2 = [log(c)/( + log(1 1
for i
N₂ + ...., N do
for j = R₁ S ................................ R do Simulate & ~ p(210) based on proposal Ó ~ N((,,
Compule p( 7)
Compute Pt min 1,.
pí p(o)
I(p < <1)
With probability pt set 0, 0 and p. p; otherwise, retain the current values of 0 and p
end for end for pecc
ENNN Li=N-Net No+1 EBn / -Na))
S.+1 2
Sort [0:] by {p:li such that P1 S P2 E ... PN
Set cx41 = PN-N.: 6 PN
t= ++1
end while
N return Samples (0,)== from ABC posterior
21

---

Published in Transactions on Machine Learning Research (09/2024)
B Experimental Details
We use the adaptive SMC ABC algorithm Drovandi & Pettitt (2011) in the preconditioning step for all experiments. We set the tuning parameters as a 0.5, 0.01 and use 1K particles for the algorithm. As the stopping rule, we set the target MCMC acceptance rate at 10%, unless otherwise specified. For the unconditional density estimator, employ unconditional normalizing flows using the Pyro package Bingham et al. 2019 with spline coupling layer using the transformation:
Yı:d g6 (Xi:a) Y(o+1):D = hap ( (d+1))DiX1.dd
(7) (8)
where X are the inputs, Y are the outputs, e.g., Xndd represents the first d elements of the inputs, ga is either the identity function or an elementwise rational monotonic spline with parameters 6, and has where & is element-wise bijection parameter, is a conditional elementwise spline, conditioning on the first d clements. Regarding the neural networks, we use four fully-connoeted layers and set the count bins to 16. Furthermore, the dimensions of the parameter space are less than three, indicating a low-dimensionall case, we also consider kernel density estimation with Gaussian kernel as the unconditional density estimator, implemented in the Scikit-learn package Pedregosa ct al.
20011). For APT and TSNPE, we use the implementation of the sbi package Tejero-Cantero et al. 2020 with default settings.
For SNPE and TSNPE, we use conditional neural spline flows Durkan et al. (2019). We use five coupling layers, with each coupling layer using multilayer perceptron of two layers with 50 hidden units. The flow is trained using the Adam optimizer with a learning rate of 5 x 10-4 and a batch size of 256. Flow training is stopped when either the validation loss, calculated on 10% of the samples, has not improved over 50 epochs or when the limit of 500 epochs reached.
B.1 Benchmarking examples
A number of popular benchmarking models in the SBI literature exist where the ground truth posteriors are available. For the two-moon model and the SLCP model, we follow the model specifications in Lueckmann let al. (2021). For the rest of the models, we refer to Lintusaari et al. (2018). After the preconditioning step, we generate 10k simulations from the unconditional normalizing flows and pair them with their simulated data to process SNPE.
B.1.1 Two moons model
The two moons model exhibits both global (bimodality) and local (crescent shape) structures in the posterior. For the proconditioning step, the total number of model simulations is around 30k. Since 0 is low-dimensional in this model, use KDE the unconditional density estimator and achieve results similar to those obtained with unconditional normalizing flows.
r croo(()) 0.25 rsin(a)
- + 02l/22
Simulator
0
where x//m ~~ (((0.,0000113
(0, 02
Prior
(01 02), 0, U-1,1) for 1,2
Dimensionality
R² RR
References
Greenberg al. (2019)
B.1.2 SLCP with Distractors model
The SLCP model is designed to have simple likelihood and complex posterior, with uninformative dimen- sions (distractors) added to the observations. The preconditioning step uses around 35k model simulations.
22

---

Published in Transactions on Machine Learning Research (09/2024)
Simulator
== (#1..... 1000); # = ply):
where re-orders the dimensions of with fixed random permutation;
1
U[1:8] (mmo, So),
V[0:100)
te E"
20
where mg
01
S
s
p8152
81 6
82 0j p tanh((;))
02]
821
$2
#² (N, to is student t-distribution with degree of freedom 2,
-j j.k NN(,,)) j k. Mj = 3e where a NN(0,) 1), Ej = otherwise
Prior
0 (01, 101,00,03.. 0s),
0, ~ ( for i = ,,...,
Dimensionality
0 RR, IR 100
References
Greenberg et al. (2019)
B.1.3 MA(2) model
The moving average model of order 2 (MA(2)) is a univariate time series model often used as a toy example in the ABC literature. We assign the true parameter values 0 (0.6, 0.2) and simulate zo 10000 In this model, To and are high-dimensional, so we use the first two lags of the autocovariance function and the variance (lag 0) as the summary statistics. The preconditioning step required 18k simulations, as the larget MCMC acceptance rale is set at 12% Given that the dimension of the parameter space is relatively low, we employ KDE as an additional unconditional density estimator, utilizing a Gaussian kernel with a bandwidth selected by Silverman method. This method achieves results similar to those obtained using an unconditional normalizing flow.
Simulator
2 (II, 1000), = e, + 01e;-1 + 02er-2 where e, 1)
Prior
0 = (01, 02). 0: ~ U(-1,1) for i = 1,2
Dimensionality
0 R2, TRRR 100
1 T
Summary Statistics
Sample autocovariances Sy(x) T
Σ
Ft for j = 0,1,2
t=1+)
References
Marin et al. (2012)
B.1.4 Univariate g-and-k model
The univariate g-and-k model is a popular benchmark in the statistical SBI literature. The four param- eters = (A, k, k) control the location, scale, skewness, and kurtosis, respectively. We assign the true parameter values as 0 = (3, 2, 0.5) and use them to simulate H, E g independent observations. For the preconditioning step, take the set of octiles as summary statistics as they are robust measure of skewness and kurtosis (Drovandi & Pettitt 2011). Setting the target MCMC acceptance rate to 20% results in 18k model simulations for the preconditioning step. Here we used larger acceptance rate to ensure that the preconditioning step does not use too many model simulations.
Simulator
(# ... 250), where x; = a + + e- tanh[2/2](1 2, N(O, 1) and 0.8
Prior
0 b,,b, and each parameter has an U (0, 10) prior 0 E R 50
Dimensionality
References
A Allingham al. (2009); Drovandi & Pettitt 12011b)
23

---

Published in Transactions on Machine Learning Research (09/2024)
B.1.5
Toad model
The toad movement model proposed in (Marchand et al.) 2017p is an individual-basee model for simulating the dispersal of Fowler's toads (Ananyrus Jowleri). It. has been considered as a lest example several times the SBI literature (e.g. Frazier et al. (2022a Drovandi Frazior (2022)).
This model captures two known behaviours of amphibians: the tendeney to return to previously visited loca- tions and small chance of Long-distance movement. Dispersal distance is modeled with Lévy alpha-stable distribution characterized by a stability factor and scale factor 7. The Lévy alpha-stable is symmetrically centred around zero but has heavy tails, allowing the simulation of both frequent short-distance and rare Long-distance movements. Toads are modeled to remain at their refuge site during the day and move to forage independently at night. After foraging, the toad may either remain at their current location or return previously visited refuge site. We follow "model 1 Marchand al., (2017), where each former refuge site has equal probability of being returned to. The probability of toad returning to former refuge site is constant probability, denoted po. Hence, there are three parameters to infer, that is 0 (a, (Pp).
The simulated data generates an observation matrix (representing Euclidean distance in metres from the origin) of 66 Loads across 63 days. We assign the true parameter values as = (1.7, 35.0, 0.6) and use them to simulate Zo E 163366 Since the data can be considered as high-dimensional, to summarise the data we first construct four sets of displacement vectors with time lags of 1, 2, and days. For each sct, we take the log difference between the 0, 0,0.1,.... quantiles, the number of absolute displacements less than 10m, and the median of the absolute displacements greater than 10m, resulting in total of 48 summary statistics (12 for each time lag).
Prior
U0,,1), 7~ 7~ UO, 100), Po 0.99
Dimensionality
0 E R I E 63x66 S(x) e gs
References
Marchand et al. (2017)
B.2 High-dimensional SVAR model
We consider our illustrative example in high- nigh-dimensional setting with A = 20, which leads to 21 parameters that need to be estimated. We set the true paramoter values as follows:
0 (-0.2764, -0.7765, 0.8231, - 0. 1972 - 0.22544 0.6334,
0.4465, - -0.961,, 0.0647, -0.1 0.0795, -0.5464, -0.9354,
0.4639. 0.7851, 0..8333, -0.1408, 0.7032, 0.8321, 0.1000),
and use these values to simulate the observation dataset. We employ the same summary statistics as in the low-dimensional example and use a uniform distribution as the prior, constrained between -1 and for the l parameters and between 0 and for a.
Since the BSL method the gold standard for this example, we use the standard BSL method proposed by Price al. (2018), with the total number of simulations set to 20 million.
B.3 BVCBM
The BVCBM simulation begins by initializing a square domain with cells arranged in a hexagonal lattice. The cell at the center the domain is identified as cancer cell, while the others are designated as healthy cells. The simulation procoeds until the tumor reaches a volume of 100 mm? accordancce with experimental measurements (Wade 2019 Kim al., 2011). When this volume is attained, the distribution of healthy and cancerous cells within the lattice recorded. This configuration then serves the starting point to
24

---

Published in Transactions on Machine Learning Research (09/2024)
simmlate tumor growth over the desired mumber of days. The model progresses by determining whether cancer cell proliferates at each timestep, as described by the equation:
Pd Po (1-2a d
(9)
dmmx
where Pa is the probability of cell division, Po is the initial division ratc, is the current cell density, and dmax the maximun densily. For cancer cells that do not proliferate, the model assesses their potential transition into invasive cells, governed by the probability Ppsc. Subsequently, the positions of all cells, both healthy and cancerous, updated using Hooke's law:
r:(1 + At) = r() + - F.())/ - -(0+ --(0)//=+(0+++
Pts(t) [r)(()]
(SL(() I|rs, (1))..
(10)
Here, (t + At) denotes the updated position of cell i, j is the cell motility coefficient, (t) is the force on cell i, À is a mechanical interaction coefficient, rı.j (t) the vector between cells and j, and st, is the natural length of the spring connecting the two cells. The parameters for the mechanical interactions, such X and μl, sourced from prior studies in the literature Meineke et al. [2001). See Jenner et al. (2020) Wang et al. 2024) for more detailed model simulation.
Four parameters 0 = (po, Ppse; dmax, Pago) control the tumor growth during a single phase, which is a period when the tumor grows based on fixed values for these four parameters. For the biphasic model, an additional parameter т introduced, representing the time at which the tumor growth pattern changes, that is, the values for 9 change. Therefore, for BVCBM, we need estimate nine parameters for two pancreatic cancer datasets, denoted as 0 (pb- pae d Dage); 02 (pō Ppacr dmax Page); and T, that 0 = (0,,0,,77 r).
The parameter Ppsc, which is the probability of tumor cell invasion into healthy cells, significantly affects the simulation time. The value of Ppsc should be around 10-5, indicating that an increase in probability will require more cells to be simulated. Moreover, a smaller value of Pasc results in simulation time. For PNPE, the total simulation time for 15k simulations (26-day datasct) and 17k simulations (32-day datasct) for the preconditioning step is approximalely 13 and 15 minutes, respectively, whereas SNPE requires around hour for the first round (i.e. based on samples from the prior) of 10k simulations. This is because the preconditioning step is effective at quickly climinating values of Prsc that lead to longer model simulation times.
25

---

Published in Transactions on Machine Learning Research (09/2024)
C Further experimental results
C.1 SVAR
We generate 10 different datasets for the SVAR model in 6 dimensions based on predefined parameter values shown in Table 3 Then, we run SNPE and PNPE with 10 different, random seed values for each dataset to investigate reproducibility. For preconditioning, we use a 10% stopping rule and train unconditional normalizing flows based on ABC posterior samples for different random seed values for cach dataset. To quantify the results, use MMD as the metric. For each dataset, compute the MMD values wetw the approximate distribution and the reference distribution, allowing to fit. By computing the MMD for the approximate posterior based each set of random seed values, we produce a boxplot corresponding to each method.
value to summarize the
DATASET
0
02
0g
04
0g
O.
1
0.9749
-0.2564
0.4230
0.8749
0.6453
-0.1772
2
0.2777
-0.2932
0.5483
-0.4679
0.1549
-0.6365
3
0.3193
-0.2509
0.7801
0.5103
-0.4358
0.0317
4
-0.9295
-0.3267
0.2930
-0.9326
-0.6113
-0.5015
5
0.5146
0.2805
-0.6695
0.7410
-0.4518
-0.4366
6 0.6401
-0.6037
0.3730
-0.6739
-0.7237
-0.5494
7
-0.9378
-0.6438
-0.8434
-0.9219
-0.4032
0.5577
8
-0.2663
0.9098
0.7383
0.1438
-0.0968
0.7992
9
0.9185
-0.1454
0.5147
-0.2921
0.0595
0.7227
10
-0.4766
0.7152
0.9845
-0.7105
-0.2426
-0.4923
Table 3: True values to generate observational datasets: Each row refers to the true parameters that we used to generate the observation datasets I, with the additional parameter o set to a constant value of 0.1.
26

---

Published in Transactions on Machine Learning Research (09/2024)
Dataset
Dataset
Dataset
0.18
0.5
0.16
0.25
0.14
0.20
0.4
0.12
0.10
3.15
0.08
0.06
0.10
0.04
ABC
PNPE
ABC
NPE
ABC
PNPE
Dataset
Dataset 5
Dataset 6
0.20
0.30
0.30
0.18
0.25
0.16
0.25
0.14
0.20
0.12
3.15
0.10
0.15
0.10
0.08
NPE
ABC
PNPE
NPE
ABC
PNPE
NPE
ABC
PNPE
Dataset
Dataset 8
Dataset 9
0.40
0.9
0.40
0.8
0.35
0.35
0.7
0.30
0.30
0.25
0.25
0.20
0.3
0.20
0.15
0.10
0.2
0.15
NPE
ABC
PNPE
NPE
ABC
PNPE
NPE
ABC
PNPE
Dataset
0.30
0.25
0.20
0.15
0.10
0.05
ABC
PNPE
Figure 10:
Performance for SVAR with 10 different datasets: The box-plots refer to the MMD
values computed from approximate posterior for each algorithm with 10 different random seed values.
27

---

Published in Transactions on Machine Learning Research (09/2024)
C.2 BVCBM
In this section, we present prior predictive distributions of tumour volumes in (a) and (b) for two pancreatic datasets (in (c) and (d) we show the samn plots but on the log scale). We also show the posterior predictive rint on the log scale obtained with different methods in Figure 12 Then we present the bivariate posterior density and marginal posterior density plots for the BVCBM additional results. It evident from Figures 13 and 14 that PNPE provides more precise estimation than SNPE for both panercatic cancer datasets.
(a) 30000
(b) 30000
25000
25000
Tumour size (mm3)
Tumour size (mm>)
20000
20000
15000
15000
10000
10000
5000
5000
0
5
10
15
20
2
0
10
15
20
25
30
Time
(days)
Time
(days)
(c) 12
(d)
Tumour size (log scale)
11
Tumour size (log scale)
10.
10
9
8-
8
7
7-
6
6
5
5
4
4
0
5
10
15
20
25
30
0
5
10
15
20
25
30
Time (days)
Time (days)
Figure 11: Prior predictive distribution for BVCBM. We sample 10k parameter values from prior
distribution and plot the prior prodictive distribution for two pancreatic datascts. In (a) and (b), the plots are in regular scale, and in (c) and (d), the plots are in log scale.
28

---

Published in Transactions on Machine Learning Research (09/2024)
(a) 12
(b) 12
Tumour size (log scale)
Tumour size (log scale)
Tumour size (log scale)
11
10
10
9
9-
8
7-
7
7-
5-
5.
10
15
20
25
0
10
15
20
25
10
15
20
25
Time (days)
Time (days)
Time (days)
(2)
(f)
Tumour size (log scale)
11-
11
11
10-
Tumour size (mm)
Tumour size (log
9
9-
8-
8
8
6-
6-
6-
5
10
15
20
25
30
0
5
20 Time (days)
10
15
25
30
4
10
Time (days)
Time (days)
Figure 12: Posterior predictive distributions for two pancreatic cancer datasets in log scale. 90% Posterior predictive interval plots for SNPE (a), preconditioning step (b) and PNPE (c) for the 26-day dataset; 90% Posterior prodictive interval plots for SNPE (d), preconditioning step (c) and PNPE (f) for the 32-day dataset.
29

---

Published in Transactions on Machine Learning Research (09/2024)
(a)
(b)
600
1000
1000
P
20
100
400
500
0.0005
10
8.00
0.01
8.00
0.01
200
0.00
0.01
0.00
0.01
0.00
0.01
0.00
0.01
0.0010
0.004
0.004
0.0010
grge
20
0.0005
0.002
0.002
0.0005
10
0.0000
0.000
0.000
250500
250500
250500
250500
0.0010
600
0.06
400
0.2
0.0010
0.0005
200
0.02
0.0005
0.0000
10
20
20
10
20
10
20
0.0000
10 20
10 20
0.0010
600
1000
1000
0.0005
200
500
500
0.00
0.01
0.00
0.01
0.00
0.01
8.00
0.01
8.00
0.00
0.01
0.0010
0.0010
0.004
400
20
0.02
0.0005
0.0005
0.002
200
10
0.01
0.0000
250
500
0.0000
0.000
250
500
0.00
250500
P/
glge
9
g
Figure 13: Bivariate density plots for the pancreatic dataset with 26-day measurements for (a) SNPE and
(b) PNPE. The diagonal entries represent the marginal posterior densities for (phacs Page T, Pf Iage).
30

---

Published in Transactions on Machine Learning Research (09/2024)
(a)
2000
1500
1000
29
400
0'000
500
0.000
0.005 0.010
200
0.000
0.005
0.010
0.000
0.005
0.010
0,000
0.005 0.030
200
0,000
0.603
5.010
6.00160
0.0000
0.0100
0.00075
600
0.0075
0,008D
0.00020
400-
0. 00022
6,0035
142
750
750
0.0035
600
0.03
0.0034
0.02
0.02
0.00010
0. 0.00025
10
20
200
10 20
35
0.00
50
70
35
0.00000
200
20
30
0 0.0005
500
500
400
0.0002
10
200
0,010
0.003
8000
0.005 0.010
0.000 0.005
0.00180
0.0000
30
0,00075
0.0006
600-
0,00050
0.034
0.004
0.0002
0.00025
01022
0.002
2,00000
Pasc
Olae
Ppsc
9
Figure 14: Bivariate density plots for the pancreatic dataset with 32-day measurements for (a) SNPE and (b) PNPE. The diagonal entries represent the marginal posterior densities for (ppecr Page T, Plee- Page).
31


Reviews:
[{"review_id": "ogBl65Y8HS", "reviewer": "Reviewer_bPPq", "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\nI think this is a pragmatic approach to setting the prior when using NPE.\n\nThe authors use a model system that has been characterized previously and is relatively tractable.\n\n“Finally, we consider a challenging real-world problem in cancer biology: calibrating the biphasic Voronoi cell-based model (BVCBM) (Wang et al., 2024) that models tumor growth.” I think this is a very interesting problem!\n\nWeaknesses\n\n“ NPE uses N training pairs of simulator parameter values and simulated datasets, {θi , xi} N i=1, to estimate the posterior distribution p(θ|x) (Papamakarios et al., 2021)” I’m confused why Papamakarios 2021 is the reference for this work, as it does not explicitly define Neural Posterior Estimation, nor do I readily see this formula. Generally, where does this simulated data and parameters come from?\n\nBSL needs to be defined. Since it is used as the reference for this work, it’s impossible to compare advances described herein.\n\nI don’t know how to interpret what is “good” from FIgures 1-6. Moreover, the results are entirely subjective. Is there any summary of the different approaches? What is “good” for posterior predictive distributions and marginal posterior distributions? In particular, Figures 5 and 6 are not interpretable. \n\nIn the ABC preconditioning step, what is the distribution of the initial parameters? How does it differ than the priors used previously?\n\n“Finally, we consider a challenging real-world problem in cancer biology: calibrating the biphasic Voronoi cell-based model (BVCBM) (Wang et al., 2024) that models tumor growth.” What is the literature in optimizing BVCBM? Is it standard Bayesian inference? Why is your method a good and reasonable approach here?", "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 describe an approach to set good initial values to Neural Posterior Estimation approaches.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "Requested Changes\n\n“... we run three rounds of SNPE and compare the results with those from BSL, considering BSL results as the gold standard for this example…” What is BSL? Can you please describe this method and associated works?\n\n“...We use summary statistics to reduce the dimension of the data. Following Thomas et al. (2020); Drovandi et al. (2023), we use the lag 1 autocovariance 1 T PT t=2 y i t y j t−1 as the summary statistics,...” What is the role of the summary statistics in NPE? I only see the importance defined for ABC.\n\nI generally find Figures 1-6 to be weak and lacking critical interpretation. Ideally, all of these figures would be in the supplement, and a more parsimonious summary and interpretation is provided by the authors to the reader.\n\nAre there more complex examples that create multimodal data distributions in one mode is very rare. In the ABC sampling approach, would these parameter values be found initially? If not, would this approach capture it?\n\n“However, APT can suffer from a ‘leakage’ issue….” and “We observed a leakage problem occurring in rounds 8 and 5 for the 26-day and 32-day datasets, respectively” What do you mean by a leakage problem? Mathematically, why does this solve the leakage problem? This seems like a pretty important choice for the cancer dataset analysis.\n\n“Figure 7: Posterior predictive distributions for two pancreatic cancer datasets. 90% Posterior predictive interval plots for SNPE (a), preconditioning step (b) and PNPE (c) for the 26-day dataset; 90% Posterior predictive interval plots for SNPE (d), preconditioning step (e) and PNPE (f) for the 32-day dataset.” Where is the comparison with a standard reference approach? What about BSL (which hasn’t been described).\n\nIs there any biological interpretation to your results? Is the modeling good or bad? What can we do with these results?", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "HXesT5Ltn9", "reviewer": "Reviewer_MRjJ", "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": "**Summary:** The paper proposes a method for Likelihood free inference (LFI) combining ABC and SNPE. Typically, the first round of SNPE involves running NPE, and the resulting approximation is used as proposal distribution for SNPE’s second round. The approach proposed here builds a different first approximation, consisting of fitting a normalizing flow to the samples produced by some ABC method.\n\nThe method’s practical usefulness comes down to whether the first approximation provided by ABC is better than the first approximation provided by NPE (in the context of SNPE).\n\nWhy may this be the case? In principle, ABC methods are known to suffer from extremely slow convergence, often slower than NPE. However, in this case, since the ABC approximation will be refined by a subsequent application of SNPE, the ABC method does not need to operate on the challenging regime of epsilon->0 (which has a big effect on the method’s efficiency).\n\nAll in all, if my understanding is correct, the claim in the paper, that the proposed method is useful in practice, can be re-interpreted as follows: The first approximation obtained by fitting a normalizing flow to ABC produced samples (in a “relaxed” configuration with epsilon>>0) is better than the first approximation provided by NPE (i.e. fitting an amortized normalizing flow with samples from the prior and the simulator). This needs to be verified empirically over a range of models and settings (see more below on requested changes, specifically empirical evaluation).\n\n**Method benefits:** The approach seems sensible when we have a single observed dataset (i.e. we are not interested in amortization, as the method cannot handle that scenario). Running NPE as a first step of SNPE may be tackling a harder task than necessary, as NPE yields an amortized approximation that can be used for any potential observation. ABC, on the other hand, focuses on producing an approximation for the set of observations at hand, a potentially simpler problem to solve. Additionally, ABC used as proposed in this work does not need to yield a perfect approximation (which is hard), as its results will be refined with SNPE. It should be useful as long as it provides a better approximation than plain NPE. As mentioned above, ABC does not need to operate in the challenging epsilon->0 regime, as its output will be refined by further rounds of SNPE.\n\n**Method drawbacks:** Reliance on convergence of ABC (for a fixed level of epsilon). While convergence may be easier to achieve for epsilon>>0, if ABC leaves important regions of the parameter space unexplored, the method will likely fail to yield good results (same holds for normal SNPE. If the first approximation provided by NPE is not a good one).", "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 a method for Likelihood free inference (LFI) combining ABC and SNPE. Typically, the first round of SNPE involves running NPE, and the resulting approximation is used as proposal distribution for SNPE’s second round. The approach proposed here builds a different first approximation, consisting of fitting a normalizing flow to the samples produced by some ABC method. Overall, I find the technical contribution of the paper limited, though it proposes a clever way of using ABC in concert with SNPE.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Overall, I find the technical contribution of the paper limited, though it proposes a clever way of using ABC in concert with SNPE. A few comments / questions.\n\n- Typo in section 2.2. Very similar paragraphs repeated.\n\n“…utilizing a neural network F and its adjustable network weights ψ, is often used\nas an NPE. In order to train qF (x,ψ)(θ), the following loss is minimized:\n\nA neural density estimator qF (x,φ)(θ) used neural network F with adjustable network weights φ is used for NPE that trained by minimize the loss:...\"\n\n- It would be good to explain the SVAR example more thoroughly. For instance, X is a 6 by 6 matrix, what are the elements of X that are non-zero and exact mapping to theta? It is also not totally clear to me what’s the size of the generated dataset to train NPE. It is also mentioned that the lag-1 autocorrelation is used as summary statistic. It is also not clear to me what the size of this summary statistic is. Since we have the lag 1 autocorrelation for each pair (i, j), are the observations 36-dimensional?\n\n- NPE-based methods do not really need the summary statistic, as they can “learn their own” with the trained network, which may be beneficial. The summary statistic is needed by ABC methods. So I’d think a comparison against NPE that takes in the raw data would be a good addition.\n\n- Overall I find that the empirical evaluation could be extended. Including a few more models would be good. Specially, since NPE has been observed to perform well in a range of settings (for instance, a few examples from the “benchmarking SBI methods” paper could be used). Additionally, I’d be interested in seeing comparisons against plain NPE in all models, where NPE is trained using the same number of model simulations as the total amount used by other methods (adding across rounds).", "broader_impact_concerns": "No", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Lv7jWJlkBA", "reviewer": "Reviewer_s2up", "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\nOverall, the paper is well-written and easy to follow. The algorithms are explained clearly and with sufficient details (especially in the appendix). The two experiments show promising effect that the proposed preconditioning method has on likelihood-free parameter inference.\n\n__Weakness:__\n\n- The method that is proposed, by combining ABC and NPE, is really quite simple and difficult to say whether it is completely novel. However, with TMLR's emphasis on technical correctness over novelty, this is not a major issue.\n- In the experiments, it would be better to have further explorations of the cost-accuracy trade-off. For example by showing e.g. exactly how many iterations of ABC should be run before observing improvements over vanilla (S)NPE and how much cost that would incur. Also, the experiments only show qualitative results on a select few examples, which puts into question the robustness of the method. By using a metric to quantify the resulting posterior estimate when compared to the output of BSL (e.g. using the Wasserstein distance, maximum mean discrepancy, KL divergence, etc...), and displaying the scores with error bars (e.g. arising from different random seeds for MCMC sampling or neural network initialisation), we would have a better understanding of how robust the method is. Also you can quantify the results in Figure 7 rather easily using e.g. negative log-predictive density, and present the scores across multiple random seeds.\n- I am also not sure if \"preconditioning\" is the right phrase to describe the proposed method. To my knowledge, \"preconditioning\" refers specifically to methods for improving the condition number of a matrix in a linear system and I'm not sure how generally this term could be applied.", "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": "In this work, the authors propose a preconditioning method for neural posterior estimators (NPE) that works by (1) running several iterations of approximate Bayesian computation (ABC) to get a rough estimate of the parameters used to generate data, (2) fitting an unconditional density estimator on the ABC output, and (3) using this as the prior in the computation of the NPE. They demonstrate two examples where standard NPE might fail and show that by preconditioning, this can be fixed.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The sentence \"A neural density estimator ... used for NPE that trained by minimise the loss\" in page 3 is incoherent and should be fixed.\n- The acronym \"BSL\" in page 4 and thereafter is not introduced.\n- The acronym \"TSNPDE\" in Section 4 is only introduced briefly in an earlier section, which is easy to miss. It might be worth reintroducing what the \"T\" in TSNPDE stands for again in Section 4 just to be clear.\n- In some places, the authors mention difficulties due to \"leakage issues\". What does this mean exactly?", "broader_impact_concerns": "Not applicable.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper unfortunately could use considerably more polish. Please address the following in the revised manuscript:
* Proofread the paper for typos and correct them, there is a fair number of them.
 * Correct "neural posterior estimator" to "neural posterior estimation" in the abstract.
 * The paragraph describing the leakage problem in Section 2.2 is very unclear. Please rewrite it to explain clearly what exactly is being leaked and why, avoiding vague phrases like "untrained areas of neural network".
* The description of the transition matrix in Section 2.2.1 is either incomplete or incorrect, as the claim that there is one non-zero off-diagonal element per column does not follow from conditions (1) and (2). Perhaps there's an unstated assumption that each variable is coupled to only one other variable? Also, condition (1) seems redundant given the sentence preceding it.
 * Make it clear how many rounds were used for (T)SNPE for the BVCBM.
 * Clarify that the claim about q recovering the true posterior after Eq. 2 is true only in the limit of infinite data.
 * Make the description of the SMC ABC algorithm more precise and explain what parameter c is supposed to be.
Decision: Accept (poster)