Learning Robust Statistics for Simulation-based Inference under Model Misspecification
Authors: Daolang Huang
Ayush Bharti
Amauri H Souza
Luigi Acerbi
Samuel Kaski
Venue: NeurIPS 2023 poster (2023)
Keywords: Simulation-based inference
model misspecification
likelihood-free inference
approximate Bayesian computation
neural posterior estimation
TL;DR: We introduce the first general approach for tackling model misspecification in simulation-based inference that relies on learning robust statistics.

Abstract: Simulation-based inference (SBI) methods such as approximate Bayesian computation (ABC),  synthetic likelihood, and neural posterior estimation (NPE) rely on simulating statistics to infer parameters of intractable likelihood models. However, such methods are known to yield untrustworthy and misleading inference outcomes under model misspecification, thus hindering their widespread applicability. In this work, we propose the first general approach to handle model misspecification that works across different classes of SBI methods. Leveraging the fact that the choice of statistics determines the degree of misspecification in SBI, we introduce a regularized loss function that penalizes those statistics that increase the mismatch between the data and the model. Taking NPE and ABC as use cases, we demonstrate the superior performance of our method on high-dimensional time-series models that are artificially misspecified. We also apply our method to real data from the field of radio propagation where the model is known to be misspecified. We show empirically that the method yields robust inference in misspecified scenarios, whilst still being accurate when the model is well-specified.

Learning Robust Statisties for Simulation-basedd Inference under Model Misspecification
Daolang Huang"
Ayush Bharti' Aalto University
Aalto University
daolang huang@aalto.f1
ayush bharti@aaltor
Amauri H. Souza
Luigi Acerbi University of Helsinki luigi. ialgii.acerbi@eelsinkkk..
Samuel Kaski
Aalto University
Aalto University University of Manchester samuel kaski @aalto. fi
Federal Institute of Ceará amauri. souzaeaalto......o. fi
Abstract
Simulation-base inference (SBI) methods such as approximate Bayesian compu- tation (ABC), synthetic likelihood, and neural posterior estimation (NPE) rely on simulating statisties to infer parameters of intractable likelihood models. However, such methods are known to yield untrustworthy and misleading inference outcomes under model misspecification, thus hindering their widespread applicability. In this work, we propose the first general approach handle model misspecification that works across different classes of SBI methods. Leveraging the fact that the choice of statisties determines the degree of misspecification regularized loss function that penalizes those statistics that increase the mismatch between the data and the model. Taking NPE and ABC as use cases, we demon- strate the superior performance of our method on nigh-dimensional time- models that are artificially misspecified. We also apply our method real data from the field of radio propagation where the model is known to be misspecified. We show empirically that the method yields robust inference in misspecified scenarios, whilst still being accurate when the model is well-specified.
SBI, we introduce
1 Introduction
Bayesian inference traditionally entails characterizing the posterior distribution of parameters assum- ing that the observed data came from the chosen model family [7]. However, in practice, the true data-generating process rarely lies within the family of distributions defined by the model, leading to model misspecification. This can be caused by, among other things, measurement errors or contami- nation in the observed data that are not included in the model, or when the model fails to capture the true nature of the physical phenomenon under study. Misspecified scenarios are especially likely for simulator-basee models, where the goal is to describe some complex real-world phenomenon. Such simulators, also known as implicit generative models [25], have become prevalent many domains of science and engineering such as genetics [71], ecology [3], astrophysics [40], economics [29]. telecommunications [8], cognitive science [79], and agent-based modeling [87]. The growing field of simulation based inference (SBI) [21, 51] tackles inference for such intractable likelihood models, where the approach relies on forward simulations from the model, instead likelihood evaluations, to estimate the posterior distribution.
Traditional SBI methods include approximate Bayesian computation (ABC) [4. 51, 76], synthetic likelihood [67, 85] and minimum distance estimation [15]. More recently, the use of neural networks
Equal contribution.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

NPE
Our mothod
3000
Simulated Observed
Thue
Statistic 2
-200
Statistic 2
Ours
2000
20
NPE
-400
10
1000
-600
Simulated
*
Observed
-800
-2000 Statistic 1
-200
2.5
5.0
7.5
Statistic
01
Inference on misspecified Ricker model with two parameters, see Section 4 for details. Left: Learned statistics obtained using NPE; each point (in black) corresponds to a parameter value sampled from the prior. The observed statistic (in red) is outside the set of statisties the model can simulate. Middle: Statistics learned using our robust method; the observed statistic is inside the distribution of simulated statistics. Right: The resulting posteriors obtained from NPE and our method. The posterior obtained from our method is close to the true parameter value, while that from NPE is completely off, going outside the prior range (denoted by dashed grey lines).
Figure I:
as conditional density estimators [72] has fuelled the development of new SBI methods that target either the likelihood (neural likelihood estimation) [52, 61, 75, 83], the likelihood Akkelihood-to-evidene to-evidence ratio (neural ratio estimation) [22, 23, 28, 30, 39, 46, 77]. or the posterior (neural posterior estimation) [38, 41, 53, 58, 59, 70, 75, 83]. A common feature of all these SBI methods is the choice of summary statistics of the data, which readily impacts their performance [14]. The summary statistics are either manually handcrafted by domain experts, or learned automatically from simulated data using neural networks [1. 17, 18, 26, 31, 47, 84].
The unreliability of Bayesian posteriors under misspecification is exacerbated for SBI methods as the inference relies on simulations from the misspecified model. Moreover, for SBI methods using conditional density estimators [5, 13, 41, 61, 69], the resulting posteriors can be wildly inaccurate and may even go outside the prior range [9, 16, 35]; see [27] for an instance of this problem in the study of molecules. This behavior is exemplified in Figure for the posterior-targeting neural SBI method called neural posterior estimation (NPH). We see that when the model misspecified, the NPE posterior goes outside the prior range. This is because the model is unable to match the observed statistic for any value of the parameters, as shown in Figure (left), and hence, NPE forced to generalize outside its training distribution, yielding posterior that is outside the prior range.
We argue that the choice of statistics is crucial the notion of misspecification in SBI. Suppose we wish to fit Gaussian model to data from some skewed distribution. If we choose, for instance, the sample mean and the skewness as statistics, then the Gaussian model would fail to match these statistics for any choice of parameters, and we would end up in the same situation as in Figure I where the SBI method is forced to generalize outside its training distribution. However, choosing the sample mean and sample variance as statistics instead may solve this issue, as the model may be able to replicate the chosen statistics for some value of parameters, even though the Gaussian model itself is misspecified for any skewed dataset. This problem of identifying low-dimensional statistic of a high-dimensionnl dataset for which the model of interest is correctly specified is termed the data selection problem [82]. We refer to the statistics that solve the data selection problem as being robust.
Contributions. In this paper, we learn robust statistics using neural networks to solve the data selection problem. To that end, we introduce a regularized loss function that balances between learning statistics that are informative about the parameters, and penalizing those choices of statistics or features of the data that the model is unable replicate. By doing so, the observed statistic does not go outside the set of statistics that the model can simulate, see Figure 1mmiddle), thereby yielding a posterior around the true parameter even under misspecification (Figure (rright)). As our method relies on learning appropriate statistics and not on the subsequent inference procedure, it is applicable to all the statistics-based SBI methods, unlike other robust methods the literature [24, 48, 81]. We substantiate our claim in Section 4 through extensive numerical studies on NPE and ABC two fundamentally different SBI methods. Additionally in Section 5, we apply our method to a real case of model misspecification in the field of radio propagation [8].
2

---

