Addressing Misspecification in Simulation-based Inference through Data-driven Calibration

## Abstract

Abstract Driven by steady progress in generative modeling, simulation-based inference (SBI) has enabled inference over stochastic simulators. However, recent work has demonstrated that model misspecification can harm SBI’s reliability.
This work introduces robust posterior estimation (ROPE), a framework that overcomes model misspecification with a small real-world calibration set of ground truth parameter measurements.
We formalize the misspecification gap as the solution of an optimal transport problem between learned representations of real-world and simulated observations.
Assuming the prior distribution over the parameters of interest is known and well-specified, our method offers a controllable balance between calibrated uncertainty and informative inference under all possible misspecifications of the simulator. Our empirical results on four synthetic tasks and two real-world problems demonstrate that ROPE outperforms baselines and consistently returns informative and calibrated credible intervals.

## 1 Introduction

Many fields of science and engineering have shifted in recent years from modeling real-world phenomena through a few equations to relying instead on highly complex computer simulations. While this shift has increased model versatility and the ability to explain complex phenomena,
it has also necessitated the development of new statistical inference methods. In particular, state-of-the-art simulation-based inference (SBI, Cranmer et al.,, ) algorithms leverage neural networks to learn surrogate models of the likelihood (Papamakarios et al.,, ), likelihood ratio (Hermans et al.,, ), or posterior distribution (Papamakarios and Murray,, ), from which one can extract confidence or credible intervals over the parameters of interest given an observation.
While SBI has proven helpful when the simulator is a faithful description of the studied phenomenon, e.g., for scientific applications (Delaunoy et al.,, ; Brehmer,, ; Lückmann,, ; Linhart et al.,, ; Hashemi et al.,, ; Tolley et al.,, ; Avecilla et al.,, ),
recent work has highlighted the unreliability of SBI methods under model misspecification (Cannon et al.,, ; Rodrigues et al.,, ).

##### Misspecification in SBI.

A model is a simplified description of a real-world phenomenon that allows reasoning about its properties. Thus, model misspecification arises in a context where the validity of simplifying assumptions cannot be totally verified. In the Bayesian inference literature (Walker,, ), a statistical model $p(\mathbf{x}_{s}\mid\theta)$ that relates a parameter of interest $\theta\in\Theta$ to a conditional distribution over simulated observations $\mathbf{x}_{s}$ is said to be misspecified if the true data generating process $p^{\star}(\mathbf{x}_{o})$ of the real observations $\mathbf{x}_{o}\sim p^{\star}(\mathbf{x}_{o})$ does not fall within the family of distributions defined by the model, i.e. $\nexists\theta\in\Theta:p^{\star}(\mathbf{x}_{o})=p(\mathbf{x}_{o}\mid\theta)$. Based on this definition, Ward et al., (); Huang et al., () have investigated solutions to improve the robustness of existing SBI algorithms
to misspecified models.
Here, we depart slightly from this definition of misspecification and instead tailor it to Bayesian uncertainty quantification, in particular, the prediction of credible regions on the parameter $\theta$ given an observation $\mathbf{x}_{o}$. We assume real-world data are sampled i.i.d. from the distribution $p^{\star}(\mathbf{x}_{o}):=\int p(\theta)p^{\star}(\mathbf{x}_{o}\mid\theta)d\theta$, where $p(\theta)$ is a known prior over the parameter we aim to infer and $p^{\star}(\mathbf{x}_{o}\mid\theta)$ is an unknown process approximated by a simulator that generates simulations $\mathbf{x}_{s}\sim p(\mathbf{x}_{s}\mid\theta)$ given a parameter $\theta$. We say the simulator is misspecified if there exists at least an observation $\mathbf{x}_{o}$ for which the simulator posterior $p(\theta\mid\mathbf{x}_{o})\neq p^{\star}(\theta\mid\mathbf{x}_{o})$, where $p^{\star}(\theta\mid\mathbf{x}_{o}):=\frac{p(\theta)p^{\star}(\mathbf{x}_{o}
\mid\theta)}{p^{\star}(\mathbf{x}_{o})}$. Notably, our definition can be interpreted as a generalized definition of miss-calibration and encompasses the definition of misspecification used in the Bayesian inference literature as a special case. While such miss-calibration often arises in practical settings where the simulator does not faithfully model some aspects of the true generating process, current SBI algorithms and their robust version may fail drastically.

##### Addressing Misspecification with a Calibration Set

In this work, we address model misspecification through a calibration set consisting of only a few pairs of real-world observations and their corresponding ground-truth labels. At first glance, having access to a calibration set may seem like an assumption that limits the applicability of such an approach. On the contrary, we argue that for any real-world task where the output of the inference is itself to be trusted, e.g., to monitor a parameter in a critical process, we will likely have access to a calibration set as obtained, for example, through a more costly procedure or measurement device. Indeed, unless we can confidently place additional assumptions on the form of the misspecification–a difficult task, since it is by definition unknown–the only way to validate the inference directly is through a validation set consisting of real-world observations and their ground-truth labels. If such a set exists, we can take a few observations to form our calibration set.

Figure: Figure 1: (left) Problem setup: we consider a real-world process which depends on some physical parameters $\theta$. Given real observations $\mathbf{x}_{o}$ of the process, our goal is to provide uncertainty quantification on the underlying parameters $\theta$. To help us, we have access to a misspecified simulator that takes parameters $\theta$ as input and produces simulated observations $\mathbf{x}_{s}$. (right) A visualization of ROPE. The training consists of two steps: (1) We approximate the posterior given the simulated data using NPE, resulting in the NSE $\mathbf{h}_{\omega^{\star}}$. (2) Using the calibration set, we fine-tune $\mathbf{h}_{\omega^{\star}}$ into $\mathbf{g}_{\varphi^{\star}}$ using the objective (7). At test time, we solve the OT problem between the representations $\{\mathbf{h}_{\omega^{\star}}(\mathbf{x}_{s}^{j})\}_{j=1}^{N_{s}}$ and $\{\mathbf{g}_{\varphi^{\star}}(\mathbf{x}_{o}^{i})\}_{i=1}^{N_{o}}$, resulting in our estimated posterior (6): the average of simulations’ posteriors weighted by the OT solution $P^{\star}$. See 21 for more details.
Refer to caption: x1.png

##### Our Contributions.

We introduce robust posterior estimation (ROPE), an algorithm that overcomes model misspecification to provide accurate uncertainty quantification for the parameters of stochastic and non-differentiable simulators.
To achieve this, ROPE relies on (1) a correctly defined prior, and (2) having access to a small, real-world calibration set of paired parameters and observations. The algorithm extends neural posterior estimation Papamakarios and Murray, () with optimal transport (OT, Peyré et al.,, ; Villani et al.,, ) to model the misspecification as an OT coupling between simulated and real-world data. We evaluate the performance of the algorithm on existing benchmarks from the SBI literature, and introduce four new benchmarks, of which two are synthetic and two come from real physical systems for which both labeled data and simulators are available. To the best of our knowledge, the latter constitute the first real-world benchmarks for SBI under misspecified models. We also perform experiments to explore the effect different calibration set sizes have on the performance of the algorithm, together with ablation studies to understand the impact of each of its components.

## 2 Background & Notations

In this section, we provide a short review of SBI and OT, as our method is at the intersection of these two fields. We start with some more fundamental definitions. We consider a simulator, implemented as a computer program $S:\mathbb{R}^{K}\times[0,1]\rightarrow\mathbb{R}^{D}$, that takes in physical parameters $\theta\in\Theta\subseteq\mathbb{R}^{K}$ and a random seed $\varepsilon\in[0,1]$ to generate measurements $\mathbf{x}_{s}\in\mathcal{X}_{s}\subseteq\mathbb{R}^{D}$.
The simulator is a simplified version of a real and unknown generative process $p^{\star}(\mathbf{x}_{o})$ that produces real-world observations $\mathbf{x}_{o}\in\mathcal{X}_{o}\subseteq\mathbb{R}^{D}$. We assume this process depends on parameters with the same physical meaning as the ones of the simulator and thus use the same notation $\theta$. Our goal is to estimate a well-calibrated and informative posterior distribution $p(\theta\mid\mathbf{x}^{i}_{o})$ for each observation in the test set $\mathbf{x}^{i}_{o}\in\mathcal{D}$, which reduces uncertainty compared to the prior distribution. As a remark, the most informative and calibrated posterior is the Bayesian posterior $p^{\star}(\theta\mid\mathbf{x}_{o})$ that corresponds to true generative process $p^{\star}(\mathbf{x}_{o}):=\int p^{\star}(\mathbf{x}_{o}\mid\theta)p(\theta)
\text{d}\theta$. To achieve our goal, we have access to (a) the misspecified simulator $S$ that embeds domain knowledge and approximates $p^{\star}(\mathbf{x}_{o}\mid\theta)$, (b) a small calibration set of labeled real-world observations $\mathcal{C}:=\{(\theta^{i},\mathbf{x}_{o}^{i})\}_{i=1}^{N_{c}}$, which enables data-driven correction of the simulator’s misspecification, and (c) a test set $\mathcal{D}:=\{\mathbf{x}_{o}^{i}\}_{i=1}^{N_{o}}$ of real-world observations arising from $p^{\star}(\mathbf{x}_{o})$ for which we want to estimate the posterior.

### 2.1 Simulation-based Inference (SBI)

Applying statistical inference to simulators is challenged by the absence of a tractable likelihood function (Cranmer et al.,, ). As a solution, SBI algorithms leverage modern machine learning methods to tackle inference in this likelihood-free setting (Lueckmann et al.,, ; Delaunoy et al.,, ; Glöckler et al.,, ). Among SBI algorithms, neural posterior estimation NPE (Papamakarios and Murray,, ; Lueckmann et al.,, ) is a broadly applicable method that trains a conditional density estimator of $p(\theta\mid\mathbf{x}_{s})$ from a dataset of parameter-simulation pairs. In this paper, we focus on making NPE robust to model misspecification.

