Robust Amortized Bayesian Inference with Self-Consistency Losses on Unlabeled Data
Authors: Aayush Mishra
Daniel Habermann
Marvin Schmitt
Stefan T. Radev
Paul-Christian Bürkner
Venue: ICLR 2026 Poster (2026)
Keywords: Bayesian models
amortized inference
robust inference
self-consistency
semi-supervised learning
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We make amortized Bayesian inference more robust with only a small amount of unlabeled real data.

Abstract: Amortized Bayesian inference (ABI) with neural networks can solve probabilistic inverse problems orders of magnitude faster than classical methods. However, ABI is not yet sufficiently robust for widespread and safe application. When performing inference on observations outside the scope of the simulated training data, posterior approximations are likely to become highly biased, which cannot be corrected by additional simulations due to the bad pre-asymptotic behavior of current neural posterior estimators. In this paper, we propose a semi-supervised approach that enables training not only on labeled simulated data generated from the model, but also on unlabeled data originating from any source, including real data. To achieve this, we leverage Bayesian self-consistency properties that can be transformed into strictly proper losses that do not require knowledge of ground-truth parameters. We test our approach on several real-world case studies, including applications to high-dimensional time-series and image data. Our results show that semi-supervised learning with unlabeled data drastically improves the robustness of ABI in the out-of-simulation regime. Notably, inference remains accurate even when evaluated on observations far away from the labeled and unlabeled data seen during training.

Published as conference paper ICLR 2026
ROBUST AMORTIZED BAYESIAN INFERENCE WITH SELF- CONSISTENCY LOSSES ON UNLABELED DATA
Aayush Mishra"
Daniel Habermann"
Marvin Schmitt Independent Scientist
Department of Statistics TU Dortmund University, Germany
Department of Statistics TU Dortmund University, Germany
Stefan T. Radev
Paul-Christian Bürkner
Department of Cognitive Science Rensselaer Polytechnic Institute, USA radevs@rpi.es
Department of Statistics
TU Dortmund University, Germany paul www.rknerrgmaaa..........e.. com
ABSTRACT
Amortized Bayesian inference (ABI) with neural networks can solve probabilistic inverse problems orders of magnitude faster than classical methods. However, ABI not yet sufficiently robust for widespread and safe application. When performing inference observations outside the scope of the simulated training data, posterior approximations are likely to become highly biased, which cannot be corrected by additional simulations due the bad pre-asymptotic behavior of current neural posterior estimators. In this paper, we propose a semi-supervised approach that enables training not only on labeled simulated data generated from the model, but also on unlabeled data originating from any source, including real data. To achieve this, we leverage Bayesian self-consistency properties that can be transformed into strictly proper losses that do not require knowledge of ground-truth pa- rameters. We test our approach on several real-world case studies, including applications high-dimensional time-series and image data. Our results show that semi-supervised lear- ing with unlabeled data drastically improves the robustness of ABI in the out-of-simulation regime. Notably, inference remains accurate even when evaluated on observations far away from the labeled and unlabeled data seen during training. The code is available at https:// wtt.b. www./oaaeeeeeeeo.oococceeeeaaeoe f floo org sel conons wttp:/elllloosseeeeeeeeoeoe istnsssttteerrre
INTRODUCTION
Theory-driven computational models (mechanistic models) are highly influential across numerous branches of science (Lavin et al., 2021). The utility of computational models largely stems from their ability to fit real data æ and extract information about hidden parameters 0. Bayesian methods have been instrumental for this task, providing principled framework for uncertainty quantification and inference (Gelman et al., 2013). However, cold-standard Bayesian methods, such as Gibbs or Hamiltonian Monte Carlo samplers (Brooks al., 2011). remain notoriously slow. Moreover, these methods are rarely feasible for fitting complex models (Dax al., 2021) or even simpler models big data settings with many thousands data points single dataset (Blei et al., 2017), when thousands of independent datasets require repeated model re-fits (von Krause et al., 2022).
"Equal contribution

---

Published as a conference paper at ICLR 2026
In recent years, deep learning methods have helped address some of these efficiency challenges (Cranmer et al., 2020). In particular, amortized Bayesian inference (ABI; Gershman & Goodman, 2014; Le et al., 2017; Gonçalves et al., 2020; Radev et al., 2023; Gloeckler et al., 2023; Elsemüller et al., 2024; Eammit-Mangion et al., 2024) has received considerable attention for its potential to automate Bayesian workflows by replacing per-dataset inference with a learned inference model. By training generative neural networks on simulations generated from a probabilistic model, ABI learns a mapping from observations to posterior distributions, allowing subsequent inference on real data to be performed in near-instant time. This amortization trades an upfront computational cost during training for substantial gains in scalability and reusability, making Bayesian inference feasible in settings involving large datasets, intractable likelihood, or real-time constraints.
However, due to the reliance on pre-trained neural networks, ABI methods can become unreliable when applied to data that is unseen or sparsely encountered during training. In particular, posterior samples from amortized methods may deviate significantly from samples obtained with gold-standard MCMC samplers when there is a mismatch between the simulated training data and the real data (Ward et al., 2022; Schmitt al., 2023; Siahkoohi et al., 2023; Gloeckler et al., 2023; Frazier al., 2024). This lack of robustness limits the widespread and safe upplicability of ABI methods.
In this work, we propose a new robust semi-supervised approach to ABI. The supervised part learns from a "labeled" set of parameters and corresponding synthetic (simulated) observations, {0,r), while the unsuper- vised part leverages an "unlabeled" data set of real observations {x*} without parameters. In contrast to other methods aiming to enhance the robustness of ABI, our approach does not require ground-truth parameters l (Wehenkel et al., 2024), post hoc corrections (Ward et al., 2022; Siahkoohi et al., 2023), or specific adversarial defenses (Gloeckler et al., 2023), nor does it entail loss of amortization (Ward et al., 2022; Huang et al., 2023) or generalized Bayesian inference (Gao et al., 2023; Pacchiardi et al., 2024).
While Schmitt et al. (2024) introduced self-consistency (SC) losses to improve simulation efficiency in fully supervised, simulation-only setting, we extend self-consistency losses to unsupervised training on real observations, enabling robust inference without access to ground-truth parameters. We prove that self- consistency losses are strictly proper, target the true analytic posterior, and are agnostic to the data distribution. These properties allow them to be seamlessly integrated with any simulation-based loss for semi semi-supervised training without introducing a trade-off in the training objective. Empirical results on diverse tasks, including high-dimensional time series and images, show that our method preserves ABI's characteristic speed while achieving strong robustness, even with as few as four unlabeled real-world observations. Posterior estimates remain accurate and well-calibrated, generalizing beyond both labeled and unlabeled training data.
2 METHODS
2.1 BAYESIAN SELF-CONSISTENCY
Self-consistency leverages a simple symmetry in Bayes' rule to enforce more accurate posterior estimation even in regions with sparse data (Schmitt et al., 2024; Ivanova et al., 2024). Crucially, it incorporates the likelihood (when available) surrogate likelihood during training, thereby providing the networks with additional information beyond the standard simulation baasattiiio---ddddddee loss typically employed in ABI (see below).
Under exact inference, the marginal likelihood is independent of the parameters 0. That is, the Bayesian celf-connistency ratio of likelihood -rrior product and posterior is constant across any set of parameter values ((1), O(1),
p(x (()) p(0(1)
p( ((.)) p(O(L))
p (z)
(I)
p(((1) x)
p(((L) x)
However, replacing p(0 I) with a neural estimator q(0 | x) (and analogously for the likelihood) leads to undesired variance in the marginal likelihood estimates across different parameter values on the right-hand hand side.
2

---

Published as conference paper ICLR 2026
Hobs =0
Hobs = 1
Hobs=2 =
Mobs =3
Mobs =5
Hobs =8
Hobs = -11
dimension 2
5
0
5
0
5
0
5
0
5
0
5
dimension
NPE SC
NPE only
reference
Figure 1: Contour plot of the normal means problem using standard NPE (red) or our semi-supervised approach (NPE + SC, blue), with the analytic posterior in gray. Symbols indicate posterior mean estimates (red cross: NPE only; blue square: NPE SC; gray triangle: reference). Each subplot shows posterior inference on observed data that are increasingly distant from the labeled training data (Mprior 0). Only the first two dimensions of the 10-dimensional posterior are shown. While standard NPE collapses to zero variance for Pobs 2. adding the self-consistent loss preserves accurate posterior estimates even far beyond both training spaces (Jobs > 3). Training was performed using the default configuration (see Section 4.
Since this variance is proxy for approximation error, we can directly minimize via backpropagation along with any other ABI loss to provide further training signal and reduce errors guided by density information. Our proposed semi suni-supervise approach builds on these advantageous properties.
2.2 SEMI-SUPERVISED AMORTIZED BAYESIAN INFERENCE
The marginal likelihood identity in Eq. (1) is straightforward, but practically sampling-base methods (e.g.. MCMC) because they do not provide a closed-form for the approximate posterior density q() xx) In contrast, can readily evaluate q(0 ()) in ABI when using neural density estimator that allows efficient density computation (e.g., normalizing flows, Kobyzev et al. (2020)). Thus, we formulate family of semi-supervised losses of the form:
used in traditional
(q", h²) argmin ((o.))-n(0..) [S(g(0 h(x)), 0)] + X (pt [c (c(rrr [paa* p(0)) (00))))))))))) h(x)))
(2)
9.h
where a strictly proper score (Gneiting& & Raftery, 2007) and C' self-consistency score. The neural net- works to be optimized are a generative model and (optionally) summary network extracting lower dimen- sional sufficient statistics from the data. will call the first loss component, (( [S(g(0 h(x)), 0),, the (standard) simulation-basee learning. E.g., this is the maximum likelihood loss for normalizing flows (Kobyzev et al., 2020; Papamakarios et al., 2021) or vector-field loss for flow matching (Liu et al., 2023; Lipman et al., 2023). We will refer to the second loss component as the (Bayesian) self-consistency loss.
loss, it forms the basis for standard ABI approaches using simulation-based
In practice, we approximate the expectations in Eq. (2) with finite amounts of simulated and real training data. That is, for instances (0n, In) ~ p(6, x) and M instances xn ~ p" (x), we employ
N
M
C
p(x" (platn
em p(0m)
)
(q",h") = argmin 9.h
[S((00 (x)))),,
Σ m=1
(3)
N
M
q(0m h(x*))
where À is the self-consistency weight described in Appendix B. Asymptotically for - co, that is, for infinite training data generated from the simulator p(0,z)). I.), universal density estimator (Draxler et al., 2024)
3

---

Published as a conference paper at ICLR 2026
minimizing strictly proper simulation- loss (Gneiting & Raftery, 2007) is sufficient ensure perfect posterior approximation for any data. By this, we mean that the posterior approximation becomes identical to the posterior we would obtain if we could analytically compute p(0 2) = p(x 0)p(0)/p(z)) This analytic posterior is sometimes also referred to as "true" or "correct" posterior (of the specified statistical model). In practice, the posterior is rarely analytic, but we can still verify the accuracy of an approximation by comparing it with the results of a gold-standard approach (if available), such as a sufficiently long, converged MCMC run (Magnusson et al., 2024).
While neural posterior approximation is perfect asymptotically, its pre-asymptotic performance, that is, when training (() h(x)) only on a finite amount of simulated data, can become arbitrarily bad (Frazier et al., 2024; Schmitt al., 2023). For any data x" that is outside the data space implied by p(0, for instance, when the model is misspecified (see Appendix B for definition), the posterior approximation q( h(x")) may be arbitrarity far away from the analytic posterior p(0 x²)). As a result, a simulation-based loss is insufficient achieve robust ABI in practice. This is where the self-consistency loss comes in: As we will show, the latter greatly improves generalization to atypical data at inference time.
One particular choice for C' is the variance (with respect to 0 ~ pO)) of the log Bayesian self-consistency ratio (Schmitt et al., 2024):
plx= 0) p()
C
(arg=po(0) p(x" 10) lo) log {(0 h(x")),,
(4)
4(0 M(x*))
which is motivated by theoretical results by Köthe (2023), who show that minimizing the variance of the Bayesian self coffisonsitt ratio also minimizes the KL divergence between true and approximate posterior. (0) can be any proposal distribution over the parameter space. practice, two natural choices for Pc are either the prior p(0) or the current approximate posterior q.(O h(x")) at training iteration f. Because the largest contributions to C' are expected to be in high-density regions of the approximate posterior, we choose (:(0) h((")) for all of our experiments, see Appendix B for further details on practical considerations.
2.3 SELF-CONSISTENCY LOSSES ARE STRICTLY PROPER
Below, we discuss the strict properness of Bayesian senf-coonistency losses, showing that they can be used alongside simulation-basel losses without altering the target posterior. To simplify the notation, we denote posterior approximators simply as q(O) r) without considering architectural details such as the use of summary networks h(v). All theoretical results and their proofs remain valid if x is replaced by h(z), provided that the summary network is expressive enough to learn sufficient statistics from x.
Proposition 1. Let C be a score that is globally minimized if and only if its functional argument is constant across the support of the posterior p(0 x) almost everywhere. Then, C applied to the Bayesian self- consistency ratio with known likelihood
C (ple) p() (p(x)0)p(0)
(5)
q(O x))
is a strictly proper loss: It is globally minimized if and only if | x)) = p(0 ()) almost everywhere.
In particular, the variance loss (4) fulfills the assumptions of Proposition 1. Proposition 2. The loss (4) based on the variance of the log Bayesian self-consistency ratio is strictly proper if the support of (ff encompasses the support of p(0 x)-
The proofs of Propositions 1 and are provided in Appendix A. The strict properness extends to semi- supervised losses of the form (2), which combine standard simulation-banl and self self-connistency losses.
Proposition 3. Under the assumptions of Proposition 1, the sepi-supervised loss (2) is strictly proper for any choice of p² (z).
4

---

Published as a conference paper at ICLR 2026
Proposition 3 holds for every instantiation of C that satisfies Proposition 1. The proof of Proposition 3 follows immediately from the fact that the sum of strictly proper losses is strictly proper. Importantly, since Proposition 3 holds independently of p" (ar, holds both in the case of a well- sell--secified model, where p" (ax) = p(æ), and also in case of any model misspecification or domain shift where p" (æ) A p(a). That is, there is no trude-off' in the semi-supervised loss (2), since both loss components are globally minimized for the same target, that is, the analytic posterior.
other words, adding the SC loss does not alter the underlying statistical model and continues to target the analytic posterior of that model. As far as we aware, all other methods aimed at improving the robustness of ABI work by either explicitly or implicitly adjusting the statistical model or the target distribution (moving away from the analytic posterior), thus introducing some form of regularization. In contrast, SC losses should not be interpreted as regularizers: they are strictly proper losses that directly target the analytic posterior of the specified statistical model (i.e., Target in the taxonomy of Elsemüller et al., 2025).
Our proposed semi-supervised loss combines standard simulation- neural posterior estimation (NPE) and unsupervised self-consistency. Both target the analytic posterior, but in complementary ways: The standard NPE loss facilitates network training on simulated (labeled) data where ground-truth parameters are known. However, it does not ensure posterior accuracy on the empirical data distribution, particularly in regions where simulated data may be sparse or unrepresentative due to model misspecification. The self-consistency loss complements NPE by enforcing Bayes' rule on observed data and explicitly rewards correct marginal likelihood estimation. This ensures alignment between prior, posterior and likelihood on the empirical data distribution-a property not guaranteed by simulation- baaeed losses alone.
For completeness, we also define strictly proper self-consistency losses for likelihood estimators, rather than posterior approximators.
Proposition 4. Suppose the posterior p(0 x)) is known and the likelihood is estimated by g(a: 0). Then, under the assumptions of Proposition 1, Bayesian celf-coonsstency ratio losses of the form
q(x () p(0)
(6)
p(0 x)
are strictly proper: They are globally minimized if and only if q(x 10) = p(r | 0) almost everywhere
The proof of Proposition 4 proceeds in the same manner as for Proposition just exchanging likelihood and posterior. Clearly, strict properness does not necessarily hold if both posterior and likelihood are unknown or approximate because any pair of approximators q(0 x)) and q(x () that satisfy q(O x)) x q(z 0) p(0)
minimize the SC loss. For example, the choices q(0 x)) p(0) and q(x: 0) x minimize the SC loss, but may be arbitrarily far from their actual target distributions p(0 x)) and p(a 0). respectively.
In other words, if both likelihood and posterior are unknown, the SC loss has to be coupled with another loss component, such as the maximum likelihood loss, to enable joint learning of both approximators g(O x)) and <(x 0). Nevertheless, SC still yields notable improvements: in our experiments, the semi-supervised loss (2) considerably enhanced the robustness of ABI even when both the posterior and likelihood are unknown.
3 RELATED WORK
The robustness of ABI and simulation- inference methods more generally has been the focus of multiple
recent studies (e.g., Frazier al., 2020: Frazier & Drovandi, 2021; Frazier et al., 2024; Dellaporta et al., 2022; Ward et al., 2022; Cemgil ct al., 2020; Gloeckler et al., 2023; Huang et al., 2023; Gao et al., 2023; Siahkoohi et al., 2023; Kelly al., 2024; Wehenkel et al., 2024; Pacchiardi et al., 2024; Schmitt al., 2023). These efforts can be broadly classified into two categories: (a) analyzing or detecting simulation gaps and (b) mitigating the impact of simulation gaps posterior estimates. Since our work falls into the latter category,
5