[57]. In [24], the authors propose a robust SBI method based on Bayesian nonparametric learning and the posterior bootstrap. More recently, robustness to model misspecification has been tackled for neural likelihood [48] and posterior estimators [81]. We use the latter method, called robust neural posterior estimation (RNPE) [81], baseline for comparing our proposed approach. Our method Related works. Model misspecification has been studied in the context of different SBI methods such as ABC [9, 33, 35, 36. 74J. synthetic likelihood [32. 34, 57], minimum distance estimation [24], and neural estimators [48, 81]. It was first analyzed for ABC methods such as the regression- adjustment ABC [5] in [35], and a robust approach was proposed [33]. In [9], the authors proposed to handle model misspecification in ABC by incorporating the domain expert while selecting the
statistics. Generalized Bayesian inference [12, 43, 49], which is a popular framework for designing methods that are robust to misspecification, has also been linked to ABC [74] and synthetic likelihood
differs from these previous robust approaches in that it is applicable across distinct SBI methods.
2 Preliminaries
Simulation- inference (SBI). Consider a simulator (P& 0 E θ}. which is a family of distributions parameterized by 0 on data space X S Rd, We assume that Fo is intractable, but sampling independent and identically distributed (iid) data X1m (x(1) (....,(()))
P9 is straightforward for x(x) E X, S i S n. Given iid observed data yi.. Q, y()) E X, from some true data-generating process Q and prior distribution of 0, we are interested in estimating a 0" e O such that Pg- = Q or alternatively. in characterizing the posterior distribution p(0), yin).. Since X' is a high- -ighndimensional space in most practical cases, it is common practice to project both the simulated and the observed data onto a wiw-dmmnniioall space of summary statistics S via a mapping X - S, such that s r)(x) is the vector of simulated statistics, and Sobs = 7(yim) the observed statistic. In the following, we assume that this deterministic map 1j is not known a priori, but needs to be learned using a neural network, henceforth called the summary network T7, with trainable parameters 4/. We now introduce two different paradigms for learning the statistics.
Learning statisties in neural posterior estimation (NPE) framework. NPEs are a class of SBI methods which are becoming popular in astrophysics [40] and nuclear fusion [37, 80), among other fields. They map each statistic vector s to an estimate of the posterior p(f]s) using conditional density estimators such as normalizing flows [41, 53, 59, 69]. The posterior is assumed member of a family of distributions qu. parameterized by v. NPEs map the statistics to distribution parameters via an inference network hig where & constitutes the weights and biases of h. such that (s) (0) s p())s). The inference network is trained on dataset {0,, ((),,s))mm, simulated from the model using prior samples { p(0), by minimising the loss (()) = - Ep(0,,) [log (s) ( Once trained, the posterior estimate is then obtained by simply passing the observed statistic Sobs through the trained network. Hence, NPEs enjoy the benefit of amortization, wherein inference on new observed dataset is straightforward after computationally expensive training phase. In cases where the statistics are not known a priori, the NPE framework allows for joint training of both the summary and the inference networks by minimizing the loss function
be
CNPE(O, v) = -Ep(0 -...).....) [log (ha(oo
(1)
on the training dataset {(0,, X xi:m,)) ) [69].
Even though NPEs are flexible and efficient SBI methods, they are known to perform poorly when the model is misspecified [16, 81J. This is because the observed statistic Sobs becomes an out-of-distribution sample under misspecification (see Figure/ for an example), forcing the inference network in NPE to generalize outside its training distribution.
Learning statistics using autoencoders. Unlike NPEs, statistics for other SBI methods are learned prior to carrying out the inference procedure. This is achieved, for instance, by training an autoencoder with the reconstruction loss [1]
CAE(M, ba) = Ep(o.x...) [(x1m - Nea (no (x1 ))
(2)
where so and tha are the parameters of the encoder and the decoder ij. respectively. The trained encoder Th is then taken as the summarizing function in the SBI method. In this paper, we will use the encoder my to perform inference in an ABC framework, which now recall.
3

---

Approximate Bayesian computation (ABC), ABC is arguably the most popular SBI method, and widely used in many research domains [64]. ABC relies on computing the distance between the simulated and the observed statistics to obtain samples from the approximate posterior of a simulator. Given a discrepancy :(::)) and a tolerance threshold 6, the basic rejection-ABC method involves repeating the following algorithm: (i) sample 6' ~ p(f), (ii) simulate x Jmm. ~ Pe., (iii) compute 1 = n7(x),)). and (iv) if p(s', Sobba) < 6, accept 0', until a sample {0) is obtained from the approximate posterior. Conditional density estimators have also been used in ABC to improve the
posterior approximation of the rejection-ABC method [5, 11, 13, 14]. These so-called regression adjustment ABC methods involve fitting function 9 between the accepted parameters and statistics pairs {68) V/ as 8 = g($)) + W/, 1 S i S ngg where { are the residuals. Once fitted, the accepted parameter samples are then adjusted as 0; = (Sobba) - ai where g(s) is the estimate of Ei0 s] and @; is the empirical residual. Despite the popularity of the regression adjustment ABC methods, they have also been shown to behave wildly under misspecification [9, 35], similar to NPEs.
3 Methodology
Model misspecification in SBI. Bayesian inference assumes that there exists a 6* E 2 such that Po. = Q Under model misspecification, Pg * Q for any 0 E O, which can lead to unreliable predictive distributions. However, this definition does not apply for models whose inference is carried out in light of certain summary statistics. For instance, a Gaussian model is clearly misspecified for a bimodal data sample. However, if inference is based on, say, the sample mean and the sample variance of the data, the model may still match the observed statisties accurately. Hence, even when a simulator is misspecified with respect to the true data-generating process, it may still be well-specified with respect to the observed statistic. The choice of statistics therefore dictates the notion of misspecificatioon in SBI. We propose an analogous definition for misspecification in SBI:
Definition 3.1 (Misseecification of model-summarizer pair in SBI). Let q4 P be a pushforward of the probability measure P2 Ps x x Po X under q : X - S, meaning that for A C S. x ... x Po(T (o())-((A)n)) Then, {π 0 E } is a set of distributions on the statistics space induced by the model for a given summarizer T3, and n Q" is the distribution of the statistics from the true data- senta-generating process. Then, the model- summarzer pair is misspecified if. AO E G T=2 + T4
n (A) 12 (m- (A)) = (P(nn
It is trivial see that if the model is well-specified in the data space X', it will also be well-specified in the statistic space S for any choice of n. The converse, however, is not true, as previously mentioned. Using this definition, we quantify the level of misspecification of model summarrizrr pair follows:
Definition 3.2 (Misspecification margin in SBI). Lei D be distance defined on the set of all Borel probability measures on S. Then, for a given 17. we define the misspecification margin En in SBI as
En
inf D (n 14 Q*)) 0
(3)
Note that the margin equal to zero for the well-specified case. Moreover, the larger the margin, the bigger the mismatch between the simulated and the observed statistics w.r.t rj. Our aim is to learn an 1) such that T is represented well whilst the momml-suumerrier pair no longer misspecified as per Definition 3. or alternatively, has zero margin as per Definition 3.2. We formulate this as constrained optimization problem in Equation (6).
Learning robust statistics for SBI. Suppose we have the flexibility to choose the summarizer y from a family of functions parameterized by v. Then, a potential approach for tackling model misspecification in SBI is to select the summarizer v) that minimizes the misspecification margin En, However, that involves computing the infimum for all possible choices of t. To avoid that, we minimize an upper bound on the margin instead, given
Enp < upper = Ep(o) [D (1 no Q*)]
(4)
Note that Equation (4) is valid for any choice of p(0), and we call upper the margin upper bound.
Minimizing the margin upper bound alone would lead to 1 being constant function, which when used for inference, would yield the prior distribution. Hence, there is trade-off between choosing an 17 that is informative about the parameters and an 77 that minimizes the margin upper bound. We

---

propose to navigate this trade-off by including the margin upper bound as a regularization term in the loss function for learning n7. To that end, let L(w, w)) be the loss function used learn By (along with additional parameters w) which, for instance, can be either LNPE from Equation (1) or CAE from Equation (2). Then, our proposed loss with robust statistics (RS) is written as
Lus(w, xb) C(w, .b) + By(a [D (n ne# (2"),
(5)
regularization
where 2
is the weight we place on the regularization term relative to the standard loss L. Minimizing the loss function in Equation (5) corresponds to solving the Lagrangian relaxation of the optimization problem:
min w.y
L(w, v))
(6)
s.t.
Ep(o) (D To++P"n+=**))) <> 0,
with À 2 0 being the Lagrangian multiplier fixed to constant value.
Estimating the margin upper bound. In order to implement the proposed loss from Equation (5),
need to estimate the margin upper bound using the training dataset {(0,, X1:n.5) ) _= ~ p(0. X1:a) sampled from the prior and the model. However, we only have one sample each from the distributions
70 7o i
I,...: : m, and one sample of the observed statistic from 70 Taking D to be the
Euclidean distance ||- Il. we can estimate the margin upper bound as 1 E Although easy to compute, we found this choice of 'D to yield overly conservative posteriors even in the well-specified case, see Appendix B.3 for the results. This is because the Euclidean distance is large even if only a handful of simulated statistics in the training data - are far from the observed statistic. As a result, the regularizer dominates the loss function while the standard loss term L related to 0 is not minimized, leading to underconfident posteriors. Hence, we need a distance DD that is robust to outliers and can be computed between the set (X1m,,)))) 1 and ny (YI:),, instead of computing the distance point-wise.
(y ).l
corresponding to different parameter values
To that end, we take to be the maximum mean discrepancy (MMD) [42], which is a notion of distance between probability distributions or datasets. MMD has number of attractive properties that make it suitable for our method: (i) it can be estimated efficiently using samples [15], (ii) is robust against few outliers (unlike the KL divergence), and (iii) it can be computed for unequal sizes of datasets. The MMD has been used in a number of frameworks such as ABC [8. 10, 50, 54, 62], minimum distance estimation [2, 15, 19, 56], generalized Bayesian inference [20, 57J, and Bayesian nonparametric learning [24]. Similar to our method, the MMD has previously been used as a regularization term to train Bayesian neural networks [66], and in the NPE framework detect model misspecification [73]. While we include the mmD-based regularizer to ensure that the observed statistic is not an out-of-distribut computing the MMD between the simulated statistics and samples from a standard Gaussian, thus ensuring that the learned statistics are jointly Gaussian. They then conduct goodness-of-fh test [42] to detect if the model is misspecified. Their method is complementary to ours, such that our method can be used once misspecification has been detected using [73].
sample in the summary space, the regularizer in [73] involves
Maximum mean discrepancy (MMD) as D. The MMD is kernel- distance between probability distributions, computed by mapping the distributions to a function space called the reproducing kernel Hilbert space (RKHS). Let He be the RKHS associated with the symmetric and positive definite function & S x S - RR called reproducing kernel [6], defined on the space of statistics S such that f(s) = (f, ((, ()),,, Vf E Hs (called the reproducing property). We denote the norm and inner product of H4 by 11H, and (., (.,))Hu, respectively. Any distribution IP can be mapped to Hs via its kernel-mean embedding pa,7 = J k(-, (ds) The MMD between two arbitrary probability measures P and Q on S is defined as the distance between their embeddings in Hy, i.c., MMD, (IP, Q)
| Mk,lll, [55]. Using the reproducing property, can express
the MqUared-MMD as
MM. Q E...- s')] - 2Es-Ps' s')] + Esss'~(k(ss,
which computationally convenient as the expectations of the kernel can be estimated using iid samples from I' and Q. Given samples of simulated statistics {np (X)m and the observed
5

---

statistic Th (V1:)) we estimate the squared-MMD between them using the V-statistic estimator [42]: MMD P0 Now Q*] =
(VI:n)]
(8)
2 Σ knn (x1:m,)) No (y1:)))
12 Σ knn (X1::s)) no (X1m..))
Note that the last term in Equation (7) is always constant for the estimator above as we only have one data-point of the observed statistic, hence we disregard it. Equation (8) therefore corresponds to estimating the MMD between the distribution of the simulated statisties for given v and Dirac
measure on Th (Yimn). ). The computational cost of estimating the squared-MMD is (1²) and its rate of convergence is (1- The NPE loss with robust statistics (NPE-RS) can then be estimated as
m
(sPF.-R((2, v) m
Σ
log (n (0:) + XMMD [ (9)
Similarly, the autoencoder loss with robust statistics (AE-RS) reads
CAE-RS(V- tha
Σ
(n (M) XMMD [In ((mo()))))) no (yimm)]
(10)
We use subset of the training dataset of size 1 < to compute the MMD instead of to avoid incurring additional computational cost in case m is large. The MMD-based regularizer can also be
used score in classification task to detect model misspecification, as shown in Appendix A.
The regularizer penalizes learning those statistics for which the observed statistic nh (Y is far from the set of statistics that the model can simulate given prior distribution. When tends to zero, maximization of the likelihood dictates learning of both á and eb in Equation (9), and our method converges to the NPE method. On the other hand, the regularization term minimized when the summary network outputs the same statistics for both simulated and observed data, i.e., when D is zero. In this case, the inference network can only rely on information from the prior p(0). As a result, for large values of X, we expect the regularization term to dominate the loss and the resulting posterior to converge to the prior distribution. Similar argument holds for the autoencoder loss in Equation (10). We empirically observe this behavior in Section 4. Hence, encodes the trade-off between efficiency and robustness of our inference method. Choosing À can be cast as a hyperparameter selection problem, for which we can leverage additional data (if available) as validation dataset, use post-hoe qualitative analysis of the posterior predictive distribution.
Role of the regularizer
4 Numerical experiments
We apply our method of learning robust statistics to two different SBI frameworks NPE [41] and ABC [5] (see Appendix B.6 for results on applying our method to neural likelihood estimator [52]). We also compare the performance of our method against RNPE [81]. by using the output of the trained summary network in NPE as statistics for RNPE. Experiments are conducted on synthetic data from two time-series models, namely the Ricker model from population ecology [67, 85] and the -hheenbe process [18]. Real data experiment
radio propagation model is presented in Section 5. Analysis of the computational cost of our method and results on a 10-dimensional Gaussian model is presented in Appendix B.4 and B.5, respectively. The code to reproduce experiments is available at https: / gi thub. com/ wwt.ccooocccceeee http//aaannnoaoo ang robust-sbi.
Ricker model simulates the evolution of population size N over the course of time i as No+1 exp(0:) N,, exp(-N, + et), ...... T, where exp(f;) the growth rate parameter, et are zero- mean iid Gaussian noise terms with variance o² and No- = 1. The observations are assumed to be Poisson random variables such that x, ~ Poiss (02N2 N.). For simplicity, we fix o = 0.09 and estimate 0 00 0,] using the prior distribution (12, 8] [0, 20)), and T 100 time-steps.
Orastein-Uhlenbeck process (OUP) is a stochastic differential equation model widely used in financial mathematics and evolutionary biology. The OU process 2; is defined as:
1t+1 = F; + A.re 1. 1= 1, 1,..., T,
Az,- = 8 (cxp((2) At + 0.5w.

---

Rocker (NPE)
OUP (NP))
Ricker (AUIG)
OUP (ABG)
15.0
Method
Method
NPS
ASC
ADC-RS jount)
7.5
5.0
2.5
1.2
1.22
1.0
NMD
0.75
0.6
0.8
0.75
0.30
0.4
0.4
0.50
0.23
0.25
0.00
0.0
0.0
0.00
064
10%
0%
10% 20%
6 10% 21%
0%
10% 200%
Figure 2: Performance of the SBI methods terms of RMSE and MMD for both the Ricker model and OUP. Each box represents the median and interquartile range (IQR), while the whiskers extend to the furthest points within 1.5 times the IQR from the edges of the box. For the well-specified case 0%), the proposed NPE-RS and ABC-RS methods perform similar their counterpart NPE and ABC, respectively. Under misspecification (e > 0%). NPE-RS RS and ABC-RS achieve lower RMSE and MMD values, demonstrating robustness model misspecification.
where 25, At 0.2, an 10, and w ~ N(O,, At). uniform prior U((0 x [-2, placed on the parameters 0 = 02]
Contamination model. We test for robustness relative to the e-contamination model [65], similar to [24]. For both the Ricker model and the OUP, the observed data comes from the distribution Q = (1 e) + elPe.. where a large proportion of the data comes from the model with Orres while proportion of the data comes from the distribution Po... Hence, [0, 1] denotes the level of misspecification within both the models, with = 0 resulting in the well-specified case. We take fru (0.5, 1.0] 0-0.5, 1.0] for OUP, and ftue [4, 10] A. [4, 100] for the Ricker model.
Implementation details. We take 1000 samples for the training data and n 100 realizations both the observed and simulated data for each O. We set À using an additional dataset simulated from the models using Oerue- For the MMD, sss,s') = exp(-Ils s' 1/1/// 3²), and its lengthscale using the median heuristic 8 / /med/2 [42], where med denotes the median of the of squared two-norm distances ||s; sy l1 for all pairs - We use / = 200 samples of the simulated statistics to estimate the MMD. For the Ricker model, the summary network 7 is composed of ID convolutional layers, whereas for the OUP, 77 is combination of bidirectional long short-term memory (LSTM) recurrent modules and ID convolutional layers. The dimension of the statistic space models. We take to be a conditional normalizing flow for all the three NPE methods. We take the tolerance 6 in ABC be the top 5% samples that yield the smallest distance.
take the kernel to be the exponentiated-quadaall kernel
of distinct data points in {
four for both the
Performance metrics. We evaluate the accuracy of the posterior, as well as the posterior predictive distribution of the SBI methods. For the posterior distribution, emoor (RMSE) as (1/N) E. (0, Otrua) )) where (0, ) are posterior samples. For the predictive accuracy, compute the MMD between the observed data and samples from the posterior predictive distribution each method. As the models simulate high-dimensional data, we use the kernel spe- cialized for time-series from [8]. The lengthscale of this kernel is set to 3 1 for all misspecification levels to facilitate fair comparison.
compute the root mean squared
Figure 2 presents the results for both the Ricker model and the OUP across 100 runs. We observe that both NPE and ABC with our robust statistics (RS) outperform their counterparts, including RNPE, under misspecification (e 10% and 20%). Moreover, performance is similar
Results.

---

25
25
15
RNPE
20
20
NPE RS (OuIP)
15
15
6 10
5
6
01
(a) Well-specified (e 0)
(b) Misspocified (e == 10%)
(c) Misspecified (e = 20%)
Figure 3: Posteriors obtained from NPE, RNPE, and our NPE-RS method for the Ricker model. We perform similar to NPE in the well-specified case, unlike RNPE. Under misspecification, NPE and RNPE posteriors drift away from fuuee going even beyond the prior range (denoted by dashed gray lines) in the case of NPE. Our method is robust model misspecification.
to NPE in the well-specified case (e = 0%). An instance of the NPE posteriors is shown in Figure for the Ricker model. Our NPE-RS posterior is very similar to the NPE posterior for the well-specified case, whereas the RNPE posterior is underconfident, as noted in [81]. Although RNPE more robust than NPE under misspecification, its posterior still drifts away from Ourue. while NPE-RS posterior stays around ftree even for = 20%. This is because our method has the flexibility to choose appropriate statistics based on misspecification level, while RNPE is bound to fixed choice of pre-selected statistics. Posterior plots for ABC and OUP are given in Appendix B.2.
Prior misspecification. So far our discussion has per- tained to the case likelihood misspecification. How- ever, in Bayesian inference, the prior is also a part of the model. Hence, a potential form of misspecification is prior misspecification, in which the prior gives low or even zero probability the "true" data generating pa- rameters. We investigate the issue of prior misspecifi- cation using the Ricker model as an example. To that end, we modify the prior distribution of 02 by setting it to Normarr (0.5, 1.0) while keeping the prior 0, un- changed. We use frrue = [4, 25] T to generate the observed data, as 02 = 25 has very low density under the lognor- mal prior, and the probability that value of 02 2 25 will be sampled is 0.00372. The result in Figure 4 shows that NPE and RNPE yield posteriors that do not include the true value of 02, whereas our NPE-RS posterior is still around 0 Ourue- This highlights the effectiveness of our method to handle cases of prior misspecification well.
80
Brs
60
RNPE
NPE-RS
3 4.4
20
0:
Figure 4: Posteriors from NPE, RNPE, and NPE-RS for Ricker model. NPE-RS is robust to prior misspecification.
Varying regularizer X. As mentioned in Section 3, our method converges NPE as À tends to zero, and to the prior distribution as X becomes high. To investigate this effect, we vary À and measure the distance between the NPE-RS and the NPE posteriors, and between the NPE-RS posterior and the prior. We use the MMD from Equation (7) to measure this distance on the Ricker model in the misspecified case (e 20%). We see in Figure 5(left) that the MMD values between the prior and the NPE-RS posteriors go close to zero for X = 10 indicating that our method essentially returns the prior for high values of X. On the other hand, the MMD values between NPE and NPE-RS posteriors in Figure 5(middle) are smallest for X = 0.01 and À = 0.1. For non-extreme values of À (i.e. or 10), we observe maximum difference between NPE-RS and NPE. This because the NPE posteriors tend to move out of the prior support under misspecification, whereas the NPE-RS posteriors remain around Oerue, as shown in Figure 3. Finally, the NPE-RS posteriors are also close to the NPE posteriors for small values of À in the well-specifice case, as shown in Figure 5(right).
8

---

Dist. prior (misepeiifeed)
Dist. NPE postetior (misspocifaed)
Dist. NPE posterior (well-specified)
1.5
1.5
1.3
1.0
1.0
1.0
MMD
MMD
0.5
0.5
0.01
0.1
1.03
10.0
100.0
1000.0
0.01
0.1
1.0
30.0
300.0
1000.0
0.01
0.1
10
10.0
100.0
1000.0
Figure 5: MMD between the NPE-RS posteriors and (left) the prior in the misspecified setting, (middle) the NPE posteriors in the misspecified setting, (right) the NPE posteriors in the well- specified setting, for different values of X.
Application to real data: A radio propagation example
Our final experiment involves stochastic radio channel model, namely the Turin model [45, 63, 78), which is used to simulate radio propagation phenomena in order to test and design wireless communication systems. The model is driven by an underlying point process which is unobservable, thereby rendering its likelihood function intractable. The model simulates high- high-dimenniooall complex- valued time-series data, and has four parameters. The parameters govern the starting point, the slope of the decay, the rate of the point-process, and the noise floor of the time-series, as shown in Figure 6.
We attempt to fit this model to a real dataset from measurement campaign [44] for which the model is known to be misspecified [8], on account of the data samples being non-iid. The averaged power of the data (square of the abso- lute value) as function of time is shown by the
-80
starting point
Roal data
NPE
RNPE
NPE-RS (ours)
-100
slope
black curve in Figure 6. The data dimension is 801 and we have the power of the complex data in decibels" as input to the summary network, which consists of ID convolutional layers, and we set uniform priors for all the four parameters. Descriptions of the model, the data, and the prior ranges are provided in Appendix C.
Power
-110
100 realizations. We take
-120
-130
noise floor
50
100
150
200
Time Delay [ns]
Figure 6: Predictive performance of NPE, RNPE and our NPE-RS method on the Turin model under misspecification. NPE-RS leads to the best fit. The shaded regions denote = std. deviation of the posterior predictive distribution of each method.
We fit the model to the real data using NPE, RNPE and our NPE-RS methods. As there is no notion of ground truth in this case, we plot the resulting posterior predictive distributions in Figure 6. Note that the multiple peaks present in the data not replicated by the Turin model for any method, which is due to the model being misspecified for this dataset. Despite that, our NPE-RS method appears to fit the data well, while NPE performs the worst. Moreover, our method is better than RNPE at matching the starting point of the data, the slope, and the noise floor, with RNPE underestimating all three aspects average. The MMD between the observed data and the posterior predictive distribution of NPE, RNPE and NPE-RS is 0.11, 0.09, and 0.03, respectively. Hence, NPE-RS provides reasonable fit of the model even under misspecification.
Conclusion
We proposed a simple and elegant solution for tackling misspecification of simulator-basse models. As our method relies on learning robust statistics and not on the subsequent inference procedure, it is applicable to any SBI method that utilizes summary statistics. Apart from achieving robust inference under misspecified scenarios, the method performs reasonably well even when the model is well-specified. The proposed method only has one hyperparameter that encodes the trade-off between efficiency and robustness, which can be selected like other neural network hyperparameters, for instance, via validation set.
In wireless communications, the power of the signal is measured in decibels (dB). A power P in linear scale corresponds to 10 lo (P) dB.

---

Limitations and future work. A limitation of our method is the increased computational com- plexity due to the cost of estimating the MMD, which can be alleviated using the sample-efficiees MMD estimator from [10] or quasi-Monte Carlo points [56]. Moreover, as our method utilizes the observed statistic during the training procedure, the corresponding NPE is not amortized anymore a limitation we share with RNPE. Thus, working on robust NPE. methods which are still amortized (to some extent) is an interesting direction for future research. An obvious extension of the work could be to investigate if the ideas translate to likelihood suffer from similar problems as SBI if all (or some) of the candidate models are misspecified.
model selection methods [68, 76] which can
Acknowledgements
The authors would like thank Dr. Carl Gustafson and Prof. Fredrik Tufvesson (Lund Univer- sity) for providing the measurement data. We also thank Masha Naslidnyk, Dr. Baançois-Xvvies Briol, and Dr. Markus Heinonen for their useful comments and discussion. We acknowledge the computational resources provided by the Aalto IT Project from Computer Science IT. This work was supported by the Academy of Finland Flagship programme: Finnish Center for Artificial Intelligence FCAL SK was supported by the UKRI Turing Al World-Leading Researcher Fellowship, [EP/W002973/1).
References
[1] Albert, C., Ulzega, S., Ozdemir, F., Perez-Cruz, F., and Mira, A. (2022). Learning Summary Statistics for Bayesian Inference with Autoencoders. SciPost Phys. Core, 5:043.
[2] Alquier, P. and Gerber, M. (2021). Universal robust regression via maximum mean discrepancy. 200002000.008400
[3] Beaumont, M. A. (2010). Approximate Bayesian computation in evolution and ecology. Annual Review of Ecology, Evolution, and Systematics, 41(1)3399446.
[4] Beaumont, M. A. (2019). Approximate Bayesian computation. Annual Review of Statistics and Its Application, 6(1):379-400) 403.
[5] Beaumont, M. A., Zhang, W., and Balding, D. J. (2002). Approximate Bayesian computation in population genetics. Genetics, 162(4))2025-200) 2035.
[6] Berlinet, A. and Thomas-Agnan, C. (2004). Reproducing Kernel Hilbert Spaces in Probability and Statistics. Springer Ocience+Business Media, New York.
[7] Bernardo, S. (2000). Bayesian Theory. John Wiley Sons.
[8] Bharti, A., Briol, F.-X., and Pedersen, Т. (2022a). A general method for calibrating stochastic radio channel models with kernels. IEEE Transactions on Antennas and Propagation, 70(6):3986- 4001.
[9] Bharti, A., Filstroff, L., and Kaski, S. (2022b). Approximate Bayesian computation with domain expert in the loop. In International Conference on Machine Learning, volume 162, pages 1893- 1905.
[10] Bharti, A., Naslidnyk, M., Key, O., Kaski, S., and Briol, F.-.. (2023). Optimally-weighted estimators of the maximum mean discrepancy for likelihood-frre inference. 1117:
[11] Bi, J.. Shen, W., and Zhu, W. (2021). Random forest adjustment for approximate Bayesian computation. Journal of Computational and Graphical Statistics, 000):1-10.
[12] Bissiri, P. G., Holmes, C. C., and Walker, S. G. (2016). A general framework for updating belief distributions. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 78(5): 103-
[13] Blum, M. G. B. and François, O. (2010). Non-linear regression models for approximate Bayesian computation. Statistics and Computing, 20(1):63-73.
10