NPE usually parametrizes the posterior with a neural conditional density estimator (NCDE), which is composed of (1) a neural statistic estimator (NSE), denoted by $\mathbf{h}_{\omega}:\mathcal{X}_{s}\rightarrow\mathbb{R}^{l}$, that compresses observations into $l\text{-dimensional}$ representations and, (2) a normalizing flow (NF, Papamakarios et al.,, ; Tabak and Vanden-Eijnden,, ) that parameterizes the posterior density as $p_{\phi}(\theta\mid\mathbf{h}_{\omega}(\mathbf{x}_{s}))$. The parameters $\phi$ and $\omega$ of the NCDE are trained with stochastic gradient ascent on the expected log-posterior probability, solving the following optimization problem

$$ $\displaystyle\phi^{\star},\omega^{\star}=$ $\displaystyle\arg\max_{\phi,\omega}\mathbb{E}_{\begin{subarray}{c}\theta\sim p (\theta)\\ \varepsilon\sim\mathcal{U}[0,1]\end{subarray}}\left[\log p_{\phi}(\theta\mid \mathbf{h}_{\omega}(S(\theta,\varepsilon)))\right],$ (1) $$

where $p(\theta)$ denotes a prior distribution over the parameters $\theta$.

Under the assumption that the class of functions represented by the NCDE contains the true posterior, solving (1) leads to a perfect surrogate $p_{\phi^{\star}}(\theta\mid\mathbf{h}_{\omega^{\star}}(\mathbf{x}_{s}))$ of the true posterior $p(\theta\mid\mathbf{x}_{s})$.
In that case, $\theta\perp\mathbf{x}_{s}\mid\mathbf{h}_{\omega^{\star}}(\mathbf{x}_{s})$, that is, the NSE $\mathbf{h}_{\omega^{\star}}$ is a sufficient statistic of $\mathbf{x}_{s}$ for the parameter $\theta$. In practice, we approach perfect training by generating a sufficiently large number of pairs $(\theta,\mathbf{x}_{s})$ and doing a search on the NCDE’s architecture and training hyperparameters. As we will see in Section 3, NPE and the sufficiency of the NSE are at the root of our algorithm. To simplify notation, we denote the NCDE learned with NPE as $\tilde{p}(\theta\mid\mathbf{x}_{s})$

### 2.2 Optimal Transport (OT)

As detailed in Section 3, our algorithm models the misspecification between simulations and real-world observations as an OT coupling. In particular, let $p(\mathbf{x}_{o})$ and $p(\mathbf{x}_{s})$ be two continuous probability measures in $\mathcal{X}_{o}$ and $\mathcal{X}_{s}$, respectively, and consider a cost function $c:\mathcal{X}_{o}\times\mathcal{X}_{s}\rightarrow\mathbb{R}$ that assigns a cost to each pair $(\mathbf{x}_{o},\mathbf{x}_{s})\in\mathcal{X}_{o}\times\mathcal{X}_{s}$. Then, we wish to find the OT coupling

$$ $\displaystyle\pi^{\star}=\arg\min_{\pi\in\Pi}\mathbb{E}_{\pi}[c(\mathbf{x}_{o} ,\mathbf{x}_{s})],$ (2) $$

where $\Pi$ is the set of joint probability measures on $\mathcal{X}_{o}\times\mathcal{X}_{s}$ whose marginals are $p(\mathbf{x}_{o})$ and $p(\mathbf{x}_{s})$.

In our setting, we only have access to a limited number $N_{o}$ of real-world observations $\{\mathbf{x}_{o}^{i}\}_{i=1}^{N_{o}}$, which we assume result from an unknown generative $p^{\star}(\mathbf{x}_{o})=\int p^{\star}(\mathbf{x}_{o}\mid\theta)p(\theta)
\text{d}\theta$. Thus, we solve the discrete and entropy-regularized version of (2). Namely, given a set $\{\mathbf{x}_{s}^{j}\}_{j=1}^{N_{s}}$ of simulated observations, we search for the doubly stochastic transport matrix $P^{\star}$, in the polytope of transport couplings
$\mathcal{B}\!:=\{P\in\mathbb{R}^{N_{o}\times N_{s}}_{+}:P\mathbf{1}_{N_{s}}=
\mathbf{1}_{N_{o}}/N_{o},P^{T}\mathbf{1}_{N_{o}}=\mathbf{1}_{N_{s}}/N_{s}\},$
that solves

$$ $\displaystyle P^{\star}=\arg\min_{P\in\mathcal{B}}\sum_{ij}^{N_{o}\times N_{s} }P_{ij}\left(c(\mathbf{x}_{o}^{i},\mathbf{x}_{s}^{j})+\gamma\log P_{ij}\right),$ (3) $$

where $\gamma>0$ can be seen as a hyperparameter that encourages entropic transport matrices. The entropy-regularized transport problem can be solved with the Sinkhorn algorithm (Cuturi,, ), which is a fixed-point iteration algorithm with an efficient GPU implementation. In our experiments, we rely on the OTT (Cuturi et al.,, ) implementation of Sinkhorn, which returns the transport coupling $P^{\star}$ given the cost matrix $C:C_{ij}=c(\mathbf{x}_{o}^{i},\mathbf{x}_{s}^{j})$ and the entropic regularization factor $\gamma$.

## 3 Modeling Misspecification with OT

We approach the problem of misspecification as a modeling task rather than seeing it as an issue of the inference algorithm. As our main modeling assumption, we assume that

$$ $\displaystyle\mathbf{x}_{o}\perp\theta\mid\mathbf{x}_{s},$ (4) $$

that is, given the simulated observations $\mathbf{x}_{s}$, the real observations $\mathbf{x}_{o}$ contain no additional information about the parameters $\theta$.

While this assumption allows us to express the misspecification independently from $\theta$, it might be violated in practice, as is the case for the two real physical systems we study in section 4. Due to the construction of our algorithm, violations do not prevent obtaining calibrated and informative posterior distributions. However, the assumption may limit uncertainty reduction by preventing our method from exploiting any signal in the real data beyond what is specified by the simulator. This can be a limiting factor for highly misspecified simulators, but if the simulator encodes phenomena that the practitioner believes are invariant across different application environments, the assumption can also prevent shortcut learning from the calibration data and benefit the generalization of the method. In Appendix A, we evaluate the method on real out-of-distribution data and demonstrate this property.

The OT coupling in (2) fulfills the definition of a joint distribution in $\mathcal{X}_{o}\times\mathcal{X}_{s}$, and together with our modelling assumption (4), we can express the posterior distribution for real-world observations as

$$ $\displaystyle p(\theta\mid\mathbf{x}_{o})=\int p(\theta\mid\mathbf{x}_{s})\pi^ {\star}(\mathbf{x}_{s}\mid\mathbf{x}_{o})d\mathbf{x}_{s},$ (5) $$

where the simulation posterior $p(\theta\mid\mathbf{x}_{s})$ can be approximated arbitrarily well with NPE (Papamakarios and Murray,, ), as NFs are universal density estimators of continuous distributions (Wehenkel and Louppe,, ).

Motivated by the factorization in (5), our algorithm computes a transport matrix $P^{\star}$ between the test set $\mathcal{D}$ and a set $\{\mathbf{x}^{j}_{s}\}_{j=1}^{N_{s}}$ of $N_{s}$ simulations generated by running the simulator on parameters from the prior $\theta^{j}\sim p(\theta)$. Thus, approximating (5), we estimate the posterior for real-world observations as a mixture of posteriors $\tilde{p}$ obtained with NPE, that is,

$$ $\displaystyle\tilde{p}(\theta\mid\mathbf{x}_{o}^{i}):=\sum_{j}^{N_{s}}\alpha_{ ij}\tilde{p}(\theta\mid\mathbf{x}^{j}_{s}),\text{ where }\alpha_{ij}=N_{o}P^{ \star}_{ij}.$ (6) $$

An interesting property of defining ROPE as in $\eqref{eq:real_posterior_def}$ is that, by design, the marginal posterior distribution over the test set, i.e., $\tilde{p}(\theta):=\int\tilde{p}(\theta\mid\mathbf{x}_{o})p^{\star}(\mathbf{x}
_{o})\text{d}\mathbf{x}_{o}$, converges to the prior distribution as the number of simulated observations $N_{s}$ approaches infinity, as expected from a well-estimated posterior distribution. A proof and further discussion of this self-calibration property is given in Appendix B.

### 3.1 Defining the OT Cost Function

In our context, an ideal coupling would assign simulations to real-world observations generated from the same parameter. Hence, we can express the corresponding ideal cost as $c(\mathbf{x}_{o},\mathbf{x}_{s})=c(\mathbf{h}_{o}(\mathbf{x}_{o}),\mathbf{h}_{
s}(\mathbf{x}_{s}))$, where $\mathbf{h}_{o}$ and $\mathbf{h}_{s}$ are any sufficient statistics for $\theta$ given $\mathbf{x}_{o}$ and $\mathbf{x}_{s}$, respectively.

As discussed in Appendix G, we can learn an approximated minimal sufficient statistic $\mathbf{h}_{\omega^{\star}}$ for the simulated observations with NPE.
Furthermore, as the simulator carries information about the true generative process and the calibration set is too small to learn a representation from real-world data only, it is reasonable to learn a sufficient statistic $\mathbf{h}_{o}$ for the real observations by fine-tuning $\mathbf{h}_{\omega^{\star}}$.
Denoting this new neural network as $\mathbf{g}_{\varphi}:\mathcal{X}_{o}\rightarrow\mathbb{R}^{l}$, the fine-tuning objective reads:

$$ $\displaystyle\mathcal{L}(\mathcal{C}):=\sum_{i=1}^{N_{c}}\lvert\mathbf{g}_{ \varphi}(\mathbf{x}_{o}^{i})-\mathbb{E}_{\varepsilon\sim\mathcal{U}[0,1]}[ \mathbf{h}_{\omega^{\star}}\left(S(\theta^{i},\varepsilon)\right)]\rvert_{2},$ (7) $$

