Pseudo-Likelihood Inference
Authors: Theo Gruner
Boris Belousov
Fabio Muratore
Daniel Palenicek
Jan Peters
Venue: NeurIPS 2023 poster (2023)
Keywords: simulation-based inference
approximate Bayesian computation

Abstract: Simulation-Based Inference (SBI) is a common name for an emerging family of approaches that infer the model parameters when the likelihood is intractable. Existing SBI methods either approximate the likelihood, such as Approximate Bayesian Computation (ABC) or directly model the posterior, such as Sequential Neural Posterior Estimation (SNPE). While ABC is efficient on low-dimensional problems, on higher-dimensional tasks, it is generally outperformed by SNPE, which leverages function approximation. In this paper, we propose Pseudo-Likelihood Inference (PLI), a new method that brings neural approximation into ABC, making it competitive on challenging Bayesian system identification tasks. By utilizing integral probability metrics, we introduce a smooth likelihood kernel with an adaptive bandwidth that is updated based on information-theoretic trust regions. Thanks to this formulation, our method (i) allows for optimizing neural posteriors via gradient descent, (ii) does not rely on summary statistics, and (iii) enables multiple observations as input. In comparison to SNPE, it leads to improved performance when more data is available. The effectiveness of PLI is evaluated on four classical SBI benchmark tasks and on a highly dynamic physical system, showing particular advantages on stochastic simulations and multi-modal posterior landscapes.

Pseudo-Likelihood Inference
Theo Gruner 1,2
Boris Belousov 3
Fabio Muratore
Daniel Palenicek 1,2 Jan Peters 1,2,3,5
Intelligent Autonomous Systems Group, Technical University of Darmstadt Aessian.Al 3 Germmnn Research Center for AI (DFKI)
Bosch Center for Artificial Intelligence Centre for Cognitive Science
theo sunoosssnoo. caanerrtu-ddrmstadde
Abstract
Simulation-Bann -assed Inference (SBI) is common name for emerging family of approaches that infer the model parameters when the likelihood is intractable. Exist- ing SBI methods either approximate the likelihood, such Approximate Bayesian Computation (ABC), or directly model the posterior, such as Sequential Neural Posterior Estimation (SNPE). While ABC efficient on low-dimensional prob- lems, on nigher-dimensionn
tasks, it is generally outperformed by SNPE which leverages function approximation. In this paper, we propose Pseudo-Likelihen Inference (PLI), a new method that brings neural approximation into ABC, making it competitive on challenging Bayesian system identification tasks. By utilizing integral probability metrics, we introduce a smooth likelihood kernel with an adap- tive bandwidth that is updated based on information-theoreiic trust regions. Thanks to this formulation, our method (i) allows for optimizing neural posteriors via gradient descent, (ii) does not rely on summary statistics, and (iii) enables multiple observations input. In comparison to SNPE, it leads to improved performance when more data is available. The effectiveness of PLI evaluated on four classical SBI benchmark tasks and highly dynamic physical system, showing particular advantages on stochastic simulations and multi-moda modal posterior landscapes.
Introduction
Parametric stochastic simulators are a well-established tool for predicting the behavior of real- world phenomena. These statistical models find widespread use in various scientific fields such as physics, economics, biology, ecology, computer science, and robotics, where they help to gain knowledge about the underlying stochastic processes [27, 34] or generate additional data for subsequent downstream tasks [38]. In both cases, the practitioner seeks to explain the observations as accurately as possible while incorporating all available information. The output of such simulator is largely determined by its parameters and their values. When estimating these parameters using Bayesian inference, given observations from physical system, which inevitably subject to measurement noise, obtain distribution over values instead of a point estimate. Additionally, there might be several parameter configurations yielding the same observation, hence rendering the resulting distribution to be multi- modal. Moreover, the likelihood function might be unknown or too expensive to evaluate for many practical use cases. The combination of these difficulties makes obtaining posterior distribution over simulator parameters challenging for state-of-the-an inference methods, both regarding effectiveness and efficiency.
SBI approaches address the issue of intractable likelihoods by using (stochastic) simulators as forward models to generate observations from proposal distributions over parameters. The approaches also often called likelihood-free which can be easily misunderstood since some of them directly approximate the likelihood [7]. ABC is a family of SBI methods that approximate the posterior with
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