---

[14] Blum, M. G. B., Nunes, M. A., Prangle, D., and Sisson, S. A. (2013). A comparative review of dimension reduction methods in approximate Bayesian computation. Statistical Science. 28(2).
[15] Briol, F.-X., Barp, A., Duncan, B., and Girolami, M. (2019). Statistical inference for generative models with maximum mean discrepancy. arXiv:1906.0523
[16] Cannon, P., Ward, D., and Schmon, S. M. (2022). Investigating the impact of model misspecifi- cation in neural simulation-bare inference. arXiv preprint arXiv:2209. 01845.
[17] Chan, J., Perrone, V., Spence, J., Jenkins, P., Mathicson, S., and Song, Y. (2018). A likelihood- free inference framework for population genetic data using exchangeable neural networks. In Advances in Neural Information Processing Systems, volume 31.
[18] Chen, Y., Zhang, D., Gutmann, M. U., Courville, A., and Zhu, Z. (2021). Neural approximate sufficient statistics for implicit models. In International Conference on Learning Representations.
[19] Ahérief-Abdellatie B.-E. and Alquier, P. (2020). MMD-Bayes: Robust Bayesian estimation via maximum mean discrepancy. In Proceesings of the 2nd Symposium on Advances in Approximate Bayesian Inference, pages 1-21.
[20] Chérief-Ahdellatif,, B.-E. and Alquier, P. (2021). Finite sample properties of parametric MMD) estimation: robustness to misspecification and dependence. arXivv1912.05737.
[21] Cranmer, K., Brehmer, J., and Louppe. G. (2020). The frontier of simulation- inference. Proceedings of the National Academy of Sciences, 117(48):30055-30062
[22] Dax, M., Green, S. R., Gair, J., Deistler, M., Schölkopf, B., and Macke, J. H. (2022). Group equivariant neural posterior estimation. In International Conference on Learning Representations.
[23] Delaunoy, A., Hermans, J., Rozet, F., Wehenkel, A., and Louppe, G. (2022). Towards reliable simulation-bas based inference with balanced neural ratio estimation. arXiv preprint arXiv:2208. 13624.
(24) Dellaporta, C., Knoblauch, J., Damoulas, T., and Briol, F.-. (2022). Robust bayesian inference for simulator-based models via the mmd posterior bootstrap. In International Conference on Artificial Intelligence and Statistics, volume 151, pages 943-970.
[25] Diggle, P. J. and Gratton, R. J. (1984), Monte carlo methods of inference for implicit statistical models. Journal of the Royal Statistical Society: Series (Methodological), 46(2):193-212.
[26] Dinev, T. and Gutmann, M. U. (2018). Dynamic likelihood inference via ratio estimation (DIRE), arXiv preprint ar Xrx:v: 1810.0989.
[27] Dingeldein, L., Cossio, P., and Covino, R. (2023). Simulation based inference of single- molecule force spectroscopy. Machine Learning: Science and Technology, 4(2):025009.
[28] Durkan, C., M.,ra,, I., and Papamakarios, G. (2020). On contrastive learning for likelihood-free inference. In International Conference Machine Learning, volume 119, pages 2771-2781.
[29] Dyer, J., Cannon, P., Farmer, J. D., and Schmon, (2022a). Black-box bayesian inference for economic agent-based models. arXiv preprint arXiv:2202.006225
[30] Dyer, J., Cannon, P. W., and Schmon, S. M. (2022b). Amortised likelihood-free inference for expensive time-series simulators with signatured ratio estimation. In International Conference on Artificial Intelligence and Statistics, volume 151, pages 11131-11144.
[31] Fearnhead, P. and Prangle, D. (2012). Constructing summary statistics for approximate Bayesian computation: automatic approximate bayesian computation. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 74(3):419-4744
[32] Frazier, D. T. and Drovandi, C. (2021). Robust approximate bayesian inference with synthetic likelihood. Journal of Computational and Graphical Statistics, 30(4):958 97669696-.
[33] Frazier, D. T., Drovandi, C., and Loaiza-Maya, (2020a). Robust approximate bayesian computation: An adjustment approach. arXiv preprint wrXiv:2008.04099