where the expectation is approximated via a Monte-Carlo approximation. The training of $\mathbf{g}$ starts from the weights $\omega^{\star}$ and optimizes (7) with gradient descent. Optimizing (7) enforces, at least on the calibration set, that $\mathbf{g}$ and $\mathbf{h}$ are close in L2 norm when they correspond to the same parameter. Thus, we define the OT cost as
$c(\mathbf{x}_{o},\mathbf{x}_{s}):=\lvert\mathbf{g}_{\varphi^{\star}}(\mathbf{x
}_{o}^{i})-\mathbf{h}_{\omega^{\star}}(\mathbf{x}_{s})\rvert_{2}$

, where $\mathbf{g}_{\varphi^{\star}}$ is the NSE obtained after fine-tuning (7). Figure 1 depicts the main training and inference steps of ROPE.

### 3.2 Entropy’s Magic

The entropic regularization of OT not only enables fast computation of the transport coupling but also provides an effective control mechanism to balance the calibration of the posterior with its informativeness.
Indeed, for small entropic regularization, the estimated posteriors have low entropy and may be overconfident, as they are sparse mixtures of a few simulation posteriors $\tilde{p}(\theta\mid\mathbf{x}_{s}^{j})$. In contrast, for large values of $\gamma$ in (3), the coupling matrix becomes uniform and the corresponding posteriors tend to the prior, as $p(\theta\mid\mathbf{x}_{o})\approx\frac{1}{N_{s}}\sum_{j}^{N_{s}}\tilde{p}(
\theta\mid\mathbf{x}_{s}^{j})$ is a Monte-Carlo approximation of $\mathbb{E}_{p(\mathbf{x}_{s})}[\tilde{p}(\theta\mid\mathbf{x}_{s})]\approx p(\theta)$. Thus, the practitioner should optimize the hyper-parameter $\gamma$ to find the right trade-off between calibration, favored by higher $\gamma$, and informativeness, favored by lower $\gamma$, of the estimated posterior distributions.

## 4 Experiments

Figure: Figure 2: Results for our method (ROPE) and the competing baselines on six benchmark tasks. For each task, we show an example of the real observations ($\mathbf{x}_{o}$) and the observations produced by the misspecified simulator ($\mathbf{x}_{s}$). We show each method’s LPP and ACAUC metrics, as computed on a labeled test set of size 2000. We report the average metrics and $\pm 1$ std. deviation over three random draws of the test set and additional sources of randomness. In some instances, e.g., tasks A and E, the likelihood of J-NPE can be $-\infty$ and is not plotted. For readability of the LPP metric, we use a linear scale between the SBI and the Prior and a logarithmic scale for values below that.
Refer to caption: x2.png

Our experiments aim to (1) empirically validate the discussion in Section 3.2, and (2) illustrate settings in which our algorithm enables uncertainty quantification under model misspecification and small calibration datasets. The experiments comprise two existing benchmarks from the SBI literature, two synthetic benchmarks, and two new benchmarks from real physical systems for which both labeled data and simulators are available. To the best of our knowledge, the latter constitute the first real-world benchmarks for SBI under misspecified models. Altogether, the benchmarks represent various types of misspecification and parameter and observation space. We briefly describe each task and provide examples of real vs. simulated observations in Figure 2. Further details about the complete setup can be found in Appendix D.

##### Task A (synthetic): CS.

We reproduce the cancer and stromal cell development benchmark from Ward et al., (). The simulator emulates the development of cancer and stromal cells in a 2D environment as a function of three Poisson rate parameters $(\lambda_{c},\lambda_{p},\lambda_{d})$. The observations are vectors composed of the number of cancer and stromal cells and the mean and maximum distance between stromal cells and their nearest cancer cell. Synthetic misspecification is introduced by removing cancer cells that are too close to their generating parent.

##### Task B (synthetic): SIR.

We also use the stochastic epidemic model from Ward et al., (), which describes epidemic dynamics through the infection rate $\beta$ and recovery rate $\gamma$. Each observation is a vector composed of the mean, median, and maximum number of infections, the day of occurrence of the maximum number of infections, the day at which half the total number of infections was reached, and the mean auto-correlation (lag 1) of the infections. Misspecification is a delay in weekend infection counts, of which $5\%$ are added to the count of the following Monday.

##### Task C (synthetic): Pendulum.

The damped pendulum is a common benchmark to assess hybrid learning algorithms (Wehenkel et al.,, ), which jointly exploit domain knowledge and real-world data. The simulator generates the horizontal position of a friction-less pendulum given its fundamental frequency $\omega_{0}\in\mathbb{R}^{+}$ and amplitude $A\in\mathbb{R}^{+}$. Randomness enters the simulator through a random phase shift and white measurement noise. As misspecified “real-world” data, we simulate observations from a damped pendulum that takes friction into account.

##### Task D (synthetic): Hemodynamics.

Following Wehenkel et al., (), we define the task of inferring the stroke volume (SV) and the left ventricular ejection time (LVET) from normalized arterial pressure waveforms. The simulator is a PDE solver (Melis,, ) that produces an $8$-second time-series $\mathbf{x}_{s}$ sampled at $125$Hz. As synthetic misspecification, the simulator assumes all arteries have constant length, whereas this parameter varies in the “real-world” data.

##### Task E (real): Light Tunnel.

We employ one of the light tunnel datasets from Gamella et al., (). The tunnel is an elongated chamber with a controllable light source at one end, two linear polarizers mounted on rotating frames, and a camera. Our task consists of predicting the color setting of the light source ($(R,G,B)\in[0,255]^{3}$) and the dimming effect of the polarizers $\alpha\in[0,1]$ from the captured images. The simulator takes the parameters $\theta:=[R,G,B,\alpha]$ and produces an image consisting of a hexagon roughly the size of the light source, with a color equal to
$[\alpha R,\alpha G,\alpha B]$.

##### Task F (real): Wind Tunnel.

We employ one of the wind tunnel datasets from Gamella et al., (). The tunnel is a chamber with two controllable fans that push air through it, and barometers that measure air pressure at different locations. A hatch controls the area of an additional opening to the outside. The dataset is a collection of pressure curves that result from applying a short impulse to the intake fan power and measuring the change in air pressure inside the tunnel. Our inference task consists of predicting the hatch position, $\theta:=H\in[0,45]$ given a pressure curve. As a simulator model, we adapt the physical model given in Gamella et al., (, Model C3, Appendix D).

##### Metrics

We consider two different metrics to assess whether ROPE provides reliable and useful uncertainty quantification.
First, given a labeled test set $\{(\theta^{i},\mathbf{x}^{i}_{o})\}_{i=1}^{N}$, we compute the log-posterior probability (LPP) as

$$ $\displaystyle\text{LPP}:=\frac{1}{N}\sum_{i=1}^{N}\log\tilde{p}(\theta^{i}\mid \mathbf{x}_{o}^{i})\approx\mathbb{E}_{\begin{subarray}{c}p(\theta)\\ p(\mathbf{x}_{o}\mid\theta)\end{subarray}}\left[\log\tilde{p}(\mathbf{x}_{o} \mid\theta)\right].$ (8) $$

The LPP is an empirical estimation of the expected cross entropy between the true and estimated posterior; thus, for an infinite test set, it is only maximized by the true posterior. LPP characterizes the entropy reduction on the estimation of $\theta$ achieved by a posterior estimator $\tilde{p}$ when given one observation, on average, over the test set. Second, the average coverage AUC (ACAUC) indicates the average calibration of $K$ 1D credible intervals extracted from the estimated posteriors, i.e.,

$$ $\displaystyle\text{ACAUC}:=\frac{1}{KN}\sum_{j=1}^{K}\sum_{i=1}^{N}\int_{0}^{1 }\alpha-\mathbf{1}[\theta^{i}_{j}\in\Theta_{\tilde{p}(\theta_{j}\mid\mathbf{x} _{o}^{i})}(\alpha)]\text{d}\alpha,$ (9) $$

where $\Theta_{\tilde{p}(\theta_{j}\mid\mathbf{x}_{o}^{i})}(\alpha)$ denotes the credible interval for the $\text{j}^{\text{th}}$ dimension of the parameter $\theta$ at level $\alpha$.
Its value is positive (negative) if, on average over different credible levels, parameter dimensionality, and observations, the corresponding credible intervals are overconfident (underconfident). The ACAUC of a perfectly specified prior distribution is zero. The integral can be efficiently approximated, as described in Appendix F. For all experiments, we compute the LPP and ACAUC on labeled test set containing $2000$ pairs $(\theta,\mathbf{x}_{o})$.

##### Baselines

As a sanity check, we compare the performance of ROPE against four reference baselines: the prior $p(\theta)$, which amounts to the lower bound on the LPP for any calibrated posterior estimator; the SBI posterior on the simulated examples which, under the independence assumption $\mathbf{x}_{o}\perp\theta\mid\mathbf{x}_{s}$, is an upper bound on the LPP for ROPE; (NPE) a posterior estimator fitted to the simulated data and applied to the real data; and (J-NPE) a posterior estimator trained jointly on the pooled simulated and real observations. The latter two baselines represent some first approaches that a practitioner may consider. Additionally, we compare the performance of ROPE to MLP, which trains a neural network that predicts the mean and log-variance of a Gaussian posterior distribution by maximizing the calibration set log-likelihood. For the CS and SIR benchmarks, we additionally run Noisy NPE (NNPE, Ward et al.,, ), which improves the robustness of NPE by introducing a Spike and Slab error model on simulated data statistics.
We do not run NNPE on tasks C-F, as this would be unfair for NNPE whose noise model is designed for low-dimensional statistics, as opposed to the temporal and spatial data considered in tasks C-F. Extending NNPE to such settings is out of the scope of our experiments.
We also run the hybrid learning method HVAE (Takeishi and Kalousis,, ), which constitutes a strong baseline when the simulator can be made deterministic (tasks C and E) but is not directly applicable when the simulator is not differentiable. More details about each method and the experimental setup can be found in Appendix D.

### 4.1 Results

##### ROPE achieves robust posterior estimation for all tasks.