---

Published as a conference paper ICLR 2026
we briefly discuss methods aimed at increasing the robustness of fully amortized approaches. E.g., Gloeckler et al. (2023) explore efficient regularization techniques that trade off some posterior accuracy to enhance the robustness of posterior estimators against adversarial attacks. Ward al. (2022) and Siahkoohi al. (2023) apply post hoc corrections based on real data, utilizing MCMC and the reverse KL-divergence, respectively,
Differently, Gao et al. (2023) propose a departure from standard Bayesian inference by minimizing the expected distance between simulations and observed data, akin to generalized Bayesian inference with scoring rules (Pacchiardi et al., 2024). Perhaps the closest work in spirit to ours is Wehenkel et al. (2024), which introduces the use of additional training information in the form of (labeled) calibration set (x* 00) that contains observables from the real data distribution as well as the corresponding ground-truth parameters.
In contrast to the methods above, our approach (a) avoids trade-off's between accuracy and robustness, (b) requires no modifications to the neural estimator after training, therefore fully maintaining inference speed, (c) affords proper Bayesian inference, and (d) does not assume known ground truth parameters for calibration set. Thus, it can be viewed as one of the first instantiations of semi-supervised ABI.
4 CASE STUDIES
4.1 MULTIVARIATE NORMAL MODEL
We first illustrate the usefulness of our proposed self-consistency loss on a controllable toy problem (Schmitt et al., 2023). The prior and likelihood are given by
~ Norrio o
prior ID),
x (k) ~ Normal (0, oik Ip)
(7)
The parameters A E IDD are sampled from a D-dimensional multivariate normal distribution with mean vector Uprior and diagonal covariance matrix o prior Ip. Here, we fix jpprior = 0 and prior 1. On this basis, K independent, synthetic data points (*()) E RR. are sampled from a D-dimensional multivariate normal distribution with mean vector 0 and diagonal covariance matrix Th ID. We set oi = such that the total information in remains constant across values of K, which simplifies comparisons across observations of varying numbers of data points. Refer to Appendix C for more details on training and neural architectures.
In our numerical experiments, we study the influence of several aspects of the normal model on the perfor- mance of NPE. To prevent combinatorial explosion, we vary the factors below separately, with all other factors fixed to their default configuration (highlighted in bold): (1) parameter dimensionality (D 2, 10, 100). (2) number of unlabeled observations for the self-consistency loss {>/ ) (M = 1, 4, 32), (3) mean j of the unlabeled observations =i, (n² 0, 2, 5). (4) inclusion of summary network (K 10) not 1). (5) likelihood function (known, estimated),
m==1
(K =
Results In Figure we depict the results obtained from (a) standard NPE (trained on the simulation- based loss only), (b) our semi-supervised NPE + SC (with self-consistency loss on known likelihood), and (c) the stad-standard (analytic) reference. We see that standard NPE already completely fails for Tobs (Mobs = 2, 0.01 and subsequently also for any larger values Hobs > 2. In contrast, our semi-supervised approach achieves almost perfect posterior estimation. This holds true even in cases where Tobs is multiple standard deviations away from all the training data, that is, from both the labeled dataset {(0,, xx) > and the unlabeled dataset {x", 1M m= These results indicate that the self-consistency criterion can provide strong robustness gains even far outside the typical space of training data.
In Figure 2, we report the maximum mean discrepancy (MMD) between the approximate and true posterior the factors parameter dimensionality and number of unlabeled observations. When varying the parameter dimensionality (Figure 2a), including the SC loss yields nearly perfect posterior approximation up to 10 dimensions, even with extreme deviations from the initial training data. It also significantly improves accuracy
6

---

Published as a conference paper ICLR 2026
parameter dimensionality
size of unlabeled training data
15-
15
MMD
10-
MMD
10-
5-
5-
10
10
observed data location
observed data location
- 2
10
100
SC) +- (no SC)
-+- 100 SC)
(no SC)
32
(a): Posterior distance between approximate and true poste- rior for varying parameter dimensionality.
(b): Posterior distance hetween approximate and true poste- rior for varying unlaheled training data sizes.
Figure 2: Posterior distance quantified by maximum mean discrepancy (MMD) to the analytic posterior for variations of the default configuration. Errorbars show +1 SDs over 10 model refits.
in the 100 parameter scenario. The dataset size factor (Figure 2b) shows robust gains, with clear improvements over NPE even when using as few as four unlabeled observations (versus 1024 labeled ones). In Figure 6 (Appendix D), we report posterior mean and standard deviation bias as well maximum mean discrepancy for all the above factors. Varying the mean Pobs of the new observations shows that, as long as the data used for evaluating the SC loss not identical to the training data (i.e., as long as Hobs * 0), including the SC loss component enables accurate posterior approximation far outside the typical space of the training data.
Figure in Appendix D, we see that the benefits of self- coffissonsttncc persist when the posterior is conditioned on more than one data point per observation (K 10), that is, in the presence of a summary network. Further, we still see clear benefits of adding the seff-consistency loss even when the likelihood is estimated by neural likelihood approximator (x:) 0), trained jointly with the posterior approximator q(0 x) on the same training data. However, with an estimated likelihood, posterior bias, especially bias the posterior standard deviation, and MMD distance to the true posterior are larger than in the known likelihood case. We also compared our semi- supervipe approach with other robustness methods- NPE- DANN, NPE-MMD, and NNPE. The results described in Appendix D.1 demonstrate the benefits of our approach over other methods.
4.2 FORECASTING AIR PASSENGER TRAFFIC: AN AUTOREGRESSIV MODEL WITH PREDICTORS
We apply our semi-supervised loss to analyze trends in European air passenger traffic data provided by Eurostat (2022a;b;c). This case study highlights that the strong robustness gains also occur in real-world scenarios and model classes that are challenging to estimate in a simulation-based inference setting. We retrieved time series of annual air passenger counts between 15 European countries (departures) and the USA (destination) from 2004 to 2019 and fit the following autoregressive process of order 1:
Ujatt Normal(a, + Vj.+,3 + U, + ",,)),
(8)
where the target quantity Mj,4+1 is the difference in air passenger traffic for country j between time t + and t. To predict yj,++1 we use two additional predictors: uz,a is the annual household debt of country j time t, measured in % of gross domestic product (GDP) and Wy, is the real GDP per capita. Training relies on a small simulation budget of N = 1024, with the self-consistency loss evaluated on real data from M (4, 8, 15) countries. Further details model and training are provided in Appendix E.
Results In Figure 3, we show exemplary results from standard NPE, our semi-supervised NPE + SC (M = 8), and Stan (Carpenter et al., 2017) reference. We see that NPE is highly inaccurate for many
7

---

Published as conference paper ICLR 2026
NPE only NPE SC reforence
0,4
0.2
0.0
0.2
-0.4
PL CZ NL PT AT FR BE SE DE
DK IE EL
ES
Country code
Figure 3: Comparison of posterior estimates for 15 countries (ISO 3166 alpha-2 codes) among standard NPE (red), NPE + SC (blue), and Stan (reference in gray; Carpenter et al., 2017). Central 50% (thick lines) and 95% (thin lines) posterior intervals of the autoregressive component are shown, sorted by lower 5% quantile per Stan (i.e., established benchmark). The SC loss was evaluated on data from M = 8 countries during training, greatly enhancing ABI's robustness in both mo-misppecifiatioo scenarios and real-data evaluations.
countries, whereas NPE SC is in strong agreement with the reference for but one country. As shown Table I, adding the self-consistency loss (M 8) strongly improves posterior estimates for all five parameters across all metrics, on average across countries. The complete results can be found Appendix
4.3 HODGKIN--UXLER MODEL OF NEURON ACTIVATION
To investigate the effect of the self-consistency loss on a model involving nigh-dimensional data, we evaluate our approach on the (2023). The Hodgkin-Hux model is classical model in neuroscience to describe neuron activation via a set ordinary differential equations. In brief, the model has parameters (electrical conductances of different ion channels, membrane capacitance and reversal potentials), and the output with observation index is 200-dimensional time series of the membrane potential. full definition of the model, as well as the training setup and description of the neural architectures, provided in Appendix G.
model, which was previously used in an ABI setting by Gloeckler
To facilitate network training, all parameters are transformed to follow standard normal distributions through appropriate transformations. For example, the parameter 9Na with marginal prior 9Na (ogNorma((log(110) 0.1) is transformed via Zyna = (log(gNu) log(110))/O.1 We denote the full set of transformed model parameters by 0. Training is performed with simulation budget of N 32,768. For each loss evaluation, the self-consistency loss is computed on a random subset of 32 samples drawn from a pool of M 1,024 unlabeled observations. These are generated by first sampling 0 ~ Normal 2). applying the inverse transformations to recover the original parameter scale, and then simulating time series of the membrane potential as above.
Results To assess the benefits of our approach in the out-of-distribution setting, Figure 4a shows posterior predictive samples inferred from data simulated with 0 ~ Normal(-2,1) 1). This contrasts with training data from 0 ~ Normal(0, 1) and self-consistency evaluation data from 0 Normal 2). When training without the self-consistency loss, the neural posterior density estimator produces samples inconsistent with the observed data, while incorporating the loss yields accurate predictions (see Appendix for further results). To quantify this, Figure 4b reports mean absolute bias differences between the estimators. The self-consistency loss consistently and strongly improves predictions across the majority of time series. Even in the worst cases, is least competitive with the estimator trained without the self- self-consistency loss.
8

---

Published as conference paper ICLR 2026
NPE only
membrane potential [mV]
80
count
60
NPE+SC SC
40
50
-50
10
20
30
40
50
60
0
-5
time [ms]
MAB(NPE)) MAB(NPE+SCC
(a): Posterior predictive samples without (top row) and with (bottom row) self-consistency loss.
(b): Quantitative evaluation of predictive bias.
Figure 4: (a) Posterior predictive samples (gray) inferred from an out-of-simulation dataset (black). NPE only produces highly biased predictions while NPE+SC yields accurate results. (b) Histogram of the mean absolute bias (MAB) difference of posterior predictions computed for 1000 out-of-simulation datasets. NPE+SC has lower bias than NPE for almost all datasets.
4.4 BAYESIAN DENOISING OF MNIST IMAGES
Finally, we demonstrate the utility of our self-consistency loss based approach in a high-dimensional (for ABI) image denoising task, following a set-up similar to Elsemüller et (2025). The parameter vector 0 E 784 is the flattened image, and the observation x E 787 is blurred version generated by simulated noisy camera. We assume an implicit prior p(0) defined by a generative model trained on blurred MNIST images of the digit "0"" (LeCun et al., 1998), and an implicit likelihood p(x[0) implemented by reapplying the same blur to each 0. This creates a challenging neural posterior posterior-lkkelhood likelihood estimation (NPLE) problem.
To construct the training set, blur all digit "0" images in the MNIST training set with fixed Gaussian filter and train a neural network to (i) sample blurred simulated images {0' J ~ p(0) and (ii) evaluate their log- probabilities. For each sample 8', then generate observation x²' p(x](") by reapplying the same Gaussian blur, yielding N 12000 pairs (#", x²) to train posterior and likelihood networks. For the SC loss. use 400 held-out MNIST test images blurred only by the likelihood model, (i.e., additional prior blur), inducing deliberate prior misspecification and enabling us to asses the robustness of SC in NPLE. More details about the training set-up and model architecture can be found in the Appendix I.
Results We perform inference on separate subset of 580 MNIST test images. Figure and Figure 16 in Appendix J show randomly chosen examples alongside the posterior mean and standard deviation maps computed from 500 posterior samples. Reconstructions using the self-consisenncy approach (NPLE+SC) smoother and more faithful to the ground-truth, whereas NPLE reconstructions are pixelated and blurry. Moreover, NPLE+SC yields coherent uncertainty maps, with elevated variance localized along digit contours where edge ambiguity is expected. In contrast, NPLE estimates exhibit scattered, patchy uncertainty across the digit and background. Additional examples are reported in Appendix J.
5 DISCUSSION
We demonstrated that Bayesian self-consistency losses significantly increase the robustness neural amor- tized Bayesian inference (ABI) out-oo-simulation data. Accurate inference outside the training distribution,
9

---

Published as conference paper ICLR 2026
Parameter (6)
D
Observation
Mean (NPLE SC)
Std Dev (NPLE SC)
(NPLE)
Figure 5: Example of denoising results for MNIST images of digit "0" in the held-out test set. The first row shows ten randomly selected MNIST images (0), the second row depicts images after applying the Gaussian blur (I), third and fourth rows depict mean and standard deviation (SD) of 500 posterior samples estimated from the corresponding blurry observations using NPLE+SC, and the fifth and sixth rows depict the mean and SD of 500 posterior samples using NPLE only. Incorporating SC loss significantly improves denoising.
such as in the presence of model misspecification, has long posed a major challenge for ABI. While self- consistency was originally introduced to improve training efficiency with slow simulators (Schmitt et al., 2024; Ivanova et al., 2024), it had not been previously explored as a remedy to simulation gaps. Existing supervised ABI approaches have been known to dramatically fail in such cases (Schmitt et al., 2023; Gloeckler et al., 2023; Huang et al., 2023), as we also illustrated in our experiments. In contrast, when optimizing for self-consistency on unlabeled aut-of-simulation data, we obtained nearly unbiased posterior estimation far beyond the training distribution. The strong robustness gains persisted even in models with several hundred parameters. However, using an approximate (neural) likelihood has not yet matched the robustness of the known-likelihood case. Finally, as self-consistency losses do not require data labels (i.e., true parameter values), we can use any amount of real data during training to improve the robustness of ABI.
Limitations and future directions Our criaace--bssed self-consistency loss relies fast density eval- uations during training, keeping times competitive. This makes free-form methods such as flow matching (Lipman et al., 2023) score-based diffusion (Song et al., 2020) less practical due their need for numerical integration. As a result, efficient self-consistency losses for free-form flows, along with joint learning of posteriors and very high-dimensional likelihoods, remains an open avenue for future research.
Our semi-supervised approach can also be combined with post-hoc correction methods such as Pareto- smoothed importance sampling (Vehtari al., 2024). When applied to NPE alone, importance sampling often fails because the learned posterior can deviate substantially from the analytic target (Li et al., 2024). By contrast, our approach produces posterior approximations that are much closer to the analytic posterior, making subsequent importance sampling more effective. This combination therefore has the potential to substantially expand the range of scenarios which accurate estimation of the analytic posterior feasible.
10

---

Published as conference paper ICLR 2026
ACKNOWLEDGMENTS
Daniel Habermann and Paul Bürkner acknowledge support of the Deutsche Forschungsgemcinscaff (DFG, German Research Foundation) Projects 508399956 and 528702768. Paul Bürkner further acknowledges support of the DFG Collaborative Research Center 391 (Spatio-Temporal Statisties for the Transition of Energy and Transport) 520388526. Stefan T. Radev is supported by NSF under Grant No. 2448380. Aayush Mishra acknowledges the computing time provided on the Linux HPC cluster Technical University Dortmund (LiDO3).
REPRODUCIBILITY STATEMENT
We provide all details necessary to reproduce our results. Hyperparameter settings and network architectures are described in Section and Appendices C, D. E, G, and The code to run all experiments is provided in the abstract. Hardware specifications are reported in Appendix K. All theoretical results are accompanied by complete proofs in Section and Appendix A.
REFERENCES
Simon Alexanderson and Gustav Eje Henter. Robust model training and generalisation with studentising flows. arXiv preprint arXiv:2006.06599, 2020.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians.
Journal of the American statistical Association, 112(518)8559 207,
Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Handbook of markov chain monte carlo. CRC press, 2011.
Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Mareus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of statistical software, 7-:1-32, 2017.
Taylan Cemgil, Sumedh Ghaisas, Krishnamurthy Dvijotham, Sven Gowal, and Pushmeet Kohli. The autoencoding variational autoencoder. Advances in Neural Information Processing Systems, 33:15077- 15087, 2020.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117445)30055-30066, 2020.
Maximilian Dax, Stephen R Green, Jonathan Gair, Jakob Macke, Alessandra Buonanno, and Bernhard Schölkopf. Real-time gravitational wave science with neural posterior estimation. Physical review letters, 127(24):241103 2021.
Charita Dellaporta, Jeremias Knoblauch, Theodoros Damoulas, and François-Xavier Briol. Robust bayesian inference for simulator-based models via the mmd posterior bootstrap. In International Conference on Artificial Intelligence and Statistics, pp. 943-970. PMLR, 2022.
Felix Draxler, Stefan Wahl, Christoph Schnörr, and Ullrich Köthe. On the universality of coupling-based normalizing flows. arXiv preprint arXiv:2402.06578. 2024.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Advances in neural information processing systems, 32, 2019.

---

Published as conference paper ICLR 2026
Lasse Elsemüiller, Hans Olischläger, Marvin Schmitt, Paul-Christian Bürkner, Ullrich Köthe, and Stefan Radev. Sensitivity amortized Bayesian inference. Transactions on Machine Learning Research (TMLR), 2024.
Lasse Elsemüller, Valentin Pratz, Mischa von Krause, Andreas Voss, Paul-Christian Bürkner, and Stefan T Radev. Does unsupervised domain adaptation improve the robustness of amortized bayesian inference? a systematic evaluation. arXiv preprint arXiv:2502.04945 2025.
Eurostat. International extra-eu air passenger transport by reporting country and partner world regions and countries, doi: 10.2990//viiaa...../a.aaaaaaaaa pacellcapppeeeeeee 2022a.
Eurostat. Household debt, consolidated including Non-profit institutions serving households % GDP, doi:10..908/TIPSDDD. 2022b.
Eurostat. Real gdp per capita, doi: (0.2908/SDG_08_10 2022c.
David T. Frazier and Christopher Drovandi. Robust Approximate Bayesian Inference With Synthetic Likelihood. Journal of Computational and Graphical Statistics, 30(4):958-976, October 2021. doi: (0000/00616000202211
David T. Frazier, Christian P. Robert, and Judith Rousseau. Model misspecification in approximate Bayesian computation: consequences and diagnostics. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(2):421-444, April 2020. doi: 10.1 10.111 111111/ta 1/rssb. 12356.
David T. Frazier, Ryan Kelly, Christopher Drovandi, and David Warne. The Statistical Aceuracy of Neural Posterior and Likelihood Estimation. arXiv preprint, 2024. doi: 10.48550/arXiv.44 1.12688....
Richard Gao, Michael Deistler, and Jakob Macke. Generalized bayesian inference scientific simulators via amortized cost estimation. Advances in Neural Information Processing Systems, 36:80191-80219, 2023.
Andrew Gelman, John B Carlin, Hal S Stern, David Dunson, Aki Vchtari, and Donald B Rubin. Bayesian Data Analysis. CRC Press, 2013.
Samuel Gershman and Noah Goodman. Amortized inference in probabilistic reasoning. Proceedings of the annual meeting of the cognitive science society, 36, 2014.
Manuel Gloeckler, Michael Deistler, and Jakob H Macke. Adversarial robustness of amortized bayesian inference. In Proceedings of the 40th International Conference on Machine Learning, 11493-11524, 2023.
Tilmann Gneiting and Adrian E Raftery. Strietly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359-378 2007.
Pedro J Gonçalves, Jan-Matthi Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri. William Podlaski, Sara A Haddad, Tim P Vogels, et al. Training deep neural density estimators to identify mechanistic models of neural dynamics. Elife, 9re56261, 2020.
Gretton, K. Borgwardt, Malte Rasch, Bernhard Schölkopf, and AJ Smola. Kernel Two-Sample Test. The Journal of Machine Learning Research, 2012.
Daniel Habermann, Marvin Schmitt, Lars Kühmichel, Andreas Bulling, Stefan T. Radev, and Bürkner. Amortized Bayesian Multilevel Models. arXiv preprint, 2024. doi: 10.48550/arXiv.2444 240. 13230.
12

---

Published as a conference paper ICLR 2026
Daolang Huang, Ayush Bharti, Amauri Souza, Luigi Acerbi, and Samuel Kaski. Learning robust statistics for simulation-ball based inference under model misspecification. Advances in Neural Information Processing Systems, 36:7289 -310, 2023.
Desi R Ivanova, Marvin Schmitt, and Stefan T Radev. Data-efficient variational mutual information estimation via Bayesian self-consistency, In NeurIPS 2024 Workshop on Bayesian Decision-making and Uncertainty, 2024.
Ryan Kelly, David Nott, David T Frazier, David Warne, and Chris Drovandi. Misspecification robuust sequential neural likelihood for simulation-based inference. Transactions Machine Learning Research, (024(uun).Arriille-mmbbe 2024.
Diederik Kingma and Max Welling. Auto-encoding variational bayes. arXiv preprint arXiv: 1.Xi::3312.6114 2013.
Ivan Kobyzev, Simon JD Prince, and Mareus A Brubaker. Normalizing flows: An introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, (3(11):3964-397 2020.
Ullrich Köthe. A review of change of variable formulas for generative modeling. arXiv preprint arXiv:2308.0266 2023.
Alexander Lavin, David Krakauer, Hector Zenil, Justin Gottschlich, Tim Mattson, Johann Brehmer, Anima Anandkumar, Sanjay Choudry, Kamil Rocki, Atulmm Güneş Baydin, et al. Simulation intelligence: Towards generation of scientific methods. arXiv preprint a 2021.
Tuan Anh Le, Atilim Gunes Baydin, and Frank Wood, Inference compilation and universal probabilistic programming. In Artificial Intelligence and Statistics, pp. 1338-1348. PMLR, 2017.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient- based learning applied to document recognition. Proceedings of the IEEE, 86( 1):2278-2324 1998.
Chengkun Li, Aki Vehtari, Paul- ChrisChist Bürkner, Stefan T Radev, Luigi Acerbi, and Marvin Schmitt. Amortized bayesian workflow. arXiv preprint arXiv:2409.04332. 2024.
Yaron Lipman. Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In The Eleventh International Conference on Learning Representations. 2023.
Xingchao Liu, Chengyue Gong, and Qiang Liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. In The Eleventh International Conference Learning Representations (ICLR), 2023.
Mâns Magnusson, Jakob Torgander, Paul-Christian Bürkner, Lu Zhang, Bob Carpenter, and Aki Vehtari. posteriordb: Testing. Benchmarking and Developing Bayesian Inference Algorithms. arXiv preprint, 2024. doi: 10.48550/azXiv.4000 .04967.
Lorenzo Pacchiardi, Sherman Khoo, and Ritabrata Dutta. Generalized bayesian likelihood- inference. Electronic Journal of Statistics, (8(2):362833688 2024.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshmi- narayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1-64, 2021.
Stefan T Radev, Marvin Schmitt, Valentin Pratz, Umberto Picchini, Ulirich Köthe, and Paul-Christian Bürkner. Jana: Jointly amortized neural approximation of complex bayesian models. In Uncertainty in Artificial Intelligence, pp. 1695-1706. PMLR, 2023.
13

---

Published as a conference paper ICLR 2026
Marvin Schmitt, Paul-Christian Bürkner, Ullrich Köthe, and Stefan T Radev. Detecting model misspecifi- cation in amortized bayesian inference with neural networks. In DAGM German Conference on Pattem Recognition, pp. 541-54155 1557.55..00 Springer, 2023.
Marvin Schmitt, Desi Ivanova, Daniel Habermann, Paul-Christian Bürkner, Ullrich Köthe, and Stefan Radev. Leveraging self-consistency for data-efficient amortized Bayesian inference. In Proceedings of the 41st International Conference on Machine Learning, pp. 03721143744. PMLR, 2024.
Ali Siahkoohi, Gabrio Rizzuti, Rafael Orozco, and Felix Herrmann. Reliable amortized variational inference with physies-based latent distribution correction. Geophysics, 8833) 2023.
Yang Song. Jascha Sohk-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.134500 13456 2020.
Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. Pareto smoothed importance sampling. Journal of Machine Learning Research, 25(72):1-58, 2024. URL http:// fm1r org/
papers/v25/19 556.556666 html.
Mischa von Krause, Stefan T Radev, and Andreas Voss. Mental speed is high until age 60 as revealed by analysis of over million participants. Nature human behaviour, 6(5):700-708, 2022.
Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian Schmon. Robust neural posterior estimation and statistical model criticism. Advances in Neural Information Processing Systems, 35:33845-33859, 2022.
Antoine Wehenkel, Juan L Gamella, Ozan Sener, Jens Behrmann, Guillermo Sapiro, Marco Cuturi, and Jörn-Henrik Jacobsen. Addressing misspecification in simulation-based inference through data-driven calibration. arXiv preprint arXiv:2405.08719 2024.
Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexan- der J Smola. Deep sets. In Advances in Neural Information Processing Systems, 2017.
Andrew Zammit-Mangion, Matthew Sainsbury-Dake and Raphaël Huser. Neural methods for amortized inference. Annual Review of Statistics and Its Application, 12, 2024.
14

---

Published as conference paper ICLR 2026
APPENDIX
A PROOFS
Proof of Proposition By assumption, is globally minimized and only if
p( 0) p() <(0 x
A
(9)
for some constant A (independent of 0) almost everywhere over the posterior's support. Accordingly, any approximate posterior solution q(0 x) that attains this global minimum has to be the form
(0 x) p(æ 0) p(0) /A
(10)
By construction, q( æ) is a proper probability density function, so it integrates to 1. It follows that
1= |
f 4(0 x) do = p( ())(() d0 A p(x) / A.
(11)
Rearranging the equation yields A p(z) and thus
q x) - p(= 0)p(0)p((a) p(0 xx)
(12)
almost everywhere.
Proof of Proposition 2. The variance distribution PC (0) reaches its global minimum (i.e., zero), if and only its argument is constant across the support of (0). Because the log is a strictly monotonic transform,
logp(>" () log ((0) log q(0 xx))
log
(13)
for some constant A implies
p(a 0) p()) A,
(14)
( x))
which is sufficient satisfy the assumptions of Proposition 1.
FREQUENTLY ASKED QUESTIONS
Q: How do you define model-misppciffcationn, simulation-gap, and domain shift?
We define model misspecification as a mismatch between an unknown (true) distribution of data .= ~ p" (x) coming from complex, real-world system, and the marginal likelihood of data p(x[M) under some model M. Concretely, we call a model misspecified
p'(x) + p( IM) f( p(x 10M) p(0 IM) d0, Va E X.
Within this framework, a simulation gap refers to the divergence between the empirical distribution of observed data p² (Ir) and the empirical distribution of simulated data p(z), where p(z) is obtained via finite samples x p(= 0, MM) for 0 ~ p(0 M). That is, even when the model not misspecified, limitations in create gaps in support.
simulation coverage
A domain shift refers to any distributional discrepancy between the simulated training data (x ~ p(x MM) and the observed data x= (x)). This includes, but is not limited to, model misspecification and
15

---

Published as a conference paper at ICLR 2026
simulation gaps, and encompasses broader generalization challenges faced when deploying amortized inference methods outside the simulation domain.
Q: What is the intuition behind using self-consistency make the inference robust?
In well-specified setting, the simulator accurately captures the true data- generating process and the learned posterior matches the true posterior. In this case, neural networks only suffer from approximation error due finite data or model capacity. Under model misspecification, the observed data are drawn from distribution p" (x) F p(z), and hence the model-implied posterior p(0 x) x p(æ (0) does not correspond to the posterior under p" (I:) but under p(x)..
such settings, neural networks face two types of error: approximation error, as before, due to imperfect learning, and extrapolation error, because the network encounters observations p" (2) that may lie outside the simulation support, and thus lacks guidance on how to sample from the posterior in those regions. In this case, minimizing the self-consistency loss encourages the learned posterior q(0 x) to remain internally consistent with the assumed model p(a, O) when evaluated on empirical samples x ~ p" (z). It enforces that the learned posterior q() I) remains consistent with the model assumptions, effectively guiding the network to extrapolate in a model-coherent way. That is, even when the data are out-of-distributie relative to the simulator, the network learns to sample in way that is self-consistent under the (misspecified) model posterior, improving robustness despite model mismatch.
Q: Why is the self-consistency loss coupled with standard simulation-based loss when self-consistency is strictly proper?
In principle, minimizing the self-consistency loss sufficient to recover the true posterior, as guaranteed by Proposition 1. However, in practice, we observe that using the SC loss as the sole training objective leads to initialization and convergence issues. This is primarily due to the density-based nature of the SC loss. At initialization, when the approximate posterior and/or likelihood are poorly calibrated, the self-consistency loss can take on extreme values, often resulting in numerical instabilities. such as exploding gradients. Consequently, optimization tends to diverge or fail to make meaningful progress. We observe that the SC loss in tandem with the NPE loss provides a more stable and effective learning signal. The amortized loss anchors the model early by using supervised simulated data, while the SC loss complements it by enforcing consistency on observed data, especially in sut-ol-simllation regimes.
Q: What is the motivation behind using the variance over the log Bayesian self-consistency ratio for C?
The choice of using the variance of the logarithm of the Bayesian self-consistency ratio, evaluated over p (0) (Eq. 4), is motivated as follows: Köthe (2023) (Eq. 33) shows that second-order Taylor expamsion of the regular NPE loss function yields KL(p(O x) ll) =)) R 2p(x)² Var (p(0)p(z10)) ). From this, it follows that minimizing Var(...) . also minimizes the KL divergence between true and approximate posterior. Because the denominator in Var(...) can cause numerical instabilities due to vanishingly small values, we instead target Var(log(... )). Minimizing this quantity via a loss function still minimizes the KL. divergence between the true and approximate posterior, which follows from Proposition 1.
(()la))
Q: What is the motivation behind using the approximate posterior 4:(0 h(**)) at training iteration / for (0))
Rearranging Bayes' rule shows that the marginal likelihood p(x) 2)0))(2 =( is constant for all 0. If we
replace the true posterior p(0 I:) by an approximate posterior q(0 | x), large fluctuations in the marginal likelihood landscape are expected to be in high-density regions of the approximate posterior. We therefore
choose the approximate posterior as the proposal distribution.
The choice of using the variance of the logarithm of the Bayesian self-consistency ratio, evaluated over p ~ (0) (Eq. 4), is motivated as follows: Köthe (2023) (Eq. 33) shows that second-order Taylor expansion
16