---

[34] Frazier, D. T., Drovandi, C., and Nott, J. (2021). Synthetic likelihood in misspecified models: Consequences and corrections. arXiv preprint arXiv:210..0346.
[35] Frazier, D. T., Robert, C. P., and Rousseau, J. (2020b). Model misspecification in approximate bayesian computation: consequences and diagnostics. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(2):421
[36] Fujisawa, M., Teshima, T., Sato, I., and Sugiyama, M. (2021). y-abe: Outlier-robust approximate bayesian computation based on a robust divergence estimator. In International Conference on Artificial Intelligence and Statistics, volume 130, pages 1783-1791.
[37] Furia, C. S. and Churchill, R. M. (2022). Normalizing flows for likelihood free inference with fusion simulations. Plasma Physics and Controlled Fusion, 64(10)):04003.
[38] Geffner, T., Papamakarios, G., and Mnih, A. (2022). Score modeling for simulation-base inference. arXiv preprint arXiv:2209. 14249.
[39] Glöckler, M., Deistler, M., and Macke, J. (2022). Variational methods for simulation-basee inference. In International Conference on Learning Representations.
[40] Green, R. and Gair, J. (2021). Complete parameter inference for GW150914 using deep learning. Machine Learning: Science and Technology, 2(3):03LT0I.
[41] Greenberg, D., Nonnenmacher, M., and Macke, J. (2019). Automatic posterior transformation for likelihood--rre inference. In International Conference on Machine Learning, volume 97, pages 2404 2414.
[42] Gretton, A., Borgwardt, K., Rasch, M. J., and Scholkopf, B. (2012). A kernel two-sample test. Journal of Machine Learning Research, 13:723-773.
[43] Grünwald, (2012). The Safe Bayesian. In Lecture Notes Computer Science, pages 169-183. Springer Berlin Heidelberg.
[44] Gustafson, C., Bolin, D., and Tufvesson, F. (2016). Modeling the polarimetric mm-wave propagation channel using censored measurements. In 2016 Global Communications Conference. IEEE.
[45] Haneda, K., Järveläänen, J., Karttunen, A., Kyrö, M., and Putkonen, J. (2015). A statistical spatio-temporal radio channel model for large indoor environments 60 and 70 ghz. IEEE Transactions on Antennas and Propagation, 6366):269422744.
[46] Hermans, J., Begy, V., and Louppe, G. (2020). Likelihood-free MCMC with amortized approximate ratio estimators. International Conference on Machine Learning, volume 119, pages 4239 4248.
[47] Jiang, B., yu Wu, T., Zheng, C., and Wong, W. (2017). Learning summary statistic for approximate Bayesian computation via deep neural network. Statistica Sinica, page 1595-1618.
[48] Kelly, R. P., Nott, D. Frazier, D. T., Warne, D. J., and Drovandi, C. (2023). Misspecification- robust sequential neural likelihood. arXiv preprint arXiv:2301.13366
[49] Knoblauch, Jewson, J., and Damoulas, T. (2019). Generalized variational inference: Three arguments for deriving new posteriors. arXiv:1904.0200
[50] Legramanti, S., Durante, D., and Alquier, (2022). Concentration and robustness of discrepancy-based ABC via Rademacher complexity. arXiv:2206.06991
[51] Lintusaari, J., Gutmann, M. U., Dutta, R., Kaski, S., and Corander, J. (2017). Fundamentals and recent developments in approximate Bayesian computation. Systematic Biology, 66:66-82. 82.
[52] Lueckmann, J.-M., Bassetto, G., Karaletsos, T., and Macke, J. H. (2019). Likelihood-free infer- ence with emulator networks. Proceedings of The 1st Symposium on Advances in Approximate Bayesian Inference, volume 96, pages 32-53.
12

---