In Figure 2, we compare the performance of ROPE (with $\gamma=0.5$) against baselines for the $6$ tasks we consider.
For all tasks, even with minimal calibration budgets, ROPE is the only method that consistently returns well-calibrated, or sometimes slightly under-confident, posterior estimation while significantly reducing uncertainty compared to the prior distribution. As the size of the calibration set increases, we see the adaptability of J-NPE and MLP as their performance improves and aligns with or outperforms ROPE. This adaptability is an expected behavior in IID settings, where real-world data eventually allows finding the minimizer of empirical risk among a class of predictors. However, on task E, where posteriors are complex conditional distributions—whose entropy increases with darker images and contain non-trivial dependencies between parameters—ROPE remains the best approach, even with a calibration set containing more than $1000$ examples. As an outlier, we observe that NPE trained on simulated data achieves the best results for the SIR benchmark (Task B), indicating that the misspecification of this benchmark is not a challenging test case for existing SBI methods.
Finally, as interpreting a numerical gap in LPP metrics can be difficult, we complement these numerical results with visual corner plots for the two real-world experiments in Figure 3 and for all tasks in subsection E.1.

Figure: Figure 3: (left) Credible intervals of the posterior estimates at levels $65\%$ and $90\%$, for a single test sample from the light-tunnel task. The black stars denote the true value of the parameter. (center) Posterior estimates for a single test sample from the wind-tunnel task, where the true parameter is denoted by a vertical black line. (right) Effect of $\gamma$ on the LPP and ACAUC scores of ROPE on the light-tunnel task for different sizes of the calibration set. The value of $\gamma$ is shown by each curve. For reference, we plot the metrics achieved by the SBI posterior and prior distribution on simulated data.
Refer to caption: x3.png

##### Ablation study.

Our algorithm combines two steps with distinct roles: (1) a fine-tuning step, which improves the generalization domain of the NSE, and (2) an OT step, aiming to model the misspecification as a stochastic mapping between simulations and observations. To better understand their respective contribution to the performance of ROPE, we look at two ablated versions of our algorithm: tuning-only which appends the fine-tuned NSE to the NF trained on simulated data $p_{\phi^{\star}}$ and directly applies it to the real observations without an OT step; and OT-only, which directly performs OT with L2-norm in the original NSE space $c(\mathbf{x}_{o},\mathbf{x}_{s})=\lvert\mathbf{h}_{\omega^{\star}}(\mathbf{x}_
{o})-\mathbf{h}_{\omega^{\star}}(\mathbf{x}_{s})\rvert_{2}$.
In Figure 2, we observe that tuning-only’s result are poor except for Task B, for which misspecification is negligible. In contrast, for tasks A, D, and F, OT-only exhibits performance on par with ROPE. Nevertheless, ROPE can significantly outperform OT-only, such as in tasks C and E where the misspecification is significant. We conclude that the OT step is crucial and fine-tuning is sometimes necessary—we recommend that practitioners first evaluate OT-only’s performance and optimize the value of $\gamma$ before using a subset of the real-world data for fine-tuning.

##### Effect of entropic regularization.

In Figure 3, we study the effect of entropic regularization by varying the regularization parameter $\gamma$.
For all values of $\gamma$, excluding $\gamma\geq 5$, we observe that both LPP and calibration consistently improve with the calibration set size.
For large values of $\gamma$, the entropic regularization dominates and pushes toward a uniform mapping, resulting in posteriors that approximate the prior distribution, as highlighted in subsection 3.2. In these cases, the LPP barely changes with the calibration set size.
These empirical results are consistent with the theoretical discussion in subsection 3.2. Furthermore, it is notable that the value of $\gamma$ does not excessively influence ROPE’s performance. As a recommendation for practitioners, our empirical evaluation suggests that values between $0.1$ and $1$ provide well-calibrated and precise credible intervals.

## 5 Discussion

While our experiments have demonstrated that ROPE enables reliable uncertainty quantification from misspecified simulator, our method and setup has certain shortcomings that we now discuss.

##### Prior misspecification.

While the importance of correct prior specification vanishes as the number of observations or sharpness of the likelihood function increases for well-specified likelihood models, the implications of violating it may severely damage ROPE’s performance because of the OT matching between simulated and real-world observations.
As a possible solution, unbalanced OT (Séjourné et al.,, ; Fatras et al.,, ) relaxes the mass conservation constraint and could be an elegant solution to couple a subset of the simulations to the real-world data.

##### Curse of dimensionality.

The dimensionality of $\theta$ may impact two critical parts of ROPE. First, with each additional parameter $\theta[K+1]$, the NSE must encode up to $K$ dependence between $\theta[K+1]$ and other dimensions $\theta[1:K]$ given $\mathbf{x}_{o}$.
While generating more simulations can address the curse of dimensionality in the simulation space, fine-tuning on a small calibration may not suffice anymore to cope with misspecification. Second, the dimensionality of the manifold on which the simulated and real-world observations are projected with the NSE grows and finding a meaningful matching between the two populations may require larger sample sizes. A potential solution is to focus marginal or 2D posterior distributions and ignore higher-dimensional dependencies in $p(\theta\mid\mathbf{x}_{o})$. Nevertheless, extending ROPE to such settings certainly opens new questions, e.g., on the development of better fine-tuning strategies that can leverage partial calibration sets where labels can be incomplete.

##### Other extensions.

Similarly to incomplete labels, certain applications have only access to noisy labels measured with a well-modeled, but noisy, measurement process. Further developing the fine-tuning stage to exploit such noisy labels would be necessary to make an approach similar to ROPE applicable.
As another direction, leveraging inductive bias embedded into the neural network architecture of neural OT and the ability to better cope with large test set appears as a promising direction to amortize the mapping between simulation and real-world data. We believe following ROPE’s strategy of modeling misspecification in SBI as an OT coupling opens up several directions to cope with more specific problem setup.

##### Conclusion

In this paper, we have argued that model misspecification in simulation-based is a modeling issue that requires real-world labeled data to be solved. Under this premise, we have introduced ROPE, an algorithm that jointly exploits a small calibration set and optimal transport to extend neural posterior estimation for misspecified simulators. Our experiments on diverse benchmarks, demonstrate ROPE’s ability to estimate well-calibrated and informative posterior distributions for various simulators and real-world examples.
As a conclusion, ROPE is a simple method that practitioners can use to predict calibrated posterior over the parameter of misspecified simulator on real-world data.

##### Broader Impact

This paper presents a framework and an algorithm to address model misspecification in simulation-based inference (SBI). SBI is predominantly applied in scientific fields where complex simulators of physical phenomena are available, such as astronomy, particle physics, or climate modeling. A priori, this circumscribes the application of our algorithm to highly specialized scientific domains in the natural sciences, precluding issues such as fairness or privacy. However, its application to the scientific domain is not exempt from societal or ethical implications, particularly when computer simulations may inform research or policy decisions. In this regard, we find some properties of the algorithm particularly promising, such as uncertainty quantification and the limitation of not drawing conclusions beyond the given expert model. However, much more work is needed to deeply understand the reliability of these properties and how they are affected by violations of the core assumptions, such as a well-specified prior. Such work should precede any sort of over-selling to practitioners about the benefits of the algorithm. Rather, we see our work as an academic contribution towards a more broad and successful application of SBI techniques; success in this endeavor, as for the establishment of any scientific tool, will require an iterative dialogue between the scientists who develop the methodology and those who use it.

## Appendix A Robustness to distribution shifts

Figure: Figure 4: Out-of-distribution performance of ROPE and some baselines. We train ROPE and other baselines on the same light-tunnel data as in task E (training distribution), but apply it to test sets originating from a target distribution where the real-world images are flipped vertically. We compare the performance on test sets from both distributions, showing the LPP and ACAUC scores for each method. For comparison, in the right plot we show again the LPP curve (light gray, dotted) attained by ROPE under the training distribution. The performance of ROPE is barely affected as it cannot exploit any signal in the real images ($\mathbf{x}_{o}$) beyond what is encoded in the simulator, and the simulator output ($\mathbf{x}_{s}$) is invariant to the transformation we consider. Because NPE is not trained on real observations, its performance, although poor, also remains virtually unchanged. On the other hand, the performance of MLP and J-NPE drops in the target distribution, as these methods are not limited in what information they can exploit from the real observations on which they are trained, potentially learning shortcuts that are not present in the target distribution. This results demonstrate that if the simulator embeds the right invariances, our modeling assumption $\mathbf{x}_{o}\perp\theta\mid\mathbf{x}_{s}$ can be favorable to out-of-distribution generalization.
Refer to caption: x4.png

## Appendix B Self-calibration Property

We say ROPE is self-calibrating because, by design, the posterior distribution marginalized over observations tends to the prior as the number of simulation increases, that is,

$$ $\displaystyle\int\tilde{p}(\theta\mid\mathbf{x}_{o})p(\mathbf{x}_{o})=p(\theta).$ (10) $$

This property is also called marginal calibration as it is a necessary condition for a posterior estimation method to be calibrated. Considering NPE, $\tilde{p}(\theta\mid\mathbf{x}_{s})$, is marginally calibrated and observations $\mathbf{x}_{o}$ are generated from the assumed prior, that is sampled from an unknown distribution $p(\mathbf{x}_{o})=\int p(\mathbf{x}_{o}\mid\theta)p(\theta)$, we can show ROPE is marginally calibrated.
Indeed, considering the Monte-Carlo approximation of the marginalized posterior distribution over the test set $\mathcal{D}_{o}:=\{\mathbf{x}_{o}\}_{i=1}^{N_{o}}$, we have,