a set of weighted particles which are obtained from Monte Carlo simulations and updated based on an empirical estimation of the intractable likelihood [48]. For an ABC upproach work well, three criteria have to be fulfilled: (i) the likelihood kernel is capable of measuring the similarity of observations meaningfully, (ii) the proposal distribution samples close to the posterior, (iii) the decision-makinn rule balances between accepting a sufficient amount of samples from the proposal and steering inference towards the posterior distribution. Constructing a suitable likelihood kernel often means tailoring summary statistics to the problem at hand. However, recent advances promise to replace this heuristic-based process by employing Integral Probability Metrics (IPMs) to measure statistical distances in observation space [3, 16, 17]. While these methods significantly increase the required number of simulations, approximations of the statistical distances [25, 8] can be computed in parallel, hence facilitating the parallelization of the whole inference pipeline. Following up on the shortcomings of ABC, the family of SNPE approaches provide Bayesian inference methods that leverage conditional neural density estimators approximate the posterior [40, 24, 41, 19]. The benchmarking study of Lueckmann et al. [32] concludes that, generally, SNPE approaches are to be preferred over ABC as they are superior in terms of expressibility and accuracy across wide range of benchmarking tasks. However, it is important to point out that the analysis of the posterior inference has solely been reported for single observations. These single-sample scenarios favor SNPE in high dimensions since ABC relies on summary statistics to evaluate the likelihood. Therefore, it remains an open question whether SNPE methods can transfer their benefits settings where the (approximated) posterior is conditioned on multiple observations at once.
Contributions. By deriving the ABC posterior from constrained variational inference objec- tive, we introduce a novel SBI method called Pseudo-Likelihood Inference (PLI). Our formulation allows for approximating the posterior with neural density estimators. PLI updates this posterior from pseudo-likelihoods which are exponentially transformed statistical distances computed using IPMs. To further remove heuristics from the inference process, we derive an adaptive bandwidth update of PLI's likelihood kernel that bounds the loss of information based on information-geomedt trust-region principles. This way, PLI can update its neural posterior solely given observations from (stochastic) black-box simulator. Moreover, the usage of IPMs enables PLI to simultane- ously condition on a variable number of observations, while SNPE methods need concatenate them and, therefore, degrade when the number of observations increases. We compare PLI against ABC and one SNPE method on two SBI benchmarking tasks as well as highly dynamical double pendulum task. For both, ABC and PLI, we investigate two IPMs: the Maximum Mean Discrep- ancy (MMD) and the Wasserstein distance. Motivated by recent benchmarking results, we chose Automatic Posterior Transformation (APT) [24] to represent SNPE approaches. Our experiments investigate the dependency of the trained density estimator's performance on the number of observa- tions, where performance is measured in observation as well as in parameter space. We show the merits and disadvantages of all methods and conclude with concrete recommendations.
2 Bayesian inference with intractable likelihoods
The objective of Bayesian inference is to find the posterior parameter distribution p(E(x].v) given a set of reference data points xi.N which are assumed to be drawn from the likelihood model (x(()). Given prior belief over the parameters p(E), the posterior is expressed via Bayes' rule
xig]xi.)) p(xi.N {) p(E).
(1)
In the following, we describe SBI methodologies that aim to approximate the posterior (1) when the likelihood is given by a simulator model. from which only sampling is possible, X1:M p(x({) but evaluating the likelihood is infeasible.
2.1 Approximate Bayesian computation
ABC methods perform Bayesian inference without explicitly computing the likelihood function (( = p(x XI:M, E ()p(x):a1(5) dx1.Ms [48]. Instead, they approximate it
(x( x / Kg(D(x x1:M)) p(X {)
using Monte Carlo samples XI:M from the simulator as the reference points and smoothening them with a kernel K (D(xi. x1:M) [28]. The kernel assesses the similarity of the reference data

---

Xi.N and the simulated data X1:M based on distance measure D, the kernel type K, and the kernel bandwidth 3. The uniform kernel 11 (D(x; xx. (see Table A.1) has emerged as the default kernel choice of many ABC methods [49, 15, 30, 31J. In this case, the bandwidth represents rejection threshold that assigns zero probability to all parameters whose simulations lie outside of the 3-ball in terms of the distance D. The uniform kernel exhibits the favorable characteristic of converging to the likelihood in the limit [28]
f P (X (E) J (xt.w) (X1:M (p(x:.///)dxx.... p(xi,w
Once the approximate likelihood (2) is obtained, ABC draws samples from the approximate posterior
(f(xx))))) ( () P(E).
(4)
There exist multiple ways of implementing this sampling procedure. rejection ABC [49], the simplest form, proposal parameters are drawn from the prior distribution p(g) and are accepted if the simulated data falls close to the true data, as measured by the kernel function. While rejection ABC yields a simple algorithm with desirable convergence properties, finding posterior samples for small bandwidths 3 in high dimensions often becomes computationally infeasible [33]. Therefore, the research on ABC focuses on three directions of improvement: (i) replacing the prior with sequentially updated proposal distribution (E( to reduce the search space during sampling, (ii) adapting the bandwidth 3 to draw samples with appropriate acceptance rate, and (iii) finding sufficient statistics to represent the simulated output in low dimensions [48]. MCMC-ABC [33] and SMC-ABC [47, 51, 15, 31] build upon sampling strategies based on Markov Chain Monte Carlo (MCMC) and Sequential Monte Carlo (SMC) to sequentially update the proposal distribution. MCMC-ABC does not allow for adaptive bandwidth, and thus, SMC sampling strategies have evolved as the leading ABC methods for these cases [15].
2.2 Sequential Monte Carlo ABC
SMC-ABC builds on SMC samplers introduced by Del Moral et al. [14]. Fundamentally, SMC- ABC approximates the posterior distribution through a sequence of intermediate target posterior distributions ps, (E)xi.N) (4) that are characterized by adaptable bandwidth parameter 3r, where t denotes the inference time. Furthermore, SMC-ABC uses importance sampling from a sequentially updated proposal distribution T (E) to improve the sample efficiency. The proposal distribution is M
represented by an empircal distribution Fe 1/M Mm 8g (E) that defined through a set
of particles {E(") Importance sampling then enables the approximation of the target posterior pa, (E]xi.N) from the proposal distribution
po, (E. xx..N) (E(1))
Pa. (E)xTON) = q: (E) = Σ W
(E); c
(E):
w!(0)
(5)
Here, We are the weights between the target posterior and the proposal distribution. SMC-ABC methods follow three steps carry out inference for the next target posterior ( (f)xt.)): (i)
A new bandwidth 3441 of the target posterior ps,., (x).vlE) is estimated. Typically, the update based on heuristics, such as the Effective Sample Size (ESS) [15] to ensure that the particle
variance does not degrade. (ii) New proposal particles <(-) are sampled from a forward Markov kernel
(1) (e < ( stay close to the target posterior of the next iteration 1 (glxi.v). ). (iii) The weights of the particles are adjusted based on approximations of (5). As the weight update is typically numerically intractable [14], different SMC-ABC methods [15, 30, 31] have been introduced which propose approximations to the optimal weight update. We refer to Appendix B for a more detailed explanation of SMC-ABC and its different approaches.
3 Pseudo-likellhood inference
The proposed PLI methodology, summarized in Figure 1, generalizes the ABC approaches by introducing exponential likelihood kernels with adaptive bandwidth updates, which are motivated from Variational Inference (VI) perspective.
3

---

prior r(E)
(x))(x))
((1)
D² (x)
100
((4)
proposal TT-1
p(x(e())
0.1
((3)
p(x(e(2)
simulate p(x(<)
0.0
0.5
1.0
p(x({(0))
p(x (()))
data xi.N
Prior
KL (Fe [ SE
Posterior
pseudo-
likelihood
E(E)
Tt (EIx)..)
<x(Elx(,n)
p(E]xi.N)
posterior
Figure 1: Schematic overview of the introduced iterative Pseudo- Secudo-Likelihood Inference (PLI) approach. (top) Based samples drawn from the simulator, the hikeldo-i (8) is evaluated based the discrepancy hetween the empirical data-generating and likelihood distributions. The har chart shows how the pseudo-likelihood evaluation changes for different bandwidths. (bottom) The evaluation of the pseudo-likelihood is used to estimate target posterior тy under trust-region constraints that from the prior distribution the final posterior.
3.1 Exponential likelihood kernels
PLI adopts the view of SMC-ABC on approximating a smoothed target posterior Pa (()x).)))) in the following denoted by Fe (E), by formulating the following constrained VI problem for each inference step I
me (<)
argmin KL (T(E) p( xi.))), (())E(())
s.t. KL (T(E) || (E) S.
The optimization is balanced between fitting the posterior distribution p(E)xi.v) and constraining the information loss between two inference steps - and (E) The loss of information incorporated as trust-region constraint with the bound 0 in the space of probability distributions
PE{) through the Kullback-Leibhe (KI.) divergence KT. (=((E) || (t-1(E)).
Theorem 1. The optimal target distribution Te(E) in the optimization problem (6) is given by
P(E)
Tx(E)
p(x
T+ TT (E)
TT)
where > 0 is a dual Lagrangian variable corresponding to the trust-region constraint.
Proof. See Appendix A.I.
The temperature parameter V/r plays the role of an adaptive step size that controls the update step from (g) to (E). In the limit of small step sizes n) - 0 at convergence, the target posterior (7) turns into the true posterior TE(E) + p())p(x (E). In the spirit of ABC-based methods, we
approximate the intractable likelihood p(x (E) with Gibbs distribution p(x".xx which call lieudd--ikelihood, and which is based on discrepancy measure between the empirical data distribution p² (x) 1/N E; &x; (x) and the simulator likelihood p(x(E)
D(p" (x)), ()),(()()))
p(xi.N ( = Z(E)
exp
(8)
28
Here, 7(5) is the normalization constant, and > 0 is a bandwidth parameter that controls the sharpness of the approximation, When the KI. divergence is used as the discrepancy measure D, we recover the true likelihood, the following lemma states.

---

Lemma 1. When D = KL and = 1/N in (8). the ngeudo-likelihon p(x1,N16) equals the true likelihood p(x
Proof. Since (L(p"(x) (x) p(x 4))
[x)p(()) [z)] -N log (()), then provided 23 = L/N, the exponential in (8) is given by exp(- (NNKL(ND (x) ((x()))) x p(xi.N (() and ( a constant. D
= = -Hpp
However, the KL divergence is intractable in SBI because we cannot evaluate the likelihood. There- fore. we propose replacing the KL divergence with Integral Probability Metrics (IPMs), such as the MMD and the Wasserstein distance, which can be evaluated on distribution samples. Although theoretical analysis less straightforward in these cases, some results have been obtained in prior works. Consistency and robustness of an MMD-based posterior estimator were shown by Chérief- Abdellatif and Alquier [5] and Wasserstein- based exponential kernels were studied by De Plaen et al. [12]. In this paper, we focus on a practical instantiation of seeudolikelihove inference, which can accommodate a variety of divergence functions and obtain superior empirical results by leveraging neural posterior approximators and adaptive step-size updates. The following subsections introduce the key components that constitute our method.
3.2 Bandwidth adaptation from trust-region principles
The Lagrangian parameter 7/ has a particularly interesting property. From (7), we see that pulling T into the pseudo- likelhood (2), yields time-dependent tempered pseudo-likelihood
(x (5) == Z(E) exp(- (xx)) p(x())) - -i(xi.N T
Here, we introduce the adaptive bandwidth 84 = (1 + ) that approaches ,3 in the limit 7/ - 0. The dual formulation of the stochastic search problem (6) leads to tractable solution for the optimal
Lagrangian parameter 17s, and hence an optimal bandwidth Bg (see Appendix A.1 for more details).
( (11m1)
g(m)) = - (1 + 17.)) log L
Pe 7x(x
(10)
While we obtain the primal optimal point in closed form (7) to obtain the optimal dual variable 7/t, we need to resort to numer-
0.3
0.3
ical optimization of the Lagrangian dual objective. The optimal bandwidth parameter 3,, that is obtained by maximizing (10), information-boundee
0.7
0.2
1.0
can be seen move the pseudo-likelihun inference stages, the proposal prior (E) (E) is typically uninfor- mative, and thus the information loss is moderate even if px (E) moves far away from (E). In the later inference steps, the proposal distribution is typically pronounced, and small devi- ations may lead to significant information loss. This intuition suggests that the bandwidth 3r should decay over iterations, and indeed Figure 2 shows that 3, quickly decays towards zero over range of values of & on the Gaussian location task (Sec. C.2). The exact decay schedule of 3g is problem-dependent. Therefore, it is convenient to set an information-los bound ε and obtain an adaptive bandwidth schedule by optimizing the dual (10) rather than pre-specifying a decay schedule by hand for each problem.
trust region update to towards the likelihood. In the early
0.1
0.0
5
10
15
20
step
Figure 2: Bandwidth 7/2 recorded for different s on the Gaussian loca- tion task. The bandwidth is mono- tonically decreasing over iterations.
3.3 Approximate Bayesian inference with pseudo likeliods
Inference (PLI) is a sequential SBI methodology based on approximating the target posterior Te (E) (7). It is closely tied to SMC-ABC by (i) sequentially approximating the target posteriors, (ii) sequentially adapting the bandwidth parameter, and (iii) sequentially updating the proposal distribution for higher sample efficiency. Instead of representing the posterior through set of weighted particles, the PLI formulation allows for various powerful neural density estimators.
A parameterized density model qs (E) is trained to approximate the PLI posterior (7) using the m- projection mingeo (L())(E) || 4 (E)). which results the Weighted Maximum Likelihood (WML)
objective with parameter samples drawn from the proposal prior (b)) ~ (g)1(((())
max Ek_ (k) log qs (E(())))
(c())) (
m(e())) (E
(e(xt.N15(k))
(1))
5

---

Thus, we derive a practical PLI algorithm by leveraging this empirical WML objective. Fur-
Algorithm Pseudo-Likelihood Inference (PLI)
I: input: reference data X prior p(E), stochastic simulator p(x[{), ), IPM D(.-))
ther details on the objective derivation are de- scribed in Appendix A.2. Our proposed PLII AI- gorithm 1 consists of four main steps. First, in lines 4-8, training pairs from the proposal and simulator are drawn, and the discrepancy measure D(x XI:M) between the observations and the simulations is evaluated for each g(*), We follow Gretton al. [25] approximate the MMD be- tween two discrete probability measures, whereas make of the entropy regularized optimal transport formulation to approximate the Wasser- stein distance [8] (see Appendix C.1). Both ver- sions facilitate parallelization on the GPU. Sec-
posterior approximator Is (E), max. iter. T
initialize proposal prior (E( = P(E)
3: for t in 1: do
sample parameters ((l:k) ot (E)
5: for each ( do
6: simulate data x XIMM ~ (k)
(R)
7: compute IPM s(k) D(x xi:N)
(k)
8:
end for
9:
update th by maximizing the dual (10)
10:
evaluate fr (x1... ((()) (9)
11:
fit qa. (E) by WML (11)
ond, in line 9, the optimal bandwidth under trust
12:
set new proposal prior TT(E) qd. (E)
region constraint is estimated, and the tempered
13:
end for
greudo-likelhoodd is evaluated. by maximizing
14: output: approximate posterior qør (E)
the dual (10). Third, in line 11, the parameterized density estimator (e(E) (E) trained to approximate the target posterior (E) via the m-projection (11). Note that the expectation w.r.t. the proposal distribution (5) enables gradient descent on the qe, estimator without requiring a differentiable simulator. Fourth, in line 12, we set the current posterior approximation go. (E) as the proposal TT(E) for the next inference step, thus leveraging bootstrapping of the density estimator. While we restrict the analysis this paper to the m-projection, we note that the i-projection can also be employed, as shown in Appendix A.3.
Normalization. The normalization term ZEE) in the definition of the pscudo-likelihood (8) requires taking an integral over the reference data xi which infeasible in practice. When the KL divergence is used in the kernel, Z (E) does not depend on &. as shown in Lemma 1. While in general, the dependence on & cannot be neglected, its influence on the weights in (10) and (11) may be negligible, provided the relative ranking of the samples is not affected significantly. In Appendix A.4, we provide an ablation study on low-dimensional problems where the integral over xi.N can be approximated by sampling. We observed that, even though the ranking correlation of the weights (k) in (1I) is different with and without estimating Z(E), the final posterior is not affected. Therefore, in the subsequent experiments, we treat it as a constant, as in the ideal case of the KL divergence. Nevertheless, this is a point where our practical implementation does not follow the theoretical derivation strictly, and this issue should be addressed in future work.
4 Experiments
We compare the PLII framework against SMC-ABC [15] and APT [24] on five diverse tasks. Our implementation is based on WBs [3], but instead of the employed r-hit kernel [30], our implementation is based on population Monte Carlo (Alogrithm 3 [31]) because we observed improved performance in preliminary studies. A summary of the different ABC methods is given in Appendix B and Table A. APT was chosen as the representative for the class of SNPE algorithms. We leverage Neural Spline Flows (NSFs) [19] as density estimators for both PLI and APT. Both neural flow configurations share the same base network architecture, but for APT, the conditional flow is augmented with an embedding network (Appendix C). All experiments are implemented in JAX [4], and each ran on single Nvidia RTX 3090. To make the experiments comparable, the simulation budgets of PLI and APT were fixed 5000 samples per inference step over 20 episodes, while ABC ran for 200 episodes on 1000 particles.
4.1 Evaluation metries
The model is compared against the reference posterior samples &* when available. We also quantify the methods' performances based on their realizations by computing the Wasserstein distance and

---

Gaussian Location
Gaussian Mixture
SLCP
SIR
Funuta
20
MMD
0.0
30DD
03
WA X
108
100
1000
100
1000
N
W-ABC
Figure Evaluation the posterior performance on five different tasks. We report the mean and 95% ci over 10 random seeds, each carried out using N data points for conditioning. We compare samples from the approximate posterior È ~ g((|xiix) against reference posterior samples using MMD and the Wasserstein distance. No posterior samples were available for the Furuta pendulum. Therefore, the performance is evaluated in the observation space. Lower values are better for all metries. PLI is the preferred method for conditioning increasing N. ABC performs better than PLI on Gaussian Location and Gaussian Mixture tasks but lags in more complex tasks. APT excels with few observations but degrades N increases.
multiple observations due to its steady improvement with
the MMD. The comparison carried out on 10000 samples each. Furthermore, we use Posterior Predictive Checks (PPCs) evaluate the predictive capabilities of the posterior models in the observation space [Eq(g(xt [D(x X1:as)] ]- Due to the computational limits, the PPCs are carried out on 1000 simulations against the reference data. Lueckmann et al. [32] also report results with classifier-basee tests and the kernelized Stein-discrepancy. However, since benchmarking is not our
focus, we restrict the analysis to comparing with the Wasserstein and MMD.
4.2 Tasks
We evaluate PLI on four common benchmarking tasks within the SBI community [32]] Gaussian Location, a Gaussian Mixture Model, system identification task Furuta pendulum representing highly dynamic continuous control system. The tasks' specifications are listed in Appendix For each task, we conduct experiments for different numbers of available reference observations N 2. 5. 10, 20, 50, 100. 1000}. The reference observations are simulated based on a pre-defined ground-trum parameter t#² Although PLI and ABC cope with varying numbers of observations N and numbers of simulations per parameter M, we choose = M for all experiments since it is required by APT. In the following paragraphs, we first discuss the results of the benchmarking tasks and present a separate discussion for the Furuta pendulum. Figure 3 gives quantitative overview of the benchmarking tasks compared to the reference posteriors, while Figure C.I in the Appendix complements the study by showing the posterior predictive performances.
Complex-Postere Posteeror,, and SIR. Further, we add
Benchmarking tasks. Each benchmarking task presents different challenges that must be addressed by the SBI methods. In Gaussian Location, the task is to infer a uni-modal 10-dimensional Gaussian distribution. Gaussian Mixture Model and SLCP feature multi-modal posteriors that require flexible density estimators. SIR is a well-known cpidemiological model that features 10-dimensional data of the dynamical system. All methods generally depict a reoccurring behavior on the different benchmarking tasks, as shown in Figure 3. For fewer observations (N S 20), APT matches the reference posterior better than the other approaches, whereas ABC and PLI match the posterior data better with increasing N. In particular, PLI consistently improves with an increasing number of reference samples. The influence of on shape the posterior is further visualized Figure C.2, which compares the posterior approximations of all methods for N = 2 and N 100 reference observations. the SL.CP task, ABC struggles to capture the multi-modality of the SL.P task. This effect is further illustrated when comparing Figures C.3 and C.4, which allow for qualitative
7

---

Step
Step 1
Step
MMD-PLI
W-PLI
MMD-ABC
W-ABC
0.12
APT
9
10
11
10
11
10
1.1
e
<00a
sin 0,
0.50
OPF
0.0
0.25
0.00
0.5
0.0
0.5 1.0
10
20
50
100
1000
N
Figure 4: Empirical and quantitative evaluation on the Furuta pendulum. (top) Snapshots of the posterior evolution on the my plane on the Furuta pendulum for N = 1000. (bottom) Predictive performance of the learned MMÍ-PPI posterior for the angular rotation sin f,. The stochasticity of the simulator is removed by synchronizing the initial state between the reference and predicted simulations. Thus, the only discrepancies between trajectories are due to the model not capturing the dynamics parameters of the system. After the inference has been completed (step 20), the predictive simulator (- MMD-PL.) can completely recover the ground truth dynamics (- Reference). (right) Evaluation of the mean accumulated error over 1000 trajectories with synchronized initial states between the simulation and the reference trajectory. All approaches improve with rising N while PLI with MMD matches the reference data best.
assessment of the posterior approximations. On SIR, PLI variants show significantly improved performance compared to the other baselines. Generally, we find that ABC and PLI perform better with MMD than with Wasserstein distance. This observation can be attributed to the Wasserstein distance not scaling well to high dimensions, which has been reported by recent studies [17, 16].
Furuta pendulum. The Furuta pendulum is an inverted double pendulum setup [21]. While the system's dynamics are inherently deterministie, small perturbations of the initial state around its unstable equilibrium point lead to highly diverse trajectories. The observation space is x dimensional, where T is the number of time steps per trajectory. We set the sampling frequency of the simulation to 100 Hz and the duration to 1 sec, resulting in 600-dimensional observations. No reference posterior is available for this task; thus, the analysis is restricted to quantifying the observed data. Given the similarity of the Wasserstein distance and the MMD in parameter and observation space the previous tasks, argue that comparison based PPCs, i.e., W ( X1:M) and MMD? (xi.N XI:M), is sufficient. In Figure 3, the Wasserstein PPC and MMD PPC exhibit divergent behaviors, with the former indicating enhanced performance as reference observations increase, in contrast the moderate improvement suggested by MMD. Therefore, we evaluate the models' predictive performances on the deterministic system by synchronizing the initial states of the reference data and the simulations in Figure This modification ensures that the similarity between two rollouts can be evaluated by the accumulated eтror err Σ, | zz/. While all approaches perform better with more reference observations, MMD-PLJ matches the reference dynamics best. Additionally, MMD is favored here, as MMDD-based PI.I and ABC outperform their Wasserstein counterparts, with both the MMD PPC plot and error plot showeasing congruent trends for N 2 20. The appended posterior plots in Figures C.5 and C.6 reveal that for N 2, all methods are widely spread over the prior region, yet converge to the ground truth. However, APT cannot recover the ground truth for N > 100, whereas PLI and ABC center around the ground truth.
Related work
In the previous sections, we have seen that PI.I algorithmically similar to ABC methods with SMC samplers. Therefore, approximating the likelihood by the empirical pscudo-likelihood (2) enables

---

drawing from the rich toolbox of existing approximate inference algorithms. This section introduces related research fields and shows how PLI fits among them.
Sequential neural density estimation. With the enriched class of neural density estimators, amor- tized SBI methods have received increasing interest in recent years. Similar to ABC, synthetic used to approximate the posterior. Sequential neural density estima- tion methods can be further classified into methods that directly train a posterior estimator [40, 24], a neural likelihood [41, 22]. or a neural ratio estimator [19, 36]. All methods have in common that they do not rely on an approximation of the posterior model but are optimized solely on pairs of parameter samples from a proposal distribution g(,) ~ pc(E) and its corresponding simulation x()) ~ (+)) We note that the original papers have only reported posteriors conditioned on a single observation x. While technically, these methods can incorporate multiple data points, this requires either stacking multiple observations or falling back to summary statistics. As noted in [50], neural likelihood estimators
samples from the simulator
sidestep these requirements by evaluating the log-likelihood of single observations and carrying out MCMC sampling on the joint log- log-likelihood. Yet, leveraging the neural likelihood restricts the evaluation of the posterior.
Summary statisties. ABC has commonly relied reducing the dimensionality of the raw ob- servations with summary statistics [48]. These summarics must be carefully chosen and are often task-specific, restricting the general applicability of ABC. Recent additions to ABC methods report on replacing summary statistics with statistical distances [17]. While direct comparison of the raw data suffers from the curse of dimensionality, comparing the observations through empirical measures sidesteps this issue [17]. Bernton et al. [3] report on augmenting the likelihood kernel with the Wasserstein distance, while Park et al. [42] leverage the kernelized approximation of the MMD [25]. Other contributions include the Cramér-von-
distance [20] and the energy distance [39]. While statistical distances are appealing due their general applicability, Drovandi and Frazier [17] con- clude that they are limited by their high computational requirements. Approaches proposed for automated summary design include ABC with indirect inference, which utilizes an auxiliary model to evaluate data summaries [23, 18].
Particle mirror descent. A posterior updating similar to ours (7) has been derived in Particle Mirror Descent (PMD) [10]. PMD tackles particle depletion by incorporating the proposal distribution of the previous round into the optimization process. Furthermore, the authors show that the proposed method converges the posterior given m posterior samples by O(1/ /m) Our version can be seen as extending their approach to the case of intractable likelihoods. We extend PMD to neural density estimators using samples from the proposal posterior (7) as a training set.
Geometric path and likelihood tempering. Rewriting the optimal posterior (7) reveals close relation of the optimal PLI posterior (4) and the geometric path formulation [6, p. 335], TE(E) x >|=>(6 (<) p(xi.v:
The optimal posterior moves from the proposal distribution (E) at inference time t to the target posterior TH T<(E)x].N) x p(x),N.& { along the geometric path that is parameterized by A. The formulation differentiates from likelihood tempering in SMC samplers [6] by leveraging the proposal instead of the prior distribution. Note, however, that for = O, the proposal mimics the prior, and thus the PLI geometric path has the same boundary values as in classical likelihood tempering. While the tempered posterior cannot be applied to SMC samplers due to its dependence on the proposal distribution, the geometric path formulation based on the prior distribution gives rise to sequential annealing ABC [1].
Generalized variational inference. Introduced by Knoblauch et al. [29], Generalized Variational Inference (GVI) is an extension of the standard variational inference framework that starts from the optimization view of Bayes' rule and generalizes it by considering different losses, divergences, and variational families. Therefore, various Bayesian inference methods can be seen as instantiations of GVI with different choices of these three parameters. In particular, ABC and PLI can be seen as KL,,), KL, P(O)) since they employ likelihood kernel as a loss. Crucially, specifying a different loss function instead of the typical log- ligeliieldoood can be shown to add robustness against model misspecification [29]. PAC-Bayes [26] can be seen as generalization of ABC as covers the whole space of possible loss functions 1 1(xi.v= E). The PAC-Bayesian theory provides a broad array of risk bounds for generalized Bayesian learning methods.
GVI with the choices P(K g ( (XI.N. X1:M).

---

Conclusion
We propose Pseudo-Likelihon Inference (PLI), a new addition to the toolbox of SBI methods. PLI targeted for Bayesian inference tasks in which the posterior is conditioned multiple observations simultaneously. For that, we derive a softened ABC posterior from a constrained variational inference problem and leverage IPMs between the empirical observations to assess the intractable likelihood. The derived posterior formulation enables the learning of flexible neural density estimators from black-box simulators, extending the range of applicability for ABC methods. Our experiments assess how well PLI, ABC, and SNPE perform based on their generative power when given varying amounts of reference observations as a condition. Given few observations, SNPE-based methods perform better than ABC and PLI, which rely on statistical distances. However, when more data is available, ABC and PLI methods perform better. When the posterior distribution is simple, ABC is efficient at reproducing it using fast particle updates. However, PLI is better option for complex posterior distributions because of its more adaptable neural density estimator. Additionally, PLI evaluates the posterior probability, which is useful in downstream tasks that require uncertainty quantification.
Limitations. In PLI, the computational cost is distributed among three main computations: the simulation, the summary statisties estimation, and the normalizing flow training. While the com- putational effort is large for every computation, PLI leverages parallelization on GPUs for all three computations. In the provided experiments, the training process of the neural model takes the main computational budget, while simulation and summary statistics are negligible. On Nvidia RTX 3090, ABC typically runs 2-10 min, while PIII and SNPE take 60-90 min, depending on the task. Simpler models however, such as multivariate Gaussian or GMM, reduce the computation time to the simulation. Furthermore, we would like to express that high-fidelity simulators might increase the simulation time significantly, making the simulation the most costly operation within the inference pipeline. Instead utilizing IPMs, one could exploit adversarial strategies to approximate the KL divergence, explored by Mescheder et al. [35] and Santana and Hernández-Lobato [45].
Acknowledgements
This work was funded by the Hessian Ministry of Science and the Arts (HMWK) through the projects "The Third Wave of Artificial Intelligence- 3AI", hessian.Al, and the grant "Einrichtung eines Labors des Deutschen Forschungszentrum für Künstliche Intelligenz (DFKI) an der Technischen Universität Darmstadt". Fabio Muratore was employed by the Robert Bosch GmbH during the time of this collaboration.
References
[1] Carlo Albert, Hans R Künsch, and Andreas Scheidegger. A simulated annealing approach to approximate bayes computations. Statistics and Computing, 2516):12 2015.
[2] Mark A. Beaumont. Approximate bayesian computation in evolution and ecology. Annual Review of Ecology, Evolution, and Systematics, 41(1):379-406 2010.
[3] Espen Bernton, Pierre E. Jacob, Mathieu Gerber, and Christian P. Robert. Approximate bayesian computation with the wasserstein distance. Journal of the Royal Statistical Society: Series (Statistical Methodology), 81 81(2):335-26, 2019.
[4] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Minn and Qiao Zhang. JAX: composable transformations of Python+NumPP programs, 2018.
[5] Madr-Eddine Chérief-Abdellatn and Pierre Alquier. band-bayes: Robust bayesian estimation via maximum mean discrepancy. In Symposium on Advances in Approximate Bayesian Inference, 2019.
[6] Nicolas Chopin, Omiros Papaspiliopoulos, et al. An introduction to sequential Monte Carlo, volume 4. Springer, 2020.
10

---

[7] Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117(43):30055-30062, 2020.
[8] Marco Cuturi. Sinkhorn distances: Lightspeed computation optimal transport. In Advances in Neural Information Processing Systems, 2013.
[9] Marco Cuturi, Laetitia Peng-Paapxxnsths, Yingtao Tian, Charlotte Bunne, Geoff Davis, and Olivier Teboul. Optimal transport tools (ott): A jax toolbox for all things wasserstein. arXiv preprint arXiv:2201.12333 2022.
[10] Bo Dai, Niao He, Hanjun Dai, and Le Song. Provable bayesian inference via particle mirror descent. In Artificial Intelligence and Statistics, 2016.
[11] Christian Daniel, Gerhard Neumann, Oliver Kroemer, and Jan Peters. Hierarchical relative entropy policy search. Journal of Machine Learning Research, 17:93:1-93:500 2016.
[12] Henri De Plaen, Michaël Fanuel, and Johan AK Suykens. Wasserstein exponential kemels. In 2020 International Joint Conference on Neural Networks. IEEE, 2020.
[13] Mare Peter Deisenroth, Gerhard Neumann. Jan Peters, et al. A survey policy search for robotics. Foundations and Trends in Robotics, 2(1-2):1-142, 2013.
[14] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. Sequential monte carlo samplers. Journal of the Royal Statistical Sociery. Series B (Statistical Methodology), 68(3):411 436,, 2006.
[15] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. An adaptive sequential monte carlo method for approximate bayesian computation. Statistics and Computing, 1005):000-1022 2012.
[16] Charita Dellaporta, Jeremias Knoblauch, Theodoros Damoulas, and Srançois-Xavier Briol. Robust bayesian inference for simulator-based models via the MMD posterior bootstrap. In Artificial Intelligence and Statistics, 2022.
[17] Christopher C. Drovandi and David T. Frazier. comparison of likelihood free methods with and without summary statistics. Statistics and Computing, 32(3):42, 2022.
[18] Christopher C Drovandi, Anthony N Pettitt, and Anthony Lec. Bayesian indirect inference using a parametric auxiliary model. Statistical Science, 30(1):72-95, 2015.
[19] Conor Durkan, lain Murray, and George Papamakarios. On contrastive learning for likelihood- free inference. International Conference Machine Learning, 2020.
[20] David T Frazier. Robust and efficient approximate bayesian computation: A minimum distance approach. arXiv preprint arXiv:2006. 14126, 2020.
[21] Katsuhisa Furuta, M Yamakita, and S Kobayashi. Swing-up control of inverted pendulum using pseudo-state feedback. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 206(4):263-269, 1992.
[22] Pierre Glaser, Michael Arbel, Arnaud Doucet, and Arthur Gretton. Maximum likelihood learning of energy-basee models for simulation-based inference. arXiv preprint arXiv:2210. 11f::22200.4766, 2022.
[23] Alexander Gleim and Christian Pigorsch. Approximate bayesian computation with indirect summary statistics. Draft paper: hutp://lett-igoosch. mee. uni-bonn. de/datadreseari/ppperss 2013.
[24] David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior trans- formation for likelihood-free inference. In International Conference on Machine Learning, 2019.
[25) Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. kernel two-sample test. Journal of Machine Learning Research, 13(1):723-773, 2012.
[26] Benjamin Guedj. A primer pac-bayesian learning, arXiv preprint
2019.

---

[27] Florian Hartig, Justin M Calabrese, Björn Reineking, Thorsten Wiegand, and Andreas Huth. Statistical inference for stochastic simulation mooersy and application. Ecology letters, 14(8):816 11--877, 2011.
[28] George Karabatsos and Fabrizio Leisen. An approximate likelihood perspective ABC methods. Statistics Surveys, 12:66 104, 2018.
[29] Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas. An optimization-cen view on bayes' rule: Reviewing and generalizing variational inference. Journal of Machine Learning Research, 23(132)11-109 2022.
[30] Anthony Lee. On the choice of MCMC kernels for approximate bayesian computation with SMC samplers. In Winter Simulation Conference, 2012.
[31] Maxime Lenormand, Franck Jabot, and Guillaume Deffuant. Adaptive approximate bayesian computation for complex models. Computational Statistics, 28()):2777-2766, 2013.
[32] Jan-Man Van-Matthis Lucckmann, Jan Boelts, David S. Greenberg, Pedro J. Gonçalves, and Jakob H. Macke. Benchmarking simulation-based inference. In Artificial Intelligence and Statistics, 2021.
[33] Paul Marjoram, John Molitor, Vincent Plagnol, and Simon Tavaré. Markov chain monte carlo without likelihoods. Proceedings of the National Academy of Sciences, (00(26):1532222 53324-15333 2003.
[34] Kevin McGoff, Sayan Mukherjee, and Natesh Pillai. Statistical inference for dynamical systems: review. Statistics Surveys, 9:209-252, 2015.
[35] Lars M. Mescheder, Sebastian Nowozin, and Andreas Geiger. Adversarial variational bayes: Unifying variational autoencoders and generative adversarial networks. ICML, 2017.
[36] Benjamin Kurt Miller, Christoph Weniger, and Patrick Forré. Contrastive neural ratio estimation. arXiv preprint arXiv:2210.06170, 2022.
[37] Fabio Muratore, Theo Gruner, Florian Wiese, Boris Belousov, Michael Gienger, and Jan Peters. Neural posterior domain randomization. In Conference on Robot Learning, 2021.
[38] Fabio Muratore, Fabio Ramos, Greg Turk, Wenhao Yu, Michael Gienger, and Jan Peters. Robot learning from randomized simulations: A review. Frontiers in Robotics and AI, 9, 2022.
[39] Hien Duy Nguyen, Julyan Arbel, Hongliang Lü, and Florence Forbes. Approximate bayesian computation via the energy statistic. IEEE Access, 1:131683-33168 2020.
[40] George Papamakarios and lain Murray. Fast e-free inference of simulation models with bayesian conditional density estimation. In Advances in Neural Information Processing Systems, 2016.
[41] George Papamakarios, David C. Sterratt, and lain Murray. Sequential neural likelihood: Fast likelihood-free free inference with autoregressive flows. In Artificial Intelligence and Statistics, 2019.
[42] Mijung Park, Wittawat Jitkrittum, and Dino Sejdinovic. K2-ABC: approximate bayesian computation with kernel embeddings. In Artificial Intelligence and Statistics, 2016.
[43] Jan Peters, Katharina Mulling, and Yasemin Altun. Relative entropy policy search. In Conference on Artificial Intelligence, 2010.
[44] Gabriel Peyré and Marco Cuturi. Computational optimal transport. Found. Trends Mach. Learn., 11(-6))3355-677 2019.
[45] Simón Rodríguez Santana and Daniel Hernández-Lobato. Adversarial a-divergence minimiza- tion for bayesian approximate inference. Neurocomputiing 471:260-274, 2022.
[46] Richard Sinkhorn. A Relationship Between Arbitrary Positive Matrices and Doubly Stochastic Matrices. The Annals of Mathematical Statistics, 35(2):876 879, 1964.
12

---

[47] S. A. Sisson, Y. Fan, and Mark M. Tanaka. Sequential monte carlo without likelihoods. Proceedings of the National Academy of Sciences, 104(6):160-1765, 2007.
[48] Scott A Sisson, Yanan Fan, and Mark Beaumont. Handbook of approximate Bayesian computa- lion (Ist ed.). CRC Press, 2018.
[49] Simon Tavaré, David J Balding. R C Griffiths, and Peter Donnelly. Inferring Coalescence Times From DNA Sequence Data. Genetics, 145(2):505-513 1997.
[50] Alvaro Tejero-Cantero, Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pedro J. Gonçalves, David S. Greenberg, and Jakob H. Macke. sbi: A toolkit for simulation- based inference. Journal of Open Source Software, 5(52):2505, 2020.
[51] Tina Toni, David Welch, Natalja Strelkowa, Andreas Ipsen, and Michael P.H. Stumpf. Approxi- mate bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of The Royal Society Interface, 6(31):187-2)) 2009.
13

---

A Algorithmic details: pseudo-likelihood inference
A.1 Deriving the optimal PLI parameter distribution
Theorem 1. The optimal target distribution (E(E) in the optimization problem (6) is given by
P(E)
1+78
Te(E)
p(x),Nl()) 1im TT-1(E))
(12)
IT-1(E))
where 7/: > a dual Lagrangian variable corresponding to the trust-region constraint.
Proof. The solution to the stochastic search problem (6) can be obtained from Lagrangian optimiza- tion. The optimization problem restated here for readability
( = arg min KL (m( || p(()xi,)))
T(E)
s.t. KL (TT(E || TT)1 S E,
/ ((E) de
We decompose the KL objective into two terms by applying Bayes' rule [logp(xi.v [E)] + KL (T)(E) p(E))
KL (T(E) p(Eix),))) x(E)
(13)
The constrained optimization problem (6) can be reformulated with Lagrange multipliers as
((T) (E) log logp(xT.N kE)g
+ / (E) (E)) dg
P(E)
+ (f T(E) log
x x(E)
(T)) de
(ʃ (EE)dE -
+
- T(E)
[- -1og p(xi.((x + log
TT())
log
TT()
+ + dg - me x.
(14)
p(E)
1e)
Here, leveraged the assumption that the likelihood p(x is fixed for all joint distributions, and thus, the joint distributions can be split into the likelihood p(x],x1() and their associated prior/proposal distributions. The gradient of the Lagrangian vanishes for the optimal parameter distribution
ac ŪT
log logglxxiii k)) log TE( p(E)
+ 1 +n log TT(E) TT-1(E)
+ 1 + (15)
x=
Reformulation yields
logp(xi,N 1+n
1+
TE p 1. (E) m,- (E)
exp
1 1+
(16)
= (n) p(E)
log logp(xt::Ve
exp
TH. TT-1(E) (E).
TT-1(E))
1+
The normalization constant Q(n) exp((1 + + 2)/(1 11))
follows by marginalization of (16)
p(E)
+
log ((
]
Q
E
exp
(e)
Tt-1(E)
1
We further obtain the dual of the Lagrangian by reinserting (16) into the Lagrangian (14)
g(n) 7E - (1 log (())).)..
(17)
14

---

A.2 Weighted maximum likelihood optimization (m-projection)
The m-projection of the optimal posterior onto the approximation family results in weighted maximum likelihood formulation
min KL (Tx(E) (E))
(18)
max /
log (E)) (E( de
P(E)
I+m
f(x7.N115 ite x0- log q (E) dg
max
11
p())
E
p(x).N Time log (E)
(19)
(()
Pt-1(E)
The weighting term w is independent of Φ. and as such, the formulation facilitates optimizing neural density estimators with gradient descent. This optimization resolves to weighted maximum likelihood where the weights are obtained from the eseudo-likelihood (8) and the importance weights. The weighted maximum likelihood formulation can be optimized in closed form for linear Gaussian models [43, 13], with Expectation- Expeetation-Mxximiaatooo (EM) using Gaussian Mixture Models (GMMs) [11] with gradient descent done in this paper.
A.3 Optimizing with the i-projection
We can reformulate the minimization problem for the i-projection in the following way
min KL (E) || (e(E))
(20)
(E)
min
log
(21)
1/(1+m))
@ (e(
p(E)
p(xi.N (5)
Fime Te- (E)
(()
1/(1+t)
p(g)
p(x ()) T-me
(E) IE)
Pt-
min KL (go(E) || (e))
E
log
Q
Te-1
The equation above alleviates the issue of back poock-popaggating through the simulator by using importance sampling. The optimization problem is fitting the posterior estimator the proposal p(E) while having a regularization term that forces the distribution to fit the reference data. The temperature parameter 1) can thus be interpreted weighting the regularization term. Small values of 1) put more emphasis on the regularization term, while large values concentrate on containing the information of the proposal. For linear Gaussian models, closed-form expression for the KL term in (21) exists. Therefore, it can be directly optimized using any non-linear optimization method.
Table A.1: Kernels used for the ABC and PLI variants during the experiments in Section 4. An IPM denoted by D plays the role of distance measure between the reference data xi and simulated samples XI;M- Parameter 3r controls the kernel bandwidth (see Section 3.1). Effective Sample Size (ESS) is defined as the inverse of the normalized weight's variance.
Kg. a, (D(x xi.N)) Algorithm
3 estimation
Update 8i
SMC ABC
ESS
arg ming, ESS( (1)
1 (D(MIMM))
aESS (1-1,, (ESS(W)-1,,1-1) Br-1) Bi QD(x (x)
PMC ABC
a-Quantile
exp
D(x 2
PLI
Trust-region
B7 3(1 + arg max. g(Te)). see (10)
15

---

A.4 Analysis of the partition function
We shed some light on the intractable log-partition function ({)) introduced in the pseudo- likelihood (8). The partition function ((E) of (8) is an integral over sample space of x X
/
N.
Z
exp
dxi:N.
(22)
3
To approximate the intractable quantity, we approximate the integral through Monte Carlo simulations with uniform distribution 4
D((XT-N) ((( )
Z(E)
RJ
exp
(23)
B
Gaussian Location
Ganssian Mixture
SLCP
0.5
0.5
0.0
Rade
0.0
-0.0
Bohr
00.5
-0.3
Rabe
10
step
5tcp
step
Furuta
0.5
0.5
-10
-1.0
Figure A.I: Spearman correlation coefficient r(w, u) weights úi and the uncorrected weights wi High values of correspond high correlation of the weight rankings, thus meaning that the weights preserve relative ordering.
[- between the partition corrected
Gaussiin Linear
Gaussian Mixture
SLCP
STR
MMDIE:-5
0.01
0.5
0.002
wig" 6
0.02
0.00
10
3.5
5
301
2.5
1.0
20
30
20
step
MMD-PLL with Z(E)
W-PLI Z(E)
Figure A.2: Training plots comparing PLI trained with and without the partition function Z(E) on N = 100 samples. Top row: MMD between posterior samples and model samples. Bottom row: Wasserstein distance between posterior samples and model samples. On all tasks but SLCP, the inclusion of the partition function does not change the posterior inference.
As we cannot sample over the whole space X, we choose to sample over the 5o interval of the exponential kernel, where x represents the mean of 100 prior simulations. We use 10000 samples from to approximate the partition function and evaluate the Z({) based on 100 samples from the target posterior (E)
We evaluate the influence of the partition function Z(E) on the performance our PIII algorithm by comparing the weights W; and w, as defined in (11) with and without Z(E), respectively. Despite the
16

---

weights being numerically different, we hypothesize that they do not change the relative ordering of the samples and, therefore, do not affect the Weighted Maximum Likelihood (WML) update (11) significantly. Therefore, we employ the Spearman correlation coefficient, nonparametric measure of rank correlation, to capture such dependencies, which we report in Figure A. Contrary to our hypothesis, the weights with and without normalization are not perfectly rank-correlated in four out of five experiments when there more than one observed data point. Potentially, this could lead to a different WML update, but when we evaluate the whole PLI algorithm with and without the partition function (E), we find that its influence is marginal for both MMD-PLI and W-PLI in Figure A.2. Only the SLCP task, an improvement of W-PLI can be seen when using the partition function, whereas MMD-PLI performs even better without it.
We have motivated PLI from VI perspective. As mentioned in the related work section (see Sec. 5), this similar to the General Variational Inference (GVI) approach [29], which also defines the posterior as a solution to an optimization problem, namely
π (E) argmin ON IE [Cx( K (= 1l p()))
(24)
TT))
When the loss is defined as the log-likelihood, Lx(E) log n(x with = L/N, the PLI objective (6) is recovered. However, when we substitute the pseudo-likelihood (8) into the GVI objective, obtain additional term in the loss which is given expectation over the log- partition function E, - [ log Z(E)]. Therefore, PLI can be seen as GVI with an additional loss term, or GVI can be equated with PLI using unnormalized pseudo-likelihon likelihood.
B Algorithmic details: sequential Monte Carlo ABC
The foundations of SMC-ABC have been laid by Del Moral et al. [14], who introduced SMC samplers. These samplers describe an approximate inference routine in which the posterior is approximated through a sequence of intermediate target posteriors. sequence of intermediate posteriors is defined by an adaptive bandwidth 3 of the approximate posterior pa, (Elxi,)) (4). Additionally, the sample efficiency of ABC is improved by replacing the prior as the sampling distribution with a proposal distribution N. (E). The proposal distribution is represented by set of partieles Te (E) 1/M Σ (E and through importance sampling an approximation of the target posterior Pp. (E)xi.N) can be obtained, (see (5)),
the context of ABC, the
M
:.O (S. X(,N) (E(
Pa. (EIxL.) = (E)
= Σ w (E); 1=1
(25)
where ( denote the importance weights. The proposal distribution (E( should ideally stay close to the target posterior Pa, (fixi,)) to improve the sample efficiency. Therefore, the proposal distribution is updated based on a Markov kernel K (Et, E(+1) which the transition probability
from & to ε1+1- The update of the proposal distribution is typically numerically intractable as it requires marginalization, i.e., integration over &, for each inference step 0:
(E++(E++1) = f -K (e)))
(26)
To alleviate the computational burden, Del Moral et al. [14] show that the joint representation of the proposal distribution n, (Eorr) can be efficiently calculated as it only requires solving the product over t transitions
Te(En:) = To) II K4t1 E+41). -=0
(27)
We define the joint proposal distribution as the empirical distribution #x(Ea))
M M- EM 8 8g (Eor) defined by a set of joint particles <0)) Thus, the joint posterior approxi-
17

---

mation of Pa. (i.N) based on the importance weights reads as
Ps. (Coit xT.N)
(E(E)))))
(Cor)):
mm
(28)
Sois
Tel (Su:t)
i=1
w! P... (Co: xi.N)
(c(E))
wi
$: (Er):
(29)
The marginal target posterior approximation (Ex) can be directly recovered from the joint approx- imation qe (foxt) Furthermore, both distributions share their weights which means that only required estimate the weights w; in order approximate the target posteriors pp, (E: x.:)). In general, the probability of the target joint posterior pa, (S ()) is intractable. Therefore, the authors introduce an auxiliary backward Markov kernel Le (E Et) simplify the computation
pa,
t-1
Pg. (Eo:((xxx)) Pa, (Er xi,)) II L. (Er+1,))))
(30)
Assuming that posterior approximation of the target posterior pa, (Er(x).v) is available through the set of weighted particles {(w, (t) , and the particles of the proposal distribution
M (E) M- E= &gte) (E) are updated based kernel transition </21 K(41(()) -541). then
e(t)
the importance weights w/+1 are updated based on the following recursion
Pa, (50a+1)
P8. (5r+1) II Le(E,+1,
W1+1
NT4 (Eor++1)
To(En) II== Ke(Er, Er+1)
(31)
Ps. n (Ex+ Le Er) pe (Ex) II 1-|
1= L1 (E Er)
(32)
Pp. (Er) (Ees Ex+1) (En) I Kr(Er.Ert
Thus, the sequential update is performed by updating the current weights we with the marginal weights W1+1. Up to now, the choice of the backward kernel has been neglected. As it is an auxiliary quantity, several approximations can be made to model the backward kerne.. Del Moral et [14] refer to the optimal backward kernel as the Markov kernel that minimizes the variance of the particles
T (E4) (E,,4(EE,,,H1 E++1)
LFDE (E.+1,
T T (E.+1)
They further show that the optimal backward kernel recovers the marginal weights from (25)
PBeca (E(xi:N) W
opt w
(33)
T+)
In general, the optimal backward kernel is numerically intractable and has led to several other approximations summarized in Table B.2. Depending on choice of approximation, a number of different SMC-ABC methods have evolved, namely the classical SMC-ABC approach by Moral et al. [15], Population Monte Carlo (PMC)-ABC [51, 2, 31], and Metropolis-Hastingg ABC [30]. Please refer to those references as well as Algorithms and for implementation details of these approaches.
18

---

Algorithm 2 Sequential Monte Carlo ABC [15]
Algorithm 3 PMC-ABC [31]
1: input: reference data N. N+ stochastic simulator p(x(E), IPM D(.-), max. iteration count T. forward kernel
prior p(E),
I: input: reference data Xi.N. prior p(E), stochastic simulator p(x(E). IPM D(.,), max. iteration count T, forward kemeel (e+1(E,,<+11) a-Quantile
((+)) resampling threshold V
2: initialize particles $0 ~ p(.)
2: imitialize particles 46 ~ p(-)
3: initialize particle weights wh 1/K
3: initialize particle weights w/ (k) 1/K
4: for in 1:T do
4: store number of best particles K., ak,
5: for each { do
5: for in 1:7 do
6: elect Ka best particles ((k)
6: simulate x (k)
(xm p(xx1xx( ((x((_-)) 1 }}
(k)
1:M
7: sample K Ka proposal particles 2(0)
7: compute IPM 5 D(xt.N (1:1))
(k)
Kr(E (k)) ((1)
8: end for
1
9: update the bandwidth Be by solving
8:
for each St do
2(1)
FIS(((W) }.33 = a FSS( (.,5(11)..1111
B-11)
9:
simulate x I.M = {xm ~ p(x({)
(I)
(1)
compute IPM 8, (1) D(x (
10:
(0 50.)
(k) wi x w
11:
end for
(1) SB
12:
update bandwidth based the empirical
a- Quantile 3 = 2 (1 ( 10 (a)
if ESS( m ). < then resample K particles E, from {
13:
update weights
11:
p(E) (E))
Wi
12:
set weights wi (k) = 1/KK
Σ
K2 (g()( (5 ((R)
Kn
end if
14: sample K particles {( ~ K(5
14: set K new particles с (E__) ; 2(0))
15:
end for
16: output: posterior particles (A)
update forward kernel Kt+1 (E., Er+1)
16: end for
17: output: posterior particles <//
C Experimental details
Here we detail the experimental configurations to reproduce the results covered in Figure and C.I. small grid search has been carried out over the learning rate, the trust-region parameter e, the batch size, and the number of training samples on the SLCP and Furuta task for N = 50. The best-fitting hyperparameters over the two tasks are reported Table C.3 and used throughout the experiments. The remaining parameters are taken from [32] to make the results comparable.
We complement Figure with Figure C.1 to compare the results in the observation space and include the Furuta pendulum.
Table B.2: Approximations of the optimal backward kernel (E:+1, E:) lead to different SMC- ABC approaches.
Algorithm
Assumption
L.
d
Optimal
me(E.) (E. (c+1))+)
Ne (E:11)
Po, (S. (x. 1 (E
Pa (x++1lx)
PMC-ABC
ग. RJ Pa.
Po. ( Ke (1)
TP,, [x) ( (E..E)11)))))
pa, (5. bi. (())
Pd (E: (xL.))
SMC-ABC
Pre+a R3 pa.
Pe, (6r41 (x(.y)
(2:
Paees (E, +1 Xi 1 (E
MH-ABC
Le(E,+1. Ex) (441(E+111
Po, (E: (e4
19

---

Gausiia Lucation
Gaussian Mixcure
SLCP
SIR
Puruta
MMD
18
06
11/12- MODP
MMDD-FLI
MMDD-BB
W-ABC
Figure C.I: Evaluation of the posterior performance on three different tasks (displayed along the columns). A node represents the mean and standard deviation 10 experiments with different random seeds, each carried out using N data points for conditioning. The samples from the approximate posterior ~ q(E x are compared against the reference posterior samples with Wasserstein distance and the MMD when available. Additionally, the log probability of the ground-truth parameter şat is evaluated and posterior predictive checks are carried out on all tasks. The ground-truth parameters are described in Appendix C. Lower values are better for all metrics.
C.1 Approximation of the integral probability metrics
In the context of this paper, we consider two instances IPMs D(p* (x), p(x{)) between the data generating distribution p" (xx and the likelihood p(x(g) the maximum mean discrepancy MMD and the squared 2- Wasserstein distance W2. To simplify the notation, we formulate the discrepancy between the pdfs p(x) and q(x) whose empirical probability distributions denoted by p(x) 1/N ΣM &x, (x) and q(y) 1/M Σ M individual samples by c(x,y) y).
E (y). Furthermore, we denote the cost between
Maximum mean discrepancy The MMD [25] can be formulated with respect evaluation kernel k(x,y) as the sum of three terms
MMD (p.D?(p.g)
E x-p(x)
[<(x,))] -
[k(x,)) y)]
[k(x, y)].
(34)
x-<(())
x--((x)) yooo(y)
yr-p())
y~uly)
As the expectations generally intractable, unbiased estimate of MMD based on samples drawn from p and q is used [25]
N.A
M
MMD (p.
Σ k(x.
Σ
k(x..ys)
Σ
x(x,,x)) (35)
N(N
NM
M(M 1)
20

---

Here, and y, represent samples drawn from the sampling distributions X1:V p(x) and VI:M ~ q(y). In this paper, a Gaussian kernel exp(-1/(2E)c(x, y)) with bandwidth is employed. The bandwidth is known to be very sensitive, which is why the kernel is evaluated over variety of bandwidths = (1, 10, 20, 40, 80, 100, 130, 200, 400, 800, 1000} by summing over the bandwidths k,),,,,) = Σ3 (,4,,, M5) [16].
Wasserstein distance. In the experiments, we consider the squared Wasserstein, which can be formulated as the solution to the optimal transport problem
W (p. 1) inf ()) dxdy.
(36)
The problem is also known as the Kantorovich problem, searching for the optimal coupling E [[p q) in the set of joint distributions that admit p and q through marginalization. For the empirical measures p and G. the Kantorovich problem can be formulated as a linear program
Wi
min Σ Pr.Cus. PEU
(37)
Here, Cig c(x;; Yj) is the cost matrix containing the pairwise comparisons between the samples drawn from p(x) and j(y). The linear program searches for the optimal coupling matrix P among the set of doubly stochastic matrices U = {P E R RRNM PMM 1/N1N, P11N = 1/M1M (1MMM} Peyré and Cuturi [44] show that introducing an entropy regularization term (P) to the objective (37) leads to an iterative scheme that can be solved with the Sinkhorn algorithm [46]. The iterative procedure enables parallelization on hardware accelerators to efficiently solve the optimal transport problem. Furthermore, it can be seen that the 2- 2-Wasserstein distance can be recovered in the limit & O. We leverage the JAX library, OTT [9], approximate W2 with Sinkhorn iterations for the computations.
C.2 Gaussian location
The Gaussian location model is 10-dimensional Gaussian model. The ten dimensional parameters ş E 1| define the means of the model N (x|[p - = & Σ = 0.1I). We choose a Gaussian prior P(E) (V{jo, 00 0.1I) for which the posterior can be recovered in closed form. The ground-truth parameter is sampled uniformly within the posterior support est ~ 1)
1.0
Bolverscn
S0813-PLL
APT
MMDDARI
3
0.0
0.0
-0.3
-0.5
-0.5
00.5
-1.0
-1.0
f1
1.0
1.0
ML
APT
MMDDABC
47
0.0
0.0
-0.2
-0.5
-0.5
-0.5
100
100
100
100
-10
-10
6
Figure C.2: Slice of the posterior through the E - E2 plane. The upper row shows experiments conducted on N 2 reference observations. The lower shows the approximate posteriors for N 100 reference observations. The dotted line represents the ground truth parameter g(gt) that was used to generate x All approaches show that the posterior becomes denser when conditioned more data.
21

---

C.3
Gaussian mixture model
The task to infer the mean parameters of two-dimensional multivariate Gaussian mixture model [47]
p(x(E) (x &. + +N(xxlu E. 0.01I)
(38)
The two-dimensional observation space represents samples from the Gaussian mixture model. We assume uniform prior p = 41-10, 10).
C.4 Simple-likelihood complex-posterior
The Simple- Likelihod Complex- Posterior (SLCP) task consists of 5-dimensional parameter space E [-3, 3] with uniform prior (1-3,, 3). The ground-truth parameter, from which the reference
Complex-Posterior
observations are generated, is {(tt) (0.7, 1.5, -1.0, 0.9, 0.6)T. The observations represent four samples from 2-dimensional Gaussian distribution
X1, X2 X4, x, Vxx,, (LE), E(E)).
(39)
For further information, we refer the SBI benchmarking paper from Lueckmann et al. [32].
C.5
SIR
The SIR model is an epidemiological time-series model that models the spreading of a disease. The name derives from the three states, (i) susceptible, (ii) infectious, (iii) recovered, that an individual can be in. The parameters the dynamics model are the contact rate 3 and the mean recovery rate ?
[ [
(0. (0,
E
(40)
The prior is log-normal distribution over and
Log (logo (0.4), 0.5)
(41)
(log(0000
(42)
We rollout the dynamics over 160 timesteps and evaluate the simulation at 20 equidistant time-steps by taking a sample from the binomial, 2; Binom (1000, I..N). Here denotes the total population at the start of the simulation. See Lueckmann et al. [32] for details on the dynamics of the SIR model.
C.6 The Furuta pendulum
This inverted double pendulum can be described by the angular deflection (Or, 0p] of the rods w.r.t. their equilibrium position [0, 0]. The equations motion can be derived by formulating the Euler-Lagrange equation [37]:
O, +
1 12 1
+ mpl? + 1 mp sin? 0p
mmplph cos 0p]
O 0p
- dpO,
=
2 mppplr cos 0p
m 13 sin 20, 8,0,
2 m m sin 0, 021
(43)
my sin 20, 0 + pg sin
Here, the mild assumption is made that the pole length is significantly greater than its diameter for which the moments of inertia of the poles around their pivot are Jj 1/3 m.l? mass matrix contains entries from the translatory and rotational movement of the two poles. As the reference coordinate systems constantly rotating w.E.t. the basis coordinate system, Coriolis forces occur. They are complemented by gravitation which works on the rotational pole. The left-hand side considers damping in the joints, represented by the damping coefficients d, and dpp and the torque which is applied from a servo motor. For this paper, we omit external forces, i.e., T = N The Furuta pendulum is set into motion by perturbing the initial state around its unstable equilibrium.
fr,p) p). The
For the system identification tasks we select the five system parameters
g
[9, 11]
9.81
1
[0.08, 0.09]
0.085
=
m,
E
[0.08, 0.1
0.095
Ip
[0.12, 0.15
0.129
[mp]
[0.02, 0.03]
0.024
with uniform prior on the predefined ranges.
22

---

Table C.3: Hyper-parameter settings of the SBI methods as used for the experiments in Section 4.
Forward slashes symbolize layers of a neural network.
Parameter
Value
PLI (Ours)
Likelihood kernel
Exponential Kernel
Trust-region threshold
0.5
Model
Neural Spline Flow (NSF)
Bijector
Rational Quadratic Spline with param size D
# Bins
10
Conditioning MLP
input dim / 50 50/ 50 D
# Bijectors/ / Transforms
5
Base distribution
Learning rate
10 10-5
Epochs
20
Train samples per iteration
5000
Batch size
125
PMC-ABC [31]
Likelihood kernel
Uniform Kernel a-Quantile, a = 0.1 (see Table A.1)
Likelihood update
0.1
Reverse transition kernel L,
reverse PMC kernel (see Table B.2)
Particles
1000
Epochs
200
Perturbation kernel
GMM with components
APT [24]
Model
Conditional Neural Spline Flow (NSF)
Bijector
Rational Quadratic Spline with param size D
# Bins
10
Bijectors/ / Transforms
5
Conditioning MLP
input dim /22 /22/32/1
32/ D
Base distribution
A((0, 1)
# Atoms
10
Learning rate
x 10
Epochs
20
Train samples per epoch
5000
Batch size
500
23

---

-2.5
-2.6
2.0
0.6
81
81
00
0.0
E1
E1
6
Figure C.3: Results of posterior inference on the SLCP task with reference observations. The unimodal distribution of the parameters & and E2) are depicted well by all approaches. On the contrary, the multi-modall modality only represented properly by the APT posterior (- Reference.
MMD-PLI, - APT. - MMD-ABC).
24

---

-2.5
-25
0.0
22
00
25
0.0
2.5
0.0
-2.5.0.0
-26
0.0
-2
00
28
00
0.0
61
55
0.0
-2.5
0.0
n.0
0.0
0.0
E1
52
83
E4
51
E1
C
Figure C.4: Results of posterior inference on the SLCP task with = 100 reference observations. Compared to the posterior given N = 2 observations (Figure C.3). the posterior (- Reference) is distributed tightly around distinct points. Here, MMD-PILI captures all modes of the posterior, - -MMD-ABC centers around a uni-mode, while - APT cannot represent the multi-modality.
25

---

0.00
0.06
0.08
0.005
0.080
0.080
0.030
0.000
0.33
0.12
0.12
110.00 05
1.10
0 005
0.13
0.18
10
11005
10
0.80000
0.19
0.5
0.00
0.065
0.025
0.020
0.12
10
0.10.05
0. G.M
0.12
0.15
0
Figure C.5: Results of posterior inference on the Furuta pendulum with N = reference observations. All methods center around the ground truth parameter. APT finds the expected correlations among the parameters while MMD-PLI and - MMD-ABC remain more widespread.
26

---

3. 0.00
0,00
2 0.083
0.005
0.060
<5 0.023
J 0.009
0.020
0.020
0.12
0.12
00.00
10:00
04000200 005
0.11
0.18
1:00.0
0010.00
OAR0020 025
0.12
0.18
0.30
0.000
03185
0.1025
8
70p
Figure C.6: Results of posterior inference on the Furuta pendulum with N 100 reference ob- servations. All models capture the ground-truth parameter well. contrast to the N 2 setting (Figure C.5) - MMD-PLI reveals pairwise correlations between the domain parameters, and - MMD- ABC is less densely distributed. Note, that APT clusters outside the ground-truth parameter.
27


Reviews:
[{"review_id": "qki8nqZlPB", "reviewer": "Reviewer_XuP4", "summary": "This work aims at proposing a sequential Monte Carlo scheme to perform approximate Bayesian computation in the context of simulation based inference. Previous works focused on directly modeling posterior distributions whereas this approach presents a novel method to approximate likelihoods utilizing kernel density methods along with an information theoretic kernel bandwidth update scheme. The authors claim that this direction allows for improvement over previous works in the case of having more reference data to learn from.", "questions": "Is there a clear explanation as for why the APT method starts to perform worse when given more reference data? Naively I would think that it would plateau in performance but it appears to severely degrade. We see somewhat similar behavior in some of the settings with both PLI and ABC methods, although to a much lesser degree (and more infrequent). I would be interested in learning if there is any insight into why this could be happening to all of these methods.  ", "limitations": "The authors did adequately discuss the limitations in my opinion.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "I think that employing variational inference to help make SMC-ABC more tractable is an interesting direction. Further leveraging this setting to derive a trust-region compliant kernel bandwidth is significant and a good direction to explore. A number of different non-trivial steps are performed in order to derive the final pseudo-likelihood inference scheme and I believe this deserves commendation.", "weaknesses": "Almost all of the derivations are based on the classic variational inference objective of reducing the KL divergence between a true distribution and the variational distribution being learned; however, the actual method proposed substitutes this for alternative divergences/\"integrated probability metrics\" such as Wasserstein distance and MMD. I understand the necessity of this due to the likelihood being unavailable in SBI; however, this does undercut a lot of the impact that the various derivations have on the actual problem in my opinion. \n\nI do not find the experimental results as convincing as they could be. For the first two benchmarks (Gaussian Location and Gaussian Mixture), it appears that MMD-ABC is either comparable to both versions this paper is proposing or better across all values of N and in both metrics. Additionally, in the appendix it can be seen in Figure C.1 that using sample based metrics (referred to as PPCs in the paper) the PLI methods appear to perform about the same or worse in several cases in all 5 tasks. I wouldn't go as far to say that they are always comparable or worse, but the story does not seem as clear cut as it is made out to be in the main set of results in Figure 3 in the main paper. Note that in line 285 it is justified to use PPCs instead of comparing posterior distributions, so I believe it should be valid to look at these results for the benchmarks as well.\n\nPlease feel free to correct me if I have misunderstood something as I could have easily missed some key details.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Answer to Reviewer XuP4\nWe thank the reviewer for the valuable feedback regarding our proposed method.\n\n**Remarks and Questions:**\n> I understand the necessity of [introducing IPMs over the KL] due to the likelihood being unavailable in SBI; however, this does undercut a lot of the impact that the various derivations have on the actual problem in my opinion.\n\nIt is indeed a crucial idea of our paper to use IPMs to compare samples when the KL divergence cannot be directly applied. Please see our detailed discussion on the impact of the different divergences in our general response.\n\n> For the first two benchmarks (Gaussian Location and Gaussian Mixture), it appears that MMD-ABC is either comparable to both versions this paper is proposing or better across all values of N and in both metrics.\n\nPlease observe that the Gaussian Location and Gaussian Mixture tasks are simple benchmarks for which we can obtain the exact posterior and where the algorithms are expected to perform well. It is not surprising that ABC even slightly outperforms PLI on these tasks because it is not training a model of the posterior. On the more challenging tasks, such as APT and SIR, there is a significant drop in the performance of ABC while the PLI performance increases with increasing $N$. This is a key result of the paper since we want to use many observations instead of only one.\n\n> Additionally, in the appendix it can be seen in Figure C.1 that using sample based metrics (referred to as PPCs in the paper) the PLI methods appear to perform about the same or worse in several cases in all 5 tasks. Note that in line 285 it is justified to use PPCs instead of comparing posterior distributions, so I believe it should be valid to look at these results for the benchmarks as well.\n\nPlease let us clarify how to interpret the additional PPC metrics in the appendix. Due to the intractability of the likelihood, the evaluation of SBI methods in the parameter space is in general infeasible. Therefore, we provided posterior predictive checks (PPCs) to evaluate the posterior predictive performance in the observation space  $\\text{PPC}^D = \\sum_i D(p^\\star(x), p(x|\\xi_i)),\\,\\xi_i \\sim\\pi(\\xi_i)$. Your observation that the PPC scores are not matching the parameter space metrics is partially correct, although we restrict this observation to the SIR task. After submission, we found a small bug where the evaluation data was drawn from an incorrect ground truth parameter. We fixed the bug and reran the evaluations. The corrected results align with the paper's claims and show that PPC can be used as a valid performance measure. We checked all other experiments, and there the results are correct. Thank you for pointing out this inconsistency.\n\n> Is there a clear explanation as for why the APT method starts to perform worse when given more reference data?\n\nWe hypothesize that it is due to the high dimensional observation space, as all observations are concatenated. As the simulation budget was fixed to make the methods comparable, this might explain the significant drop in performance. In our general rebuttal response above, we discuss the performance drop in further detail.\n\n> We see somewhat similar behavior in some of the settings with both PLI and ABC methods\n\nHere we would like to stress that this behavior is distinct from the behavior seen in APT. The slight decay of performance for large $N$, i.e., $N=1000$ can be attributed to the fact that the reference posterior samples are only approximations obtained from MCMC sampling for the SLCP and SIR tasks. For large $N$, obtaining reliable posterior samples with MCMC becomes a challenging task. Note that this phenomenon is not visible in the Gaussian location and Gaussian mixture tasks as these posteriors can be evaluated analytically. \n"}, {"review_id": "M0Aj0vQRzw", "reviewer": "Reviewer_AzJv", "summary": "This paper introduces Pseudo-Likelihood Inference (PLI) as a novel method inside the simulation-based inference (SBI) family. SBI methods are used to estimate model parameters when the likelihood function is unknown or intractable. Existing SBI methods, such as Approximate Bayesian Computation (ABC) and Sequential Neural Posterior Estimation (SNPE), either approximate the likelihood or directly model the posterior. The contribution by the authors, PLI, combines neural approximation with ABC to achieve competitive performance in challenging Bayesian system identification tasks. The authors make use of integral probability metrics and introduce a smooth likelihood kernel with an adaptive bandwidth that is updated based on information-theoretic trust regions. This formulation allows for optimizing neural posteriors via gradient descent and does not rely on summary statistics. PLI performance improves the more data is available, particularly in scenarios involving stochastic simulations and multi-modal posterior landscapes. The effectiveness of PLI is evaluated through experiments conducted on four classical SBI benchmark tasks as well as a highly dynamic physical system. The results demonstrated particular advantages of PLI in handling stochastic simulations and multi-modal posterior landscapes.", "questions": "* What are the limitations and potential challenges associated with implementing PLI? The paper presents PLI as an effective method for challenging Bayesian system identification tasks, but it is important to address potential limitations or practical considerations when applying this approach in real-world scenarios. For example, what computational resources are required by PLI? Are there specific problem domains where its performance may be less favorable?", "limitations": "* The paper could improve by comparing PLI with more SBI methods.\n  \n* The evaluation focuses on a limited set of benchmark tasks and one physical system, potentially limiting the generalizability of PLI to diverse problem domains.\n\n* There is insufficient discussion about the computational efficiency and resource requirements of implementing PLI.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "* PLI is a novel approach addressing the limitations of existing methods that achieves competitive performance in challenging Bayesian system identification tasks. In fact, the proposed method outperforms SNPEs and other likelihood-free models, particularly when more data is available.\n\n* The usage of the adaptive bandwidth for the smooth likelihood kernel allows for better adaptation to different problem settings and improves the accuracy of parameter estimation.\n\n* The method enables to directly use gradient descent optimizing neural posteriors using gradient descent directly on raw observations or data points without requiring any intermediate summarization step.\n\n* The effectiveness of PLI is evaluated in multiple experiments conducted on four classical Simulation-Based Inference benchmark tasks and in a highly dynamic physical system, demonstrating its applicability across various domains.", "weaknesses": "* The changes induced by replacing the KL divergence with the MMD and the Wasserstein distance should further explored. Consequences on changing the KL divergence for other metrics has been studied before in works such as [1], where an adversarial approach is used to evaluate the KL divergence and other divergences using only samples. I strongly suggest discussing this possibility in the article, since it seems that this alternative can be useful here.\n\n* To provide more robust experimental results, when comparing different models I would suggest also using other proper scoring rules [2].\n\n* Since efficiency is an important aspect when dealing with high-dimensional tasks, there is limited discussion about computational efficiency throughout the paper. Considering that likelihood-free inference methods can be computationally demanding, it would be valuable to discuss any potential limitations or considerations regarding computational resources required by PLI compared to other SBI approaches\n\n* While the paper highlights PLI's advantages over SNPEs, it would be valuable to understand how PLI performs against other popular approaches s.a. Markov Chain Monte Carlo methods. A comparative analysis would provide a more comprehensive understanding of PLI's strengths and weaknesses relative to existing techniques.\n\n* (_minor_) Due to the formulation of several function-space Bayesian inference models, if you deemed it fitting it could be interesting discussion on the function-space approaches such as [3-6]. I think most of these models can be encompassed as (S)NPE methods and that could enrich the discussion in the paper.\n\n* (_minor_) Please correct typos, s.a. \"posteriror\" in line 143, \"divergnce\" in line 156.\n\n* (_minor_) The red hyperlink boxes to equations appear with an extra line inside, which difficults reading. Please check this.\n\n---\n\n**References**:\n\n[1] Santana, S. R., & Hernández-Lobato, D. (2022). Adversarial α-divergence minimization for Bayesian approximate inference. Neurocomputing, 471, 260-274.\n\n[2] Gneiting, T., & Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477), 359-378.\n\n[3] Ma, C., & Hernández-Lobato, J. M. (2021). Functional variational inference based on stochastic process generators. Advances in Neural Information Processing Systems, 34, 21795-21807.\n\n[4] Rodrı́guez Santana, S., Zaldivar, B., & Hernandez-Lobato, D. (2022, June). Function-space Inference with Sparse Implicit Processes. In International Conference on Machine Learning (pp. 18723-18740). PMLR.\n\n[5] Ma, C., Li, Y., and Hernández-Lobato, J. M. (2019). “Variational implicit processes”. In: International Conference on Machine Learning, pp. 4222–4233.\n\n[6] Sun, S., Zhang, G., Shi, J., and Grosse, R. (2019). “Functional variational Bayesian neural networks”. In: International Conference on Learning Representations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Answer to Reviewer AzJv\nWe thank the reviewer for the valuable feedback and the relevant references that we will include in the final version of the paper.\n\n**Remarks and Questions:**\n> The changes induced by replacing the KL divergence with the MMD and the Wasserstein distance should further explored.\n\nThank you for this suggestion. We provided further insights in our general rebuttal response above regarding the choice of the divergence measure. Indeed, we agree that the effect of the divergence measure deserves further study and that there are other choices that one could consider, e.g., Kernelized Stein Discrepancy (KSD), the adversarial estimation of the KL, or the alpha-divergences. We believe that it is a fruitful direction for future work, and we hope that our paper stimulates this research by showing strong empirical results using MMD and Wasserstein distances. We will include a paragraph in Section 3.3 that discusses adversarial approaches as a way to approximate the KL and provide a reference to the suggested paper:\n*Alternatively, one could exploit adversarial strategies to approximate the KL divergence, as explored in the studies by Mescheder et al. (2017) and Santana and Hernández-Lobato (2022). Incorporating these methods into the PLI framework requires training an additional adversarial model. We leave this investigation to future work.*\n\n> Using other proper scoring rules\n\nThanks for this suggestion. Indeed, the choice of the evaluation metric is important, as was also noted in the SBI benchmark [1]. Currently, we employed MMD and Wasserstein distance, the former being an instance of a kernelized scoring rule and the latter having a clear geometric interpretation. Since both metrics yield similar results with growing $N$, we did not further investigate this choice. Are there any particular proper scoring rules you suggest we include in the final version of the paper? \n\n> Considering that likelihood-free inference methods can be computationally demanding, it would be valuable to discuss any potential limitations or considerations regarding computational resources required by PLI compared to other SBI approaches\n\nPLI has the same computational requirements as other SBI methods that train neural models, such as SNPE. Only compared to ABC does PLI take longer, but note that ABC does not have a training phase, whereas PLI can be queried very cheaply after it is trained. Please see our overall response for more details. We will add these clarifications to the final version of the paper.\n\n> While the paper highlights PLI's advantages over SNPEs, it would be valuable to understand how PLI performs against other popular approaches s.a. Markov Chain Monte Carlo methods.\n\nWe have compared against SMC-ABC --- the strongest method amongst ABC-based approaches. MCMC-ABC is not so popular as it requires a fixed rejection threshold, making the problem significantly harder when the prior distribution is not close to the target distribution. We highlight that the reference posterior samples for the SLCP and SIR tasks are drawn from MCMC sampling with the true likelihood. Thus, we are already implicitly comparing against MCMC sampling and see that Wasserstein and MMD converge to the MCMC samples.\n\n> Due to the formulation of several function-space Bayesian inference models, if you deemed it fitting it could be interesting discussion on the function-space approaches\n\nAs far as we know, Functional Variational Inference (FVI) has not been considered in the Likelihood-Free Inference (LFI) context so far. It sounds interesting and we will mention it as a potential future work, together with the references to the suggested papers.\n\n> What are the limitations and potential challenges associated with implementing PLI? Are there specific problem domains where its performance may be less favorable?\n\nWe want to highlight two scenarios in which PLI's *\"performance may be less favorable\"*. First, PLI is not amortized and won't perform well when reference data of different parameter configurations is present. Second, PLI relies on the MMD or Wasserstein distance to assess the similarity of the data, which can be a hard task when only a few observations are available. Additionally, all SBI methods assume that the simulation matches the data-generating distribution closely. Accounting for modeling errors is an ongoing large field that is far from solved.\n\n[1] Lueckmann, J. M., Boelts, J., Greenberg, D., Goncalves, P., & Macke, J. (2021, March). Benchmarking simulation-based inference. In International conference on artificial intelligence and statistics (pp. 343-351). PMLR."}, {"review_id": "6SnnUq8ggB", "reviewer": "Reviewer_9BA6", "summary": "This paper proposes an approch to deal with posterior inference when the likelihood is intractable, which is often the case for problems involving simulators of physical processes.\nThe paper investigates Approximate Bayesian Computations (ABC) methods, and it proposes incorporating neural approximations through variational principles.\nThe main strenghts of the proposed approach include the possibility to carry out optimization of neural posteriors via gradient descent, the lack of reliance on summary statistics, and the possibility to handle multiple observations as input. \n", "questions": "As a non-expert on ABC I found the paper interesting and well written, although I find it difficult to assess the novelty of this work.\nReplacing the KL with some objective that can be estimated through samples is a good idea in general, and this paper does it in eq. 6; not sure how novel the ABC community would consider this.\nSimilarly, moving away from modeling the posterior through particles in favor of neural density estimators yields clear advantages, and again I'm not sure how novel this would be considered in this literature.\nHaving said that, I found the motivation behind the work convincing and the idea well executed with experiments showing improvements over alternative approaches from the state-of-the-art. \nAs a result, I have an overall positive opinion of the paper which combines sensible ideas to produce an approach which seems to work well across a number of experiments.\n", "limitations": "I haven't seen any specific text on the limitations of the approach and I think this could be included in the revision. \n", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is well written and easy to follow.\n\nThe methodological developments seem solid and well motivated.\n\nExperiments seem to support the main message of the paper.\n", "weaknesses": "Experiments could have included a larger variety of scenarios involving other types of data and simulators.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Answer to Reviewer 9BA6\nWe thank the reviewer for the valuable feedback. We appreciate the *\"overall positive opinion\"* on the methodology but would like to summarize the main contributions of PLI again as the reviewer *\"find[s] it difficult to assess the novelty of this work\"*\n\nWe want to motivate the novelty of PLI from the perspective of ABC and APT. The idea of replacing summary statistics with statistical measures has been extensively studied in ABC and has been shown to improve ABC methods significantly. The PLI formulation allows leveraging statistical measures to train neural posterior models, thus making the training more robust compared to APT. In the SBI community, training neural posterior models has been typically restricted to generating a training dataset using the simulator and optimizing the neural model based on the maximum likelihood principle. The introduction of weighted maximum likelihood can be considered novel from the perspective of SBI. In short, individually, the two ideas - using statistical measures and training neural posteriors - are already known and popular in their respective fields, but our PLI algorithm brings them together in the context of ABC and thereby achieves state-of-the-art performance.\n\n**Remarks:**\n> Experiments could have included a larger variety of scenarios involving other types of data and simulators.\n\nWe evaluate the PLI methodology based on common benchmarking tasks in the SBI community and extend the benchmark to scenarios where $N>1$. We agree that it is interesting to evaluate how well SBI methods scale to other types of data, in particular scaling to high dimensional problems. However, note that we are already going beyond what is commonly done in SBI papers regarding problem dimensionality. In particular, we introduced the Furuta pendulum task as a new benchmark featuring multi-modal high-dimensional observation space. As experiments with high-dimensional observation and parameter spaces are quite compute-intensive, we believe it is reasonable to consider them in future work. For additional details on the choice of evaluation tasks, please refer to our overall rebuttal response.\n"}, {"review_id": "cgkkNtIKx3", "reviewer": "Reviewer_dHap", "summary": "This paper proposes a Pseudo-Likelihood Inference approach that brings neural approximation into ABC. This new approach can condition on a variable number of observations unlike previous neural-based SBI approaches. ", "questions": "* The hyperparameters seem to be set without mention of any hyperparameter optimization. What effect does the number of particles have and the number of steps have on the overall performance? What do bins mean on the Supplementary materials?", "limitations": "* The idea appears novel and is the main reason for the slightly positive score. The experimentation could definitely be improved, as well as the explanation behind various hyperparameters and what the limitation of the approach is when it comes to computational cost.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "\n* The incorporation of Sequential Monte Carlo into SBI is a strength. SMC seems like it should be related to SBI as it follows the structure (as shown in Fig. 1) of annealing to the desired distribution.\n* The literature review is a strength and highlights the authors' knowledge on the topic.\n* The paper is well-written.\n", "weaknesses": "* The experiments section could be strengthened to highlight the advantage of the approach compared to SNPE. For example, is the a better way to highlight that $N \\not= M$ for PLI? It seems a shame not to highlight this strength of the new approach. Also, does $N=M$ mean that for APT, when $N=1$, only one training point is used to train the normalizing flow?\n* In Lueckmann et al. (2021) there are actually 10 benchmark experiements. Is there a reason that only four were used here? Additionally, if “benchmarking is not the focus” (line 247) then what is the focus? This motivation could be clearer.\n* The computational cost is not highlighted in the paper. Is the cost of retraining the density estimator at each step equivalent to APT?\n* It is a weakness that the code was not provided. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Answer to Reviewer dHap\nWe thank the reviewer for the insightful questions and feedback. \n\n**Remarks and Questions:**\n>  Is there a better way to highlight that $N\\neq M$ for PLI?\n\nIndeed, choosing $M > N$ is a strength of PLI and ABC that we could have highlighted more. We expect both methods to be more competitive against APT when using $M>N$, particularly when only a few reference observations $N$ are available. Our reasoning for restricting the evaluation to $M=N$ was to keep the simulation budget comparable.\n\n> In Lueckmann et al. (2021) there are actually 10 benchmark experiments. Is there a reason that only four were used here? Additionally, if “benchmarking is not the focus” (line 247) then what is the focus?\n\nWe intend to guide the experiments as closely as possible to Lueckmann et al.'s (2021) benchmarking paper—note, however, that there the number of observations is fixed to a single sample N=1. As the generation of the reference posterior samples with MCMC methods proved to become more unstable for larger N, we only showed the results for environments in which we found good posterior samples for up to N=1000 reference observations. For further details on the choice of environments, we refer to our general response.\n\n> The computational cost is not highlighted in the paper. Is the cost of retraining the density estimator at each step equivalent to APT?\n\nYes, the training and simulation times of PLI and APT are equivalent. The only computational overhead of PLI is the calculation of the pseudo-likelihood, which is negligible. For a more thorough discussion of the computational effort, we refer to our general rebuttal response above.\n\n> It is a weakness that the code was not provided.\n\nWe, unfortunately, did not have the code cleaned and ready to publish for the initial submission. As of now the code is clean and will be published with the final paper.\n\nQuestions:\n> The hyperparameters seem to be set without mention of any hyperparameter optimization.\n\nWe did perform a grid search for the learning rate, the number of training samples, the trust-region parameter $\\varepsilon$, and the batch size on the SLCP and Furuta tasks for $N=50$. The remaining model-specific parameters of APT were taken from the SBI benchmark and adopted to PLI to make the methods as comparable as possible. We added a note to the experimental details section in the appendix.\n\n> What effect does the number of particles have and the number of steps have on the overall performance?\n\nOne of the main difficulties of simulation-based inference is the exploration-exploitation trade-off in the parameter space. Increasing the number of particles, in general, helps exploration as the search space is better covered. The number of steps is typically coupled to the perturbation kernel. Increasing the number of steps facilitates smaller perturbations, thus enabling a more fine-grained search in the parameter space. In SMC-ABC, this quantity can be tuned by the choice of either the $\\alpha$-Quantile in PMC-ABC or the effective sample-size parameter $\\alpha$ and resampling threshold $V$ in SMC-ABC. Further information can be found in [1].\n\n> What do bins mean on the Supplementary materials?\n\nThe number of bins is a particular parameter of neural spline flows [2], denoting the number of rational-quadratic polynomial segments. The introduction of monotonic polynomial splines has been popularized in normalizing flows that use coupling transforms [2]. Coupling transforms is a common technique to make the Jacobian determinant of the change of variables tractable.\n\n[1] Lenormand, M., Jabot, F., & Deffuant, G. (2013). Adaptive approximate Bayesian computation for complex models.\n\n[2] Durkan, C., Bekasov, A., Murray, I., & Papamakarios, G. (2019). Neural spline flows."}, {"review_id": "VhsIK1nbwu", "reviewer": "Reviewer_7dVy", "summary": "This method perform simulation-based inference with a hybrid of previous methods, in this case integrating a neural likelihood-type step within a SMC-ABC algorithm for performing posterior inference. It considers use of distance metrics such as the Wasserstein and the Maximum Mean Discrepancy to compare simulated and observed data, while performing particle-based updates in a fairly standard way for SMC-type methods. It performs a comparison on five different test examples, with comparisons with different ABC methods and a neural likelihood method.", "questions": "The Furuta pendulum is a nice example. Are there other examples of using chaotic systems rather than stochastic systems as test cases for simulator based inference? Do you treat the variation between two adjacent deterministic chaotic systems as being basically random, and is there a downside to making that assumption?\n\nIs “APT” in the results section basically targeting the posterior directly without the SMC-type exploration? Why does it seem to do well up to a given value of N then fall off quite extremely? Is a more concentrated posterior (larger N) more difficult to explore, and the SMC-type method helps with exploration of a complex posterior? Constructing a neural likelihood approximation at every iteration must add quite a large computational overhead that should be acknowledged.\n\nWhat is the relevance of model misspecification here? This is acknowledged in the section describing General Variational Inference, and it would be nice to see at least a small discussion of the possible advantages of using the MMD or Wasserstein compared to the KL for the “true” likelihood.\n\nWhat is the expected scaling with dimensionality of the parameter space? You have constrained yourself to low-dimensional examples in this paper. The density-ratio type methods tend to do quite well in higher-dimensional spaces as neural nets can fit in this spaces with enough data, but SMC suffers quite severely with the curse of dimensionality. Does the use of neural networks in PLI help overcome this curse? If yes, great! If no, is it worth the overhead of training a neural net in a lower-dimensional space when there might exist other methods that do the job adequately? Would the straight APT method do better than any of the SMC-type methods in a higher-dimensional space?\n", "limitations": "No clear immediate social impact of this work. It really depends on the model on which inference is being performed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The method is fairly well presented, at least the steps setting up the framework of the algorithm used.\n\nThe examples used are appropriate, albeit limited to low-dimensional parameter spaces. The Furuta pendulum, which I had not seen before, is particularly nice.\n\nYou have chosen a few competing methods which, by no means exhaustive, at least make similar assumptions to PLI and therefore the comparison is revealing.\n\n“Related Work” is fairly thorough.\n", "weaknesses": "There are a few little typoes in the language that I noticed: \"Posteriror” “numerica”. I imagine there are more that a standard spellcheck would help with.\n\nThe actual algorithm and computational steps used to estimate the relevant divergences could be much more thoroughly explained, i.e. “Compute IPM” without saying how to for MMD vs Wasserstein vs other is unclear. You provide a reference for this estimation but I think it should be better explained in this article itself.\n\nI’m not sure that I really buy the analogy with variational inference: the proposal distribution formulation built with weights is fairly standard from the SMC literature, and I don’t think the variational angle is particularly clear. I imagine that you do have a legitimate point to make in this respect, so if you could expand your reasoning on this line then that would be a good thing.\n\nA few specific claims that I didn't really agree with:\n\n“There might be several parameter configurations yielding the same observation, hence rendering the resulting distribution to be multi-modal.” - Multi-modality is not always a result of non-identifiability, just different parameter values having comparable likelihoods, not necessarily the same observed data.\n\n“KL divergence is intractable in SBI because we cannot evaluate the likelihood” - can be evaluated without the likelihood with density ratio estimation and classifiers as in classifier ABC and f-divergences. FIND CITATIONS.\n \nThe term “likelihood kernel” is introduced without explanation.\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Answer to Reviewer 7dVy\nWe thank the reviewer for the extensive feedback.\n\n**Remarks and Questions:**\n> [explain] computational steps used to estimate the relevant divergences.\n\nWe will add:\n*We follow [1] to approximate the MMD for two discrete probability measures, whereas entropy-regularized optimal transport is used to approximate the Wasserstein distance [2]. Both versions facilitate parallelization on the GPU.*\n\nAnd we added a subsection to Appendix C describing the exact implementation of MMD and the Wasserstein distance.\n\n> [...] the variational angle is [not] clear.\n\nViewing PLI from either the SMC or the VI perspective makes the method interesting. At each inference iteration, we solve the variational inference problem in eq. (6). The closed-form solution of the constrained optimization problem results in a weighted formulation which shares many characteristics with SMC methods. Thus, our method is similar to Particle Mirror Descent [3] and Particle-Based Variational Inference [4] in that they also represent the variational distribution using particles. Therefore, we also consider the particle-based representation as a type of variational family.\n\n> Multi-modality is not always a result of non-identifiability [...]\n\nThough our statement does not exclude the case of having different observations with the same likelihood, we agree that we should rephrase it for completeness:\n*Additionally, several parameter configurations may yield observations with comparable likelihood, rendering the posterior distribution multi-modal.*\n\n> [KL] can be evaluated without the likelihood [...]\n\nTo be precise, we indeed cannot \"evaluate\" the likelihood, but of course, if the KL divergence is approximated with another estimator, then we could use it. This might be an interesting idea to explore in future work. As we intended to outline that the KL cannot be estimated without specifying an adversarial model, we update the formulation accordingly:\n*Without specifying additional estimators, the KL is intractable in SBI, ...*\n\n> The term “likelihood kernel” is introduced without explanation.\n\nThe likelihood kernels are presented in Table A.1. in the appendix. In the ABC setting, the likelihood kernel represents the uniform kernel with rejection threshold $\\beta$, whereas the likelihood kernel in PLI refers to the exponential kernel in eq.(6). We will ensure it is well defined in the final version.\n\n> Are there other examples of using chaotic systems [...]\n\nWe are unaware of examples of using chaotic systems in the SBI literature. This seems like an exciting direction for investigation.\n\n> Do you treat the variation between two adjacent deterministic chaotic systems as being basically random, and is there a downside to making that assumption?\n\nWe did not make this assumption explicitly. But it is a very insightful question how and if such assumptions influence the applicability of SBI methods for deterministic chaotic systems.\n\n> Is “APT” in the results section basically targeting the posterior directly without the SMC-type exploration?\n\nIndeed \"APT\" is directly targeting the posterior distribution during training. The method leverages a proposal distribution and the simulation data to construct a training dataset based on maximum likelihood estimation to train the posterior. The APT formulation allows bootstrapping the posterior of previous iterations as a new sampling distribution, similar to PLI, and accounting for it through importance reweighting.\n\n> Why does [APT] fall off quite extremely? [...]\n\nAs APT is solely trained on tuples of parameters and simulations, higher dimensional problems that include larger N might require more data or an inductive bias, as introduced in PLI through the weights.\n\n> What is the relevance of model misspecification here?\n\nModel misspecification denotes the discrepancy between the data generating distribution $p^\\star(x)$ of the data generating process and the simulated samples $p(x|\\xi)$. If the data generating model contains outliers, evaluation of $\\text{KL}(p^\\star(x)|| p(x|\\xi))$ for such an outlier would lead to large weights, which heavily influence inference [5]. [6] and [7] show that MMD is more robust to outliers as the MMD posterior concentrates tightly around the data generating distribution at an $\\sqrt{N}$-rate\n\n$$\\lim_{N\\rightarrow \\infty}E_{\\xi \\sim \\pi(\\xi)}[\\text{MMD}(p^\\star(x), p(x|\\xi))] \\leq \\text{inf}_\\xi \\text{MMD}(p^\\star(x), p(x|\\xi)).$$\n\nThe robustness property is interesting in real-world scenarios where the data generation is often accompanied by sensory noise, and models of physical systems can only be accurate to a certain degree.\n\n> What is the expected scaling with dimensionality of the parameter space? You have constrained yourself to low-dimensional examples in this paper.\n\nWe chose those experiments because they are prevalent benchmarking tasks in the related literature. Further, we did introduce the Furuta benchmark to show the scalability. This paper aims to introduce and justify the presented PLI approach. For future work, we agree that scaling PLI to higher-dimensional problems is exciting and should outline the strength of the methodology even more.\n\n> Would the straight APT method do better than any of the SMC-type methods in a higher-dimensional space?\n\nWe expect PLI and APT to scale better to high-dimensional problems than the ABC approaches. We believe that the weighted maximum likelihood approach adopted in PLI should make the optimization of the posterior more robust than APT, especially when scaling to higher dimensions. As APT is only trained on pairs of parameters and observations, we suspect that the additional inductive bias of PLI should be beneficial in guiding the training.\n\n[1] Cuturi (2013)\n\n[2] Gretton et al. (2012)\n\n[3] Dai et al. (2016)\n\n[4] Liu et al. (2019)\n\n[5] Knoblauch et al. (2019)\n\n[6] Chérief-Abdellatif et al. (2020)\n\n[7] Dellaporta et al. (2022)"}]

Author rebuttal:
### Answer to all Reviewers

We thank the reviewers for their insightful feedback. We are pleased that the **novelty and significance** of the paper are recognized by the reviewers.
- **XuP4:** A number of different non-trivial steps are performed [...] and I believe this deserves commendation
- **AzJv:** novel approach addressing the limitations of existing methods
- **9BA6:** methodological developments seem solid and well motivated

The reviewers point out multiple times that the **experiments support the paper’s claims**, that PLI is **competitive** compared to state-of-the-art baselines, and that the **extension to the novel Furuta pendulum** experiment is a strength:
- **AzJv:** The effectiveness [of PLI] is evaluated in multiple experiments [...], *demonstrating its applicability across various domains*.
- **9BA6:** Experiments seem to support the main message of the paper.
- **7dVy:** The Furuta pendulum, which I had not seen before, is particularly nice

We are further pleased about the positive comments regarding **presentation** and **clarity** of the paper:
- **9BA6:** the paper is well written and easy to follow
- **dHap:** The literature review is a strength and highlights the authors' knowledge on the topic.
- **7dVy:** method is fairly well presented

The reviewers raised several good questions. We want to answer these here:

**Choice of the discrepancy measure $D$.**
Some reviewers expressed interest in the theoretical justification for employing MMD or Wasserstein divergence in the pseudo-likelihood kernel. They also questioned whether Lemma 1's convergence result and Appendix A1's derivation for the KL are still valid in the context of other divergences. We aim to answer this question empirically by comparing our posterior to the reference posterior samples. We demonstrate that PLI posteriors, indeed, concentrate on the true posterior. Though rigorous theoretical convergence guarantees remain a work in progress, we cite MMD-Bayes [1] (l. 156-161) for concentration bounds around the posterior predictive distribution for a posterior based on an exponential MMD kernel. In contrast, there is no such result for W-PLI, despite its similarity to Wasserstein-ABC [2], which has convergence guarantees for the uniform kernel. We hope that our empirical findings will stimulate research into the concentration properties of PLI with Wasserstein and other divergences.

**Computational cost of PLI.**
As some reviewers asked about the computational cost, we will include the following summary in the final paper:

*In PLI, computational costs are spread across three components: simulation, summary statistic estimation, and normalizing flow training, all benefiting from GPU parallelization. In our experiments, neural model training accounted for the primary computational effort, while simulation and summary statistic computation were minor. PLI and SNPE required 60-90 min per task, compared to ABC's 2-10 min on an Nvidia RTX 3090. If a simpler model is used, e.g.,  GMMs, PLI's computation time matches ABC's. Commonly, high-fidelity simulators are employed in SBI methods, which are expensive to query. Therefore, simulations will likely constitute a substantial proportion of the total compute in more realistic situations.*

**Performance of APT drops for a large number of observations $N$.**
Several reviewers noticed that the performance of APT drops for large $N$. Since the increase in the number of reference observations means that the dimensionality of the input to the conditional normalizing flow (CNF) increases with $N$, the search space increases significantly for APT, which is trained on pairs of parameters and simulations. In contrast, PLI does not suffer from the increase of dimensionality because the density estimator is not conditioned on the data and has an additional inductive bias through the pseudo-likelihood evaluation.

**Updated experimental results on the SIR task.**
We investigated the issue raised by Reviewer XuP4 that the evaluation metrics in the observation space and the parameter space were disagreeing on the SIR task in Figure C.1. After the initial submission, we made the same observation and identified the source of that issue: the evaluation data (observations) for the SIR task had been generated from a different ground truth parameter yielding incorrect evaluation data. This bug explains the contradicting results between the evaluation in the observation space (upper two rows) against the metrics in the parameter space (lower three rows) in column 4. We have updated the evaluation data in the meanwhile and rerun the evaluation script for the SIR task (column 4). We further checked all other experiments, and there the data was correct. We provide an updated Figure C.1 in the rebuttal PDF.

**Choice of the evaluation tasks.**
Some reviewers suggested including more tasks from the SBI benchmark [3]. We want to clarify why we focused only on a subset of tasks. Namely, we did not use the same setup as in the SBI benchmark, where there is only $N=1$ observation, but we extended it to $N>1$ observations. Therefore, we aimed, on the one hand, at verifying our approach on tasks that have an analytical posterior (Gaussian Location Model and Gaussian Mixture Model) and on which all baseline methods are expected to work well, and on the other hand, at tackling the most challenging problems, SLCP and SIR, that require MCMC sampling to obtain the posterior samples on which the baselines would struggle. As getting accurate posterior samples with MCMC for increasing $N$ proved difficult, we restricted our evaluation to these four benchmark tasks.

[1] Chérief-Abdellatif & Alquier (2020). MMD-Bayes: Robust Bayesian estimation via maximum mean discrepancy.

[2] Bernton et al. (2019). Approximate Bayesian computation with the Wasserstein distance.

[3] Lueckmann et al. (2021). Benchmarking simulation-based inference.

Meta-review: The reviewers all recommend acceptance (6-6-6-6-6). The author-reviewer discussion has led to a number of remarks, questions and suggestions that can be used to improve the paper. We request the authors to address as much as possible the reviewers' comments in the final version.
Decision: Accept (poster)