[53] Lueckmann, J.-M., Goncalves, P. J., Bassetto, G., Ócal, K., Nonnenmacher, M., and Macke, J. H. (2017). Flexible statistical inference for mechanistic models neural dynamics. In Advances in Neural Information Processing Systems, volume 30.
[54] Mitrovic, J., Sejdinovic, D., and Teh, Y.-W. (2016). DR-ABC: Approximate Bayesian computa- tion with kernel-based distribution regression. In Proceedings of the International Conference on Machine Learning, pages 1482-1491.
[55] Muandet, K., Fukumizu, K., Sriperumbudur, B., and Schölkopf, B. (2017). Kernel mean embedding of distributions: A review and beyond. Foundations and Trends@ in Machine Learning, 1011-
[56] Niu, Z., Meier, J., and Briol, F.-X. (2021). Discrepancy-based inference for intractable generative models using quasi-Monte Carlo. arXiv:2106.11561.
[57] Pacchiardi, L. and Dutta, R. (2021). Generalized bayesian likelihood-free inference using scoring rules estimators. arXiv preprint arXiv:2104.03889.
[58] Pacchiardi, L. and Dutta, R. (2022). Likelihood-free inference with generative neural networks via scoring rule minimization. arXiv preprint arXiv:2205.15784
[59] Papamakarios, G. and Murray, (2016). Fast e-free inference of simulation models with bayesian conditional density estimation. In International Conference on Neural Information Processing Systems, page 1036-1044.
[60] Papamakarios, G., Pavlakou, T., and Murray, I. (2017). Masked autoregressive flow for density estimation. Advances in neural information processing systems, 30.
[61] Papamakarios, G., Sterratt, D., and Murray, (2019). Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows. In International Conference on Artificial Intelligence and Statisties, volume 89, pages 837-848.
[62] Park, M., Jitkrittum, W., and Sejdinovic, D. (2015). K2-ABC: approximate Bayesian computa- tion with kernel embeddings. International Conference on Artificial Intelligence and Statistics,
[63] Pedersen, T. (2019). Stochastic multipath model for the in-room radio channel based on room electromagnetics. IEEE Transactions on Antennas and Propagation, 67(4):2591-2603.
[64] Pesonen, H., Simola, U., Köhn-Luque, A., Vuollekoski, H., Lai, X., Frigessi, A., Kaski, S., Frazier, D. T., Maneesoonthom,, W., Martin, G. M., and Corander, J. (2022). ABC of the future. International Statistical Review.
[65] Peter J. Huber, E. M. R. (2009). Robust Statistics. WILEY.
[66] Pomponi, J., Scardapane, S., and Uncini, A. (2021). Bayesian neural networks with maximum mean discrepancy regularization. Neurocomputing, 453:428-437,
[67] Price, L. F., Drovandi, C. C., Lee, A., and Nott, D. J. (2018). Bayesian synthetic likelihood. Journal of Computational and Graphical Statistics, 27(1):1-11
[68] Pudlo, P., Marin, J.-M., Estoup, A., Cornuet, J.-M., Gautier, M.. and Robert, C. P. (2015). Reliable ABC model choice via random forests. Bioinformatics, 32(6)885--866
[69] Radev, S. T., Mertens, U. K., Voss, A., Ardizzone, I.., and Köthe, U. (2022). Bayesflow: Learning complex stochastic models with invertible neural networks. IEEE Transactions on Neural Networks and Learning Systems, 33(4):1452-146)
[70] Ramesh, P., Lueckmann, J.- Boelts, J., Tejero-Canterd A., Greenberg, D. S., Goncalves. P. J., and Macke, J. H. (2022). GATSBI: Generative adversarial training for simulation-ban inference. International Conference on Learning Representations.
[71] Riesselman, A. J., Ingraham, B., and Marks, S. (2018). Deep generative models of genetic variation capture the effects of mutations. Nature Methods, 16-0022
13

---

[72] Rothfuss, J., Ferreira, F., Walther, S., and Ulrich, M. (2019). Conditional density estimation with neural networks: Best practices and benchmarks. arXiv preprint arXiv:1903.0095: 00954.
[73] Schmitt, M., Bürkner, P.-C., Köthe, U., and Radev, S. (2021). Detecting model misspecifica- tion in amortized bayesian inference with neural networks. arXiv e-prints.
[74] Schmon, S. M., Cannon, P. W., and Knoblauch, (2020). Generalized posteriors in approximate bayesian computation. arXiv preprint arXiv:201 L.0664.
[75] Sharrock, L.,, Simons, J., Liu, S., and Beaumont, M. (2022). Sequential neural score estima- tion: Likelihood-free inference with conditional score based diffusion models. arXiv preprint arXiv:2210.0487..
[76] Sisson, S. A. (2018). Handbook of Approximate Bayesian Computation. Hall/CRC.
Chapman and
[77] Thomas, O., Dutta, R., Corander, J., Kaski, S., and Gutmann, M. U. (2022). Likelihood-free inference by ratio estimation. Bayesian Analysis, 17(1).
[78] Turin, G. L.,, Clapp, F. D., Johnston, T. L.,, Fine, S. B., and Lavry, D. (1972). A statistical model urban multipath propagation. IEEE Transactions on Vehicular Technology, 21(1):1-9.
[79] van Opheusden, B., Acerbi, L., and Ma, W. J. (2020). Unbiased and efficient log- -ieeeikl estimation with inverse binomial sampling. PLoS Computational Biology, 10012)x
[80] Vasist, M., Rozet, F., Absil, O., Mollière, P., Nasedkin, E., and Louppe, G. (2023). Neural posterior estimation for exoplanetary atmospheric retrieval. Astronomy & AAstrophysis, 672:A147.
[81] Ward, D., Cannon, P., Beaumont, M., Fasiolo, M., and Schmon, S. M. (2022). Robust neural posterior estimation and statistical model criticism. In Advances in Neural Information Processing Systems.
[82] Weinstein, N. and Miller, J. W. (2023). Bayesian data selection. Journal of Machine Learning Research, 24(23):1-72.
[83] Wiqvist, S., Frellsen, J., and Picchini, U. (2021). Sequential neural posterior and likelihood approximation. arXiv preprint arXiv:2102.06522.
[84] Wiqvist, S., Mattei, P.-A., Picchini, U., and Frellsen, (2019). Partially exchangeable networks and architectures for learning summary statistics in approximate Bayesian computation. In International Conference on Machine Learning, volume 97, pages 6798-6807,
[85] Wood, S. N. (2010). Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466(7310): 100 102-1104.
[86] Zaheer, M., Kottur, S., Ravanbakhsh, S., Poczos, B., Salakhutdinov, R. R., and Smola, A. J. (2017). Deep sets. Advances in neural information processing systems, 30.
[87] Zbair, M., Qaffou, A., and Hilal, K. (2022). Approximate bayesian estimation of parameters of an agent-based model in epidemiology. In Lecture Notes in Networks and Systems, pages 302-31 314. Springer International Publishing.
14

---