$$ $\displaystyle\int\tilde{p}(\theta\mid\mathbf{x}_{o})p(\mathbf{x}_{o})$ $\displaystyle\approx\frac{1}{N_{o}}\sum_{i=1}^{N_{o}}\tilde{p}(\theta\mid \mathbf{x}_{o})$ (11) $\displaystyle=\frac{1}{N_{o}}\sum_{i=1}^{N_{o}}\sum_{j=1}^{N_{s}}N_{o}P^{\star }_{ij}\tilde{p}(\theta\mid\mathbf{x}_{s}^{j})$ (12) $\displaystyle=\sum_{j=1}^{N_{s}}\left[\sum_{i=1}^{N_{o}}P^{\star}_{ij}\right] \tilde{p}(\theta\mid\mathbf{x}_{s}^{j})$ (13) $\displaystyle=\frac{1}{N_{s}}\sum_{j=1}^{N_{s}}\tilde{p}(\theta\mid\mathbf{x}_ {s}^{j})$ (14) $\displaystyle\approx p(\theta),$ (15) $$

where we use the definition of the transport matrix to get $\sum_{i=1}^{N_{o}}P^{\star}_{ij}=\frac{1}{N_{s}}$. The last approximation tends to be exact as the number of simulations increases, if the NPE is marginally calibrated.

## Appendix C ROPE Algorithm

Figure: Algorithm 1 Posterior Inference using Robust Neural Posterior Estimation (ROPE)

## Appendix D Experimental Setup

In this section, we provide more details on our experiments. For completeness, we provide details on the neural architectures and training hyperparameters. However, we encourage the reader interested in reproducing our experiments to examine our code directly (a link to the code will be made available in the public version of the paper).

For all methods training on calibration set we keep always keep $20\%$ of the calibration to monitor validation performance and we select the best model based on this metric.

For the MLP we use the same architecture as the NSE for all our experiments and optimize its parameters on the calibration set with Adam and a learning rate equal to $0.0003$, we select the best model based on the LPP attributed to the validation subset of the calibration set.

### D.1 Task A: CS & Task B: SIR

We refer the reader to Ward et al., () for more details about the simulator and prior distribution. We use the exact same setting as theirs.

#### Neural architecture & Training Hyperparameters

For all methods we use the same backbone MLP as the NSE with ReLU activations and layers composed of $[4K,16K,16K,12K,3K]$ neurons, where $K$ is the dimensionality of $\theta$. The NF is a 1-step UMNN-MAF (Wehenkel and Louppe,, ) with $[100,100,100]$ neurons for both the autoregressive conditioner and normalizer. For NNPE, we train the UMNN-MAF on simulations poluted by Spike and Slab errors. We train models with Adam and a learning rate equal to $0.0005$ and all other parameters set to default. We optimize the SBI model for $10^{6}$ gradient steps and select the best model on random validation sets containing $10^{5}$ simulations.

### D.2 Task C: Pendulum

#### Description

The first task is inspired from the damped pendulum benchmark commonly used to assess hybrid learning algorithms. Given a $2$D physical parameter $\theta:=[\omega_{0},A]$, where $\omega_{0}\in\mathbb{R}^{+}$ denotes the fundamental frequency and $A\in\mathbb{R}^{+}$ the amplitude of a friction-less pendulum, the simulator generates the horizontal position of the pendulum at $200$ discrete times during uniformly sampled in a $10$ seconds interval as

$$ $\displaystyle\mathbf{x}_{s}:=[\theta(t=0),\dots,\theta(t=10s)]\in\mathbb{R}^{2 00}$ $\displaystyle\text{ where }\theta(t)=A\cos(\omega_{0}t+\varphi)\quad\varphi \sim\mathbb{U}(-\pi,\pi).$ (16) $$

The relationship between the parameters and the simulation is thus stochastic as $\varphi$ accounts for an unknown phase shift when the measurements start.
We generate real-world observations synthetically by replacing $\theta(t)$ from (16) by

$$ $\displaystyle\tilde{\theta}(t)=e^{\alpha t}A\cos(\omega_{0}t+\varphi)\quad \varphi\sim\mathbb{U}[-\pi,\pi]\quad\alpha\sim\mathbb{U}[0,1],$ $$

where $\alpha$ represents the effect of friction. We also add Gaussian noise on both simulated and real-world data to represent the inaccuracy of a sensor measuring the pendulum’s position. The prior distribution is a product of uniform distribution, $p(\theta:=[\omega_{0},A])=\mathcal{U}[0,3]\times\mathcal{U}[0.5,10]$.

#### Neural architecture & Training Hyperparameters

##### Neural Posterior Estimator.

The NSE is a 1D convolutional neural network, with the architecture described in Algorithm 2.

Figure: Algorithm 2 Convolutional Neural Network for Tasks A and D.

The NCDE is a one-step discrete normalizing flow with an autoregressive conditioner and a UMNN (Wehenkel and Louppe,, ) as the normalizer.
The autoregressive conditioner is a MADE with ReLU activation and $3$ layers of $100$ neurons that output a $10$ dimensional vector to the UMNN. The UMNN has an integrand net with $3$ layers of $100$ neurons with ReLU activations. For training the NPE, we use a batch size of $100$ and a learning factor equal to 1e-4. NPE is trained until convergence. Other parameters are set to default values and should marginally impact the NPE obtained.

##### ROPE NSE.

We have selected the best NPE based on the validation set with $10000$ examples generated with the simulator. The NPE is fixed to one best-of-all model.
We fine-tune the NCDE with a learning rate equal to 1e-5 for $5000$ gradient steps on $80\%$ the full calibration set. We use a 1-sample Monte Carlo estimate of the expectation in (7).

##### J-NPE

To train J-NPE, we simply randomly use a batch composed of $50\%$ of simulated pairs $(\theta,\mathbf{x}_{s})$ and of $50\%$ $(\theta,\mathbf{x}_{o})$ from the calibration set. We use the same architecture and hyper-parameters as the SBI NPE. The best model is selected based on the best training set performance. We do $50$ epochs with $50000$ simulated examples for each epoch. The batch size is $100$.

##### HVAE

For the HVAE, we re-use the NPE model as the physics encoder and replace the decoder with a deterministic version of the simulator, thus removing the Gaussian noise on a random phase shift. In addition, we follow the approach of Takeishi and Kalousis, () and have 1) a real-world encoder that maps $\mathbf{x}_{o}$ to $\mathbf{z}_{a}$, 2) a reality-to-physics encoder, and 3) a physics-to-reality decoder. The real-world encoder has the same architecture as the NSE of the NPE and outputs the mean and log-variance of a $5$D latent vector $\mathbf{z}_{a}$. The reality-to-physics and physics-to-reality also have the same architectures and are two conditional 1D U-Net with neural network architecture described in Algorithm 3.

Figure: Algorithm 3 UNet1D Architecture

To train the HVAE, we freeze the parameters of the NPE and optimizes the ELBO as well as a calibration loss that evaluates the likelihood assigned to the true physical parameters. All distributions are parameterized by Gaussian with mean and log-variance predicted by the neural networks. We do not use any additional losses as we expect constraining NPE and using the calibration set should already provide the necessary support to use the physics in a meaningful way. The HVAE is trained on the $2000$ test examples as it is the only real-world data, calibration set aside, that we have access to. We use a batch size equal to $100$ and a learning rate equal to 1e-3. We believe obtaining a better HVAE is possible. However, we emphasize the complexity of setting up a good HVAE for the only purpose of statistical inference over parameters.

#### Datasets

For this task, we can generate samples $(\theta,\mathbf{x}_{s})$ on the fly to train the NPE. The calibration and test sets are also generated randomly by sampling from the prior distribution and using the damped pendulum simulator.

### D.3 Task D: Hemodynamics

#### Description

Inspired by Wehenkel et al., (), we define the task of inferring important cardiovascular parameters from normalized arterial pressure waveforms measured at the radial artery. The simulator
uses many physiological parameters that modulates the heart function, physical properties of the $116$ main arterial segments, and behavior of the vascular beds. Our inference concerns two parameters of the heart function, $\theta:=[\text{SV},\text{LVET}]$, the stroke volume (SV) is the amount pumped out from the left ventricle over the heart beat modeled, and the left ventricular ejection time (LVET) is the time interval between opening and closure of the aortic valve. Other parameters, such as the heart rate or arteries’ stiffness, are considered as nuisance effects and are randomly sampled from a realistic population distribution. An additional source of randomness is added by modeling measurement errors with a white Gaussian noise and randomizing the starting recording time with respect to the cardiac cycle. The simulator produces $8$-second timeseries $\mathbf{x}_{t}\in\mathbb{R}^{1000}$ sampled at $125Hz$. As synthetic misspecification, the simulator assumes all arteries have the same length over the population considered, whereas "real-world" data are artificially generated by also varying the length of arteries and account for the effect of human’s height. The simulator is based on the openBF PDE solver (Melis,, ) specialized for hemodynamics, which is not differentiable and takes approximately one minute to simulate one sample on a standard CPU. This synthetic tasks represent a common scenario in which a simulator, although faithful to the effect of certain parameters, misses additional degrees of freedom that exists for the real-world data.

#### Neural architecture & Training Hyperparameters

Figure: Algorithm 6 CNN Architecture for Task C.

##### Neural Posterior Estimator.

The NSE is the 1D convolutional neural network described in Algorithm 6.
The NCDE is a $5$-step discrete normalizing flow with an autoregressive conditioner and affine normalizers.
Each of the $5$ autoregressive conditioners is a MADE with ReLU activations and $4$ layers of $300$ neurons that output $4$ dimensional vectors used to parameterize the affine transformations. For training the NPE, we use a batch size of $100$ and a learning factor equal to 5e-4. NPE is trained until convergence. Other parameters are set to default values and should marginally impact the NPE obtained.

##### ROPE NSE.

We have selected the best NPE based on the validation set with $2000$ examples generated with the simulator. The NPE is fixed to one best-of-all model.
We fine-tune the NCDE with a learning rate equal to 1e-5 for $2000$ gradient steps on $80\%$ of calibration set. We use a 1-sample Monte Carlo estimate of the expectation in (7).

##### J-NPE

To train J-NPE, we simply randomly use a batch composed of $50\%$ of simulated pairs $(\theta,\mathbf{x}_{s})$ and of $50\%$ $(\theta,\mathbf{x}_{o})$ from the calibration set. We use the same architecture and hyper-parameters as the SBI NPE. The best model is selected based on the best training set performance. We do $50$ epochs with $6000$ simulated examples for each epoch. The batch size is $100$.