---

Published as conference paper ICLR 2026
of the regular NPE loss function yields (p((
From this, follows that minimizing Var(...) also minimizes the KL divergence between true and approximate posterior. Because the denominator in Var(... can cause numerical instabilities due to vanishingly small values, we instead target Var(log(.... }). Minimizing this quantity via a loss function still minimizes the KL divergence between the true and approximate posterior, which follows from Proposition
x)llq(0 x)) 2p(z)² (0
p(0)p(a
Q: How to choose the self-consistency weight X?
From the prope-pro of self-consistency losses, follows that À does not control trade-off between the standard NPE loss and the SC loss, as both optimize for the same training objective. Instead, its main role is to stabilize training by preventing initialization and convergence issues. For experiments with low-dimensional parameter spaces, a default choice of À typically works well. However in high-dimensional settings, such as the image-denoising experiment, the approximate posterior can be particularly bad in the early epochs, leading unstable training. In such cases, we found it to be effective to use linear schedule for À that keeps it at for the first few epochs and then gradually increases to 1.
C DETAILED SETUP OF THE MULTIVARIATE NORMAL CASE STUDY
From the multivariate normal model described in Section 4.1, simulate a labeled training dataset with budget of N = 1024, that is, N independent instances of 0n (the "labels") with corresponding observations ( In = { / each consisting of K data points. This labeled training dataset {(On, ] is used for optimizing the standard simulation-bol
loss component. The self-consistency loss component is optimized
on an additional unlabeled dataset {= JM of M 32 independent sequences x {am =k=1 which, for the purpose of this case study, are simulated from
=(k)
(15) Since the self-consistency loss does need labels (i.e., the true parameters having generated wr))). we could have also chosen any other source for x², for example, real-world data. Within each training iteration t, the variance term within the self-consistency loss was computed from L = 32 samples g(l) qe(( xx)) from the current posterior approximation.
[(k] ~ Normal(p", Ip).
To evaluate the accuracy and robustness of the NPEs, we perform posterior inference on completely new (k) 1K observations Tobs (z (oobs/ae 1== each consisting of K independent data points sampled from
(16) The mean values Mobs E 10. 11} are progressively farther away from the training data. While conceptually simple and synthetic, this setting is already extremely challenging for simulation- algorithms because of the large simulation gap (Schmitt et al., 2023): standard NPEs are only trained on (labeled) training data that are several standard deviations away from the observed data the model sees at inference time.
(k) Fobs Normal(Hobe do = 0.01/p)
inference
The faithfulness of the approximated posteriors q( Tobe) are assessed by computing the bias in posterior mean and standard deviation as well as the maximum mean discrepancy (MMD) with a Gaussian kernel (Gretton et al., 2012) between the approximate and true (analytie) posterior.
The analytic posterior for the normal means problem is a conjugate normal distribution
p(0 Tobs) Normal (Apost Ip),
(17)
where Mpost is a D-dimensional posterior mean vector with elements
Uprior ② prior
K ((Fobo)d
(Hpost)d = 0 post
+
(18)
on
17

---

Published as conference paper ICLR 2026
Post is the posterior variance (constant across dimensions) given by
-1
K
(19)
post
prior
d lik
and (Nobb)d is the mean over the Dth dimension of the K new data points ff 1K Ik=1-
For the NPEs g(0 I we use neural spline flow (Durkan et al., 2019) with 5 coupling layers of 128 units each utilizing ReLU activation functions, L2 weight regularization with factor 7 10 5% dropout and a multivariate unit Gaussian latent space. The network is trained using the Adam optimizer for 100 epochs with batch size of 32 and learning rate of 10-1 These settings were the same both the standard simulation-based loss and our proposed superviee loss. For the conditions involving estimated likelihood g(x) f). we use the same configuration for the likelihood network as for the posterior network. For the summary network (xx) (if included),
deep set architecture (Zaheer et al., 2017) with 30 summary dimensions and mean pooling, equivariant layers each consisting of dense layers with 64 units and ReLU activation function. The inner and outer pooling functions also use 2 dense layers with the same configuration. The likelihood network as well as the summary network are jointly trained with the inference network using the Adam optimizer for 100 epochs with a batch size of 32 and learning rate of 10-4
18

---

Published as conference paper ICLR 2026
D COMPREHENSIVE RESULTS FOR THE MULTIVARIATE NORMAL CASE STUDY
Parameter dimensionality
posterior mean blas
posterior SD blas
posterior distance
15
Oet. Otpe
0.0
-0.2
10
Mest
-0.6
20
location
location
data
Parameter dimensionality
10
-4-
100
(no SCI
10 (no SC)
---
100 (no SC)
Size of the unlabeled training data
posteriar mean
posterior SD bias
posterior distance
0.2
Mouue
Oone
MMD
-0.4
5-
-0.6
observed data location
observed data location
observed data location
dataset size
+
0 (no SC)
- 32
Mean of the unlabeled training data
posterior mean bias
posterior SD bias
posterior distance
0.50
15
Met. Mous
0.25
-0.25
-0.5
10
10
0
8
8
10
observed data location
abserved data Iocation
observed data location
dataset
Figure 6: Bias of posterior mean, bias of posterior standard deviation and posterior distance quantified by maximum mean discrepancy the analytic posterior for variations of the default configuration outlined in Section 4. NPE approximators with the added self ceff-consistency loss component are shown in blue. NPE approximators using just the standard dimulation-based loss are shown in red. Irrespective of the varied factor and for all metrics, adding the self-consistency loss component is always a drastic improvement over the standard simulation-basee based loss alone. The plots show that adding the self-consistency loss component provides strong robustness gains even in high-dimensional spaces (top row) or when the self-consistent loss is evaluated on little data (center row). Variation of the mean of the unlabeled training data show that adding the self-consistency loss drastically improves posterior estimation long as data used for evaluating the celf-sonsitt loss is at least slightly out-of- disstribbution compared to the original training data (μ² 2 1). Errorbars show +1 standard deviations over 10 model refits training data.
19

---

Published as conference paper ICLR 2026
Likelihood function
posterior mean bias
posterior SD bias
posterior distance
15
Omee
1.0
MMD
10
0.5
0.0
0.5
PAF
10
10
10
observed data location
observed data location
observed data location
estimated
analytic
likelihood
estimated (no SC)
analytic (no SC)
Summary network
posterior mean bias
posterior SD bias
posterior distance
0.25
15-
Must Mooue
MMD
10
-0.25
-0.00
observed data location
observed data location
abserved data location
with with SC)
without
summary network
SC)
Figure 7: Bias of posterior mean, bias of posterior standard deviation and posterior distance quantified by maximum mean discrepancy to the analytic posterior when the likelihood is estimated (top row) and in presence of summary network (K 10 data points; bottom row). In the setting where the likelihood function is estimated, we observe lower bias of the posterior mean and lower maximum mean discrepancy the true posterior when the self-consistency loss component added compared the standard simulation- based loss alone. However, we do see some bias of the posterior standard deviation, although with reversed signed compared to the standard loss. The self-consistency loss provides strong robustness gains in the presence of summary network (and known likelihood) in terms of all metrics, Errorbars show 11 standard deviations over 10 model refits on new training data.
D.1 COMPARISON TO OTHER ROBUSTNESS METHODS
Existing approaches for enhancing the robustness of amortized Bayesian inference typically fall into one of three categories, each with its own limitations. First, regularization-base methods (e.g., Gloeckler et al., 2023; Huang et al., 2023) introduce additional constraints inevitably altering the original training objective, thereby deviating from the true Bayesian posterior that the self-consisenncy loss seeks to approximate.
Second, some methods employ supervised corrections using small amount of real-world, labeled data (e.g.. Wehenkel et al., 2024) which requires access to ground-truth parameters, shifting the setting from semi-supervised to fully supervised learning.
Third, unsupervised domain adaptation (UDA) techniques such as NPE-DANN and NPE-MMD have been proposed to improve robustness under domain shift by aligning the distributions of simulated and observed data in the summary feature space (Elsemüller et al., 2025). These losses govern the alignment the summary space between simulated and observed data, effectively adjusting the observed data resulting in posteriors
20

---

Published as a conference paper ICLR 2026
posterior mean blas
posterior SD bias
MMD
1.2
14
10
12
Hest - -Abuee
0.8
0.6
-6
0.4
NNPE
NPE-MMD
0.2
NPE-DANN
NFE+SC
2 4 observed data location
10
0 2
8
10
10
observed data location
observed data location
Figure 8: Bias of posterior mean, bias of standard deviation, and posterior distance quantified by maximum mean discrepancy, relative to the analytic posterior. NPE SC provides far more robustmess gains compared to NPE-MMD, NPE-DANN, and NNPE for the case of prior misspecification. Errorbars show =1 standard deviations over 10 model refits on new training data.
conditioned on transformed version of the observation (Target 2, Section 3, Elsemüller et al., 2025), rather than the analytic posterior. Moreover, as shown by Elsemüller et al. (2025), these UDA methods and NNPE (Ward et al., 2022) exhibit notably poor performance under prior misspecification.
That said, we conducted additional experiments with NPE-MMD (Huang et al., 2023), NPE-DANN (Elsemüller et al., 2025), and NNPE (Ward et al., 2022) on the two-dimensional multivariate normal model, replicating the prior misspecification setting described in Section 4.1. Since both UDA methods rely on a summary network for feature alignment, we employed a shared DeepSet architecture (2-layer network with 30- dimensional output) across all methods, including NPE+SC, to ensure a fair comparison. For self-consistenge loss, we used 32 unlabeled observed data Hlos = during training. We used Adam optimizer with a learning rate of 0.0005 trained for 40 epochs using 1,028 simulations. All the hyperparameter configurations were kept the same as described in Section 4.1 and Appendix C. For NPE-MMD, NPE-DANN, and NNPE, we utilized the implementation by Elsemüllerr et al. (2025).
The results clearly show that NPE-MMD, NPE-DANN, and NNPE fail under prior misspecification. Also, with increase in misspecification, the results for these methods became even worse while NPE+SC maintains the robustness gains. UDA methods also require large unlabeled data during training while we see strong robustness gains in NPE+SC with as few 4 unlabeled data points. As these methods have different training objective (Target 2, Section 3, Elsemüller et al., 2025) rather than the analytic posterior, we do not include them in other case study comparisons. Instead, we focus standard simulation-based NPE and state-of-the-art MCMC via Stan, which target the analytic posterior and serve as appropriate baseline and benchmark, respectively, for evaluating our semi-supervised NPE+SC approach.
DETAILED SETUP OF THE AIR TRAFFIC CASE STUDY
The parameters a, country-level intercepts, 8, the autoregressive coefficients, Y are the regression coefficients of houschold debt and ô, are the regression coefficients of GDP per capita, and o, is the standard deviation of the noise term. This model was previously used within ABI in Habermann et al. (2024). As commonly done for autoregressive models, we regress on time period differences to mitigate non-stationarity. This is critical for simulation-based inference because when B3 1, exponential growth quickly produces
21

---

Published as conference paper ICLR 2026
unrealistic air traffic volumes. Moreover, amortizing over covariate spaces, such as varying GDP per capita between countries, can lead to model misspecification such fluctuations are underrepresented in training.
For the air traffic model defined in Section 4.2, we set independent priors on the parameters as follows:
~ Normal 0.5)
Bj ~ Normal (0, 0.2)
7 Normal(O, 0.5) 8, Normal(O, 0.5)
(20)
lo)( ~ (11,all-1,,0)))
For the NPEs g( x)), we use neural spline flow (Durkan al., 2019) with 6 coupling layers of 128 units each utilizing exponential linear unit activation functions, L2 weight regularization with factor y = 10-3 5% dropout and multivariate unit Gaussian latent space. These settings were the same for both the standard simulation-bail based loss and our proposed semi-supervised loss. The simulation budget was set to N 1024. For the summary network, we use a long short-term memory layer with 64 output dimensions followed by two dense layers with output dimensions of 256 and 64. The inference and summary networks jointly trained using the Adam optimizer for 100 epochs with batch size of 32 and a learning rate of 5 x 10-4
COMPREHENSIVE RESULTS FOR THE AIR TRAFFIC CASE STUDY
Table I: Posterior metrics NPE and NPE augmented with self-consistency loss (NPE+SC) relative to Stan. For each parameter, absolute bias in posterior means and standard deviations, as well as Wasserstein distance, are reported. Values are shown as mean (SE) across 15 countries. Self-consistency loss was evaluated using M {4, 8, 15) countries during training.
Parameter
NPE
NPE+SC (M 4)
NPE+SC (M 8)
NPE+SC (M 15)
Mean bias (0x
MStanl )
a
0.079 0.001
0.003 0.009
0.014 0.013
0.002 + 0.001
0.153+0.00
0.012 + 0.026
0.031 + 0.023
0.001 + 0.002
y
0.087 0.044
0.048 0.037
0.006 0.026
0.002 + 0.003
0.053 0.033
0.046 + 0.033
0.042 + 0.023
0.003 + 0.004
log(a)
0.207 0.076
0.148 0.058
0.002 + 0.002
Standard deviation bias (oo Ostan
a
0.033 0.011
0.020 0.008
0.020 0.007
0.002 + 0.001
B
0.055 0.010
0.011 + 0.007
0.004 + 0.003
0.001 + 0.002
0.058 0.015
0.022 0.016
0.035 0.013
0.005 + 0.003
]
0.038 0.015
0.018 + 0.016
0.031 + 0.012
0.005 + 0.002
log(ar
1.09
0.058 0.017
0.011 + 0.008
0.004 + 0.002
Wasserstein distance
a
0.086 0.019
0.033 0.009
0.035 0.013
0.006 + 0.001
8
0.161 + 0.024
0.070 + 0.026
0.054 + 0.023
0.009 + 0.002
0.154 0.045
0.112 0.037
0.068 0.026
0.014 + 0.003
+ 1.119
0.102 + 0.033
0.064 + 0.023
0.013 + 0.002
log(o)
0.304 0.076
0.282 0.076
0.170 + 0.058
0.011 + 0.002
22

---

Published as conference paper ICLR 2026
1.5
NPE only
1.0
NPE SC
reference
0.5
0.0
BE
SE
CZ
DE
DK
EL
AT
FI
PL
IT
NL
IE
PT
ES
FR
0.5
0.0
-0.5
PL
CZ
NL
FI
PT
AT
FR
BE
SE
DE
DK
IE
IT
EL ES
0
-1-
IE
DK
NL
IT
BE
PT
EL
ES
PL
AT
DE
SE
CZ
FR
FI
0
-1
DK
IE
IT
EL
DE
SE
BE
AT
PT
FI
CZ
NL
PL
ES
FR
0
l0g
-1
-2
CZ PL AT
EL PT SE DK ES BE FR IE NL DE
Country code
Figure 9: Comparison of posterior estimates between standard amortized NPE (red circles), NPE augmented by our self-consistency loss (NPE+ SC: blue squares) and Stan (reference; gray triangles). The plots illustrate central 50% (thick lines) and 95% (thin lines) credible intervals of all five parameters for different countries, sorted by the lower 5% quantile according to Stan. Abbreviations follow the ISO 3166 alpha-2 codes. The self-consistency loss was evaluated on data from M 4 countries during training.
23

---

Published as conference paper ICLR 2026
1.5
NPE only
1.0
NPE SC
reference
0.5
0.0
BE
SE
CZ
DE
DK
EL
AT
FI
PL
IT
NL
IE
PT
ES
FR
0.5
0.0
-0.5
PL
CZ
NL
FI
PT
AT
FR
BE
SE
DE
DK
IE
IT
EL
ES
1
0
IE
DK
NL
IT
BE
PT
EL
ES
PL
AT
DE
SE
CZ
FR
FI
-1
DK
IE
IT
EL
DE
SE
BE
AT
PT
FI
CZ
NL
PL
ES
FR
0
109
-1
-2
CZ PL AT
EL PT SE DK ES BE FR IE
NL
DE
Country code
Figure 10: Comparison of posterior estimates between standard amortized NPE (red circles), NPE augmented by our self-consistency loss (NPE + SC; blue squares), Stan (reference; gray triangles). The plots illustrate central 50% (thick lines) and 95% (thin lines) credible intervals of all five parameters for different countries, sorted by the lower 5% quantile according to Stan. Abbreviations follow the ISO 3166 alpha-2 codes. The self-consistency loss was evaluated on data from M = 8 countries during training.
24

---

Published as conference paper ICLR 2026
1.5
NPE only
1.0
NPE SC
reference
0.5
0.0
BE
SE
CZ
DE
DK
EL
AT
FI PL
IT NL
IE
PT
ES FR
0.5
0.0
-0.5
PL
CZ
NL
PT
AT
FR BE SE
DE
DK
IE
IT
ES
0
-1
DK
NL
BE
PT
ES
PL
AT
DE
SE
CZ
FR FI
1
0
-1
DK
IE
IT
DE
SE
BE
AT
PT
FI
CZ
NL
PL
ES
FR
logo o
-1
-2-
CZ PL AT
EL PT SE DK ES BE FR IE
NL
DE
Country code
Figure 11: Comparison of posterior estimates between standard amortized NPE (red circles), NPE augmented by our self-consistency loss (NPE + SC; blue squares), Stan (reference; gray triangles). The plots illustrate central 50% (thick lines) and 95% (thin lines) credible intervals of all five parameters for different countries, sorted by the lower 5% quantile according to Stan. Abbreviations follow the ISO 3166 alpha-2 codes. The self-consistency loss was evaluated on data from M 15 countries during training.
25

---

Published as conference paper ICLR 2026
F. COMPARISON TO VAE AND ONLINE TRAINING OF NPE
We conducted additional experiments using variational auto-encoder (VAE) (Kingma & Welling, 2013) baseline with a known decoder (likelihood) and a neural encoder (posterior network). The encoder a feed- forward network with three hidden layers of 128 units each, followed by two output heads that parameterize the mean and standard deviation of a standard normal base distribution. The model was trained on 40,960 simulated pairs until convergence. We trained NPE and NPE SC in online setting that leads to training data the same size as VAE. We also compare our method against Schmitt et al. (2024) where simulated data is used evaluate SC loss during training. We similarly used 40,960 simulations for training NPE and 32 simulations for evaluating self-consistency loss for 40 epochs.
Results comparing the above methods for the air-traffic experiment reported in Table and Figure 12. These results show that the performance of VAE with known likelihood is comparable to NPE (which does use explicit likelihood density) but it is much worse than our semi-supervised (NPE + SC) approach. Our results also show that the method of Schmitt et al. (2024) does not overcome model misspecification and performs poorly compared our cun--super approach.
Table Posterior Wasserstein distances, mean bias, and standard deviation bias for online trained NPE, Schmitt et al. (2024), VAE, and our proposed NPE SC method relative to Stan. The self-consistency loss was evaluated on data from M countries during training. Metrics are averaged over all 15 countries.
Parameter
Mean bias (1s MStan
NPE (online)
Schmitt et al. (2024)
VAE
NPE (online) SC
0.075 + 0.028
0.026 + 0.019
0.060 + 0.017
0.008 + 0.002
0.087 0.020
0.052 0.028
0.096 0.020
0.013 + 0.003
0.101 + 0.014
0.200 + 0.071
0.108 + 0.020
0.014 + 0.003
0.096 0.018
0.112 = 0.072
0.110 0.028
0.009 + 0.002
log(o)
0.189 + 0.051
0.042 + 0.141
0.142 + 0.024
0.009 + 0.002
Standard deviation bias (fo
Ostau I)
a
0.030 + 0.008
0.003 + 0.021
0.029 + 0.006
0 0.001
0.011 + 0.002
0.039 + 0.039
0.050 0.004
0.058 + 0.016
0.022 + 0.033
0.079 + 0.014
0.009
0.059 0.017
0.001 = 0.041
0.083 0.015
0.016
log(a)
0.170 + 0.009
0.066 + 0.034
0.026 + 0.006
0.013 = 0.002
Wasserstein distance
Q
0.080 0.028
0.068 + 0.017
0.010 0.002
B
0.088 + 0.020
0.133 + 0.024
0.102
1018
0.113 0.015
0.278 + 0.050
1018:0000
0.110 0.020
0.240 + 0.049
0.128 0.027
0.017 1 0.002
log(a)
0.227 + 0.045
0.459 + 0.073
0.144 0.024
0.016 + 0.002
26

---

Published as conference paper ICLR 2026
2
NPE (online)
Schmitt
VAB
NPE SC
Stan
0
CZ
BE
SE
FI
EL
AT
DK
PL
IT
PT
IE
NL
DE
ES FR
1
0
-1-
NL
CZ
PL
FR
PT
BE
AT
DE
SE
DK
IT IE EL
2
1
0
-1
IE
PT
DK
PL
BE
ES
IT
CZ
AT
SE
FI
DE
1
0
EL
DK
IT
IE
SE
AT
FI
DE
BE
CZ
PT
PL
NL ES FR
0
-2
FI
CZ
PL
AT
EL
PT SE DK ES BE FR NL IT DE
Country code
Figure 12: Comparison of posterior estimates among online trained NPE, Schmitt al. (2024), VAE, and proposed cumi-supervised NPE+SC method against Stan (reference). The self-consistency loss was evaluated data from M 15 countries during training. The results clearly show that our method is consistent with Stan and outperforms all other methods by large margin.
27

---

Published as conference paper ICLR 2026
G DETAILED SETUP OF THE NEURON ACTIVATION CASE STUDY
G.1 MODEL DESCRIPTION
The prior and likelihood are given by
INa Log 0.12). 9x LogNormal(log(36)..11 0.13), 9M LogN)omal(oog(22).00555 0.52)
ENA N50, 52), Ex Normal(-77,5*) (7),,²², Bleak Normal(-55,55 55,5", C. Normal 0.002
Vi. Student-t (t), 0.1², df=10),
where INa: 9k. 9M denote the maximum conductances (in mS/cm²) of sodium and two different types of potassium channels; ENa; Pk, Freak the sodium, potassium, and leak reversal potentials (in mV) for sodium, potassium, and leak currents; Cm is the membrane capacitance (in jF/cm²). The membrane voltage Vn (1) time is obtained by solving set of five ordinary differential equations. Here, X LogNormal o²) denotes lor-normal distribution with (X) Normal(H,, o²)
The membrane voltage Vm (t) evolves according to the classical Hodgkin-Huxley model, extended with an additional slow (muscarinic, M-type) potassium channel with current IM. The total current across the membrane modeled as:
dV INa + Ix IM
Cm dt
Ieak + [mn(),
(21)
where INa: IK. IM. Ieak are the sodium, potassium, M-type potassium, and leak currents respectively; and (t) is an externally applied current, which set to be pulse input of 3.248 nA between fon 10ms and forr = 50ms. The ionic currents are calculated as:
NN INam 3 hi ENa - Vm)
(22)
1K gKn"t (Ex Vm)
(23)
JM 9MP P.x Vm)
(24)
heak Dlenk (Eeak Vm),
(25)
where the leak conductance is fixed to fleak variables take the form:
0. 1mSmmm?m and m, h, n, p are gating variables. The gating
dx
Jxx (Vo)
E mn,m,
(26)
dt To (Vm)
with (Vm) (Vmn)))(a, (v(m)) and Ta (Vm)) where a₂ and 32 are voltage-dependen rate functions defined as:
1/(a + Bo(V for n, m and h.
(V())) 0.032
(0.2(.. 75))
0.28 exp(O..((v.. 100))
Bn (n
0.2
0.2
am (m())) 0.32
expppxx2p5(()) 0V. 73))
0.28 exp(0.2((V. 100))
Bmn(V
0.25
0.2
(V
77)
4
(V())) 0.128 exp(
(L(())) 1 + exp(-0.2 100))
18
28

---

Published as conference paper ICLR 2026
For the gating variable p the M-type potassium channel, sigmoidal steady-state activation and custom time constant are used:
600
Poo (Vm) 1+ exp((O.1(Vmm
0.1(V + 35)) (Vm) 3.3 exp(0.05(Vmn + 35)) exp(- + 35))
Numerical integration of the system is performed using fixed-step Euler method over time window [0. 60] with time step At 0.01. Voltage traces Vn (t) are downsampled every 30th observation, and 200-dimensional time series are simulated according to Vi, Student-i(V 0.1 df=10).
G.2 NETWORK ARCHITECTURE AND TRAINING
For the NPEs g( ( we use a neural spline flow (Durkan et al., 2019) with 10 coupling layers of 256 units each utilizing ReLU activation functions, L2 weight regularization with factor 7 10 1.-3, 5% dropout and multivariate unit Gaussian latent space. These settings were the same for both the standard simulation-baser loss and our proposed semi-supervised loss. For the summary network, layer with 100 output dimensions followed by sequence of dense layers with output dimensions of 400, 200, 100, and 50, respectively. The inference and summary network are jointly trained using the Adam optimizer with batch size of 256 for 100 epochs and a fixed learning rate of 5 x 10-1, followed by 100 epochs with a fixed learning rate of 10-6 and a final run of 100 epochs with learning rate of 5 x 10-6
long short-term memory
COMPREHENSIVE RESULTS OF THE NEURON ACTIVATION CASE STUDY
NPE only
membrane potential [mV]
100
0
80
50
count
SC
50
40
20
-50
10
20
30
60
10
15
time [ms]
MAB(NPE) MAB(NPE+SC)
(a): Posterior predictive samples without (top row) and with (bottom row) self-consistency loss.
(b): Quantitative evaluation of predictive bias.
Figure 13: (a) Posterior predictive samples (gray) inferred from an in-simulation dataset (black) with parameters 0 ~ Normal(O, 1). Both NPE only and NPE+SC produce predictions that are consistent with the observed data. However, samples from NPE+SC are much closer to the ground truth. (b) Histogram of the mean absolute bias (MAB) difference of posterior predictions computed for 1000 out-of-simulation datasets. NPE+SC has lower bias than NPE for almost all datasets. Mean absolute bias is defined as Ei= |ir. Ds.el for time series Vi, with observation index at time 1,..., T.
MAB (yr,t; (i.s)
Ut.1 = - ((())) denotes the mean of the posterior predictive distribution at time computed over posterior samples.
29

---

Published as conference paper ICLR 2026
In-distribution data (e Normal(O, 1))
NPE only
NPE+SC
50
50
membrane potential [mV]
-50
50
10
20
30
40
50
60
10
20
30
40
50
60
50
50
0
-50
10
20
30
40
50
60
50
50
-50
-50
10
20
30
40
50
60
50
50
0
-50
-50
10
20
30
40
50
60
10
20
30
40
50
60
50
-50
10
20
40
50
60
time [ms]
Figure 14: Posterior predictive samples (gray) inferred from 5 simulation datasets following the same distribution as the training data. Both NPE only and NPE+SC show predictions that are consistent with the observed data.
30

---

Published as conference paper ICLR 2026
Out-of-distribution data (0 Normal(-2, 1))
NPE only
NPE+SC
50
membrane potential [mV]
-50
-50
30
40
50
0
-50
10
20
30
40
50
60
10
20
30
40
50
60
50
-50
-50
50
-50
-50
10
20
30
40
50
60
10
20
30
40
50
60
50
0
-50
-50
10
20
40
50
60
time [ms]
Figure 15: Posterior predictive samples (gray) inferred from 5 out-of-simulatte datasets, generated from parameter draws 0 ~ Normal(-2, 1). While NPE only produces highly biased predictions, NPE+SC is consistent with the observed data.
31

---

Published as a conference paper ICLR 2026
I DETAILED SETUP OF THE MNIST IMAGE DENOISING CASE STUDY
We implement the jointly amortized posterior and likelihood networks (Radev et al., 2023) using two normalizing flows with fully connected cffine-cuullng layers that operate on the flattened 784-pixel vectors. We use the same network architecture that was used by Radev et al. (2023). As both 0 and whose intrinsic dimensionality is significantly lower than their raw pixel count, we use identical 4-laeer convolutional neural networks summary networks for both posterior and likelihood networks. These summary networks terminate in a global wverage-pooling layer to produce a 128-dimensional summary of the original or blurred image, respectively. The posterior network itself is implemented neural network (cINN) consisting of 12 conditional affine-coupling layers; each coupling layer embeds its conditional information via an internal fully connected network with single hidden layer of 512 units and ReLU activations. The likelihood network adopts exactly the same ienditiionalcoupliin architecture. In both cases, we employ a multivariate Student-T distribution in the latent space (Alexanderson & Henter, 2020; Radev et al., 2023), which enables more stable maximuu-llkellhowe training at elevated learning rates.
images
conditional invertible
The prior network that was used to generate blurred MNIST simulations followed the same architecture as the posterior network defined above. We applied the Gaussian blur with PSF = 1.0 to the 5, 923 images of digit "0" in the MNIST training dataset to train the prior network using Adam optimizer for 120 epochs with a batch size of 32, learning rate of 1 x 10-3, and a 15% dropout. After training, we generated 12000 blurred images (#) of the digit 0 to train the posterior and likelihood networks. A Gaussian blur with PSF = 1.0 was further applied to these images to generate observations (z) which represent images from noisy camera. The posterior and likelihood networks along with summary networks were jointly trained on {0' , [=1)=1 pairs for 100 epochs with batch size of 32 using learning rate of 1 x 10-4, and 15% dropout.
12000
For self-consistency loss, the MNIST test set of digit "0" was divided into two subsets comprising 400 and 580 images respectively. The subset with 400 images was used for training self-consistency loss. A Gaussian blur with PSF = 1.0 was applied to these images to generate observations (x*). No prior blur was applied to these images. This represents a prior misspecification scenario the simulated images used to train NPLE were already blurred before applying the noisy camera while the MNIST images used for inference do not have a prior blur. This misspecification scenario depicts the effectiveness of utilising self-consistency loss overcome prior misspecification. The self-consistency loss was activated at epoch 21, with its weight linearly ramped from zero to one by epoch 40. Training was performed using minibatches of 16 images, and 32 consistency samples were drawn to estimate the variance. The inference was performed on the other held-out subset comprising 580 images and all the results in Figures 5, 16 and 17 use the MNIST images from this subset.
32

---

Published as conference paper ICLR 2026
J COMPREHENSIVE RESULTS OF THE MNIST IMAGE DENOISING CASE STUDY
Parameter (8)
Observation DO
Moan INPLE SC)
Std Day INPLE SC)
Mean (NPLE)
Std Dev (NPLE)
Parameter (8)
Observation
Mean INPLE sc!
Std Dev INPLE Sc)
Mean (NPLE)
Std Dev (NPLE)
Figure 16: Examples of denoising results for MNIST images of digit "0" in the held-out test set. The first row shows ten randomly selected MNIST images (0), the second row depicts the same images after applying the Gaussian blur (z), third and fourth rows depict the mean and standard deviation of 500 posterior samples estimated from the corresponding blurry observations using NPLE + SC based model, and the fifth and sixth rows depict the mean and standard deviation of 500 posterior samples from model based on NPLE only. Incorporating self-consistency loss significantly improves denoising as the means of reconstructed unblurred image are smoother, less-pixelated and better resemble the ground truth. The darker regions in the standard deviation show the regions of higher variability in the outputs. The standard deviation maps of NPLE + SC based upproach are far more coherent showing high variability only along the edges.
33

---

Published as conference paper ICLR 2026
Parameter (8)
Posterior samples (NPLE)
Parameter (0)
samples
Figure 17: Ground-truth images and the corresponding posterior draws for seven randomly selected MNIST "0" digits from the held-out MNIST test set. In each the leftmost panel shows the true image (6), and the following panels show ten independent samples from the approximate posterior. The top-figure shows the posterior draws using the standard NPLE based model and the bottom figure shows posterior draws from combining self-consistency loss to the NPLE based model. It can clearly be seen that NPLE+SC posterior draws are better reconstruction the original image whereas NPLE based posterior samples highly pixelated.
34

---

Published as conference paper ICLR 2026
K COMPUTATIONAL RESOURCES FOR EXPERIMENTS
1. Multivariate normal model:
The experiments were run on 16-core AMD Ryzen 5950x CPU, equipped with GB system RAM.
2. Air traffic case study: The experiments were conducted on single MacBook Pro (M3, 2024) equipped with Apple's 's M3 chip and 16 GB of unified RAM, rumning macOS Sonoma 14.6. We did not use the GPU cores.
3. Neuron activation case study: The experiments equipped with 32 GB of system RAM.
16-core AMD Ryzen 5950x CPU,
4. MNIST image denoising:
The experiments were run on a -ereh--effrrmncc compute cluster using a GPU-equipped compute node featuring single NVIDIA Tesla P100-PCIE with 12 GB of dedicated HBM2 memory, paired with 16 GB of system RAM. The training for the longest experiment took 100 minutes.
35


Reviews:
[{"review_id": "fV7WZirjpl", "reviewer": "Reviewer_bLfT", "summary": "This paper introduces self-consistency loss proposed to improve training efficiency into amortized Bayesian inference (ABI) to enhance its robustness. In addition to the loss based on labeled simulation data, it includes the self-consistency loss defined with unlabeled real data. The later adds extra information to improve the robustness of ABI.", "questions": "1. Does the condition of proposition 2 always hold? Does proposition 3 depend on on proposition 2? If so, the statement of proposition 3 reads like it holds regardless. Some clarification will be appreciated.\n\n2. Figure 4 -(b) is misleading. I understand that it represents the difference between MAB (NPE) and MAB (NPE+SC). However it looks that MAB (NPE) is 0 while MAB (NPE+SC) is bigger and concentrated around 5. Why not plotting the histograms of each using different colors?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The independence of true parameter values is the most striking feature of the proposed method. The method is justified by theories. Various numerical evidences including through simulation and real-data applications are strong and convincing.", "weaknesses": "The conditions of propositions 2 and 3 are not clear to me (See question below). Figure 4 (b) is a little bit misleading. Overall, these are minor defects.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yWfI3oVL2Y", "reviewer": "Reviewer_Lkfz", "summary": "The authors propose to include a self-consistency loss regularizer to the typical neural posterior estimation (NPE) class of objectives for amortized inference (more generally, a loss with respect to any score). The regularizer, unlike the main objective, is computed over real data, not simulated draws from the model. The regularizer enforces self-consistency (Eq. 1) on the observed dataset, draws from which are assumed to be the target for inference. The motivation provided by the authors is that this results in more accurate inference, especially in distribution-shift cases where the real observations are substantively different from the simulated data with which $q$ is fit.", "questions": "- What precisely is meant by “sufficient statistics” in line 158? Sufficient in the formal statistical sense (if so, can you elaborate?), or do you mean more informally low-dimensional learned representations of $x$ that are known to have low reconstruction error, or something like that?\n- The experiments appear to utilize a fixed simulation budget (e.g., 1024, line 723). Although this is the setting explored in Schmitt et al., my understanding of “standard” NPE typically allows unlimited samples from the model, e.g., fresh samples are drawn for each minibatch used to compute each gradient. I think the authors should also vary the simulation budget in their case study, as not doing so unfairly penalizes NPE (besides have the wrong model, it also has only finitely many samples from this model, and so will overfit).", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The writing and proposed methodology are clear and easy to understand\n- The proposed regularizer is intuitive: as the true posterior demonstrates  self-consistency property, it’s a natural extension for the variational posterior to satisfy this condition (approximately)\n- The problem setting is a significant one, as distribution shift in simulation-based inference (SBI) problems continues to be a robust area of research.\n- The others take care to formalize their result more carefully with respect to strictly proper losses, proving that the regularizer does not alter the problem's solution.", "weaknesses": "- The method is combinatorial: the NPE (score-based) objective for SBI is well-established, and Schmitt et al. (2024) introduced the self-consistency loss.\n- I don’t find the degree of novelty to be adequate enough to differentiate the work from Schmitt et al. In the case where the simulation model is correct, the proposed method is exactly identical to Schmitt et al. Thus, the authors’ main contribution, to me, seems to be applying the method to the setting of a misspecified simulator. While this is undoubtedly an important setting as acknowledged above, the novelty of applying the method to a new setting is fairly limited\n- Experimental evaluation is insufficient. To me, comparing only to NPE isn’t a fair comparison. The experiments involve deliberately mispecifying the model, so no one would expect NPE to perform well. At a minimum, the authors should compare to a VAE trained with the incorrect prior/likelihood (but with the real dataset), as well as the method of Schmitt et al., which imposes self-consistency on the simulated data. Beyond these, comparisons to some of the many domain adaptation approaches to SBI from the literature would be in order.\n- I find the literature review of relevant related work to be lacking. For example, the Autoencoding Variational Autoencoder (Cemgil et al.) has a similar notion of self-consistency applied to real data, but in the VAE setting. This line of work would be a promising area for future ablation studies that I think would improve the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tdWOVokjmN", "reviewer": "Reviewer_L5ck", "summary": "This paper proposes a semi-supervised approach to amortized Bayesian inference using self-consistency losses on unlabeled data, demonstrating improved robustness in out-of-distribution settings. While the theoretical contributions are solid, key claims about high-dimensionality and robustness are not adequately supported by the experimental evidence.", "questions": "See above.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Improving robustness of SBI to data distribution shifts in an important problem.", "weaknesses": "- The abstract uses \"bad pre-asymptotic behavior\" without definition, making the core problem inaccessible.\n\n- When would practitioners have unlabeled real data but not be able to generate more simulations?\n\n- The paper claims \"high-dimensional\" capability but experiments max out at 100 parameters with significant performance degradation (Figure 2a shows MMD increasing substantially). The MNIST example (784D) is modest by modern standards. Please scale back claims.\n\n- Head-to-head comparisons on all case studies are essential. Also missing: ablations on the choice of p_C(theta), sensitivity to lambda schedule, and performance with different amounts of unlabeled data.\n\n- The paper claims to maintain ABI's speed but provides zero timing comparisons, memory usage, or analysis of overhead from computing SC loss and sampling q_t(theta|x) during training.\n\n- Propositions 1-4 establish strict properness under idealized conditions (universal approximators, infinite data), but there's no finite-sample analysis or characterization of approximation error with realistic network capacities. The gap between asymptotic guarantees and finite-sample behavior needs some comments.\n\n- The MNIST prior misspecification (double-blur vs single-blur) is artificial. How does the method perform under realistic misspecification where the error structure is unknown during training?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "OCIqDhN3Wk", "reviewer": "Reviewer_ipsv", "summary": "The authors propose a new loss term that uses unlabeled data to improve the robustness of amortized neural posterior estimation methods. They thoroughly motivate and explain their loss term and prove that it is a proper loss (in the relevant context of SBI) and present promising results for their method on a simple Gaussian mixture scenario, as well as three more challenging setups (Air traffic data, Hodgkin-Huxley model, MNIST denoising).", "questions": "I would be happy to see your response regarding my concern about the contribution of the paper. \n\nWhy should I use SBI, and this method in particular, when I have access to the likelihood under the simulator? Why not directly leverage the likelihood and rely on samples instead?\n\nAre there any further (theoretical) insights into the case where both the likelihood and posterior are unknown? \n\nI am looking forward to your response!\n\nI would be very happy to increase my score if you further justify the applicability of using NPE plus your robustness term when the likelihood is known (i.e. in what case would this be a relevant setup?) OR provide more evidence that using your method with an estimated likelihood is theoretically sound and works well experimentally.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The core idea of the paper is neat. I also appreciate that it's a simple idea!\n\nFurthermore, the approached problem of robustness is highly relevant in SBI (and beyond). It also seems like the availability of unlabelled real-world data, that is necessary for this approach, is a very reasonable assumption that is often fulfilled in practice. I see this as a major strength of this approach compared to other recent methods. \n\nThe overall presentation is good and the structure of the paper very clear and well-written. \n\nRelated work is comprehensively discussed. \n\nThe idea is original, but some parts of the theory are already present in other papers, albeit for a different purpose.", "weaknesses": "Contribution: \n\nI have one major concern with the paper: You do assume for the majority of your theory and experiments that the likelihood under the simulator is known. Doesn't this defeat the purpose of SBI to some degree? If the likelihood is know, why not use a method that directly leverages this information?\n\nI am aware that the case of unknown likelihood + unknown posterior is briefly discussed, but more theoretical insight would be very helpful. \n\nAnalogously for the experiments. I believe that the case of unknown likelihood + unknown posterior is the actually relevant one, but it is not really evaluated or thoroughly discussed. \n\nI am also unsure about the benefit of the method if used in combination with post-hoc methods such as SNIS; I believe that many people practically interested in the method might want to use post-hoc methods in combination with this approach.\n\nPresentation: \n\nPlease note that those two points are not crucial: \n\nSome parts of the methodology seem a bit repetitive and do not really present very interesting insights. For instance Proposition 4 seems a bit too redundant for this much space. I am wondering if it would help the paper to shorten this section and, for instance, include a brief introduction to amortized SBI in the beginning to help readers less familiar with this field. \n\nSimilarly for the description of the multivariate Normal case-study: I feel like it takes up a bit too much space for the practical relevance this case-study has. Instead more results/figures for the other case-studies could be presented in the main body. \n\nNit: There is a typo around line 392: \"enablung\"\n\nExperiments: \n\nPlease also report metrics for the denoising MNIST case-study. \n\nSeeing C2ST scores for all experiments would also be nice. \n\nWhile I believe that the experiments overall are similar in scope to what is usually done in SBI, but very small compared to other fields of ML. I would be interested in seeing larger-scale experiments (higher dimensionalities, larger neural networks, more extensive training); but this is not a critical point.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)