Supplementary Materials
The appendix is organized follows:
In Appendix A, we present the results detecting model misspecification.
In Appendix B, we provide further implementation details and results for the numerical experiment in Section 4.
Appendix B.I: Implementation details
Appendix B.2: Additional posterior plots
Appendix B.3: Results for being the Euclidean distance
Appendix B.4: Computational cost analysis
Appendix B.5: Adversarial training on Gaussian linear model
Appendix B.6: Experiment with neural likelihood estimator
In Appendix C, we provide details of the radio propagation experiment of Section 5.
A Detecting misspecification of simulators
Considering that existing SBI methods can yield unreliable results under misspecification and that real-world simulators are probably not able to fully replicate observed data in most cases, detecting whether the simulator is misspecified becomes necessary for generating confidence in the results given by these methods. As misspecification can lead to observed statistics or features falling outside the distribution of training statisties, detecting for essentially boils down to class of distribbitiiino detection problems known as novelly detection, where the aim is to detect the test sample Sobs come from the training distribution induced by { This two-label classification problem can potentially be solved by adapting any of the numerous novelty detection methods from the literature. We propose the following two simple novelty detection techniques for detecting misspecification:
Distance- based approach. We assign score to the observed statistic based on the value of the margin upper bound, as introduced in the main text. We use the MMD as the choice of distance D, and estimate the MMD between the of simulated statistics [ss ] and the observed statistic 8obs- This MMD-based score can be used in classification method to detect misspecification.
Density-based approach. In this method, the training samples (8: 1m used fit generative model g. and the log- ligelikelihood of the observed statistics under are used as the classification score. We a Gaussian mixture model (GMM) with &s components as g, having the distribution
((s) = Σ (,p(s)u E,),
(11)
where E. nis and Σ, are the weight, the mean and the covariance matrix associated with the 7t component, and y denotes the Gaussian pdf. The score In g(Sobe) can then be used to classify it as either being from in or of the training distribution.
Experimental set-up. We test the performance of the proposed detection methods on the Ricker model and the OUP with the same contamination model as given in the main text. For each of these simulators, we first train the NPE method on m = 1000 training data points, and fit a GMM with A 2 components them. We then generate 1000 test datasets points, half of them from the well-specified model and the other half from the misspecified model, and compute their score. The area under the receiver operating characteristic (AUROC) is used the performance metric.
Baseline. We construct baseline for comparing performance of the proposed detection methods. The baseline is based on the insight that under model misspecification, the NPE posterior moves away from the parameter value (even going outside the prior range). Therefore, we take the root mean squared error (RMSE), defined (1/N ENn (0: )2) where (0.) J are posterior samples, as the classification score.
15

---

Ricker
OUP
1.0
1.0
0.9
AUROC
0.8
0.8
0.7
0.7
0.6
RMSE
0.5
GMM
0.5
MMD
0.4
0.01
0.15
0.01
0.05
0.15
Misspecification level
Misspecification level
Figure 7: Misspecification detection experiment. AUROC of the proposed detection methods (GMM and MMD) versus misspecification level for the Ricker model and the OUP. The RMSE-based baseline is shown in blue.
Results. The AUROC of the classifiers for different levels of misspecification (c in the main text) is shown in Fig. 7 for both the models. The proposed GMM-based detection method performs the best, followed by the MMD-based method. The RMSE-based baseline performs the worst at the classification task. We conclude that it is possible detect model misspecification in the space of summary statistics using simple to use novelty detection methods.
B
Additional details and results of the numerical experiments
B.1 Implementation details
We implement our NPE-RS models based publicly available implementations from https: / //githb.. comm/aacelab///ooooocoooaaaaaaaa sbb.. We use the NPE-C model [41] with Masked Autoregressive Flow (MAF) [60] as the backbone inference network, and adopt the default configuration with 50 hidden units and 5 transforms for MAF. The batch size is set to 50, and we maintain a fixed learning rate of 5 x 10-4 The implementation for RNPE is sourced directly from the original repository https: h////gtttbbb/ com/ foom/daniellin daniellwrd227 Sarerdd77/ngggg /rpp..
Regarding the summary network in NPE tasks, for the Ricker model, employ three ID convolu- tional layers with 4 hidden channels, and we set the kernel size to 3. For the OUP model, we combine three ID convolutional layers with one bidirectional LSTM layer. The convolutional layers have 8 hidden channels and a kernel size equal to 3, while the LSTM layer has 2 hidden dimensions. We pass the data separately through the convolutional layers and the LSTM layer and then concatenate the resulting representations to obtain our summary statistics. For the Turin model in Section 5, we utilize five ID convolutional layers with hidden units set to [8, 16. 32, 64, 8], and the kernel size is set to 3. Across all three summary networks, we employ the mean operation as our aggregator to ensure permutation invariance among realizations.
In ABC tasks, we incorporate autoencoders as our summary network. For the Ricker model, the encoder consists of three ID convolutional layers with 4 hidden channels, where the kernel size is set to 3. The decoder comprises of three ID transposed convolutional layers with the same settings as the encoder's 'sse convolutional layers, allowing for data reconstruction. For the OUP model, we adopt similar summary network as the one used for the Ricker model but with a smaller stride.
In NPE tasks, we use 1000 samples for the training data, along with 100 realizations of both observed and simulated data for each value O. We also use 1000 samples for training the autoencoders. For ABC, we use 4000 samples from the prior and accept We take p to be Euclidean distance the rejection ABC and normalize the statistics by the median absolute deviation before computing the distance to account for the difference in their magnitude.
200 samples giving tolerance rate of 5%.
B.2 Additional posterior plots
We now present examples of the remaining posterior plots, apart from the one shown in the main text. The posterior plots for OUP using the NaE-based methods is shown in Figure 8. The observations
16

---

are similar to the Ricker model example in the main text: we see that our NPE-RS method yields similar posterior as NPE in the well-specified case, whereas RNPE posteriors are underconfident. When the model is misspecified, NPE posterior goes far from the true parameter value. The NPE-Rs posteriors, however, are still around Orue: demonstrating robustness to misspecification.
2
2
Orue
RNPE
-2
62
(a) Well-specified (e 0)
(b) Misspecified (e 10%)
(c) Misspecified (e 20%)
Figure 8: and NPE for different degrees of model misspecification.
process. Posteriors obtained from our method (NPE-RS), RNPE,
Similar behavior is observed in the ABC case for both the Ricker model and OUP in Figure 9 and Figure 10, respectively. The ABC posteriors go outside the prior range under misspecification, while ABC with our robust statistics yields posteriors closer to Iruee In Table 1, we report the sample mean and standard deviations for the results shown Figure of the main text.
25
20
Ours
15
15
10
10
5
0
91
(a) Well-specified (e 0)
(b) Misspecified (c 10%)
(c) Misspecified (e = 20%)
Figure 9: Ricker model. Posteriors obtained from our method (ABC-RS) and ABC for different degrees model misspecification.
B.3 Results for D being the Euclidean distance
We present results for D being the Euclidean distance in the well-specified case of the Ricker model Figure [a). As mentioned in Section of the main text, this choice leads to very underconfident posteriors. This is because the Euclidean distance is not robust distance: it becomes large even a few points are far from the observed statistic. As a result, using this as the regularization term penalizes most choices of summarizer and we lear statistics that are very concentrated around the observed statistic (orange dot). Although a good choice for being robust, Euclidean distance leads to statistics that are not informative about the model parameters, yielding posterior that is similar to the uniform prior. Hence, we used the MMD as the distance the margin upper bound, which provides better trade-off between robustness and efficiency (ún terms of learning about model parameters).
17

---

15
2
Brue
ARC
1
Dars
10
5
0
2
0
(a) Well-specified (6 = 0)
(b) Misspocified (e = 10%)
(c) Misspecified (e = 20%)
Figure 10: Ornstein-Uhlenbed process, Posteriors obtained from our method (ABC-RS) and ABC for different degrees of model misspecification.
Table Performance of the SBI methods in terms of RMSE and MMD for both Ricker and OUP, We report the average (+1 std. deviation) values across 100 runs for varying levels of misspecification.
RMSE (4)
MMD (J)
6 0%
10%
6 20%
< 0%
<= 10%
6 20%
NPE
2.16 (3.07)
7.86 (1.57)
11.2 (1.70)
0.04 (0.07)
0.74 (0.09)
1.06 (0.17)
Ricker
RNPE
3.27 (0.35)
5.51 (0.58)
7.14 (1.15)
0.06 (0.05)
0.51 (0.19)
0.79 (0.25)
NPE-RS (ours)
2.18 (2.66)
2.19 (1.01)
4.66 (4.15)
0.09 (0.14)
0.21 (0.16)
0.42 (0.37)
ABC
1.46 (0.44)
6.95 (0.25)
9.79 (0.96)
0.01 (0.01)
0.85 (0.02)
1.18 (0.04)
ABC-RS (ours)
1.20 (0.51)
3.16 (1.08)
(1.28)
0.01 (0.02)
0.17 (0.15)
0.18 (0.16)
NPE
0.79 (0.62)
1.26 (1.18)
2.59 (2.75)
0.01 (0.01)
0.34 (0.15)
0.63 (0.29)
RNPE
0.78 (0.09)
0.87 (0.10)
0.98 (0.15)
0.01 (0.01)
0.22 (0.13)
0.49 (0.26)
OUP
NPE-RS (ours)
0,74 (0.70)
0.62 (0.33)
0.63 (0.36)
0.02 (0.05)
0.09 (0.09)
0.21 (0.17)
ABC
0.50 (0.07)
1.20 (0.40)
5.16 (2.39)
0.05 (0.03)
0.88 (0.21)
0.92 (0.23)
ABC-RS (ours)
0.44 (0.06)
0.62 (0.23)
0.88 (0.48)
0.02 (0.02)
0.26 (0.17)
0.50 (0.38)
B.4
Computational cost analysis
We now present a quantitative analysis of the computational cost of training with and without our MMD regularization term. The results, presented in Table 2, are calculated on an Apple M1 Pro CPU. As expected, we observe a higher runtime for our method due to the computational cost of estimating the MMD from 200 samples of simulated data. The total runtime also depends on the number of batchsize Nuach hence, as Nuuth increases, the proportion of runtime used for estimating MMD reduces. As a result, we see that for large Noach. the increase in the computational cost of our method with robust statistics is not significant.
B.5 Adversarial training on Gaussian linear model
To verify the robustness of our method higher dimensional parameter space, we an experiment on the Gaussian linear model, where the data x E IR 10 is sampled from (0, 1.1 A uniform prior 1]) is placed on the parameters 0 e 10 We take 0uuee = [0.5, ..., 0.5] 0g = [2. ..... 2] To introduce contamination to the observed data, employ the same approach outlined in the main text of our paper. However, there is a slight divergence in our experimental setup. In this example, we employ adversarial training, meaning that the model is trained on observed data with high degree of misspecification (c = 20%). while we perform inference on data with a lower misspecification degree (c 10%). For the summary network, comprises two linear layers, each with a width of 20 hidden units, paired with the ReLU activation function. The decoder is constructed with single linear layer of 20 hidden units.
utilize the DeepSet [86] architecture. The encoder
The results are shown in Table 3. NPE-RS outperforms NPE in terms MMD between the posterior predictive distribution and the observed data, which highlights the effectiveness of our approach in
18

---

20
One NPE (Fuclidean)
15
4
2
4
6
B
(a) Posteriors
(b) Summary statistics
Figure 11: Ricker model. Posteriors and summary statistics for D being the Euclidean distance.
Table 2: Comparison of computational costs across different models on Ricker model. We report the mean value (standard deviation) derived from 20 updates. We use different batch size Nbach and generate 100 realizations for each 0.
Runtime (seconds)
Nbuch 50
Nbach 100
Nbwch = 200
NPE
0.22 (0.03)
0.46 (0.04)
0.87 (0.03)
NPE-RS (ours)
1.26 (0.05)
1.53 (0.14)
1.92 (0.10)
ABC
(0.04)
1.41 (0.04)
3,29 (0.27)
ABC-RS (ours)
1.79 (0.04)
2.71 (0.25)
4,25 (0.46)
high- -ich-dimensional parameter spaces, even though the observed data was not used in the training of NPE-RS. This points towards the possibility that by employing adversarial training, we might achieve robustness against lower levels of misspecification whilst still being amortized.
Table 3: Performance comparison between NPE and NPE-RS for Gaussian linear model. We use MMD between the posterior predictive distribution and the observed data as metric. We report the average (+1 std. deviation) values across 100 runs.
NPE
NPE-RS
A
20
50
100
MMD
0.26 (0.02)
0.19 (0.04)
0.18 (0.06)
0.21 (0.08)
B.6 Experiment with neural likelihoud estimators
In this section, explore the performance of method when paired with Neural Likelihood Estimators (NL.)). NLE are a class of methods that leverage neural density estimators to directly estimate likelihood functions, bridging the gap between simulators and statistical models.
For this experiment, we adopt NLE-A as proposed by [61]. The original implementation can be found at https //github. com/mace httm//cckelabbbeee Similar method to NLE by pre-emptively training an autoencoder with our regularization term to learn the summary statistics. We refer network and simulator are consistent with those described in Appendix B.
approach with ABC, adapt
adapted method as NLE-RS. The configurations for our summary
Figure 12 presents the posterior plots for the Ricker model using the NLE-based methods. Consis- tent with our observations in the previous experiments, NLE-RS still demonstrates robustness to misspecification, while the NLE posterior tends to deviate away from the true parameters.
19

---

30
Ours NLE
20
10
2.5
5.0
7.5
01
Figure 12: Ricker model. Posteriors obtained from our method (NLE-RS) and NLE. We set 10% for this experiment.
C Details of the radio propagation experiment
In this section, we describe the data and the Turin model used in Section 5 of the main text.
Data and model description. Let B be the frequency bandwidth used to measure radio channel data at K equidistant points, leading to a frequency separation of Aʃ = B/ - 1). The measured transfer function at kth point, Yg, is modelled
Y H Wk, k == 1,,,..... 1........ K -,,,
where Hg is the transfer function the keth frequency, and We is additive zero-mean complex circular symmetric Gaussian noise with variance oW- Taking the inverse Fourier transform, the time-domain signal (t) can be obtained
K-1 Σ Y. exxp((22kkA()). K=0
v(t)
K
The Turin model defines the transfer function H3 = expp(((22 time-delay and on is the complex gain of the 1th component. The arrival time of the delays is modeled as one-dimensional homogeneous Poisson point processes, i.e., 77 ~ PPP(R+, v). with O. The gains conditioned the delays modeled as iid zero-mean complex Gaussian random variables with conditional variance || = Go exp(-7//T) The parameters of the model are 0 = [Go, T, v, T...T. The prior ranges used for the parameters are given in Table
where тi is the
Table 4: Prior distributions for the parameters of the Turin model.
Go
T
Prior U(10- 9,000 8) u(10 9 10 80-8) u(10", 5 x 109) 4(10 10 9)
The radio channel data from [44] is collected in a small conference room of dimensions 3 4 x m² using vector network analyzer. The measurement was performed with bandwidth of B 4 GHz, and 801. Denote each complex-valued time-series by ý E R and the whole dataset by Yi:, where 100 realizations. We take the input to the summary network to be Yı:- 10 1og
Scatter-plot of learned statistics. In Figure 13 and Figure 14, we show the scatter- ocatter-plots of the learned statisties using the NPE and our NPE-RS method, respectively. We observe that the observed statistics (shown in orange) is often outside the set of simulated statistics (shown in blue) for the NPE method. Hence, the inference network forced to generalize outside its training distribution, which leads to poor fit of the model, as shown in Section 5 of the main text. On the other hand, the observed statistic is always inside the set simulated statistics (or the training distribution) for our method in Figure 14, which leads to robustness against model misspecification.
20

---

5
5
52
54
56
Pairwise scatter-plots summary statistics learned using NPE method for the Turin model. Each blue dot corresponds to simulated statistic obtained from a parameter value sampled from the prior. The orange dot represents the observed statistic.
Figure 13:
21

---

59
51
5t
52
5g
54
56
Figure 14: Pairwise scatter-plots of summary statistics learned using our NPE-RS method for the Turin model. Each blue dot corresponds to simulated statistic obtained from a parameter value sampled from the prior. The orange dot represents the observed statistic.
22


Reviews:
[{"review_id": "eBUJVaiqQX", "reviewer": "Reviewer_DQRY", "summary": "The paper addresses the \"data selection problem\", i.e. identify a low-dimensional statistic (e.g. mean and variance for Gaussian distribution) of a high-dimensional dataset for which the model can replicate even when misspecified. The paper propose using the auto-encoding framework to automatically extract statistics and a penalty on mismatched statistics, i.e. statistics that the model is unable to replicate, as a general approach to handle model misspecification. Empirical results show robust inference in misspecified scenarios whilst still being accurate in well-specified scenarios.", "questions": "1. Since the extracted statistics are changing/ unknown beforehand, how will these statistics fit into downstream applications for downstream applications that were fixed beforehand? Perhaps the downstream tasks are retrained on the new robust statistics.\n2. In the experiment implementation, should the model be the misspecified case while the observations be the true case? Or both can have some degree of misspecification while we try to extract the true parameters?", "limitations": "1. Authors mentioned the limitation of using the observed statistic during the training procedure and recommended working on it as future work.", "rating": 4, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper introduces the \"data selection problem\". Sorry, I am not a subject matter expert in the field of \"data selection\" thus I cannot discuss on the originality and significance of this work. However, I learnt about the data selection problem from the paper.", "weaknesses": "1. The ideas of auto-encoding, reducing mismatch and tuning the regularizer are not new, but they are perhaps new in the application to \"data selection\". Please discuss further.\n2. From the paper, since it is an auto-encoding framework on the model, the statistics are informative of the model. On the other hand, how to ensure the statistics are informative of the observations? Please discuss further.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "* **\"The ideas of auto-encoding, reducing mismatch and tuning the regularizer are not new, but they are perhaps new in the application to \"data selection\". Please discuss further.\":** Yes, our contribution is to show how to solve the important problem of simulation-based inference under model misspecification by framing it as a data selection problem. As models in real world are often misspecified, addressing this problem is necessary for the application of the inference methods in many fields.\n\n* **\"From the paper, since it is an auto-encoding framework on the model, the statistics are informative of the model. On the other hand, how to ensure the statistics are informative of the observations?\":** Statistics are functions of the observations, and should be informative about the model parameters to perform inference. The first term in the loss function makes sure the statistics are informative, while the regularizer term ensures they are robust.\n    \n* **\"Since the extracted statistics are changing/ unknown beforehand, how will these statistics fit into downstream applications for downstream applications that were fixed beforehand? Perhaps the downstream tasks are retrained on the new robust statistics.\":** Yes, the new robust statistics can be used for the downstream tasks. It is a good point that if the downstream application is known beforehand, the statistics can be chosen according to that task. To do that, the first term in Equations 9 and 10 can be replaced by a term that minimises the downstream loss. Thank you for pointing this out; we will add a comment about this to the Discussion.\n    \n* **\"In the experiment implementation, should the model be the misspecified case while the observations be the true case?:**  A model is misspecified or well-specified with respect to a given set of observations, which in the simulated experiments we have chosen to either match or not match the modelling assumptions. So yes, the observations are always the ’true case’."}, {"review_id": "b2OPdr7uZ1", "reviewer": "Reviewer_EXCq", "summary": "This paper describes a method to perform simulation based inference, focusing on neural posterior estimation and approximate Bayesian computation (ABC), under model miss specification when performing inference using summary statistics of the data by using an MMD loss between the simulated and observed data summaries as a means to regulate miss specification.  Experiments are performed to empirically examine the efficacy of the method, and the method is compared to standard neural posterior estimation and ABC, as well as a method called robust neural posterior estimation.\n\n\n\n\n\n\n\n", "questions": "How does this work differ from ref [72] ? Can you provide a much clearer description of the differences is needed to understand the novelty of this work. Can you compare to this work?\n\nOn line 257, is is stated that lambda is set using simulations using theta_true. What does it mean that you set lambda using a data set with known theta_true? isn't this something that needs to be estimated? Doesn't this greatly reduce the challenge of inference, and this information is not realistically available? \n\nWhy does the summarizer act on on sets of iid data, rather than on each iid data example? Acting on each example individually seems to be much closer to practical usage. Do you have experiments in this setting? Would this significantly impact the computation driven by MMD, it seems this may make the method impractical?\n\nPlease provide more details on how RPNE was used in the experiments, and how the summary statistics were determined.\n\n", "limitations": "The authors adequately discuss limitations", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The paper is tackling an important issue, the method is interesting, and the description of the method and its motivation are made clear.\n", "weaknesses": "It is not clear how this work differs from ref [72], Schmitt et. al, \"Detecting Model Misspecification in Amortized Bayesian Inference with Neural Networks\". Much clearer description of the differences is needed to understand the novelty of this work\n\nIt is not clear how to determine the amount of regularization, i.e. how to set the hyperparameter lambda. This seems to be a key missing piece of information on how to practically use this method.\n\nThe paper seems to indicate that the summarizer acts on sets of iid samples. But is it not clear why must the summarizer act on the set of x's, rather than summarizing each observation, and using the fact that the summaries will remain iid. This seems to come back as a constraint later in line 179, but does not seem well motivated, or practically how people use summary statistics for iid data. This could also significantly impact the ability to train a conditional normalizing flow due to the massive reduction in information by summarizing over a set of examples rather that summarizing per example.\n\nIn the experiments, it is not clear what summary statistic is used for RNPE? Summary statistics are not described in the RNPE paper, but rather attempting to model the p(theta|x), so how was this choice made? Does it affect the results?\n\nMoreover, in the experiments, why not compare to the method in ref [72]?\n\nIn terms of reference, there are also methods to learn robust summaries prior to inference, e.g. using pivots, such as in Louppe, et. al, \"Learning to pivot with adversarial networks\", and similar domain adaptation approaches, some specifically using MMD.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "* **\"How does this work differ from ref [72] ? Can you provide a much clearer description of the differences is needed to understand the novelty of this work. Can you compare to this work?\":** See the general response to all the reviewers.\n    \n* **``It is not clear...how to set the hyperparameter lambda'':** We propose to set the hyperparameter $\\lambda$ either by using posterior predictive checks or via inference results on a held-out validation dataset (which can be a subset of the observed data). As shown in Figure 5, we found the inference results to not be very sensitive to changes in $\\lambda$ in various settings. We further argue that setting $\\lambda$ is similar to setting other hyperparameters such as the choice of architecture, number of layers, choice of activation, learning rate etc., which is now an accepted part of fitting any deep learning-based model. \n\n* **\"Why does the summarizer act on sets of iid data, rather than on each iid data example? Acting on each example individually seems to be much closer to practical usage. Do you have experiments in this setting? Would this significantly impact the computation driven by MMD, it seems this may make the method impractical?\":** We do not specify explicitly whether the summarizer acts on sets of iid data or on each iid data example. In fact, the latter is a special case of the former. The summarizer can act on each data example as well, however, that means there are as many summaries as number of data points, thus invoking the curse of dimensionality when computing distances in ABC. Moreover, computing statistics on the whole dataset can be more informative about parameters that govern the distributional behaviour of the data. We refer to [68] that define the summarizer in the same general way as us. For the experiments shown in the paper, the summary network we used first summarizes each data example, and then aggregates them into a single statistics vector for each $\\theta$, similar to [68].\n\n* **\"Summary statistics are not described in the RNPE paper, but rather attempting to model the $p(\\theta|\\mathbf x)$, so how was this choice made?\":** We respectfully disagree. The RNPE paper does use summary statistics as $\\mathbf{x}$, see the sentence before Section 2.2 in the RNPE paper that says \"Hereafter, we perform model criticism and inference using handcrafted summary statistics...\", and also Section 4.2 where they mention the statistics used for SIR and CS tasks.\n\n* **\"Please provide more details on how RPNE was used in the experiments, and how the summary statistics were determined.\":** The summary statistics learned using the joint NPE framework were the ones used for RNPE. That is, we used the output of the trained summary network in NPE as statistics of RNPE, thus making sure that both the methods used the same statistics so that comparing them was fair.\n\n* **\"In terms of reference, there are also methods to learn robust summaries prior to inference, e.g. using pivots, such as in Louppe, et. al, \"Learning to pivot with adversarial networks\", and similar domain adaptation approaches, some specifically using MMD.\":** Good point. We have made a note of it in the related works.\n\n* **\"What does it mean that you set lambda using a data set with known $\\theta_{\\mathrm{true}}$?\":** We do not assume the method knows the true theta; $\\theta_{\\mathrm{true}}$ simply is notation for the $\\theta$ used for simulating the observed and validation data. We use the validation set, generated from true $\\theta$, to set the value of $\\lambda$."}, {"review_id": "NoSfgjgD0q", "reviewer": "Reviewer_KhFo", "summary": "The authors propose a method to make neural posterior estimation robust to misspecification. The method relies on adding a regularizer to NPE loss function. The regularizer is implemented as the MMD between the embedding of the approximate (learned) posterior and the embedded observation. They apply their method to two low-dimensional benchmark tasks and a real-world example from radio propagation.", "questions": "The authors claim that RNPE is not amortized. However, as far as I understand, RNPE does not require retraining for new data (yes, it requires MCMC, but this can be very fast, especially for low-d parameter space). Please clarify if I misunderstood this or clarify in your paper.", "limitations": "State more explicitly that the method does **not** converge to the true posterior.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 4, "contribution": 3, "strengths": "**Originality**: \nThe paper tackles an important issue: how to make neural simulation-based inference robust to misspecification. The method is novel and intuitive.\n\n**Quality**: \nThe method is applied to benchmark tasks and to a real-world problem (with iid data samples and high-dimensional data).\n\n**Clarity**: \nThe paper is well-written and easy-to-follow. Figures are intuitive and well-designed.", "weaknesses": "My main concerns with the paper are that it (1) uses only two very low-dimensional benchmark tasks and (2) overstates it contributions.\n\n**Originality**:\n- the method is a very straightforward extension to Schmitt et al 2021 (which is only cited with a passing reference). The authors should clearly discuss how the works are related.\n\n**Quality**: \n- The authors repeatedly emphasize that their method performs on par with NPE for well-specified models. I do not think that their results support this claim though. NPE-RS performs significantly worse than NPE on well-specified data on one of the (only) two benchmark tasks. In order to make this claim, the authors should use significantly more benchmark tasks, ideally with varying parameter and data-dimensionality, and demonstrate that the behaviour shown in figure 2, left is a rare exception. In addition, the authors would have to show that NPE-RS converges (at least very closely) to the true posterior with many simulations (see below).\n\n- Convergence: The method proposed by the authors does no longer converge to the true posterior distribution for well-specified models, which further emphasizes that the method does not replace detection of misspecification. I would appreciate an investigation of how regularization strenght trades-off performance on well-specified data vs robustness to misspecification.\n\n- how is the hyperparameter chosen? The authors claim that the regularizer can be selected with a validation set (L353). It is unclear to me how this would work. What would be the loss function used to assess performance on the validation set? In L231 the authors also say that this might be done based on the posterior predictive distribution, but I think this is tricky because it (1) requires (potentially many) more simulations and (2) can easily lead to pathological cases where the posterior is off but the predictive distribution is good. Please elaborate on how the hyperparameter should be set in practice.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "* **\"The method is a very straightforward extension to Schmitt et al 2021.\":** See the general response to all the reviewers.\n\n* **\"... the paper uses only two very low-dimensional benchmark tasks \":** We respectfully disagree. While the tasks may be low dimensional in the number of parameters, they are very high dimensional in terms of data (100, 50, and 801 for Ricker, OUP and Turin model, respectively). Since misspecification occurs in the data space, we argue that the dimensionality of the data is more relevant to the problem under study. We would also like to point out that we added the arguably more important real-world case of misspecification (over additional synthetic benchmark examples).\n\n* **\"NPE-RS performs significantly worse than NPE on well-specified data on one of the (only) two benchmark tasks.\":** We disagree that NPE-RS performs *significantly* worse, and present additional results on different performance metrics in Figure 1 of the attached pdf. We see that the NPE-RS posterior is close to the NPE posterior in terms of MMD, and has comparable, if not better, empirical coverage.\n\n* **\"Convergence: The method proposed by the authors does no longer converge to the true posterior distribution for well-specified models, which further emphasizes that the method does not replace detection of misspecification.\":** Correct, our method does not converge to the true posterior. That is the price to pay to be robust to misspecification. What we claim is that our method achieves posterior consistency by leveraging Theorem 1 from Frazier et al. (2018), {Asymptotic Properties of ABC}, Biometrika, that says that as long as the statistics are informative about $\\theta$ (which is the case for our method), the resulting ABC posterior concentrates on the true theta in the well-specified case. Unfortunately, such a result does not exist for NPE; however, as shown in Figure 1 of the attached pdf, the NPE-RS posterior is close to the NPE posterior in the well-specified case. Nevertheless, we agree that our results do not explicitly show that detecting misspecification is unnecessary, and we will remove that claim. \n\n* **\"The authors would have to show that NPE-RS converges (at least very closely) to the true posterior with many simulations (see below).\":** We do not claim that NPE-RS converges to the true posterior, but rather to the NPE posterior as $\\lambda$ goes to 0, as shown in Figure 5(right).\n    \n* **\"I would appreciate an investigation of how regularization strength trades-off performance on well-specified data vs robustness to misspecification.\":** This is exactly what we investigated in Figure 5 by varying the value of $\\lambda$ in different settings.\n    \n* **\"How is the hyperparameter chosen? In L231 the authors also say that this might be done based on the posterior predictive distribution, but I think this is tricky because it (1) requires (potentially many) more simulations and (2) can easily lead to pathological cases where the posterior is off but the predictive distribution is good. Please elaborate on how the hyperparameter should be set in practice.\":** \n    Our method only has one hyperparameter $\\lambda$, which we propose to set either by using posterior predictive checks or via inference results on a held-out validation dataset (which can be a subset of the observed data), as it depends on the degree of misspecification in practice. As shown in Figure 5, we found the inference results to not be very sensitive to changes in $\\lambda$ in various settings. We further argue that setting $\\lambda$ is similar to setting other hyperparameters such as the choice of architecture, number of layers, choice of activation, learning rate etc., which is now an accepted part of fitting any deep learning-based model. \n\n    We argue that under model misspecification, it becomes difficult to say if the posterior is off or not (as there is no notion of a true posterior), as standard inference techniques are not reliable. We therefore rely on checking if the predictive distribution is accurate.\n\n* **\"The authors claim that the regularizer can be selected with a validation set (L353). It is unclear to me how this would work. What would be the loss function used to assess performance on the validation set?\":** We test the performance of the method on the validation set with different $\\lambda$ values. The performance metric can be any loss function, such as MMD (which is what we used) or KL divergence, between the model's predictive distribution and the validation dataset. The same metric can also be computed on the space of statistics instead of the data. \n    \n* **\"The authors claim that RNPE is not amortized. However, as far as I understand, RNPE does not require retraining for new data (yes, it requires MCMC, but this can be very fast, especially for low-d parameter space). Please clarify if I misunderstood this or clarify in your paper.\":** Agreed. Here, we refer to the amortization of the entire inference procedure, and not just the surrogate neural network, as noted in Box 1 of Lueckmann et al., Benchmarking Simulation-based Inference, AISTATS, 2021. Due to the additional MCMC step, the inference procedure is not amortized, which is also the case with neural likelihood and ratio estimators (NLE and NRE). We will clarify the text to reflect this distinction."}, {"review_id": "7f6Y3I2D0e", "reviewer": "Reviewer_1V4M", "summary": "This paper proposes a general approach to handling model misspecification for SBI. The paper introduces a regularized loss function that penalises mismatches between the observed data and the learned model. The paper focuses on NPE and ABC as the likelihood-free approaches. ", "questions": "* Could the authors provide more details about the step in Eq. (4) to go from infimum to upper bound? I found this challenging to understand.\n* The motivation behind defining a Q and a P that are different is a bit confusing. Why were these definitions highlighted? It seems that Q is equal to P when using the $\\theta_{true}$, but is it a requirement to know $\\theta_{true}$ prior to the experiments in order to simulate from Q?\n* The last sentence of section 5 highlights that the MMD is superior for the NPE-RS method. Is this surprising given it is built into the loss function? Are there other statistics that could be used here that are specific to the application domain?\n", "limitations": "This is well captured by the paper. This is appreciated!", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "* The radio propagation example is interesting, and shows what the potential strength of using the MMD regularizer. \n* The structure of the paper is clear.\n", "weaknesses": "\n* The novelty of the work is a concern in light of reference [72] equation (9), compared to equation (9) in this paper. Additionally, equation (10) appears similar to the InfoVAE (https://arxiv.org/abs/1706.02262), which also uses an auto-encoder with a MMD loss in the latent space. \n* A minor weakness is that the paper is difficult to follow. While the structure made sense, the descriptions of the approach and the motivation for it were challenging to understand. For example, it seemed implicit that this approach did not perform amortised NPE until it was highlighted this was the case in the conclusion. This is not in itself a criticism of the approach, but it was very hard to follow what consisted of the training data, and observations. For example, line 255 mentions using 1000 samples for the training data and 100 realisations of both the observed and simulated data for each $\\theta$. Under an amortised setting we would not have access to the observed $\\theta$. It is not clear what this setup is from the description. Further questions in the next section highlight some of these confusions that likely stem from descriptions of the approach that could be improved.\n* An additional minor weakness is only comparing to ABC and NPE. Is there any limitation on applying this approach to other estimators such as Neural Likelihood Estimators? Also for ABC, the $\\rho$ was only described in the supplementary materials as the Euclidean distance. After defining the MMD as the better regularizer, why would this not also be used/incorporated into the discrepancy?\n\nTypo: line 180, the $i$ should be subscript?\n\n* In Figure 5, it would be useful to ensure the three plots share the same y axis.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "* **\"The novelty of the work is a concern in light of reference [72]\":** See the general response to all the reviewers.\n\n* **\"equation (10) appears similar to the InfoVAE (https://arxiv.org/abs/1706.02262), which also uses an auto-encoder with a MMD loss in the latent space.\":** InfoVAE tackles the issue of learning meaningful latent features in variational autoencoders by modifying the ELBO objective, whereas we tackle the problem of learning robust statistics for simulation-based inference. Moreover, the MMD in InfoVAE is computed between variational distribution and prior distribution, while we compute the MMD between the simulated and observed statistics.\n    \n* **\"Is there any limitation on applying this approach to other estimators such as Neural Likelihood Estimators?\":** No, our method can be applied to Neural Likelihood Estimators (NLEs) in the same way as we applied it to ABC: by learning a summary function up front using neural networks. See Figure 2 in the attached pdf for NLE results under misspecification.\n    \n* **\"After defining the MMD as the better regularizer, why would this not also be used/incorporated into the (ABC) discrepancy?\":** MMD is used as the discrepancy in ABC to circumvent the need to summarize data into summary statistics, as mentioned in ref [61]. However, this requires the need to select an appropriate kernel function for the data, which is not always feasible. In most practical cases, summary statistics are used in ABC along with the Euclidean distance, which is what we used in our experiments. Euclidean distance is preferred here as the data is summarized into one statistic vector for each parameter value, while MMD is useful when computing distance between datasets.\n    \n* **\"line 255 mentions using 1000 samples for the training data and 100 realisations of both the observed and simulated data for each $\\theta$. Under an amortised setting we would not have access to the observed $\\theta$.\":** We are not sure we understand. Did you mean \"in the real-world setting we would not have access to the true $\\theta$\"? In that case, yes, that is the situation in the real-data experiment of Section 5. In section 4, we simulate observed data using the true $\\theta$ to check if we are robust to misspecification in the parameter space.\n    \n* **\"Could the authors provide more details about the step in Eq. (4) to go from infimum to upper bound?\":** Here we utilize the fact that the average value of a random variable (its mean) will be greater than or equal to its minimum value (infimum). \n    \n* **\"The motivation behind defining a Q and a P that are different is a bit confusing. Why were these definitions highlighted?\":** We introduce $P_{\\theta}$ and $Q$  to define model misspecification in SBI, see Section 2 of ref [24] for similar notation. \nHere $Q$ is the unknown true data-generating process, and $P_\\theta$ is the model we are trying to fit. We use $\\theta_{\\mathrm{true}}$  in experiments to simulate the observed data. In practice, $\\theta_{\\mathrm{true}}$ is not known and we would only have samples from $Q$, which is the case in Section 5. \n\n* **\"The last sentence of section 5 highlights that the MMD is superior for the NPE-RS method. Is this surprising given it is built into the loss function? Are there other statistics that could be used here that are specific to the application domain?\":** The MMD in Section 5 is used to test the performance of the different methods based on their predictive distribution, and uses a different kernel than the MMD in the loss function that summarizes the data, thus they are unrelated. Nevertheless, by using the KL divergence estimator instead of MMD between the observed data and predictions, we get the values 4.59, 5.81 and 2.30 for NPE, RNPE and NPE-RS, respectively, showing again that our NPE-RS method fits the data the best.\n    \n* **\"In Figure 5, it would be useful to ensure the three plots share the same y axis.\"** Agreed. We will edit the figure accordingly.\n* **\"Typo: line 180, the $i$ should be subscript?\":** Yes, thanks for the careful reading."}]

Author rebuttal:
We thank the reviewers for their careful consideration. The reviewers agreed on the following strengths of the paper:
* **Relevance:** Reviewers EXCq and KhF0 agree that the ''paper is tackling an important issue''
* **Good presentation:** EXCq: ''description of the method is clear''. KhF0: ''The paper is well-written and easy-to-follow. Figures are intuitive and well-designed''. 1V4M: ''structure of the paper is clear''.
* **Methodology:** KhF0: ''method is novel and intuitive''. EXCq: ''method is interesting''.
* **Real-world experiment:** KhF0: ''method is applied to real-world problem with high-dimensional data''. 1V4M: ''radio propagation example is interesting, and shows the potential strength of using the MMD regularizer''. 

The reviewers highlight concern regarding the novelty of the work in light of reference [72], Schmitt et. al, ``Detecting Model Misspecification in Amortized Bayesian Inference with Neural Networks'', which we clarify here. 
* **The purposes of the methods are different:** while we tackle the problem of robust inference under model misspecification, [72] addresses the problem of detecting model misspecification.
* **Role of the regularizer:** we include the regularizer to learn a summarizing function that ensures that the observed statistic is not an out-of-distribution sample in the summary space. On the other hand, [72] proposes to add an MMD regularizer term between the simulated statistics and samples from a standard Gaussian, which ensures that the learned statistics are jointly Gaussian.
* **How statistics are used:** we use the learned statistics to perform inference, while [72] conducts a goodness of fit test of Gretton et al, (2012) to detect if the model is misspecified.
* **The scope is different:** while their method is only applicable for NPE, our method can be used to perform robust inference in other SBI methods as well, such as ABC, NLE, NRE.

As [72] does not provide a method for robust inference, comparing our method with theirs is infeasible. Their method is complementary to ours, such that our method can be used once misspecification has been detected using [72]. This is the reason why in the paper we compared our method to RNPE, which addresses the same problem of robust inference as we do.

We also include additional results in the attached pdf document in response to some of the questions raised by reviewers KhFo and 1V4M.
We address all the individual comments of the reviewers in separate responses to each.

Meta-review: This paper an important problem in simulation based inference: how to handle misspecification since inferences are based on the assumption that the model is correct. All reviewers agree that the proposed method addresses this issue through the use of carefully constructed summary statistics. However, it is crucial that the authors adjust their paper to address the points made during the discussion phase.
Decision: Accept (poster)