##### HVAE

There is no HVAE for this experiment as the simulator is non-differentiable.

#### Datasets

For this task, we cannot generate samples $(\theta,\mathbf{x}_{s})$ on the fly to train the NPE. For the purpose of this experiment, we have generated $10000$ simulations and real-world observations. Our fine-tuning strategy approximates (7) by finding the simulations with the closest parameter value.

### D.4 Task E: Light Tunnel

#### Description

We use one of the light-tunnel datasets from the causal chamber project (Gamella et al.,, , [causalchamber.org](https://causalchamber.org)). In particular, we use the data from the [ap_1.8_iso_500.0_ss_0.005](ap_1.8_iso_500.0_ss_0.005) experiment in the [lt_camera_v0](https://github.com/juangamella/causal-chamber/tree/main/datasets/lt_camera_walks_v0) dataset. The light tunnel is an elongated chamber with a controllable light source at one end, two linear polarizers mounted on rotating frames, and a camera that takes images of the light source through the polarizers. We refer the reader to Gamella et al., (, Figure 2) for a complete schematic. Our task consists of predicting the color setting of the light source ($(R,G,B)\in[0,255]^{3}$) and the dimming effect of the linear polarizers $\alpha\in[0,1]$ from the captured images. As a misspecified simulator of the image-generating process, we adopt the simple model described in Gamella et al., (, Model F1, Appendix D). A Python implementation is available through the [causalchamber](https://pypi.org/project/causalchamber/) package ([models.model_f1](models.model_f1)); visit [causalchamber.org](https://causalchamber.org) for more details. As input, the simulator takes the parameters $\theta:=[R,G,B,\alpha]$ and produces an image consisting of a hexagon roughly the size of the light source, with an RGB color vector equal to
$[\alpha R,\alpha G,\alpha B]$. The factor $\alpha:=\cos^{2}(\theta_{1}-\theta_{2})$, where $\theta_{1},\theta_{2}$ denote the angles of the two polarizers, corresponds to Malus’ law (e.g. , Collett,, ), which models the dimming effect of the polarizers as a function of their relative angle. Besides the obvious misspecification with respect to image realism (see Figure 2), the model ignores other important physical aspects, such as the spectral response of the camera sensor or the non-uniform effect of the polarizers on the different colors—more details can be found in Gamella et al., (, Appendix D.IV.2.2). The prior is uniform over colors and polarizer angles, which leads to a non-uniform prior over the dimming effect $\alpha$.

#### Neural architecture & Training Hyperparameters

##### Neural Posterior Estimator.

The NSE is the 2D convolutional neural network described by Algorithm 5.

The NCDE is also a one-step discrete normalizing flow with an autoregressive conditioner and a UMNN (Wehenkel and Louppe,, ) as the normalizer. The autoregressive conditioner is a MADE with ReLU activation and $3$ layers of $500$ neurons that outputs a $10$ dimensional vector to the UMNN. The UMNN has an integrand net with $4$ layers of $150$ neurons with ReLU activations.
For training the NPE, we use a batch size of $100$ and a learning factor equal to 5e-4. NPE is trained until convergence. Other parameters are set to default values and should marginally impact the NPE obtained.

##### ROPE NSE.

We have selected the best NPE based on the validation set with $10000$ examples generated with the simulator. The NPE is fixed to one best-of-all model.
We fine-tune the NCDE with a learning rate equal to 1e-4 for $2000$ gradient steps on on $80\%$ of the calibration set. We use a 1-sample Monte Carlo estimate of the expectation in (7).

##### J-NPE

To train J-NPE, we simply randomly use a batch composed of $50\%$ of simulated pairs $(\theta,\mathbf{x}_{s})$ and of $50\%$ $(\theta,\mathbf{x}_{o})$ from the calibration set. We use the same architecture and hyper-parameters as the SBI NPE. The best model is selected based on the best training set performance. We do $50$ epochs with $1000$ simulated examples for each epoch. Simulations are generated randomly for each batch by sampling the prior and simulating for the corresponding parameters. The batch size is $100$.

##### HVAE

For the HVAE, we re-use the NPE model as the physics encoder and use the simulator as is as it is differentiable without additional effort. In addition, we follow the approach of Takeishi and Kalousis, () and have 1) a real-world encoder that maps $\mathbf{x}_{o}$ to $\mathbf{z}_{a}$, 2) a reality-to-physics encoder, and 3) a physics-to-reality decoder. The real-world encoder has the same architecture as the NSE of the NPE and outputs the mean and log-variance of a $5$D latent vector $\mathbf{z}_{a}$. The reality-to-physics and physics-to-reality also have the same architectures and are two conditional 2D U-Net with the architecture described by Algorithm 7.

Figure: Algorithm 7 2D UNet

To train the HVAE, we freeze the parameters of the NPE and optimizes the ELBO as well as a calibration loss that evaluates the likelihood assigned to the true physical parameters. All distributions are parameterized by Gaussian with mean and log-variance predicted by the neural networks. We do not use any additional losses as we expect constraining NPE and using the calibration set should already provide the necessary support to use the physics in a meaningful way. The HVAE is trained on the $2000$ test examples as it is the only real-world data, calibration set aside, that we have access to. We use a batch size equal to $100$ and a learning rate equal to 1e-3. We believe obtaining a better HVAE is possible. However, we emphasize the complexity of setting up a good HVAE for the only purpose of statistical inference over parameters.

#### Datasets

For this task, we can generate samples $(\theta,\mathbf{x}_{s})$ on the fly to train the NPE. However, the calibration and test sets are real-world data. We ensure there is not overlap between calibration and test set. The is no randomization and the test set is constant for all experiments, the calibration set are also fixed for a given calibration set size.

### D.5 Task F: Wind Tunnel

#### Description

We use one of the wind-tunnel datasets from the causal chamber project (Gamella et al.,, , [causalchamber.org](https://causalchamber.org)). In particular, we use the data from the [load_out_0.5_osr_downwind_4](load_out_0.5_osr_downwind_4) experiment in the [wt_intake_impulse_v1](https://github.com/juangamella/causal-chamber/tree/main/datasets/wt_intake_impulse_v1) dataset.
The tunnel is a chamber with two controllable fans that push air through it and barometers that measure air pressure at different locations. A hatch precisely controls the area of an additional opening to the outside (see Gamella et al.,, , Figure 2). The data is a collection of pressure curves that result from applying a short impulse to the intake fan load and measuring the change in air pressure using one of the barometers inside the tunnel. Our inference task consists of predicting the hatch position, $\theta:=[H]\in[0,45]$ given a pressure curve (see Figure 2). As a simulator model, we combine the models A2 and C3 described in Gamella et al., (, Appendix D); we numerically solve the ODE in model A2, and add stochastic components to simulate the sensor noise and the unknown time point at which the impulse is applied. This results in the simulator being neither differentiable nor deterministic. A Python implementation of the complete simulator is available in the [causalchamber](https://pypi.org/project/causalchamber/) package ([models.simulator_a2_c3](models.simulator_a2_c3)); visit [causalchamber.org](https://causalchamber.org) for more details. Misspecification arises from the many simplifying assumptions needed to model the complex dynamics of the airflow inside the tunnel—more details can be found in Gamella et al., (, Appendix D.IV.1.2).

##### Neural Posterior Estimator.

The NSE and NCDE have the same 1D convolutional neural network as for Task A. For training the NPE, we use a batch size of $100$ and a learning factor equal to 5e-4. NPE is trained until convergence. Other parameters are set to default values and should marginally impact the NPE obtained.

##### ROPE NSE.

We have selected the best NPE based on the validation set with $10000$ examples generated with the simulator. The NPE is fixed to one best-of-all model.
We fine-tune the NCDE with a learning rate equal to 1e-4 for $20000$ gradient steps on on $80\%$ of the calibration set. We use a 1-sample Monte Carlo estimate of the expectation in (7).

##### J-NPE.

To train J-NPE, we simply randomly use a batch composed of $50\%$ of simulated pairs $(\theta,\mathbf{x}_{s})$ and of $50\%$ $(\theta,\mathbf{x}_{o})$ from the calibration set. We use the same architecture and hyper-parameters as the SBI NPE. The best model is selected based on the best training set performance. We do $50$ epochs with $10000$ simulated examples for each epoch. The batch size is $100$.

##### HVAE.

There is no HVAE for this experiment as the simulator is non-differentiable.

#### Datasets

For this task, although slightly slower than Task A and B, we can generate samples $(\theta,\mathbf{x}_{s})$ on the fly to train the NPE. However, the calibration and test sets are real-world data. We ensure no overlap between the two sets for all calibration set sizes. All sets are fixed for all experiments.

## Appendix E Additional results

### E.1 Corner plots

Figure: Figure 5: Three corner plots for task A with a calibration set with $50$ samples.
Refer to caption: x5.png

Figure: Figure 6: Three corner plots for task B with a calibration set with $50$ samples.
Refer to caption: x8.png

Figure: Figure 7: Three corner plots for task C with a calibration set with $50$ samples.
Refer to caption: x11.png

Figure: Figure 8: Three corner plots for task D with a calibration set with $50$ samples.
Refer to caption: x14.png

Figure: Figure 9: Three corner plots for task E with a calibration set with $50$ samples.
Refer to caption: x17.png

Figure: Figure 10: Three corner plots for task E with distribution shift with a calibration set with $50$ samples.
Refer to caption: x20.png

### E.2 Calibration plots

Figure: (a)
Refer to caption: extracted/2405.08719v1/neurips2023_otml/figures/calibration/calibration_Light_Tunnel.png

## Appendix F Computing ACAUC

Figure: Algorithm 9 Statistical Calibration of Posterior Distribution

## Appendix G Learning Minimal Sufficient Statistics with Neural Posterior Estimation

We now discuss why NPE may learn a minimal sufficient statistic. First, under a sufficiently large validation set, NPE’s objective function is only optimal on the validation set if NPE models the true posterior as defined implicitly by the prior $p(\theta)$ and the likelihood corresponding to the simulator $S$. This consistency has been proven in (Papamakarios and Murray,, ) and is the motivation to use such an objective when estimating density. Second, some normalizing flows, such as autoregressive UMNN flows (Wehenkel and Louppe,, ), are universal approximators of continuous densities. In addition, neural networks are also universal function approximators. As such, we can claim that it is always possible to parameterize the NCDE $p_{\theta}(\theta\mid\mathbf{h}_{\omega}(\mathbf{x}))$ such that the class of functions its parameters represent contains the true posterior.
We directly observe that $\mathbf{x}$ is only used by the NCDE through $\mathbf{h}_{\omega}(\mathbf{x})$. Thus, under perfect training $p_{\theta^{\star}}(\theta\mid\mathbf{h}_{\omega^{\star}}(\mathbf{x}))=p(\theta
\mid\mathbf{x})$ and $\mathbf{h}_{\omega^{\star}}(\mathbf{x})$ is a sufficient statistic for $\theta$ given $\mathbf{x}$ under the simulator’s model.

Without additional constraint, we cannot claim anything about the minimality of $\mathbf{h}_{\omega^{\star}}(\mathbf{x})$. Nevertheless, we can enforce the neural network $\mathbf{h}_{\omega^{\star}}(\mathbf{x})$ has an information bottleneck and thus reduces the information carried. In practice, we ensure the output dimension of $\mathbf{h}_{\omega^{\star}}(\mathbf{x})$ and the NCDE achieves optimal performance on the test set. Recalling that, in the context of SBI, we can generate as many samples as needed and obtain estimators that closely approach the simulation’s posterior and a minimal sufficient statistic.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: In section 3.1, the authors say that "it is reasonable to learn a sufficient statistic $\mathbf h_o$ by fine-tuning $\mathbf h_{\omega^\star}$". How do we know that such a fine-tuning will produce sufficient statistics? Are there any guarantees?

In most of the experiments, the performance of RoPE does not var as the calibration set size increases. Is there a reason why the performance of RoPE saturates in those cases even with $N_c = 10$? Are these very simple problems?
soundness: 2
strengths: The paper addresses a relevant practical problem in the field of SBI. The paper is largely well-written, and the experimental evaluation seems thorough. The idea of using the optimal transport coupling for learning statistics is neat (approach is similar to the work of Huang et al. (2023) who also learn statistics to be robust to model misspecification, without the calibration set of course).
summary: The paper presents a robust posterior estimation method (RoPE) for simulation-based inference (SBI) when the model is misspecified but the the underlying parameters of interest can be physically measured via (possibly expensive) experiments. RoPE uses such a calibration dataset to better quantify the posterior uncertainty under model misspecification scenarios.
weaknesses: I have concerns regarding the soundness of the proposed method. The main claim is that RoPE provides well-calibrated uncertainty quantification, which I am not convinced of. The self-calibration property discussed in Appendix F is under the assumption that the NPE posterio is well-calibrated, which itself is not guaranteed and is often overconfident as pointed out by [1] (and therefore making NPE well-calibrated is an active area of research in SBI [2]). Moreover, the expectation in equation 10 in the proof is estimated using $N_o$ Monte Carlo samples from the test set, which I am assuming, would be limited (unless having a large test set is one of the assumptions/requirements of RoPE, in which case it should be mentioned clearly since the other methods do not have such a requirement). Basically I do not see why RoPE solves the problem being addressed. Please let me know if I am missing something.

[1] Hermans et al. (2022). A Crisis In Simulation-Based Inference? Beware, Your Posterior Approximations Can Be Unfaithful. TMLR.
[2] Falkiewicz et al. (2023). Calibrating Neural Simulation-Based Inference with Differentiable Coverage Probability. NeurIPS.

Looking at how the authors define the model misspecification problem, it seems like the problem is of posterior calibration. Calling it "model misspecification" only makes it confusing since the definition is in terms of posteriors, and is not necessarily a property of the model (the example in Appendix A focuses on posterior bias rather than calibration). I do not see what is gained by calling this a "model misspecification problem", as a well-specified model may still produce mis-calibrated posteriors (mentioned in lines 92-93). I really think that the story becomes straightforward if the problem being addressed is of obtaining calibrated SBI posteriors (i.e. a "calibration problem"), both in cases when the model is well- or mis-specified (right?), which can be done using some additional information (in the form of calibration dataset), which is available in scenarios such as the cardiovascular system and industrial process. If the authors agree with this statement, then the more appropriate baseline for RoPE is the method of [2], instead of NPE-RS and NNPE.

The equation in line 190 is the same as equation 3 from the Ward et al. (2023) paper (NNPE method). Seems like the primary difference between the two is that RoPE uses OT to learn $p(\mathbf x_s | \mathbf x_o)$, while NNPE uses a spike and slab model. Please elaborate on the key differences between the two approaches, as it would help gauge the originality of this work. Also, it is not clear to me what is the role of the assumption in equation 3 then, as NNPE does not seem to have that assumption but still arrives at the same equation (correct me if I am wrong). I suppose lines 193-194 states why this assumption is needed, but it is not clear to me. Plus, there needs to be more discussion on how limiting is this assumption, in which cases is it satisfied and not satisfied, and the implications for when it does not hold.

The notation is unclear at times. For instance, what is the difference between $\mathbb{P}^\star$ and $p^\star(\mathbf x_o)$? Sometimes the author say $p^\star(\theta)$ is the true distribution of parameters in the real-world, but then they use $p^\star(\theta | \mathbf x_o)$ which is not introduced.

In section 2.2, it would be great if you could explain what an OT coupling is and what is the general idea (in simple words) before going into the OT details. This would help the readers not familiar with OT to get a sense of what this does without needing to understand every detail.

Minor comments:
* In equation 6, I think there should be $\varphi$ in the left-hand side instead of $\phi$, right?
* typo in line 817: "i.d.d." instead of "i.i.d".

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: Can the authors explain the insensitivity of performance to the size of the calibration set?
Can the method be extended to incorporate a less-than-perfect calibration set (e.g., as coming from a well-calibrated posterior estimator from another data set)?
soundness: 3
strengths: - The paper contributes an interesting approach to a hard problem in SBI and scientific modeling in general.
- The OT formulation and ensemble-based correction are well motivated and the assumptions are clearly stated.
- The evaluation in the setting of low-dimensional models is comprehensive and features good comparisons. The results are encouraging even with a very small calibration set.
summary: This work contributes another contender to the family of “robust” methods for neural SBI, which is known to be unreliable for out-of-distribution (OOD) data. The basic idea is to use a “calibration set” of observables and ground-truth measurements as an additional training signal to inform the neural network responsible for embedding the observables into a lower-dimensional representation about distribution shifts and reduce its susceptibility to OOD observables. During inference, an optimal transport (OT) cost matrix is used to derive weights for approximating the posterior over observables as a weighted average over an ensemble of posteriors obtained from the training simulations.
weaknesses: - The basic premise of the work is that “ground truth” parameters or latent quantities are available during training. In other words, it assumes that the simulator is parameterized in the same way as reality. This severely limits the applicability of the proposed method beyond certain domains of physics physics, as “true” latent quantities are not available or not even sought-for in a majority of applications of Bayesian inference outside of physics. Nevertheless, the authors openly acknowledge this limitation and I don’t think that it is an unsurmountable dealbreaker. On a related note, this limitation extends to the formal definition of misspecification assuming there is a “true prior”, and correspondingly, a “true posterior” that represents the updated prior which is hardly actionable - the marginal likelihood $p(x_0)$ is a model-implied quantity and de Finetti’s representation theorem about the existence of a factorization apples only to infinitely exchangeable sequences without any guarantees that $dim(\theta)$ is finite. 
- The authors are quite selective in their choice of cited related work. For instance, the authors can find extensive discussions on model misspecification in likelihood-free inference in [1, 2, 3, 4, see also references therein]; [5] also discuss an alternative definition of model misspecification through the scope of the marginal likelihood, which is also explicitly discussed in [6]. [7] propose to view misspecification as latent-space distortion and correct it during inference using the reverse KL. This goes on to show that the proposed paper can profit from a slightly more in-depth discussion on conceptual and empirical approaches to dealing with model misspecification and perhaps offer a categorization of methods designed to *detect* vs. methods designed to *fix* potential errors. Additionally, reference [8] regarding universal density approximation is misleading; see [9 and references therein] for a comprehensive treatment of the university of coupling-based normalizing flows.
- I have two concerns regarding the scalability of the method. First, it seems that the transport matrix can grow very large for settings with millions of observations and millions of simulations. This, coupled with the need to optimize the optimize the $\gamma$ trade-off parameter makes the method attractive only for small-data, small-simulation-budget settings. Second, the empirical evaluation only considers models with extremely low parameter dimensionality $\theta$ (i.e., the largest parameter space has 4 parameters), which is a serious limitation and leaves the question of generalizability to more complex models completely open. As such, I believe the claims of the paper to be a bit over the top. I do appreciate that the evaluations were performed on fairly large tests sets of 2000 simulations and that recent contenders and good ablations were considered.
- The presentation of the results is rather hard to follow. I believe some of the model details (e.g., example data sets) could have been delegated to the appendix and sushi plots summarized more succinctly. I am a bit concerned by the lack of sensitivity regarding the size of the calibration set (apart from task E).

[1] Nott, D. J., Drovandi, C., & Frazier, D. T. (2023). Bayesian inference for misspecified generative models. Annual Review of Statistics and Its Application, 11.

[2] Frazier, D. T., Robert, C. P., & Rousseau, J. (2020). Model misspecification in approximate Bayesian computation: consequences and diagnostics. Journal of the Royal Statistical Society Series B: Statistical Methodology, 82(2), 421-444.

[3] Frazier, D. T., Kohn, R., Drovandi, C., & Gunawan, D. (2023). Reliable Bayesian inference in misspecified models. arXiv preprint arXiv:2302.06031.

[4] Kelly, R., Nott, D. J., Frazier, D. T., Warne, D., & Drovandi, C. (2024). Misspecification-robust sequential neural likelihood for simulation-based inference. Transactions on Machine Learning Research.

[5] Schmitt, M., Bürkner, P. C., Köthe, U., & Radev, S. T. (2023, September). Detecting model misspecification in amortized Bayesian inference with neural networks. In DAGM German Conference on Pattern Recognition.

[6] Masegosa, A. (2020). Learning under model misspecification: Applications to variational and ensemble methods. Advances in Neural Information Processing Systems, 33, 5479-5491.

[7] Siahkoohi, A., Rizzuti, G., Orozco, R., & Herrmann, F. J. (2023). Reliable amortized variational inference with physics-based latent distribution correction. Geophysics, 88(3), R297-R322.

[8] Wehenkel, A., & Louppe, G. (2019). Unconstrained monotonic neural networks. Advances in neural information processing systems, 32.

[9] Draxler, F., Wahl, S., Schnörr, C., & Köthe, U. (2024). On the universality of coupling-based normalizing flows. arXiv preprint arXiv:2402.06578.

---

confidence: 4
contribution: 4
initial_rating: 6
presentation: 2
questions: 1) Definition of misspecification in SBI: There are several parts in the definition that
I find challenging to follow. In line 090 it says

> our simulator models the relationship between the real observations xo and the
> parameters of interest as they appear in the calibration set

and for that sentence it is followed, that the current definition is insufficient
because a model may be well-specified by not well-calibrated. I cannot follow this line
of argument. Are you defining mispecification based on posterior calibration here? Can
you please clarify? (I have seen the example in Appendix B, but it shown only that the
posterior is biased, not badly calibrated).

2) Related to my first question, I am wondering how your definition differs from the one
   given in https://arxiv.org/abs/1904.04551, who introduced a synthetic likelihood
   approach for misspecified simulators. 

3) In line 201, $\pi$ is introduced as a joint distribution. However, below in equation
   (4) it is then used as conditional distribution. Can you clarify this change please? 

4) In line 257, RoPE$^*$ is defined with $\tau < 1$, but how much smaller than 1 does it
   have to be? Further below, it is set to $\tau=0.9$. How is chosen in practice? 

5) In line 409, it is claimed that RoPE significantly reduces the uncertainty compared
   to the prior. How is this measured? Is LPP your proxy for this? If so, please clarify
   in the text. 

6) How is ACAUC related to the more common posterior calibration metrics like SBC or
   expected coverage? In line 381, it says that it is zero for a perfectly specified
   prior. Shouldn't it refer to the calibration performance of the posterior instead? 

7) The **SBI** defined in line 387 is unclear to me. How is it obtained precisely and
   how can it be the upper bound on the performance in misspecified simulators? How can
   it below the prior baseline in the ACAUC plots for some tasks? 

8) How can the MLP baseline work so well when it assumes a Gaussian posterior? Are the
   benchmarking task posteriors all close to Gaussian? 

9) The choice of $\gamma$ can have a substantial effect on the accuracy of the
   posterior, tending towards the prior for too large values of $\gamma$. This would be
   highly undesirable in practice. How can the practitioner tune $\gamma$ for real-world
   applications, i.e., outside the benchmarking tasks where it seems clear what the upper
   and lower bound of the metric are?

10) General remark on the figures: Figure 1 is great, but figure 2-4 are packed quite
    challenging to read. The dotted lines are barely visible and quite different than
    shown in the legend. Some lines have markers, others do not. Especially in the ACAUC
    plots it is almost impossible to tell apart the different lines. Smaller line width
    would be a start, but a general improvement in the choice of color code and markers
    would be even better. 

11) Some questions on the RoPE performance: Why does it perform so badly compared to NPE
    on the SIR example? You mentioned that the misspecification in SIR is weaker. How
    does RoPE perform on well-specified simulators in general then? Could this be a
    general weakness of the method?
    Also, RoPE seems to be generally underconfident in terms of ACAUC and this does not
    improve with calibration set size. This is concerning as well. I think it would be
    important to check and compare the posterior predictive distributions as well.


**Additional feedback**

- The reference to figures from the main text is a bit off. For example, figure 1
  appears on page 2, but is referenced only on page 5. Figures 3 and 4 are referenced
  before figure 2. 
- A lot of content is moved to the appendix. It would be better to move parts of it back
  to main text, e.g., algorithm 1. Or give more context in the main text. For example,
  to me it became clear only by looking at the algoritm, that at inference time, one
  needs to run additional simulations to match the test set of observed data for
  computing the weighted posterior.
soundness: 3
strengths: The paper is well-written and well-structured. The first two sections
give a good introduction to the topic, accurately cite previous work and clearly state
the contributions and the specific SBI setting the method is designed to solve. The
paper uses two existing benchmarks and introduces four new model misspecification
benchmarks, which is a valuable contribution for the community. The performance of the
proposed method is promising and the ablation and hyperparameter sensitivity studies
help understanding the method and tend to make it easier for practitioners to choose
parameter in practice.
summary: The paper introduces a method for addressing model misspecification in simulation-based
inference (SBI)  settings. Common SBI algorithms usually assume that the simulator (the
model) is well-specified, i.e., that it can simulate data close to the observed data.
However, if this is not the case, the resulting approximated posterior distributions can
be substantially biased. This paper proposes a procedure to mitigate this problem for
scenarios where multiple observed data points are available.

The authors propose to first follow a standard SBI approach of learning a neural-network
based embedding on the simulated data, followed by neural posterior estimation (NPE) on
the embedded data. Their key contributions consists of three additional steps. First,
they fine-tune the neural embedding using the set of observed data points. Second, they
then solve the optimal transport problem between the embedding simulated data and the
embedded observed data. Finally, at inference time, they use the optimal transport
solution between the embeddings to obtain the posterior given the observed data as a
weighted sum of the posteriors given simulated data.

The proposed method is evaluated on six benchmarking tasks for misspecification, two
existing tasks and four new tasks. It is compared to baselines and competing algorithms
and performs well in most settings. Additionally, the authors show results for ablations
and the method's behaviour for different hyperparameter settings.
weaknesses: The new definition of model misspecification in SBI is a bit unclear
(see questions below) and should be clarified. Some of the benchmarking results of the
proposed method are unclear and concerning and need clarification. In general, the
figures in the experiments sections are difficult to read and make the overall results
unclear (see suggestions below).

**Initial recommendation:**

Overall, I believe this could be a valuable contribution for the field and I am
tentatively recommending to accept this paper. Model misspecification is an urgent
problem in SBI and methods for addressing it have been scarce so far. The proposed
method and additionally the proposed benchmarking tasks are a promising next step.
However, there are several questions that should be clarified and several improvements
to the visualizations of the results that should be made. The performace of RoPE on
several benchmarking tasks seems problematic to me, but this might be due to a
misunderstanding of the ACAUC metric and needs clarification (see question 11 below).

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: - I also don’t really see in Section 3, where exactly the real-world labels (the parameters underlying the real data) are actually used. Can you explain this?
- In EQ 3, you make a central assumption. Can you explain why you think this assumption is usually justified? Asked differently, when would it not be justified?
- In equation 6: I don’t fully understand the logic of why 6 is a valid approach to learn a summary network of the real data. Can you explain this to me in more detail?
- I don’t quite understand the results in relation to the LPP and ACAUC metrics:
    - So 0 would be optimal for ACAUC, right?
    - And for LPP higher values would be better, right?
    - If so, SBI looks better than ROPE. I know this is probably not how these metrics should be read, so I am merely asking to clarify this point for me.
soundness: 3
strengths: - The paper investigates an important issue.
- The goals of this paper are clearly explained
- Detailed examples and ablation studies are performed.
summary: The paper details with the important issue of making neural posterior estimation methods more robust against model misspecification. The authors suggest to use a small set of labeled real-world data to calibrate the posterior inference in the face of a misspecified simulator. At the moment, I see some important weakness that make me cautios at recommending acceptance at this stage. But I am open and happy to change my mind during the discussion period.

After the discussion period, I have updated my score from 5 to 6.
weaknesses: - The main weakness I see is the following: How do we gain access to the *labeled* real-world observations? This seems like quite a strong assumption. I would need to understand if and when this assumption is justified before I would consider recommending acceptance.
- As an expert in SBI, I struggle to understand the introduction of optimal transport in Section 2.2. Adding some more details there might help readers which are not already experts in OT to understand what is going on.
- The applied metrics in the experiments are reasonable but still make it hard to see in which way we actually deviate from the targeted true posterior. The authors could perhaps also show something like bias in posterior mean and SD or some other easy to understand quantities.
- Curse of dimensionality: These scaling issues are mentioned by the authors themselves and seem to be severe in the number of model parameters(?) I would like to see the scaling laws discussed in more detail here if possible.

Meta-review: This paper proposes RoPE, a robust posterior estimation framework for simulation-based inference (SBI) that addresses the issue of model misspecification. It does so by using a small real-world calibration set of ground-truth parameter measurements and formalizing the “misspecification gap” using optimal transport (OT).  

I thank the authors and all the reviewers for their engaging discussions. It seems the reviewers agree the problem of misspecification in SBI is important and the idea of using OT in this setting is novel. Throughout discussions between authors and reviewers, some of the reviewers have increased their scores, while others remain hesitant. The main concerns about this paper relate to (1) the soundness of the approach and the definition of the problem; (2) the motivation for using OT; (3) hyper-parameter setting; (4) scalability and scope of the method and (5) a better experimental evaluation. 

There was a lot of back and forth between the authors and Reviewer YfYx about (1) to the point that the reviewer remains unconvinced and I think the paper will need much more clarity on this, and how this ties up to the methodological development. I think (2)-(5) all remain valid concerns. Furthermore, the paper does have a limited scope, as pointed out by the reviewers, there are many problems where labeled real-world parameter measurements are simply not available. Considering all this, I tend to agree with Reviewer YfYx and Reviewer 6VTN that the paper is not ready yet and will benefit from incorporating all the feedback given and addressing the most pressing concerns. Therefore, I recommend rejection.

Decision: Reject