Amortized Bayesian Workflow
Authors: Chengkun LI
Aki Vehtari
Paul-Christian Bürkner
Stefan T. Radev
Luigi Acerbi
Marvin Schmitt
Venue: Accepted by TMLR ()

Abstract: Bayesian inference often faces a trade-off between computational speed and sampling accuracy. We propose an adaptive workflow that integrates rapid amortized inference with gold-standard MCMC techniques to achieve a favorable combination of both speed and accuracy when performing inference on many observed datasets. Our approach uses principled diagnostics to guide the choice of inference method for each dataset, moving along the Pareto front from fast amortized sampling via generative neural networks to slower but guaranteed-accurate MCMC when needed. By reusing computations across steps, our workflow synergizes amortized and MCMC-based inference. We demonstrate the effectiveness of this integrated approach on several synthetic and real-world problems with tens of thousands of datasets, showing efficiency gains while maintaining high posterior quality.

Published in Transactions on Machine Learning Research (02/2026)
Amortized Bayesian Workflow
Chengkun Li University of Helsinki
Aki Vehtari
aki. dbi.vchtari@oatto..
ELLIS Institute Finland, Aalto University
Paul-Christian Bürkner
paul www.
TU Dortmund University
Stefan T. Radev Rensselaer Polytechnic Institute
cudevs@rpil edu
Luigi Acerbi
Suigg..uerrbiShhcsinn.....
University of Helsinki
Marvin Schmitt
mail.marriinschhsstt@@@maaaa com
Independent Scientist
Reviewed OpenReview: https: // openreview. net/ pormf/ id= V7dJJKKD
Abstract
Bayesian inference often faces trade-off between computational speed and sampling accuracy. We propose adaplive workflow that integrates rapid amorlized inference with gold-standard MCMC techniques to achieve a favorable combination of both speed and accuracy when performing inference on many observed datasets. Our approach uses principled diagnostics to guide the choice inference method for each dataset, moving along the Pareto front from fast amortized sampling via generative neural networks to slower but guaranteed-accure when needed. By reusing computations across steps, our workflow synergizes amortized and MCMC-based inference. We demonstrate the effectiveness of this integrated approach several synthetic and real-world problems with tens of thousands of datasets, showing efficiency gains while maintaining high posterior quality.
MCMC
1 Introduction
In many statistical modeling applications, from finance to biology and neuroscience, often aim to infer unknown parameters 0 from observables modeled as joint distribution p(0. y) (e.g., Raulo et al., 2023; Seaton et al., 2023; George et al., 2022; Landmeyer et al., 2020; Chen et al., 2019; Malén et al., 2022; Schneider et al., 2018; Tsilifis Ghosh, 2022). The posterior p(0lv) y) is the statistically optimal solution to this inverse problem, and there are different computational approaches to approximate this target distribution.
Markov chain Monte Carlo (MCMC) methods constitute the most popular family of posterior sampling algorithms and still remain the gold standard for modern Bayesian inference due their theoretical guarantees and powerful diagnostics (Gelman et al., 2013; 2020). MCMC methods yield autocorrelated draws conditional a fixed datasct Yobs- As consequence, the probabilistic model has to be re-fit for cach new datasot, which involves repeating the entire MCMC procedure from scratch. Modern implementations equip MCMC with saate-of-the-ate extensions, for example, through Hamiltonian dynamics (HMC; Neal, 2011), by minimizing the required tuning by users (NUTS; Hoffman Gelman, 2014), or by parallelizing thousands of chains on GPU hardware (ChEES-FE HMC: Hoffman al., 2021). The rell-established Boyesian workflow (Gelman et al., 2020) leverages these tools in iterative process of model specification, fitting, evaluation, and revision.

---

Published in Transactions on Machine Learning Research (02/2026)
While powerful, this approach becomes computationally burdensome when applied independently to large collections of datasets.
Differently, amortized Bayesian inference (ABI) aims to learn a direct mapping from observables to the corresponding posterior p(0 y), using flexible function approximators such as deep neural networks (Cranmer et al., 2020; Radev et al., 2020; Greenberg al., 2019; Papamakarios et al., 2021; Wildberger et al., 2023; Sharrock et al., 2024; Zammit-Mangion et al., 2025). Amortized inference typically follows a two-stage approach: (i) training stage, where neural networks learn to distill information from the probabilistic model based simulated examples of observations and parameters (0. ~ p(0) p(y|0); and (ii) an inference stage where the neural networks approximate the posterior distribution for an unseen datasct Yobs in near-instant time without repeating the training stage. In other words: The upfront training cost is amortized by negligible inference cost on arbitrary amounts of unseen test data. Owing to its reliance on simulated data, amortized inference this form overlaps with simulation-based inference (Cranmer et al., 2020), which originated from postcrior computations for models with intractable likelihood.
However, amortized inference lacks the powerful diagnostics and gold-standard guarantees associated with MCMC samplers in the standard Bayesian workflow (Gelman et al., 2020). Yet, applying a standard workflow is computationally prohibitive at scale. In modern Bayesian computation, MCMC and ABI occupy different ends of a Pareto frontier (see Figure 1): the former provides reliable accuracy at high cost, while the latter offers near-instant inference speed with limited per-dataset reli- ability (Hermans et al., 2022; Schmitt et al., 2023; Lueckmann et al., 2021).
Accuracy guarantees weak
MCMC
strong
Pareto front
+ PSIS x
Amortized inference
instant Inference speed per data set
slow
In this paper, we propose an adaptive workflow that consistently yields high-quality posterior draws while remaining computa- tionally efficient. Our proposed workflow moves along the Parelo front, enabling aat--and-accurate inference when possible, and slow-but-guarantee autt-guarranted--ccuraaa accuratee inference when necessary (sce Figure 1). It combines the strengths of ABI and MCMC by incorporating diagnostic checks to guide inference decisions and reuse compulations wherever pos- sible. The resulting amortized Bayesian workflow therefore offers a principled, scalable, and diagnostic-driven approach for efficient posterior inference on many observed datasets; see Figure 2 for a conceptual overview." 1 To summarize, our contributions are:
Figure Our workflow adaptively moves along the Pareto front and reuses previous computations.
Design of- systematic guidance through-an an adaptive Bayesian workflow for accelerating Bayesian inference, which combines the strengths of amortized inference, importance sampling, and MCMC in a theoretically motivated and modular manner.
Empirical validation of the workflow and of its inference speedup, demonstrating the applicability of the workflow on both synthetic and large-scale, real-world problems.
2 Integrating amortized inference into the Bayesian workflow
Our adaptive workflow starts with neural network training to enable subsequent amortized inference on large number of unseen datasets typically well into tens of thousands. This training phase is conceptually identical to standalone amortized inference training (e.g., Radev et al., 2020; Cranmer et al., 2020). For the inference phase, however, we develop a principled control flow that guides the analysis. Based on stato-of-the-art diagnostics that are tailored to cach step along the workflow, we propose decision critcria to select the appropriate inference algorithm for each observed dataset. In order optimize the overall efficiency, our workflow contains mechanisms to reuse previous computations along the way.
The software implementation available https://github: http://ghhhubb.oom////////annaat Baw.:iiidd.aaaeeiianaooorrtrtree
2

---

Published in Transactions on Machine Learning Research (02/2026)
Step Amortized inference
Step 2
Step 3 ChRFS-HMC with amortized inits
Pareto-smeothed importance sampling
re-use draws
re-use
draws
Neural donsily eslimator
draws 0 mus
Importance weights
Pareto distriuttion
Initalize S superchains
Wammup and sampling
6,00
Example draws
Diagnostics OK?
Example
Diagnostics
Example draws
Dingnosties OK?
nc
NUTS revise model
OK?
yes
yes
yes
Accept
Accept
Accept drows
drawo
draws
K-256 observed dalasces
PSIS on remaining 64 datasets
CIEES-HIMC ca remaining 8 datasets
Accept amortized draws for 192/256 data
Accept PSIS draws for 56/64 data
Accepi CHÉES-HMC draws for 888 dita sets
Figure 2: Our adaptive workflow leverages near-instant amortized posterior sampling when possible and gradually resorts to slower-bur more accurate sampling algorithms. As indicated by the blue dashed arrows, we reuse the draws from the amortized posterior in Step 1 for the subsequent steps in the form of PSIS proposals (Step 2) and initial values in ChEES-HMC (Step 3).
2.1 Training phase:
simulation-bart based optimization
In ABI, a neural estimator da with trainable parameters & typically minimizes strictly proper scoring rule S (Gneiting & Raftery, 2007; Pacchiardi & Dutta, 2022) in expectation over the joint model p(0, y) p(0)p(y) 0),
arg argmin min ((,,)) -p()))) [S(1o(-19) 1y). 8)].
(1)
A popular choice is the logarithmic scoring rule, S(2a(- v): 0) = log go (0) y), which amounts to the forward Beulbackk-eibber (KL) objective used for training normalizing flows in AB1 (Greenberg et al,, 2019; Radev et al., 2020). Score- based formulations that target time-dependent gradient Vo, log v) are also possible (Sharrock et al., 2024; Gloeckler et al., 2024). Since most Bayesian models are generative by design, we can readily simulate M synthetic tuples of parameters and corresponding observations from the joint probabilistic model,
(@(m), (m)) p(0.y)
gim) p(0). (mm) ply| 0) for m
which results in the training set {(((m) g ) for optimizing Eq. 1. Throughout this paper, we use coupling-based normalizing flows (Durkan et al., 2019; Papamakarios et al., 2021) as a flexible conditional density estimator dø and the forward KL divergence as the training objective. However, our proposed workflow is agnostic to the specifie choice of generative backbone used for amortization, as long the model supports efficient sampling (see Section 2.2.1) and density evaluations (see Section 2.2.2).
Diagnostics. Since the neural network training algorithm hinges on simulated data, we cannot evaluate the amortized posterior estimator on real data just yet. However, we can easily simulate a synthetic test set ((e4), y of size JJ from the joint model via Eq. 2. In this closed-world setting, we know which "truc" parameter vector a(5) generated each simulated test datasct y(u). A key diagnostic for evaluating the amortized posterior est.maator is simulation--nsed calibration checking (SBC; Talts et al., 2018; Säilynoja et al., 2022; Modrák et al., 2025; Yao & Domke, 2023). Formally, SBC involves (1) defining test quantity f f x Y + R (e.g., marginal projections 0 or the log likelihood p(y 0)), (2) computing this statistic for the true data- daata-generating parameter 05), and (3) comparing it to the empirical distribution of the same statistic derived from amortized posterior draws given y(o) (Modrák et al., 2025). The rank of the true statistic within the posterior draws should be uniformly distributed if the amortized posterior estimator is well-calibrated.
We recommend assessing uniformity using the graphical approach by Säilynoja et al. (2022), which reveals the type of miscalibration present (e.g., bias over- wven/uunderdddspersson and is therefore useful for guiding
3

---

Published in Transactions on Machine Learning Research (02/2026)
improvements to amortized training. The choice of test quantity in SBC determines the sensitivity of the check; for example, the log-likelihood test quantity is typically more sensitive at detecting discrepancies than marginal projections (Modrák et al., 2025); using expressive neural classifiers is also possible (Yao & Domko, 2023). We further note a tradc-off: imposing stricter criteria can improve the fidelity of the amortized estimator but will also tend reject otherwise practically useful amortized estimators.
By default, we use marginal projections as the test quantities for SBC and complement SBC checking with parameter recovery checking, where parameter estimates compared against known ground-truth parameters via direct visualization (Radev et al., 2020; 2023). Parameter recovery checking provides practical insight into whether the learned inverse mapping from y to A is effective and helps mitigate known failure mode of SBC with marginal projections as test quantities, in which the posterior approximation simply recovers the prior. We refer to Appendix A for further details and corresponding pseudocodes.
Note. Amortized inference lies at the intersection of Bayesian modeling and deep learning. unlocking massive potential for scalable posterior inference. However, this also comes with the practical challenges inherent to training deep neural networks. While detailed treatment of neural architecture design and optimization exceeds the scope of this paper, practitioners libraries like sbi (Boclts et al., 2025) or BayesFlow (Radev et al., 2023), which provide modern plug-and-play components as well as sensible defaults for a wide range of applications. We summarize a set of best practices and actionable recommendations for training amortized posterior estimators in Appendix B.
use established simulation-
inference
If simulation-based calibration checking and parameter recovery diagnostics pass, proceed to Step 1. Otherwise, tune the training hyperparameters (e.g., simulation budget, training epochs, learning rate, or neural network architecture) and re-train the amortized network.
Training phase:
2.2 Inference phase: posterior approximation on observed datasets
Once the amortized estimator is capable of yielding sufficiently accurate posterior draws in closed-world settings (i.c., in-distribution), we use the pre-trained neural network to achieve rapid amortized postcrior inference on a lotal of K observed datasets (8) 1 {yam))=11 Recall that a given pre-trained amorlized neural estimator may be perfectly suitable for some real datasets while it is utterly untrustworthy for others. Therefore, we want assess on per-dataset basis whether the amortized posterior draws are trustworthy and should be accepted, whether should proceed to slower algorithm with stronger accuracy guarantees. The ics in the inference phase are evaluated conditionally each observed dataset, with the ultimate goal of determining whether the set of current posterior draws is acceptable for that specific dataset.
2.2.1 Step 1: Amortized posterior draws
We want to exploit the rapid sampling capabilities of the amortized posterior estimator 46 as much as possi- ble, as long as the sampled posteriors are trustworthy according to set of principled diagnosties. Therefore, (k)
0.4
Training samples (null distribution)
Test samples
Donsity
0.3
OOD cut-off
the natural first step for each observed dataset Yobs is
0.2
to query the amortized posterior and sample S posterior
draws 0/8) ,...,0g Figure 2, first panel).
(("),...,03 ((k.) (0) (x) in near-instant time (see
0.1
0.0
5
10
15
20
25
Diagnostics. Like other neural network approaches (Yang el al., 2024), amortized inference may yield un- faithful results under distribution shifts (Schmitt et al., 2023; Ward al., 2022; Huang al., 2023). To address this, we detoct whether an observed dataset Yobs is out-of- distríbution (OOD) relative to the data-generating process
Mahalanobis distance
Figure 3: Illustration of our caapiing-based hy- pothesis test that flags OOD datasets (to the right of the OOD cut-off).
4

---

Published in Transactions on Machine Learning Research (02/2026)
p)0,y). We first compute low-dimensional summary statistic s(y) R for each dataset. ? The summary statistics from the training dataset (y(()) are used to approxinate the Mahalanobis distance by estimat- ing their empirical mean H4 and covariance Mg. Then, for any test dataset y, its Mahalanobis distance to the training set is:
- Das(y) = (c,))) T TE (s(y)
us))
(3)
We compute (D (m) for all training datasets to establish frequentist sampling distribution of distances under the mull hypothesis (i.c., of in-distribution datasets). Given new observed datasct yobs: we compare its Mahalanobis distance (M((o to the empirical distribution of training distances. We define the OOD rejection rule as:
Reject (Jbbb)) = II Dar (Wobx) > Quantile (IDM (m ))
(4)
where a is by default set to 0.05 and we flag datasets whose Mahalanobis distances fall in the right tail of the empirical training distances as out-of-distribution (see Figure 3). The type-1 error rate a (false rejection) of this test can be set relatively high to obtaín a conservative test that will flag many datasots for detailed investigation in further steps of our workflow.
In a nutshell, this a sampling- based hypothesis test for distribution shifts, similar in spirit the kernel-based test proposed by Schmitt et al. (2023). Since the amortized estimator has no guarantecs nor known error bounds for data outside of the empirical support of the joint model p(0, y) (Flsemüller el al., 2024; Schmitt et al., 2023; Frazier et al., 2024; Elsemüller et al., 2025), we propagate such out-of-distribution datasets to Step 2. It is worth noting that a smaller Mahalanobis distance does not necessarily imply better posterior quality and that this OOD test is only intended to filter out datasets that are most likely to be problematic for the amortized estimator specifically, those requiring extrapolation outside the ellipsoid defined by the training summary statisties.
Alternative diagnostics. In addition to the proposed out-of-distributtio test, more sophisticated data- conditional diagnostics can further assess the accuracy of amortized posterior draws for individual datasets and enhance the reliability of accepted amortized draws. Examples include posterior simulation-based calibration checking (posterior SBC; Säilynoja et al., 2025) or the local classifier two-sample test (L-C2ST; Linhart et al., 2023), to name few. These diagnostics each offer distinct advantages and limitations, but typically require substantially more computation than the OOD test.
Posterior SBC is conceptually straightforward and offers necessary conditions for the accuracy of amortized posterior samples by assessing consistency. However, it requires additional simulations for each test datasct and requires training the amortized estimator on inputs that effectively double the size of the original observations. T-C2ST, which trains classifiers to inguish between (a yy) p(y) and the joint distribution p(0,y) provides theoretically sufficient and necessary conditions for amortized inference accuracy. In practice, however, its effectiveness can be very sensitive to several factors, including classifier design choices (c.g., data pre-processing and optimization stralegies), classifier calibration, and the relative sizes of the simulation budgets allocated to classifier training and amortized estimator training.
The choice to apply these additional diagnostics depends on context-speeific factors, including the number of observed datasets, the relative computational cost of simulations versus likelihood evaluations," 3 and the dimensionality of the observations. Ultimately, whether amortized posterior draws are deemed acceptable hinges on the accuracy requirements of the specific application. By default, recommend the OOD test for its simplicity, efficiency, and suitability as first-line diagnostic.
TTh summary statistic can be either learned by the amortized estimator qe in the training phase be based domain knowledge.
For example, if likelihood evaluations rolatively cheap, instcad of applying sophisticated diagnostics Step 1, is often worthwhile to process directly Step 2, where Partto-smonthed importance sampling can offer more informative and powerful diagnostics.
5

---

Published in Transactions on Machine Learning Research (02/2026)
If the observed dataset passes the OOD test (i.e., Mahalanobis distance below the threshold), accept the amortized draws; otherwise, proceed to Step 2.
Step 1:
2.2.2 Step 2: srroth-s importance sampling
In this step, we use Parcto-smoothed importance sampling (PSIS) (Vehtari et al., 2024) to both improve and assess the quality of the amortized posterior draws of datasets which have previously been rejected (see Figure 2, second pancl). Based on the amortized posterior draws from Step 1, PSTS computes importance weights w = ((y(") 1d=) p()))((())/p)))))))))))))))))) p(0s)/qq ( for each observed dataset y (as in default importance sampling). Then, PSTS fits a generalized Pareto distribution to the largest importance weights, which in turn is used to smooth the tail of the weight distribution (Vehtari al., 2024). Finally, these smoothed importance weights are used for computing posterior expectations and for improving the posterior draws with the sampling importance resampling (SIR) scheme (Rubin, 1988). While the utility of standard importance sampling for improving neural posterior draws has previously been investigated (Dax et al., 2023), we specifically use the PSIS algorithm, which is self-diagnosing (see Diagnostics below) and therefore better suited for a principled workflow. Further details of PSIS are provided in Appendix A.
Diagnostics. We use the Pareto-k diagnostic to gauge the fidelity of the PSIS-refined posterior draws. Pareto-k is the estimated shape parameter of the generalized Pareto distribution and quantifies the tail heaviness of the largest importance weights. According to Vehtari et al. (2024), for moderate sample size (S > 2000), Pareto-k S 0.7 indicates that PSIS estimates are reliable;* when k >> 0.7, the minimum sample size for obtaining a reliable Monte Carlo estimate through (Pareto-smoothed) importance sampling rapidly grows infeasibly large in practice, implying that the amortized posterior is poor proposal for importance sampling correction and the correspouding dataset should be rouled to Step 3. This & threshold is consistent with the established practice of using PSIS to improve and assess the quality of posterior approximations obtained from variational inference (Yao et al., 2018; Dhaka et al., 2021; Zhang et al., 2022).
Note. The posterior estimator in ABI is typically mode-covering since it optimizes the forward KI, divergence in Eq. 1. When the neural network training is insufficient (e.g., small simulation budget or poorly optimized network), this may lead to overdispersed posteriors. Fortunately, this tends to err in the right direction, and PSTS can generally mitigate overdispersed mode-covering draws in low to moderate dimensions (Dhaka et al., 2021). In contrast, variational inference typically optimizes the reverse KL divergence (Rezende & Mohamed, 2015), which implies mode-seeking behavior that. is less favorable for importance sampling.
Step 2: If Pareto-k S 0.7, accept the importance sampling results; otherwise, proceed to Step 3.
2.2.3 Step 3: Many-chains MCMC with amortized initializations
If PSIS does not yield satisfactory results, we resort to an MCMC sampling scheme as a safe fallback option. In our amortized workflow, the MCMC step is augmented by reusing computations from the previous steps as initialization values. In principle, this step can incorporate any MCMC algorithm suited to the problem at hand. Examples include slice sampling for models with non-differentiable likelihoods (Neal, 2003), or HMC (Neal, 2011) samplers when gradients are available.
In this work, we use the ChEES HMC algorithm (Hoffman et al., 2021) as an instantiation of MCMC. Most notably, ChEES-HMC supports the execution of thousands of parallel chains on a GPU for high- throughput sampling (Sountsov et al., 2024). Amortized posterior draws from previous steps provide a natural and convenient choice for initializing MCMC chains to accelerate convergence (Figure 4). This approach is conceptually similar to using methods like parallel quasi-Newton variational inference (i.e., Pathfinder; Zhang et al., 2022) to obtain initial values for MCMC chains. However, the amortized initial values are drawn in parallel in near-instant time, while Pathfinder requires re-fitting the varia- tional approximation for each new observed dataset. For the purpose of ChEES-HMC initializations
4 For small sample size (S < 2000), the threshold of Pareto-& is min(1 (S).
6

---

Published in Transactions on Machine Learning Research (02/2026)
with multimodal posterior distributions, it is again desirable that the amortized posterior draws are typ- ically mass-covering (cf. Step 2). See Appendix A for additional details on the ChEES-HMC algorithm.
Diagnostics. In this last step, we use the nested R diagnos- tic (Margossian et al., 2024), which is specifically designed to assess the convergence of the many-but-short MCMC chains. 5 If the diagnostics this step indicate unreliable inference, we recommend resorting to the overarching Bayesian workflow (Gel- man et al., 2020) and addressing the computational issues that even persist when using the (ChEES-)HHN) algorithm. This could involve increasing the number of warmup iterations, using the established NUTS-HMC algorithm (Hoffman & Gelman, 2014; Carpenter et al., 2017), or revising the Bayesian model specification and parametrization.
*He Amortized initializations
ChEES-HMC samples
Target posterior
01
Figure 4:
We initialize many ChEES-HMC chains with amortized draws.
Step 3: If (nested) Å is below the convergence threshold (e.g., 1.01), accept the MCMC draws. Otherwise, increase warm-up or revise the model according to the standard Bayesian workflow (Gelman et al., 2020).
2.3 Related work
Both simulation-basel inference and amortized inference have seen rapid progress over the past decade (Zammit-Mangion et al., 2025; Cranmer et al., 2020; Lavin et al., 2021), driven by the need to perform Bayesian inference in complex models with intractable likelihoods (e.g., Dingeldein et al., 2024; Wehenkel et al., 2024; Zhou et al., 2024; Ghaderi- Kanggvvri et al., 2023; von Krause et al., 2022; Bieringer et al., 2021; Radev et al., 2021). These advances have been fueled by modern generative modeling, such as normalizing flows (Papamakarios et al., 2021; Radev et al., 2020; Greenberg et al., 2019), transformers (Müller et al., 2022; Chang et al., 2025; Whittle et al., 2025), diffusion models (Song et al., 2021; Sharrock et al., 2024; Linhart. et al., 2024; Geffner al., 2023; Gloeckler et al., 2024), consistency models (Song al., 2023; Schmitt et al., 2024b), and flow matching (Lipman et al., 2023; Wildberger et al., 2023). Practical software toolkits such BayesFlow (Radey et al., 2023) and sbi (Boelts et al., 2025) further make these simulation-based inference techniques accessible to practitioners in user-friendly interfaces.
To address the potential systematic errors of (amortized) neural posteriors, several works propose corrections using importance reweighting schemes (Dax et al., 2023; Starostin et al., 2025), augmented training objectives (Delaunoy et al., 2022; Mishra et al., 2025; Orozco et al., 2025; Schmitt al., 2024a), or post-hoe corrections (Siahkoohi et al., 2023). Simultaneously, hybrid approaches that combine density estimators with MCMC have gained traction (Salimans et al., 2015; Hoffman et al., 2019; Gabrié et al., 2022; Midgley et al., 2022; Arbel et al., 2021; Cabezas et al., 2024; Grenioux et al., 2023). These include using variational approximations or learned flows as preconditioners for MCMC (Hoffman et al., 2019; Cabezas & Nemeth, 2023), adaptive proposal mechanisms (Parno & Marzouk, 2018; Gabrié et al., 2022), and initialization strategies to accelerate convergence or improve diagnostics (Zhang et al., 2022; Wang et al., 2023; Starostin et al., 2025).
More broadly, automated Bayesian inference has been contral design goal of probabilistic programming systems such Stan (Carpenter et al., 2017), PyMC (Oriol et al., 2023), (Num)Pyro (Bingham et al., 2019; Phan et al., 2019). These libraries provide general-purpose inference engines typically gradient-based MCMC and variational inference that can be applied to wide range of likelihood--bsse models and are accompanied by well-developed diagnostic recommendations and workflow guidelines (Gelman et al., 2020). However, they do not natively support amortized inference across many datasets, and inference must be rerun from scratch for each dataset.
Our proposed workflow builds on and complements these lines of work by integrating amortized inference, likelihood- based correction, and many-chain MCMC into unified, modular, and diagnostic-driven pipeline for accelerating Bayesian inference. It dynamically adapts the inference strategy to the dataset at hand using
in more conventional settings involving long MCMC chains, the standard R diagnostics (Gelman & Rubin, 1992; Vehtari et. al., 2021) can be applied.
7

---

Published in Transactions on Machine Learning Research (02/2026)
amortized posterior draws when they are adequate and escalating to PSIS and MCMC improving the robustness of amortized inference and the overall efficiency of posterior computation. This modular design provides a practical foundation for principled amortized inference across diverse data regimes.
3 Experiments
In this section, we empirically evaluate the effectiveness of our proposed amortized Bayesian workflow across various synthetic and real-world problems. We also examine how reusing amortized posterior draws in subsequent steps can improve the downstream sampling performance. The source code reproduce all experiments is available in the supplementary material.
3.1 Procedure
Training settings. For each problem, begin by training the amortized posterior estimator simulated parameter-obserratton pairs (i..,, simulation- training). We verify that the model performance is satisfactory in closed-world setting, as diagnosed by simulation-based calibration and parameter recovery checking (see Section 2.1). Details on diagnostic results, simulation budgets, and training hyperparameters are provided in Appendix C.
Inference settings. For the out-of-distribution diagnostics in Step 1, we use the 0.05 as the rejection threshold. We compute Mahalanobis distances in the summary statistics using 10,000 training simulations. We draw 2,000 posterior samples from the amortized posterior at Step 1. In Step 2, we correct the amortized draws using PSIS, rejecting draws if Pareto-k > 0.7. Step 3 uses CAEES-HMC with convergence determined by nested R < 1.01. We run 2048 chains in parallel (16 superchains, each with 128 subchains), with 200 warmuu steps and single sampling step, for a total of 2048 posterior draws.
Evaluation metrics. To assess the quality of postcrior draws from workflow, compare them to reference posterior draws using Iwo evaluation metries: the 1-Wasserstein distance (W1) and the mean marginal total variation distance (MMTV). The W1 distance quantifies the overall discrepancy between full joint distributions. MMTV measures the lack of overlap between marginal distributions and takes value in the range 1]: for example, an MMTV value of 0.2 implics that, on average, the approximate posterior draws and reference draws share an 80% overlap for their marginal distributions. For both metrics, lower values indicate better posterior approximation quality. As rule of thumb, MMTV values below 0.2 indicate good posterior approximation fidelity (Acerbi, 2020; Li et al., 2025).
3.2 Applications
We apply the proposed workflow to four posterior inference problems, including both simulated benchmarks and real-world experimental datasets. These case studies were chosen to reflect range of commonly encountered statistical inference scenarios, including classical distributional parameter estimation and analyses of large- scale datasets arising in psychology and cognitive modeling. We describe each problem briefly below, with further details provided Appendix C.
Generalized extreme value distribution (GEV). We consider parameter inference for the generalized extreme value (GEV) distribution, which models the maxima of samples from a distribution family. Each observation yi is modeled as:
yi GEV(H, a,))
55)
where μ E IR is the location, a E Ry0 is the scale, and E IR is the shape parameter. We follow the prior specification from Caprani (2021). For each dataset, we collect N = 65 i.id. observations and infer the posterior distribution over the parameter vector == (M, We generate a total of K 1000 test datasets by deliberately simulating from model with 2x wider prior distribution to emulate out-o--distribution settings in real applications (see Appendix C for details).
8

---

Published in Transactions on Machine Learning Research (02/2026)
Bernoulli GLM. The Bernoulli generalized linear model (GLM) is classical model with binary outcomes, included in the SBI benchmark suite (Tueckmann et al., 2021). Each observation y/ E {0, 1 is modeled as:
y (ernoulli(a(;; 0)),
(6)
where v; _. is fixed input vector, 0 E I is the parameter vector, and (-)) denotes the logistic function. We generate K = 10, 000 in-distribution test datasets by sampling parameters from the model prior and simulating corresponding observations {y)) (Lueckmann et al., 2021).
Psychometric curve fitting. Psychometric functions are widely used in perceptual and cognitive science to characterize the relationship between stimulus intensity and the probability of specific response (Wichmann & Hill, 2001). We use the overdispersed hierarchical model from Schütt et al. (2016), where the number of correct trials yi al stimuli level x; is modeled as:
(---1) Pi,
L (pp (1 (+-1))(-7))
yi Binomial (na, Pi Beta
where is the number of trials, 17 E [0, 1] controls overdispersion, and pi (1,;; m, w. X, 7) is the expected success probability given by the psychometric function s)(2;; m, w, 7) = + (1 - X --7) S(x; m. w), where S is a sigmoid function (e-g., cumulative normal), m is the threshold, w is the width, À is the lapse rate for infinitely high stimulus levels, and y is the guess rate for infinitely low stimulus levels. In total, the model parameters are 0 (m, w, X, Y, n)). Our empirical evaluation uses 8,526 mouse behavioral datasots from the International Brain Laboratory public database (The International Brain Laboratory et al., 2021).
Decision model. The drift-diffusion model (DDM) is popular evidence accumulation model for psycho- logical models of human decision making (Rateliff& & McKoon, 2008). It describes a two-choice decision task as a stochastic process in which noisy evidence accumulates over time until it reaches one of the decision boundaries. The evolution of the decision variable z(t) is modeled as
dz(() = dll + dW),
(8)
where v is the drift rate (the average rate of evidence accumulation), a is the noise scale, and W(t) denotes a standard Wiener process. A decision is made when z(t) reaches either a positive or negative boundary, typically placed symmetrically at Ea, where a the boundary separation. The model also includes a non-decision time parameter т, capluring processes that are not part of the decision process. We adopt. the model specification from von Krause et al. (2022), which extends the standard DDM to incorporate experimental condition effects via six paramcters: = (V1, 02 01, a22 Tn). The test datascts consist of 15,000 participants from the online implicit association test. (IAT) dalabase (Xu et al., 2014; von Krause et al., 2022), providing a large-scale, real-world benchmark for Bayesian inference in cognitive modeling.
3.3 Main Results
Table 1 summarizes the performance of the proposed amortized Bayesian workflow across the four problems described in Section 3.2. Step 1 (ABI) exhibits extremely low time per accepted datasct (TPA). with most of the cost. incurred as a one-time expense during the training phase-includin training, and diagnostic evaluation. Once trained, ABT incurs negligible marginal cost (& 1sec) when applied to a new dataset. Datasets flagged as out-of-distribution in Step are forwarded to Step 2 for correction via PSIS. PSIS is highly effective, successfully correcting most rejected amortized draws and substantially reducing the number of datasets requiring full MCMC. Only small subset of datasets progresses to Step 3, where ChEES-HMC. is used for high-fidelity sampling. As the most computationally expensive component, CLEES-HMC is applied selectively, allowing the workllow lo retain both accuracy and efficiency. Overall, the amortized workflow completes inference for nearly all datasets. 6 Compared to using ChEES-HMC for all datasets, our workflow achieves substantial computational savings- apprrss-appxxxittllly over 5x, 120x, 60x, and 15x faster for the GEV, Bernoulli GLM, psychometric curve, and decision model tasks, respectively.
prior simulation, model
A small number of datasets with particularly difficult properties require extended MCMC runs converge.
9

---

Published in Transactions on Machine Learning Research (02/2026)
Table 1: Summary of our amortized Bayesian workflow across four problems. For each step, we report the number accepted datascts, wall-clock time (minutes), and time per accepted dataset (TPA) in seconds. The time for the training phase includes amortized estimator training, simulations, and diagnostics evaluations. The time for Step 1 includes amortized posterior draws and the OOD test. The TPA for Step accounts for both the training phase and Step 1. "Workflow total" aggregates the results of our method across all steps. As baseline reference, "Baseline workflow total" is an estimate of the total required runtime for ChEES-HMC on all datasets.
Problem
Step
Accepted datasets
Time (min)
TPA (s)
Training phase
0.4
Step 1: Amortized inference
523/1000
0.1
GEV
Step Amortized + PSIS
357 357/477
0.8
0.1
Step 3: ChEES-HMC w/ inits
87/120
11
7
Workflow total (ours)
967/1000
15
0.9
Baseline workflow total
85
Training phase
0.8
0.007
Step 1: Amortized inference
9519/10000
0.3
Bernoulli GLM
Step 2: Amortized + PSTS
425/481
0.4
0.06
Step 3: TMEES-HMC. w/ inits
56/56
4
4
Workflow total (ours)
10000/10000
0.03
Baseline workflow total
688
Training phase
0.06
Step 1: Amortized inference
7213/850
0.4
Psychometric curve
Step 2: Amortized + PSIS
1215/1313
4
0.2
Step 3: ChEES-HMC w/ inits
69/98
26
22
Workflow total (ours)
8497/8528
0.3
Baseline workflow total
2217
Training phase
85
0.4
Step 1: Amortized inference
13/98/15000
Decision model
Step 2: Amortized + PSIS
827 /27/1502
47
3
Step 3: ChEES HMC w/ inits
554/675
526
57
Workflow total (ours)
14879/15000
659
3
Baseline workflow total
11594
Figure 5 presents the quality of posterior draws using the WI distance (top row) and MMTV distance (bottom row), comparing draws from each step of the workflow against reference posteriors obtained via well-tuned NUTS. Rejected amortized draws (ABIX) exhibit markedly worse performance than accepted ones (ABI)), confirming the effectiveness of the OOD diagnostics." PSIS-corrected draws offer accuracy comparable to CLEES-HMC samples, with only a slight decrease in quality. While amortized draws accepted in Step 1 are less accurate than those produced by PSIS or ChEES-HMC, they still provide high-quality approximations across the majority of datascts, implicd by the WI and MMTV metrics. These results demonstrate that the proposed workflow not only scales efficiently but also consistently produces high-quality posterior estimates.
3.4 Advantage of amortized initializations for MCMC
One major goal of our workflow is minimize reliance on expensive MCMC by maximizing the reuse of computations. Even when ABI and the PSIS refinement fail to yicld acceptable posterior draws after Step 2, can still leverage the amortized outputs to accelerate MCMC in Step 3.
To evaluate whether amortized posterior estimates remain useful in such cases, we test their effectiveness as initializations for thEES-HMC chains. We conduct experiments on 20 randomly selected test datasets that progress to Step 3 of the workflow. This indicates that both the amortized posterior draws and their Pareto-
For the Bernoulli GLM, the rejected amortized draws appear good quality because the datasets are drawn directly from the same prior used during training (l.e., in-distributtion.
10

---

Published in Transactions on Machine Learning Research (02/2026)
GEV
Bernoulli GLM
Psychometric curve
Decision model
2.4
WT
1.6
0.8
0.0
0.8
MMTV
0.5
0.2
0.0
ABI()) ADI()) PS18 C-HMC ABI()) ABI() PS1S C-HMC ABIV) ADI(K) PSTS C-HMC ABI(V) ABI(0) PS1S C-HMC
Figure 5: Evaluation of posterior draws across four problems based on two metrics: W1 distance (top row) and MMTV distance (bottom row). Lower values indicate better posterior approximation. ABI(V) and ABI(x) denote accepted and rejected draws, respectively, from amortized Bayesian inference in Step 1. PSIS denotes importance-weighee draws accepted in Step 2, and C-HMC denotes draws accepted via ChEES-HMC in Step 3. Metrics are computed on up to 100 datasets for each type of draws.
GEV
Bernoulli GLM
10
Nested R
Figure 6: Using amortized posterior draws initializa- tions for ChEES-HMC re- duces the required warmup in the GEV and decision model tasks. We show modian=IQR across 20 test datasets in Step 3.
10-1
10-2
10-2
10-3
10-3
Psychometric curve
Decision model
Nested R R-1
10-1
10-1
10:
10-2
10
10
1030
100
200
300
500
1050
100
200
300
Warmup iterations
Warmup iterations
Amortized
Amortized PSIS
Random Init.
smoothed refinement are deemed unacceptable, quantified by Pareto-k 0.7 in Step 2. We compare three initialization methods for CnEFS-HMC chains: (1) amortized posterior draws, (2) PSTS-refined amortized draws, and (3) a random initialization scheme similar to Stan (Carpenter et al., 2017). We run the chains for varying numbers of warmup iterations, followed by a single sampling iteration. As described in Section 2, we use the nested R value to gauge whether the chains converged appropriately during the warmup stage, as quantified by the common R threshold of 0.01 (Vehtari al., 2021).
Figure 6 shows that amortized posterior draws (and their PSIS-refined counterparts) can significantly reduce the required number of warmup iterations to achieve ChEES-HMC chain convergence, though the draws themselves have previously been flogged as unacceptable. For the GEV problem and the docision model, chains initialized with amorlized draws converge faster than those using random initialization. In the Bernoulli GLM, all methods perform similarly. For the psychometric curve model, random initialization leads to faster convergence for the carly stage, but amortized draws still reach the convergence threshold similar speed at iteration 200, indicating competitive performance. These findings are particularly relevant in the many-short-chains regime, where computational cost is dominated by the warmup phase. For instance,
11

---

Published in Transactions on Machine Learning Research (02/2026)
with 2048 parallel chains, every single post-warmup step yields 2048 posterior samples, leading to enormous efficiency gains from shorter warmup.
Overall, these results demonstrate that amortized inference may provide suitable initializations for ChEES- HMC. However, the added benefit of initializing chains with PSIS-refined amortized draws (Step 2) instead of raw amortized draws (Step 1) remains unclear. While PSTS often accelerates convergence, it occasionally degrades worst-case performance (see upper error bounds for GEV task in Figure 6). We further study the impact of initialization for the popular NUTS sampler (Hoffman & Gelman, 2014), with similar results: amortized initializations reduce the required warmup in most cases (see Appendix E).
4 Discussion
We presented an adaptive Bayesian workflow to combine the rapid speed of amortized inference with the undisputed sampling quality of MCMC. Our amortized workflow enables fundamental shift in the scale and fcasibility of Bayesian inference. Applying traditional MCMC (e-g., ChEES-HMC) within standard Bayesian workflow to every dataset independently would require approximately 10 days of GPU computation across our experimental suite. In contrast, our amortized workflow completes inference in half a day, achieving speedups ranging from over 5x to 120x depending on the problem. Crucially, high-quality posterior draws are retained through a cascade of diagnostics and selective escalation to PSIS and MCMC. In conclusion, our workflow efficiently uses resources by (i) applying fast amortized inference when the results are accurate; (ii) refining draws with PSIS when possible; and (iii) amortized initializations of slower but accurate MCMC chains when needed.
Modularity and practical flexibility. A key strength of the proposed workflow lies in its modular structure, which allows practitioncrs tailor each component to the specific constraints and objectives of their application. In cases where preliminary analysis or low-latency decision-making is essential (e-g., real-time experimental pipelines) where likelihood evaluations are computationally expensive, the workflow can operate in a lightweight mode using amortized inference with out-of-distribution rejection alone (i.e., Step in our workflow). Conversely, in high-stakes applications where accuracy is paramount, analysts can enforee escalation of all amortized draws through PSIS and, if needed, proceed to full MCMC to guarantee statistical robustness. The choice of MCMC sampler in Step 3 is also fully interchangeable: alternative algorithms such as slice sampling, ensemble samplers (e.g., emcee; Foreman- Macaey et al., 2013), or NUTS can be substituted if the model is -in-different,,ble, multimodal, or requires richer exploration.
Furthermore, while our paper focuses the trade-off between wall-clock inference speed and posterior quality, practical deployments may also involve additional factors, such inference cost (e.g., monetary expense for GPU/CPU hours) and hardware availability. Consequently, the most suitable workflow variant can differ across settings. For example, when GPU resources are limited, launching parallel MCMC chains on CPUs offers a practical alternative, making the workflow more accessible for a broader range of users.
Applicability, limitations and future directions. Our proposed workflow targets likelihood-basee Bayesian models for which prior predictive simulation and likelihood evaluation are possible. It is most beneficial in repeated-infereet infereen-ee regimes (many datasets or frequent re-fits), with moderate effective dimen- sionality, and when a good amortized estimator can be trained once and subsequently reused. Hence, it not universally suitable and does not yield inference speedup gains for all Bayesian models. Training amortized models requires upfront investment in optimization and simulation. In our experiments, we found that default neural network hyperparameter settings, such normalizing flow architectures, summary network configurations, and optimizer settings, generally yield good performance.
However, in more challenging cases, such as the GEV problem, adjustments may be necessary, guided by training phase diagnostics. The simulation burden can be exacerbated in high-dimensional (Z 10 parameters) or weakly identifiable models, where neural estimators may struggle to approximate complex inverse maps. Alternative amortized inference approaches (sce, e.g. Mittal al., 2025) could be explored in future work to complement sinulation- based amortized inference in such scenarios. In settings where likelihood evaluations
12

---

Published in Transactions on Machine Learning Research (02/2026)
are particularly expensive, iterative refinements of the amortized estimator individual datasets (Glöckler et al., 2022) may also be practical alternative to likelihood-based corrections in Steps 2 and 3.
Moreover, while our diagnostic for the amortized posterior draws in Step 1 is effective and highly efficient in practice, it remains an imperfect proxy for the true posterior approximation error and can occasionally result. the acceptance of poor-quality amortized draws (ef. Figure 5). An additional empirical study in Appendix D shows that (1) the Mahalanobis distance and the posterior quality metrics are positively correlated when OOD datasets present; (2) however, some low-distance datasets still yield poor metries, highlighting the limitation that the OOD diagnostic camnot fully guarantee accuracy and motivating enforced escalation to PSTS (Step 2) when higher accuracy is a requirement. Future work could explore even more effective discrepancy measures, polentially tailored to the task at hand.
More broadly, the workflow supports a compelling vision training amortized models once and reusing them across tasks or studics-a a strategy well suited to applications ranging from psychology to computational biology, among others. In such settings, our layered diagnostics and selective escalation are crucial for maintaining reliability and efficiency. This positions the workflow practical bridge between amortized inference and traditional Bayesian rigor, enabling scalable yet trustworthy inference.
Acknowledgments
We thank Alexander Fengler for the help with PyMC implementation of the drift-diffusion model. CL and LA were supported by the Research Council of Finland (grants nunber 356498 and 358980 to LA). AV acknowledges the Research Council of Finland Flagship program: Finnish Center for Artificial Intelligence, and Academy of Finland project 340721. STR is supported by the National Science Foundation under Grant No. 2448380. MS and PB acknowledge support of Cyber Valley Project CyVy-RF- 2021-16, the DFG under Germany's s Excellence Strategy FXCC2075 390740016 (the Stuttgart, Cluster of Excellence SimTech).. PB additionally acknowledges the support of DFG Collaborative Research Cluster 391 "Spatio-Tempor"" Statistics for the Transition of Energy and Transport" 520388526. MS acknowledges travel support from the European Union's Horizon 2020 research and innovation programme under grant agreements No 951847 (ELISE) and No 101070617 (ELSA), and support from the Aalto Science-IT project. The authors wish to thank the Finnish Computing Competence Infrastructure (FCCI) for supporting this project with computational and data storage resources. The authors also acknowledge the research environment provided by ELLIS Institute Finland.
References
Luigi Acerbi. Variational Bayesian Monte Carlo with noisy likelihoods. Advances in Neural Information Processing Systems, 33:8211-8222, 2020.
Michael Arbel, Alex Matthews, and Arnaud Doucet. Annealed flow transport Monte Carlo. In Proceedings of the 38th International Conference on Machine Learning, pp. 318-330. PMLR, July 2021.
Sebastian Bieringer, Anja Butter, Thco Heimel, Stefan Höche, Ullrich Köthe, Tilman Plchn, and Stefan T Radev. Measuring qed splittings with invertible networks. SciPost Physics, 10(6):126, 2021.
Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman. Pyro: Deep universal probabilistic programming. Journal of Machine Learning Research, 20(28):1-6, 2019. ISSN 1533-7928.
Jan Boelts, Michael Deistler, Manuel Gloeckler, Álvaro Tejero-Cantero, Jan-Matthis Lueckmann, Guy Moss, Peter Steinbach, Thomas Moreau, Fabio Muratore, Julia Linhart, Conor Durkan, Julius Vetter, Benjamin Kurt Miller, Maternus Herold, Abolfazl Ziacemehr, Matthijs Pals, Theo Gruner, Scbastian Bischoff, Nastya Krouglova, Richard Gao, Janne K. Lappalainen, Bálint Mucsányi, Felix Pei, Auguste Schulz, Zinovia Stefanidi, Pedro Rodrigues, Cornelius Schröder, Faried Abu Zaid, Jonas Beck, Jaivardhan Kapoor, David S. Greenberg, Pedro J. Gonçalves, and Jakob H. Macke. Sbi reloaded: toolkit for simulation-bart inference workflows. Journal of Open Source Software, 10(108):7754, April 2025. ISSN 2475-9066. doi: 10.21 105/joss.00754.
13

---

Published in Transactions on Machine Learning Research (02/2026)
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Mflne, and Qiao Zhang. JAX: composable transformations of Python+-Numly programs, 2018. URL http://github. cwwmcom////o/o google/jax.
Alberto Cabezas and Christopher Nemeth. Transport elliptical slice sampling. In Proceedings of The 26th International Conference on Arlificial Intelligence and Statistics, pp.
PMLR, April 2023.
Alberto Cabezas, Louis Sharrock, and Christopher Nemeth. Markovian flow matching: Accelerating MCMC with continuous normalizing flows. Advances in Neural Information Processing Systems, 37:104883-104111 December 2024.
Colin Caprani. Generalized extreme value distribution. http:///wo..o/////o/o//////////o html, 2021. PyMC Examples. Accessod: 2025-05-13.
Bob Carpenter, Andrew Gclman, Matthew D Hoffman, Daniel Loe, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of statistical software, 76(1), 2017.
Paul E Chang, Nasrulloh Loka, Daolang Huang, Ulpu Remes, Samuel Kaski, and Luigi Acerbi. Amortized probabilistic conditioning for optimization, simulation and inference. In Proceedings of the International Conference on Arlificial Intelligence and Statistács (AISTATS). PMLR, 2025.
Gang Chen, Paul-Christian Bürkner, Paul A. Taylor, Zhihao Li, Lijun Yin, Daniel R. Glen, Joshua Kinnison, Robert W. Cox, and Luiz Pessoa. An integrative Bayesian approach to matrix based analysis in neuroimaging. Human Brain Mapping, 40(14):4072 4090, 2019. doi: 10.1002/hbb 10002/hbb.26665 24686.
François Chollet et al. Keras. bttpss//keraat io, 2015.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation based inference. Proceedings of the National Academy of Sciences, 2020.
Maximilian Dax, Stephen R. Green, Jonathan Gair, Michael Pürrer, Jonas Wildberger, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf. Neural importance sampling for rapid and reliable gravitational-wave inference. Phys. Rev. Lett., 130:171403, Apr 2023. doi: hht/PhyysRee..t....... URL https://link, aps. wwwo/oggg//....oooooooooooo.oo /100//10.110 1101 www.oooy///eeeeet...oeom 130. 1tt.::eeee................... 171403.
Arnaud Delaunoy, Joeri Hermans, François Rozet, Antoine Wehenkel, and Gilles Louppe. Towards reliable simulation-based inference with balanced neural ratio estimation. Aduances in Neural Information Processing Systems, 35:20025 20037, December 2022.
Akash Kumar Dhaka, Alejandro Catalina, Manushi Welandawe, Michael R Andersen, Jonathan Huggins, and Aki Vehtari. Challenges and opportunities in high dimensional variational inference. In Advances in Neural Information Processing Systems, volume 34, 7787 7798. Curran Associates, Inc., 2021.
Joshua V. Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan, Dave Moore, Brian Patton, Alex Alemi, Matt Hoffman, and Rif A. Saurous. TensorFlow distributions, 2017.
Lars Dingeldein, David Silva-Sánchez, Luke Evans, Edoardo D' D'Imprima, Nikolaus Grigorieff, Roberto Covino, and Pilar Cossio. Amortized template-matching of molecular conformations from cryo-clectron microscopy images using simulation-barn inference. bioRaiv, pp. 2024.07.22 604154, 2024. doi: 10. 10.1101/2024.07 604154, URL
23. 604154. abstract.
Conor Durkan, Artur Bekasov, Tain Murray, and George Papamakarios. Neural spline flows. Advances in neural information processing systems, 32, 2019.
Lasse Elsemüller, Hans Olischläger, Marvin Schmitt, Caul-Christian Bürkner, Ullric Koethe, and Stefan T. Radev. Sensitivity-aware amortized Bayesian inference. Transactions Machine Learning Research, 2024.
14

---

Published in Transactions on Machine Learning Research (02/2026)
Lasse Elsemülleer, 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, 2025.
Alexander Fengler, Yang Xu, Bera Krishn, Aisulu Omar, and Michael J. Frank. HSSM: generalized toolbox for hierarchical Bayesian estimation of computational models in cognitive neuroscience. Manuscript in preparation, 2025. URL https://github. com/ hean/Inccbrow//HHe
D. Frreman-Macker D. W. Hogg, D. Lang, and J. Goodman. emcee: The MCMC hammer. PASP, 125: 306-312, 2013. doi: 00.1086/670067
David T. Frazier, Ryan Kelly, Christopher Drovandi, and David J. Warne. The statistical accuracy of neural posterior and likelihood estimation, 2024.
Marylou Gabrié, Grant M. Rotskoff, and Eric Vanden- Eijnden. Adaptive Monte Carlo augmented with normalizing flows. Proceedings of the National Academy of Sciences, 211140),,,004220119, March 2022. doi: 10.1073 211773/pma..2111 pmas22209420119,
Tomas Geffner, George Papamakarios, and Andriy Mnih. Compositional score modeling for simulation- based inference. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), Proceedings of the 40th international conference on machine learning, volume 202 of Proceedings of machine learning research, pp. 11098 11116. PMLR, 2023. URL https: //proceedings. mlr. wwt.srrerss//////eeeeee 20.002/gofnerr23a.mm html.
Andrew Gelman and Donald B. Rubin. Inference from Iterative Simulation Using Multiple Sequences. Statistical Science, 7(4):457 472, November 1992. ISSN 0883-4237, 2168-8745. doi: 10.1214/ss/ 0ts//1177111100
Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin. Bayesian Data Analysis. Chapman & /CC, Philadelphia, PA, 3 odition, 2013.
Andrew Gelman, Aki Vehtari, Daniel Simpson, et al. Bayesian workflow. arxiv preprint, 2020.
J.-P. George, P.-. Bürkner, T. G. M. Sanders, M. Neumann, C. Cammalleri, J. V. Vogt, and M. Lang. Long-term forest monitoring reveals constant mortality rise in European forests. Plant Biology, 24(7): 1108-1119, 2022. doi: 10.1111/plb.134/
Amin Ghaderi- Kangavar,, Jamal Amani Rad, and Michael D. Nunez. A general integrative neurocognitive modeling framework to jointly describe EEG and decision-making on single trials. Computational Brain & Behavior, 6(3):317-377 376, 2023. ISSN 2522-0861, 2522-087X. doi: 11t///1..100/2111111111
URL
Manuel Glöckler, Michael Deistler, and Jakob H. Macke. Variational methods for simulation-based inference. In International Conference on Learning Representations, 2022.
Manuel Gloeckler, Michael Deistler, Christian Weilbach, Frank Wood, and Jakob H Macke. All-in-one simulation-based inference. In Proceedings the International Conference on Machine Learning (ICML), pp- 15735 -1766,, 2024.
Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359 378, 2007.
Pedro J Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F Podlaski, Sara A Haddad, Tim P Vogels, David S Greenberg, and Jakob Macke. Training deep neural density estimators to identify mechanistic models of neural dynamics. eLife, 9ce56261, September 2020. ISSN 2050-084X. doi: 10. 10.7554/eLife.560 56261.
David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood- free inference. In International Conference on Machine Learning, 2019.
15

---

Published in Transactions on Machine Learning Research (02/2026)
Louis Grenioux, Alain Oliviero Durmus, Eric Moulines, and Marylou Gabrié. On sampling with approximate transport maps. In Proceedings of the 40th International Conference on Machine Learning, pp- 11698-1111111 11733.
PMLR, July 2023.
Joeri Hermans, Arnaud Delaunoy, François Rozet, Antoine Wehenkel, Volodimir Begy, and Gilles Louppe. A crisis in simulation-based inference? Beware, your posterior approximations can be unfaithful. Transactions on Machine Learning Research, September 2022. ISSN 2835-8856.
Matthew Hoffman, Pavel Sountsov, Joshua V. Dillon, lan Langmore, Dustin Tran, and Srinivas Vasudevan.
Neu NeuTra-lizing bad geometry in Hamiltonian Monte Carlo using neural transport, March 2019.
Matthew Hoffman, Alexey Radul, and Pavel Sountsov. An adaptive-MCMC scheme for setting trajectory lengths in Hamiltonian Monte Carlo. Arindam Banerjee and Kenji Fukumizu (eds.), Proceedings of The 24th International Conference on Artificial Intelligence and Statisties, volume 130 of Proceedings of Machine Learning Research, pp. 3907 3915. PMLR, 13 15 Apr 2021.
Matthew D. Hoffman and Andrew Gelman. The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(47)11593 1623, 2014.
Daolang Huang, Ayush Bharti, Amauri Souza, Luigi Acerbi, and Samuel Kaski. Learning robust statistics for simulation-based inference under model misspecification, 2023.
Ravin Kumar, Colin Carroll, Ari Hartikainen, and Osvaldo Martin. Arviz a unified library for exploratory analysis of bayesian models in python. Journal of Open Source Software, 4(33):1143, 2019. doi: 10.21105/ 0114.. URL attps://doi... org/ 20105/ 01143.
Nils C. Landmeyer, Paul-Christian Bürkner, Heinz Wiendl, Tobias Ruck, Hlans-Peee Hartung, Heinz Holling, and Meuth. Disease- modiffyig treatments and cognition in relapsing-remitting remitting multiple sclerosis: A
meta-analysis. Neurology, 94(2)):2373-238) 2020. doi:
Alexander Lavin, Hector Zenil, Brooks Paige, et al. Simulation intelligence: Towards a new generation of scientific methods. arXiv preprint, 2021.
Juho Lee, Yoonho Lee, Jungtack Kim, Adam Kosiorek, Seungjin Choi, and Yee Whye Teh. Set transformer: A framework for attention-based permutation-invariant neural networks. In Proceedings of the 36th International Conference on Machine Learning, pp. 3744-3753, 2019.
Chengkun Li, Bobby Huggins, Petrus Mikkola, and Luigi Acerbi. Normalizing flow regression for Bayesian inference with oflline likelihood evaluations. In 7th Symposium on Advances in Ayprozimale Bayesion Inference, 2025.
Julia Linhart, Alexandre Gramfort, and Pedro L. C. Rodrigues. L-C2ST: Local diagnostics for posterior approximations in sinulation-based inference. In Thirly-seventh Conference Neural Information Processing Systems, 2023.
Julia Linhart, Gabriel Victorino Cardoso, Alexandre Gramfort, Sylvain Le Corff, and Pedro L. C. Rodrigues. Diffusion posterior sampling for simulation based inference in tall data settings, June 2024. URL http: www.xiv..rg/bb////oo 07593. arXiv:2401.07593 (es, stat].
Yaron Lipman, Ricky T. Q Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In The Eleventh International Conference Learning Representations, 2023.
llya Loshchilov and Frank Hutter. Decoupled weight decay regularization. In International Conference on Learning Representations, 2019.
Jan-Matthis Lueckmann, Jan Boclts, David Greenberg, Podro Goncalves, and Jakob Macke. Benchmarking simulation- based inference. In Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, pp. 343 351.. PMLRR March 2021.
16

---

Published in Transactions on Machine Learning Research (02/2026)
Tunlia Malén, Tomi Karjalainen, Janne Isojärvi, Aki Vehtari,
Bürkner, Vesa Putkinen, Valtteri Kaasinen, Jarmo Hietala, Pirjo Nuutila, Juha Rinne, and Lauri Nummenmaa. Atlas of type dopamine receptors in the human brain: Age and sex dependent variability in large PET cohort. Neurolmage, 255: 119149, 2022. doi:
Charles C. Margossian, Matthew D. Hoffmaan Pavel Sountsov, Lionel Riou-Durand, Aki Vehtari, and Andrew Gelman, Nested R: Assessing the convergence of Markov chain Monte Carlo when running many short chains. Bayesian Analysis, pp. 1 28, 2024. doi: 24-BA1453.
URL https://dol http:///ooo.or///...o 1214/
Laurence Illing Midgley, Vincent Stimper, Gregor N. C. Simm, Bernhard Schölkopf, and José Miguel Hernández- Lobato. Flow annealed importance sampling bootstrap. In The Fleventh International Conference on Learning Representations, September 2022.
Aayush Mishra, Daniel Habermann, Marvin Schmitt, Stefan T. Radev, and Paul-Christian Bürkner. Robust amortized Bayesian inference with self-consistency losses unlabeled data, 2025.
Sarthak Mittal, Niels Leif Bracher, Guillaume Lajoie, Priyank Jaini, and Mareus Brubaker. Amortized In-Context Bayesian Posterior Estimation, February 2025.
Martin Modrák, Angie H. Moon, Shinyoung Kim, Paul Bürkner, Niko Huurre, Kateřina Faltejsková, Andrew Gelman, and Aki Vehtari. Simulation- based calibration checking for Bayesian computation: The choice of test quantities shapes sensitivity. Bayesian Analysis, 20(2):461-48) June 2025. ISSN 1936-0975, 1931-6690.
doi: 10.1214/23-BA. 23--33-4444
Samuel Müller, Noah Hollmann, Sebastian Pineda Arango, Josif Grabocka, and Frank Hutter. Transformers can do Bayesian inference. In International Conference on Learning Representations (ICLR), 2022.
Radford M. Neal. Slice sampling. The Annals of Statistics, 31(3):705-767, June 2003. ISSN 0090-5364, 2168-8966. doi: 10. 10.1214/a0s/1056652246
Radford M. Neal. MCMC using Hamiltonian dynamics. In Handbook of Morkov Chain Monte Carlo. Chapman and Hall/CRC 2011. ISBN 978-0-429-13850-8.
Abril-Pla Oriol, Andreani Virgile, Carroll Colin, Dong Larry, Fonnesbeck Christopher J., Kochurov Maxim, Kumar Ravin, Lao Jupeng, Luhmann Christian C., Martin Osvaldo A., Osthege Michael, Vieira Ricardo, Wiecki Thomas, and Zinkov Robert. PyMC: A modern and comprehensive probabilistic programming framework in Python. PeerJ Computer Science, 9:e1516, 2023. doi: 1..7177/peerj-cs.11//
Rafael Orozco, Ali Siahkoohi, Mathias Louboutin, and Felix J Herrmann. Aspire: iterative amortized posterior inference for Bayesian inverse problems. Inverse Problems, 41(4), 2025.
Lorenzo Pacchiardi and Ritabrata Dutta. Likelihood-ffre inference with generative neural networks via scoring rule minimization, May 2022. S11i:2205511784.
George Papamakarios, Fric 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.
Matthew D. Parno and Youssef M. Marzouk. Transport map accelerated Markov chain Monte Carlo. SIAM/ASA Journal Uncertointy Quantification, 6(2):645-688) January 2018. doi: 10.1137/17M113446000
Du Phan, Neeraj Pradhan, and Martin Jankowiak. Composable effects for flexible and accelerated probabilistic programming in numpyro. arXiv preprint or 19t:1912.1155. 2019.
Stefan T Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone, and Ullrich Köthe. Bayesflow: Learning complex stochastic models with invertible neural networks. IEEE transactions on neural networks and learning systems, 2020.
17

---

Published in Transactions on Machine Learning Research (02/2026)
Stefan T Radev, Frederik Graw, Simiao Chen, Nico T Mutters, Vanessa M Eichel, Till Bärnighausen, and Ullrich Köthe. Outbreakflow: Model-based Bayesian inference of disease outbreak dynamics with invertible neural networks and its application to the covid-19 pandemics in germany. PLoS computational biology, 2021.
Stefan T. Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ulliic Köthe, and Panl-Christian Bürkner. BayesFlow: Amortized Bayesian workflows with neural networks. Journal of Open Source Software, 8(89):5702, 2023. doi: 10.21105/joss.05702. URL https: //doi. org//0.. 21 21105 05702.
Roger Ratcliff and Gail McKoon. The diffusion decision model: Theory and data for two-choice decision tasks. Neural Computation, 20(4):873 922, 2008.
Aura. Raulo, Paul-Christian Bürkner, Jarrah Dale, English, Curt Lamberth, Josh A. Firth, and Coulson. Social and environmental transmission spread different sets of gut microbes in wild mice. bioRaiv preprint, 2023. doi: 10.1101/2023.09 27.200.4994
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In Francis Bach and David Blei (eds.), Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 1530 1538, Lille, France, 07 09 Jul 2015. PMLR. URL https://proceedingg mlr.. wtt..mm...ee//////eeeeeeee wetp///dee./eeeeeeee...... htmll.
Donald B. Rubin, Using the STR algorithm to simulate posterior distributions. In Bagesion statistics 8. Proceedings of the third Valencia international meeting, 1-5 June 1987, pp. 395 402. Clarendon Press, 1988.
Teemu Säilynoja, Paul-Christian Bürkner, and Aki Vehtari. Graphical test for discrete uniformity and its applications in goodness of-heess-offftt evaluation and multiple sample comparison. Statistics and Computing, (2):1-21, 2022.
Tim Salimans, Diederik Kingma, and Max Welling. Markov chain Monte Carlo and variational inference: Bridging the gap. In Proceedings of the 32nd International Conference on Machine Learning, pp. 1218 1226. PMLR, June 2015.
Marvin Schmitt, Paul-Christian Bürkner, and Köthe. Detecting model misspecification in amortized Bayesian inference with neural networks. Proceedings of the German Conference on Pallern Recognition (GCPR), 2023.
Marvin Schmitt, Desi R. Ivanova, Daniel Habermann, Ullrich Koethe, Paul-Christiin Bürkner, and Stefan T. Radev. Leveraging self-consistency for data-efficient amortized Bayesian inference. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarleit, and Felix Berkenkamp (eds.), Proceedings of the 41st international conference on machine learning, volume 235 of Proceedings mlr. .tree//////
of machine learning research, pp. 43723 43741. PMLR, 2024a. URL v235/ 2235//chmitte wcwsitt24a.a..m.......mm html.
Marvin Schmitt, Valentin Pratz, Ullrich Köthe, Paul-Christian Bürkner, and Stefan T. Radev. Consistency modcls for scalable and fast simulation-basod inference. Proceedings of the 38th international conference on neural information processing systems, 2024b. URL http://arxiv http::///axxvvooo/aasaeese abgg222//312.... 05440.
]loaa Schneider, Harald Kugel, Ronny Redlich, Dominik Grotegerd, Christian Bürger, Paul-Christian Bürkner, Nils Opel, Katharina Dohm, Dario Zaremba, Susanne Meinert, et al. Association serotonin transporter gene AluJb methylation with major depression, amygdala responsiveness, 5-HTTLPR /rs25531 / polymorphism, and stress. Neuropsccopparnnacolog. 43(6):1308 1113300-866,,,111111111111111111 2018. doi: 10. 2t/08//ppp201..77.
Heiko H. Schütt, Stefan Harmeling, Jakob H. Macke, and Felix Wichmann. Painfree and accurate Bayesian estimation of psychometric functions for (potentially) overdispersed data. Vision Research, 122:105-122; May 2016. ISSN 0042-6989. doi: 10.1016/j.visres.. 201.
18

---

Published in Transactions on Machine Learning Research (02/2026)
Fiona M. Seaton, David A. Robinson, Don Monteith, Inma Lebron, Paul- Christanistian Bürkner, Sam Tomlinson, Bridget A. Emmett, and Simon M. Smart. Fifty years of reduction in sulphur deposition drives recovery in soil pH and plant communities. Journal of Ecology, 111(2):464 478, 2023. doi: 10.1111/13652744.//
Louis Sharrock, Jack Simons, Song Liu, and Mark Beaumont. Sequential neural score estimation: Tikelihood- free inference with conditional score based diffusion models. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (eds.), Proceedings of the 4 1st international conference on machine learning, volume 235 of Proceedings of machine learning research, pp. 44565 44602. PMLR, 2024. URL https: : /ttp///pcceeeee...///// mlr, press/v235/ 2ttaas//223/saarrck224aa ssarrocc22aaaa... htmm..
Ali Siahkoohi, Gabrio Rizzuti, Rafael Orozco, and Felix J. Herrmann. Reliable amortized variational inference with physics-based latent distribution correction. GEOPHYSICS, 88(3), 2023.
Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International conference learning representations, 2021.
Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever. Consistency models. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), Proceedings of the 40th international conference on machine learning, volume 202 of Proceedings of machine learning research, pp. 32211 32252. PMLR, 2023. URL
mlr. wwl...pres//v2///eee https///22/s/oog23ee song23a. hong23a.html. html.
Pavel Sountsov, Colin Carroll, and Matthew D. Hoffman. Running Markov chain Monte Carlo modern hardware and software, November 2024.
Vladimir Starostin, Maximilian Dax, Alexander Gerlach, Alexander Hinderhofer, Álvaro Pejero-Canter and Frank Schreiber. Fast and reliable probabilistic reflectometry inversion with prior-amortized neural posterior estimation. Science Advances, 11(11)xead19668, March 2025. ISSN 2375-2548, doi: 10.1126/ 2//////cavv.aa/a adr9668...
Teemu Säilynoja, Marvin Schmitt, Paul-Christin Bürkner, and Aki Vehtari. Posterior SBC: Simulation-based calibration checking conditional on data. arXiv:2502.032) 0320779 2025.
Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman. Validating Bayesian inference algorithms with simulation-based calibration. aiviv preprint, 2018.
The International Brain Laboratory, Valeria Aguillon-Rodrigues Dora Angelaki, Hannah Bayer, Niccolo Bonacchi, Matteo Carandini, Fanny Cazettes, Gaelle Chapuis, Anne K Churchland, Yang Dan, Eric Dewitt, Mayo Faulkner, Hamish Forrest, Laura Haetzel, Michael Häusser, Sonja B Hofer, Fei Hu, Anup Khanal, Christopher Krasniak, Thes Laranjeira, Zachary Mainen, Guido Meijer, Nathaniel J Miska, Thomas D Mrsic-Flogel, Masayoshi Murakami, Jean-Paul Noel, Alejandro Pan- Vazquez, Cyrille Rossant, Joshua Sanders, Karolina Socha, Rebecca Terry, Anne E Urai, Hernando Vergara, Miles Wells, Christian J Wilson, Ilana B Witten, Lauren E Wool, and Anthony M Zador. Standardized and reproducible measurement of decision-making in mice. eLife, 10:c63711, May 2021. ISSN 2050-084X. doi: 10.7554 71.5554/CLafcc.3711.
Panagiotis Tsilifis and Sayan Ghosh. Inverse design under uncertainty using conditional normalizing flows Tn ATAA SCITECH 2022 Forum. American Institute of Acronautics and Astronautics, January 2022. doi: 00.251//6.202--066 URL http: / /txx/xxxxxx//// doi. orgg 1010000..1.11 2514/6. 2022-0661.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. Rank- normalization, folding, and localization: An improved R for assessing convergence of MCMC (with discussion). Bayesian Analysis, 16(2), 2021. ISSN 1936-0975. doi: bo.121 URL http://dx.dof doi. http://dx..do..og/1....oo orgg//..114//0-BBBB 1214/20-BA1221.
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.
Mischa von Krause, Stefan T. Radev, and Andreas Voss. Mental speed is high until age 60 as revealed by analysis of over a million participants. Nature Human Behaviour, 6(5):700-705, May 2022. ISSN 2397-3374. doi: 10. 10110////11562-02-01133333
19

---

Published in Transactions on Machine Learning Research (02/2026)
Yu Wang, Mikolaj Kasprzak, and Jonathan H. Huggins. A targeted accuracy diagnostic for variational approximations. In Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, pp. 8351 8372. PMLR, April 2023.
Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian Schmon. Robust neural postcrior ostimation and statistical model criticism. Advances in Neural Information Processing Systems, 35:33845 33859,, 2022.
Antoine Wehenkel, Jens Behrmann, Andrew C. Miller, Guillermo Sapiro, Ozan Sener, Marco Cuturi Camcto, and Jörn-Henrik Jacobsen. Simulatioon-bssed inference for cardiovascular models. In NeurIPS workshop, 2024. URL http:///xxxvvooog/////
13918.
George Whittle, Juliusz Ziomek, Jacob Rawling, and Michael A Osborne. Distribution transformers: Fast approximate Bayesian inference with on-thc-fly prior adaptation. arXiv preprint arXiv:2502.02463 2025.
Felix A. Wichmann and N. Jeremy Hill. The psychometric function: 1. Fitting, sampling, and goodness of fit. Perception & Psychophysics, 63(8))1293 1313, November 2001. ISSN 1532-5962. doi: 1B.3558//B001194444
Jonas Wildberger, Maximilian Dax, Simon Buchholz, Stephen Green, Jakob H Macke, and Bernhard Schölkopf. Flow matching for scalable simulationn-ased inference. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in neural information processing systems, volume 36, pp. 16837 16864, 2023.
Kaiyuan Xu, Brian Nosek, and Anthony G. Greenwald. Psychology data from the race implicit association test on the project implicit demo website. Journal of Open Psychology Data, 2014. doi: 10.5334/jopd.ac.
Jingkang Yang, Kaiyang Zhou, Yixuan Li, and Ziwei Liu. Generalized out-of-distribution detection: A survey. International Journal of Computer Vision, (32 2024.
Yuling Yao and Justin Domke. Discriminative calibration: Check bayesian computation from simulations and flexible classifier. Advances in Neural Information Processing Systems, 36:36106 -6111,, 2023.
Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Yes, but did it work?: Evaluating variational inference. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 5581 5590. PMLR, 10 15 Jul 2018.
Manzil Zahecr, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola. Deep Sets. In Aduances in Neural Information Processing Systems, volume 30. Curran Associales, Inc., 2017.
Andrew Zammit-Mangion Matthew Sainsbury-Dale, and Raphaël Huser. Neural methods for amortized inference. Annual Review of Statistics and Its Application, 12(Volume 12, 2025):311-338 2025. ISSN 2326-831X. doi: 10.1146 /annr//esstatttstiss-- 1122233-3112222221111212
Lu Zhang, Bob Carpenter, Andrew Gelman, and Aki Vehtari. Pathfinder: Parallel quasi-Newton variational inference. Journal of Machine Learning Research, 23(306):1 49, 2022. URL http://jmlz: http::///mmmo..o/aaappeee prgg/apers/ v23/21-0889. html.
Lingyi Zhou, Stefan T. Radev, William H. Oliver, Aura Obreja, Zehao Jin, and Tobias Buck. Evaluating sparse galaxy simulations via out-of-distribution detection and amortized Bayesian model comparison, October 2024.
20

---

Published in Transactions on Machine Learning Research (02/2026)
Appendix
This appendix provides additional details and analyses to complement the main text, included in the following sections:
Background, Appendix A
Best practices for training amortized estimators, Appendix B
Experiment details, Appendix C
Additional experimental study of the OOD diagnostic Step 1, Appendix D
Amortized initialization for NUTS, Appendix E
A
Background
This section provides a concise overview of the diagnostics and algorithms used in our workflow, including simulation-based calibration chocking, paramcter recovery chocking, out-of-distribution diagnostic with Mahalanobis distance, Pareto-smoothed importance sampling, the ChEES-HMC algorithm, and the nested R convergence diagnostic. Pseudocodes also given for reference.
Simulation-based calibration checking. Simulation- calibration (SBC; Talts et al., 2018; Modrák et al., 2025) principled technique for assessing the calibration of posterior distributions estimated by Bayesian inference procedures, particularly useful in simulation-based amortized inference settings. SBC is based on the idea that if the posterior p(0) y) is correctly specified, then the rank of the true parameter 0. among posterior draws should follow uniform distribution. Formally, SBC defines test statistic f e x Y + IR (e.g., a component of 0, or the log-likelihood p(y 0)). For each simulated dataset y generated from the joint model p(0, y), the test statistic is evaluated the ground-truth parameter 9(3) and compared to the same statistic evaluated over posterior samples { ~ (0) The rank of f((!)) y()) among {f(e)) { recorded. Repeating this process for all simulated datasets yields distribution of rank statistios, which should be uniform under well-calibrated inference. Deviations from uniformity signal systematic bias (e.g.; over/unne ser/under-dispersion in the posterior approximation. We the graphical approach proposed by Säilynoja et al. (2022) to assess the uniformity of the rank statistics in SBC. This method provides visual diagnostics for identifying systematic biases or miscalibrations in the posterior approximation by plotting the empirical cumulative distribution function (ECDF) and confidence bands (95%). The pseudocode of SBC is provided in Algorithm Examples of SBC checking results using this approach are provided in Appendix C.
Algorithm SBC diagnostic
Require: Joint model p(0. y); amortized posterior qo y); scalar test function f; number of simulated datasets J (e.g., 200); posterior draws S (e.g., 1000)
for I,....,J do
2: Sample 8(5) ~ p(f), y(o) p(yl p(y10(1) 0(-))
Draw g(5)) (0) y(u) for 1.....
4: Compute T(() = ((((1) y(t)), T,(d) = -(0)), و)
Compute rank p() ES (T! T! + uniform(0, Es T[T!") T())
end for
7: Compare empirical ranks أنرر against the uniform(0,S) distribution using the graphical approach of Säilynoja et al. (2022) to identify miscalibration patterns. Alternatively, the uniformity test based on the scalar statistic in Eq. 7 of Modrák et al. (2025) can also be applied.
21

---

Published in Transactions on Machine Learning Research (02/2026)
Parameter recovery checking. Parameter recovery is complementary diagnostic to SBC and provides direct visualization of posterior approximation in recovering true generative parameters (Radev al., 2020; 2023). The idea is to simulate a collection of datasets {y6}} along with their corresponding ground-truth parameters {05)) } from the joint model, and assess whether the posterior distributions (o(0) y effectively recover these known values. In our workflow, we compare the posterior median extracted from each posterior to the corresponding ground-truth values, along with the median absolute deviation to indicate uncertainty. These comparisons are visualized using scatter plots, with correlation coefficients quantifying the strength of recovery. While not a direct measure of posterior calibration or correctness, parameter recovery provides important practical insight into whether the learned inverse mapping from y to is effective. The pseudocode of the parameter recovery diagnostic is provided A Algorithm Examples of parameter recovery checking results using this approach are provided in Appendix C.
Algorithm Parameter recovery diagnostic
Require: Joint model p(0, y); amortized posterior To (0 y); number of datasets J; posterior draws S for j L,....,J do
Sample o!5) p(0), y(و) ppy| 6!
3: Draw 0(4)) ~ (0) y( for s s S.....
4: for each parameter component k do
Compute posterior median 0())
median, [0()]
Optionally compute dispersion mcasure (e.g., median absolute deviation) for (00 7: end for
8: end for
9: For each k, plot 0(:0) vs. 0l1) and report correlation
OOD diagnostic with Mahalanobis distance. The out-of-distribution (OOD) diagnostic used in Step tests whether an observed dataset Yobs falls outside the support of the prior predictive distribution (i.e., the training distribution for the amortized estimator). We work with low-dimensional summary statistics s(y) E Rd, which are either learned (e.g., via a summary network hx)) or hand-crafted with domain knowledge. In the latter case, the amortized estimator 96 must be trained using these same hand-crafted statistics. The pscudocode for computing the Mahalanobis distance and checking whether an observed dataset Yobs OOD is provided in Algorithm 3.
Algorithm OOD diagnostic with Mahalanobis distance (Step 1)
Require: Training datasets ( Jm=1 (e.g., M = 10000); summary statistic function s(-): rejection level a (e.g., a = 0.05) 1: Compute 3 (m) (m)) for all m
M
M
2: Compute empirical mean M
(m) and covariance Eg- = M
Σ (s(m)) (s) Ma)) my
m=1
3: For each m, compute Mahalanobis distance DM (a (m)) = (((lm) (a) Σs (s((m) - Ms))
Let Ta be the empirical (1 au--quattle of
5: procedure TESTOOD(ychn)
6: Compute Sobe (Yobbs and Dar(Nobs)
7: return (Mane))))
> = OOD, 0 = in-distribution
8: end procedure
Pareto-smoothed importance sampling. Pareto-smoothed importance sampling (PSTS; Vehtari et al., 2024) is robust method for improving the stability and reliability of importance sampling (IS) estimates. Given target distribution ply (0) p(0) and proposal distribution (0), with samples 0, ~ (0). PSTS
22

---

Published in Transactions on Machine Learning Research (02/2026)
stabilizes the raw importance weights wa p(y (a)) p(0s)/qs by modeling the tail behavior of the importance weights. Specifically, the distribution of extreme importance weights can be approximated by a generalized Pareto distribution (GPD):
2
k exp (
k x k == 0,
u k)
99)
where & is the shape parameter, u is the location parameter and a is the scale parameter. The number of finitc fractional moments of the importance weight distribution depends k: a generalized Pareto distribution has 1/k finite moments when k > 0. To stabilize the importance sampling estimate, the extreme importance weights are replaced with well-spaced order statistics drawn from the fitted generalized Pareto distribution, leading to a more stable and efficient IS estimator. The estimated shape parameter & serves as a diagnostic for the reliability of the importance sampling estimate. The pseudocode of PSTS provided in Algorithm
Algorithm PSIS weights and Pareto-k diagnostic (Step 2)
Require: Observod data Yobs; log- ligllikl log ply 0); log-prior log p(0); amortized posterior ( y); draws 8g qo. S.....
1: for 1, ..... do
Lg = log p(yobs 10,) + log)(0.) log (0. yobs)
3: end for
4: te = t. - max, C,
5: WN = exp(f.)
6: Sort w. obtain w(1) S .... w(S)
Choose tail size M = [min(0.28,
and define tail MS) W())
8: Fit GPD to the M largest importance weights "(S-M+1): w(s) and obtain shape estimate &
9: Replace the M largest weights with smoothed values from the fitted GPD to obtain stabilized weights w.
10: Normalize: W. = w./ ////E=. w,
11: Use 0s, W. as weighted PSIS-corrected posterior draws; treat them as reliable if & &
Given the PSIS-stabilized weights w. from Algorithm 4, one can either compute weighted Monte Carlo estimates directly (Vehtari al., 2024) or apply the SIR procedure in Algorithm obtain approximately unweighted posterior draws for downstream (e.g., visualization or MCMC initialization).
Algorithm Sampling importance resampling (STR) using PSTS weights
Require: PSIS-corrected weighted sample {(0,, (.)) desired number of resampled draws S' (e.g., S' S)
1: Define a categorical distribution on indices E {1....,S} with probabilities WI....,
2: for j=jj=1 1,..,, do
3: Sample index Ij ~ Categorical( ws) & Typically with replacement; weighted sampling without replacement useful for generating unique initializations for MCMC chains
Set 8, = 01,
end for
6: Return unweighted PSIS-corrected posterior draws 10,1
ChEES-HMC algorithm. The ChEES-HMC algorithm (Hoffman et al., 2021) is massively parallel and adaptive extension of Hamiltonian Monte Carlo (HMC) designed to leverage single-instructon multiple-data (SIMD) hardware accelerators such as GPUs. This enables rapid generation of posterior draws following an initial warm-up phase. During warm-up, ChEES-HMC adaptively tunes the trajectory length T and step size by maximizing the "Change the Estimator of the Expected Square" (ChEES), a heuristic that serves as a proxy for reducing autocorrelation the second moments of the Markov chain. ChEES-HME is particularly suitable for our amortized workflow, as we can easily generate a large number of good starting points (amortized draws) to launch many short. MCMC chains. For our experiments, we used ChERS-HMC
23

---

Published in Transactions on Machine Learning Research (02/2026)
to run 2048 parallel chains, organized into 16 superchains with 128 subchains cach. This configuration is essential for computing the nested R diagnostic (Margossian et al., 2024), which assesses convergence across large number of short chains. The pseudocode explaining the use of ChEES-HMC in Step 3 is provided in Algorithm 6.
Algorithm Use of ChEES-HMC in Step
Require: Log-posterior density log p(O, Bohs) p(yos 0) + logp(0) and its gradient w.r.t. 0; number of superchains K; number subchains per superchain M; warm-up length Nwarmup (e.g., Nearmup = 200); after warm-up sampling length N (e.g.. N =)
Collect K unique amortized posterior draws or PSIS-corrected draws for chain initialization, each of these K draws must have finite log- pog-posterior density value.
For each superchain group k = L...., K, initialize M subchains at the same initial state (total K M chains)
3: Run ChEES- HC warm-up for Narmmup iterations to adapt step size and trajectory length L
4: Fix (c, L) and run N iterations to collect draws from KM parallel chains
Compute nested R for each parameter component of 0
6: If nested R is below a chosen threshold (e.g. < 1.01), accept the combined draws as Step 3 posterior draws; otherwise increase warmup length, use alternative MCMC algorithm (e.g., NUTS-HMC) or revise the model
Nested R diagnostic. The potential scale reduction factor R (Gelman & Rubin, 1992; Vehtari et al., 2021) is arguably the most popular diagnostic for assessing the convergence of MCMC chains. The basic idea is that multiple MCMC chains starting from overdispersed initial points should produce similar Monte Carlo estimalors if they have converged, i.e. the impact of initialization vanishes as the chains converge to the stationary distribution. Nested R diagnostic (Margossian et al., 2024) extends the classical Rt diagnostic for monitoring convergence many-short-chain MCMC samplers such as the ChEES-HMC algorithm. It requires organizing chains into K superchains, each consisting of M subchains that share the same initial point. Thus, one can assess convergence through comparing the variability between superchains and the variability within superchains, similar to the standard R diagnostic. We provide the pseudocode for computing the nested R diagnostic in Algorithm 7 for reference.
Algorithm Nested R diagnostic
Require: Postcrior draws {B(mmk)) after warm-up, where (1....,K) K} (superchains), m f1...., MY (subchains), n E (1,...,N) N) (draws); scalar function of interest
Compute scalar values f((mk)) &- f(0(mmnk) for all m,
Compute subchain means: f((mk)
E==1= f(amk)
3: Compute superchain mcans f((-)) +
1 ΣM f(mm)) and overall mean f( + 1 IK f(.
M
K
4: for each superchain & 1,...., ,....,K K do
5: Compute between-chain variance Br:
M > 1, B₃ + * MM (f(-mk) 7 else Ba + 0.
Compute within-chain variance W6:
N 1, We - 3 EM [1 EN (f(mm)) - F(mblyy); else We <- 0.
8:
N-I
end for 10: Compute between-suuperchain variance: B.
K-I E 1 (f( F
1
11: Compute within-superchain variance: W. [Ke
(Bg + Ws)
12: Return nested Rt
w. A W.
24

---

Published in Transactions on Machine Learning Research (02/2026)
B Best practices for training amortized estimators
Amortized inference approaches problems of Bayesian modcling with methods from deep learning. While the precise training setup naturally depends on the concrete problem hand, some general principles have proven help across wide range of amortized inference applications. We summarize these here as initial guidance for applied practitioners.
Rely on established tooling. Modern librarics for amortized inference such as sbi (Boelts et al., 2025) and BayesFlow (Radey et al., 2023) provide well-tested neural density estimators, training loops, and data pipclines. many cases, their default architectures and optimization scttings already yicld strong performance without manual luning. First and foremost, we strongly recommend starting from these defaults and only introducing additional complexity if the diagnosties indicate deficiencies.
Monitor the training process. In many amortized inference settings, data are generated on the fly by a forward simulation program (sce Eq. 2), and training does not rely on a fixed dataset. In this case, classical data splits into training and validation set are less meaningful, since each minibatch effectively constitutes fresh data from the joint model. Nonetheless, it still important monitor training progress with multiple signals, such as the training loss, calibration diagnostics, or summary statistics of posterior samples. These checks help assess whether the model continues to improve or has already reached performance platean after short period. In other settings, amortized inference may be trained on fized set of simulated data, for example due to an expensive simulator or precomputed datasets. In such cases, holding out a validation set is strongly recommended to detect overfitting and guide selection of the amortized estimator.
Track multiple performance signals. Regardless of whether data are simulated on the fly or fixed, recommend monitoring multiple indicators during training (e.g., after each epoch). Loss curves provide coarse signal of optimization progress but are often insufficient on their own. For example, normalizing flows are usually trained with negative log-likelihood loss, which does not account for mode coverage in multi-modal posterior distributions. Complementary diagnostics such as simulation- based calibration, paramcter rocovery, or posterior prodictive checks on held-out datasots offer more direct insight into postcrior quality and failure modes.
Assess and adjust model expressiveness. When training stagnates or diagnostics indicate systematic errors, the limiting factor is often model expressiveness rather than optimization details. Underexpressive models may show symptoms such as poor parameter recovery, persistent miscalibration, posterior collapse toward the prior (i.e., data insensitivity), even when the training loss decreases. A pragmatic strategy is lo begin with an overly expressive architecture (i.e., many trainable weights) to establish a performance baseline, and then gradually simplify the model. Conversely, diagnostics remain unsatisfactory, increasing model capacity often more effective than tuning hyperparameters of the optimizer.
Accept that training is iterative. Even with modern tooling, amortized inference training may require many iterations during development, especially for complex weakly identifiable models. The objective is not to find a universally optimal neural architecture, but to reach a regime where the amortized posterior is reliable enough to enter the adaptive workflow proposed in this paper, where subsequent diagnostics and correction steps can take over.
C Experiment details
In this section, we provide additional experimental details omitted from the main text for brevity.
Evaluation metrics. To assess the quality of posterior approximations produced by each step of the workflow, we compare them against reference posterior draws obtained via a well well-tuned No-U-Turn Sampler (NUTS). Specifically, we precomputed NUTS-based posterior samples for subset of 5000 test datasets, which
25

---

Published in Transactions on Machine Learning Research (02/2026)
serve a ground-truth reference for evaluation. We then evaluate the 1-Wasserstein distance (W1) and the mean marginal total variation (MMTV) distance up to 100 datasets from each inference step: Step (amortized inference), Step 2 (amortized + PSIS), and Step 3 (ChEES-HN HMC with amortized initializations). These metrics are reported in the main text (Figure 5).
Neural network architecture for amortized inference. For all experiments, we use coupling nasslin-bbased normalizing flow implemented in BayesFlow (Radev et al., 2023). The flow consists of 6 transformation layers, each comprising invertible normalization, two affine coupling transformations, and random permutation between elements. Before entering the coupling flow network as conditioning variables, the observed dataset is encoded into a dinessiimenio summary statistic hw(y) via a summary network hy. This summary network is implemented cither as a DeepSet architecture (Zahecr et al., 2017) or a Serarasfformrr (Lee et al., 2019), depending on the problem setting. Both architectures are designed to handle permutation-invaree data structures. For the Bernoulli GLM, we bypass the summary network and directly use the known 10-dimensional sufficient statistics (Lueckmann et al., 2021). The specific choice of summary network for each application is described in the respective problem descriptions below.
Training-phase optimization. For all problems, the neural network is optimized via the AdamW optimizer (Loshchilov & Hutter, 2019) with weight decay of 10 and a cosine decay learning rate schedule (initial learning rate of 2.5 x 10-4, a warmup target of x 10-4, a = 0) as implemented in Keras (Chollet et al., 2015). A global gradient clip norm of 1.5 is applied. Training is performed with a batch size of 512 for 300 epochs, 9 with cosine decay steps set to the product of batch size and epochs. A held-out validation set is used to monitor optimization and select the best-performing model checkpoint.
Space transformation. Following standard practice in Bayesian computation (e.g., PyMC; Oriol et al., 2023, Stan; Carpenter et al., 2017), we transform parameters to an unconstruined space for inference. The amortized neural estimator is trained to estimate parameters this unconstrained space. PSTS operates independently of the parameterization and thus remains unaffected by this transformation. ChEES-HMC also performs inference in the unconstrained space. All evaluation metrics (W1 and MMTV distances) are computed in this space. However, parameter recovery and simulation- calibration plots are shown in the original constrained space for better interpretability.
Computing infrastructure and software. For all applications, the full workflow including amortized training, inference, diagnostics, Pareto-smoothed importance sampling, and ChEES-HMC sampling was conducted on single NVIDIA V100 GPU (32GB), cores of an AMD EPYC 7452 processor, and 8-16GB RAM. For runtime details across experiments, refer to Table 1 in the main text. The core code base was built using BayesFlow (Radev et al., 2023) (MIT license), Py (Oriol et al., 2023) (Apache-2.0 license), ArviZ (Kumar et al., 2019) (Apache-2.0 license) and JAX (Bradbury et al., 2018) (Apache-2.0 license). We used the implementation of ChEES-HMC provided by TensorFlow Probability (Dillon et al., 2017) (Apache-2.0 license).
Below, we provide details for each problem to complement the main text.
C.1 Generalized extreme value distribution
Problem description. Following Caprani (2021), the prior distribution for the parameters of the generalized extreme value distribution (GEV) is defined as:
# (3388,0000 0.04)
HalI-Normal(0.0009)
(10)
{ ~~ Truncated-Normaal0. 0.04) with bounds [-0.6, 0.6].
For the generalized extreme value distribution problem, 1000 reference datascts were used, corresponding the total number of test. datasets. For Bernoulli GLM, we only train for 100 epochs.
26

---

Published in Transactions on Machine Learning Research (02/2020)
M
16-
0.4
Estimate
3.0
-12
922
6.5
Ground truth
Ground tnuth
Ground bruth
(a) Parameter recovery checking.
Rant ECDP
0.0
Pank ECD 95% Confidunce Bands
Tank cco!
95% Cont dence Bands
25% Contidenco Ba
ECDP Difference
0:65
0.5
M
-0.13
0.8 Fractional rank statistc
28 Fractional rank
Fractional
(b) Simulation- based calibration checking.
Figure 7: Training-phase diagnosties for the GEV problem. The paramcter recovery is strong for the parameters a, and good for the shape parameter &. Simulation caaulation-based calibration checking indicates good calibration for all parameters. Parameter recovery and simulation-based calibration checking indicate acceptable convergence of the amortized posterior estimator.
Simulation budgets. We use 10,000 simulated parameteeroobservation pairs for training the amortized estimator, 1000 for validation, and 200 for training-phase diagnostics, including parameter recovery and simulation-based calibration.
Summary network. We use DeepSet as the summary network. The dimensionality of the learned summary statisties is 16. The DeepSet has a depth of 1, uses a mish activation, max inner pooling layers, 64 units in the equivariant and invariant modules, and 5% dropout.
Training-phase diagnostics. The closed-world diagnostics (parameter recovery and simulation-based calibration checking) in Figure 7 indicate that the neural network training has successfully converged to an acceptable posterior estimator within the scope of the training set.
Test datasets. In order to emulate distribution shifts that arise in real-world applications while preserving the controlled experimental enviromment, we simulate the "observed" datasets from joint model whose prior is 2x wider (i.e., with 4x the variance) than the model used during training. More specifically, the prior is specified as:
0(3.8,0.16)
Half-Normal (0.00)))
(11)
{ Truncated-Normal(, 0.16) with bounds [-1.2, 1.2].
C.2 Bernoulli GLM
Problem description.
Following Lueckmann et al. (2021), set the prior for
0 N~N (0, | (
(12)
27

---

Published in Transactions on Machine Learning Research (02/2026)
where the matrix F is defined such that Fi,-2 = F Fu, V and Fly = 0 otherwise,
1000 where each o2 E R100 for i 9. The task duration is set to T = 100, with fixed input vectors {wi)e
Corresponding observations are denoted by { Further details can be found Lueckmann et al. (2021); Gonçalves et (2020).
Simulation budgets. We use 10,000 simulated parameter obbervationn pairs for training the amortized estimator, 1000 for validation, and 200 for moainin--phsse diagnostics, including parameter recovery and caulatton-bbsed calibration.
Summary network. For Bernoulli GLM, the 10-dimensional sufficient summary statistic for each dataset can be computed as Vy" where y = [y,..... y/20o] and V [1i,-.. U100]. We therefore use this summary statistic for amortized training directly without relying on separate summary neural network.
Training-phase diagnostics. The closed-worid diagnostics (parameter recovery and simulation-basee calibration checking) in Figure 8 indicate that the neural network training has successfully converged to an acceptable posterior estimator within the scope of the training set.
Test datasets. We generate K 10,000 in-distribution test datasets by sampling parameters from the model prior and simulating corresponding observations {yx ! from the Bernoulli distribution.
C.3 Psychometric curve fitting
Problem description. We adopt overdispersed psychometric model (Schitt et al., 2016) with the error function (erf) as the sigmoid function in the psychometric function:
(b;;; m. w, X, 7) = + (1 -1-A8 x 7) arf(,mm,u m.
(13)
where m is the threshold, w is the width, À is the lapse rate, and y is the guess rate.
The full probabilistic model is defined as follows:
m Beta 2);
w Half-Normal(O. 1),
X, Beta(1, 10),
m = 2m - 1,
(14)
P: (2 m. Wt, X, 7);
P: ~ Beta (+--)) Pii
1 17 -1) (1 p.) (1-p))
y: ~ Binomial pt),
where n; denotes the number of trials, and 2; the stimulus level. Stimuli are presented at nine fixed levels: oi (-100.0, -25.0, 122.5, -6.25, 0.0, 6.25, 12.5, 25.0, 100.0) and each value is further normalized by dividing by 100.
Simulation budgets. We use 50,000 simulated parameter- pairs for training the amortized estimator, 1000 for validation, and 200 for training-phhe phase diagnostics, including parameter recovery and simulation-based calibration.
Summary network. We use DeepSet the summary network, which maps the input dataset to 16-dimensional summary statistic. The DeepSet has depth of 2, uses a gelu activation, mean inner pooling layers, 64 units in the equivariant and invariant modules, and 5% dropout.
Training-phasc diagnostics. The closed-worid diagnostics (parameter recovery and simulation-based calibration checking) Figure 9 indicate that the neural network training has successfully converged to an acceptable posterior estimator within the scope of the training set.
28

---

Published in Transactions on Machine Learning Research (02/2026)
F- 0.258
Fm 0.947
0.271
Fmm 0.272
r- 0.72
Extmat
rom 0.453
0.047
0.937
3,000
0.920
Ecinvee
(a) Parameter recovery checking.
65% Canl Eieeiieeo Sandd
C33
witl Carfi dinanc baunds
65% Evef deice Ronds
wi's Crefidenoo Ronde
60% Confilenco Rand5
ECOP Diffeeree
0.55
6/5
6:55
My Wn
045-
045
0.00
1//4
N
6
33% Dent Eenncc tanth
denice Zenth
95% Bonds
99% Ovniibmne Bunds
59% Cenfiencc Sands
EC3P fffieence
0.0
6:0
045-
043
M1/
0.0
c.0
Tractionll FUrn statoli
Praciional Funi stataMc
Traciional Fani Stat
Fractional rurts Jtaridic
Practional runc statbuc
(b) Simulation- calibration checking.
Figure 8: Training-phar diagnostics for the Bernoulli GLM problem. The parameter recovery strong for all parameters. Simulation-bassd calibration checking indicates good calibration for all parameters. Parameter recovery and simulation- based calibration checking indicate acceptable convergence of the amortized posterior estimator.
Test datasets. Our empirical evaluation uses 8,526 mouse behavioral datascts from the International Brain Laboratory public database (The International Brain Laboralory et al., 2021). We retrieve the data using the provided API with the argument "aat=""biasdchhoiceekrr"""" which corresponds to behavioral data collected after the mice have completed training. Each datasct processed into an observation tensor of shape (9, 3), where each row contains the number of correct trials y4, the total number of trials n,, and the stimulus level z;-
29

---

Published in Transactions on Machine Learning Research (02/2026)
15
0.422
0.580
08
02
2.5
83
20
25
64
82
82
30
8.
02
03
02
(a) Parameter recovery checking.
032
Kark EEDS
CJD
Park Frite
CL12
Poie entoe
022
Kork FC-M
Condireenn Sands
aeece Bands
wits Ronds
Conlidenco Roads
Canfilence Sands
ECOR Dfffrence
0.5
c.
6/76
ore
AAM
24
20 2n
05 D.S.
26
03
30 24 08 62
(b) Simulation- based calibration checking.
Figure 9: Training-phase diagnostics for the psychometric curve fitting problem. Recovery is good for in and while the other parameters exhibit weaker recoverability. Simulation-based calibration checking indicates excellent calibration for all parameters. Parameter recovery and simulation- based calibration checking indicate acceptable convergence of the amortized posterior estimator.
C.4 Decision model
Problem description. Following von Krause et al. (2022), we specify the prior distributions for the drift-diffusion model parameters
U1, 2/2
Gamma 2, 1),
01 02 ~ Gamma 66, 0.15),
(15)
Te ~
Gamma (3, 0.15),
Tn Gamma (3, 0.5),
where all Gamma distributions use the shape scale parametrization." 10 We implement the drift-diffusion model likelihood using the hsam package (Fengler et al., 2025) and PyMC.
Simulation budgets. We use 100,000 simulated parameter observation pairs for training the amortized estimator, 1000 for validation, and 200 for training-plie phasse diagnostics, including parameter recovery and simulation-based calibration.
Summary network. We use Set'Transformer as the summary network, which maps the input dataset to 16-dimensional summary statistic. The SefTransformer has two set attention blocks, followed by a pooling multi-head atlention block and a fully connected output layer. Each multilayer perceptron (MLP) in the set. blocks has two hidden layers of width 128, with gelu activation and 5% dropout
Training-phase diagnostics. The closed-world diagnostics (parameter recovery and simulation based calibration checking) in Figure 10 indicate that the neural network training has succossfully converged to an acceptable posterior estimator within the scope of the training
Test datasets. The test datasets consist of 15,000 participants pre-processed from the online implicit association test (IAT) database (Xu et al., 2014; Krause et al., 2022). Each test dataset is tensor of
10The prior distributions for the boundary separation parameters 8o ø2 differ slightly from those Krause al. (2022) due different parametrization of boundary separation.
30

---

Published in Transactions on Machine Learning Research (02/2026)
Grandouth
from
Grandooth
Crund
Grand
Parameter recovery checking.
M// Mv
w/ New
(b) Simulation- calibration checking.
Figure 10: thiining-paase diagnostics for the decision model. Parameter recovery is strong for all parameters. Simulation-base calibration chocking indicates good calibration for all parameters except Tn, which shows mild deviations, suggesting occasional overestimation by the amortized estimator for this parameter. Parameter recovery and simulation-basel calibration checking indicate acceptable convergence of the amortized posterior estimator.
shape (120, 4), where each row corresponds to a single trial and contains the response time, missing data mask, experiment condition type, and stimulus type.
D Additional experimental study of the OOD diagnostic in Step 1
To further investigate the relationship between the Mahalanobis distance the OOD diagnostic and the quality of the amortized posterior, we visualize this relationship using scatter plots in Figure 1la for the four tasks considered in the main text. For each task, we use around 1000 test datasets and compute the Pearson correlation coefficient r. The Mahalanobis distance is positively correlated with the two posterior quality metrics (W1 and MMTV) for the GEV, psychometric curve, and decision model tasks, where out- of- di-distribution test datasets are present. For the Bernoulli GLM, the correlation is negative; here, all test datasets were generated from the same distribution (prior simulations) as the training datasets and the Mahalanobis distance is not informative. From Figure 11a, we see a key limitation the Stop-1 OOD diagnostic: the Mahalauobis distance is clearly not a perfect proxy for the posterior quality. In particular, the amortized estimator may still yield low-quality posterior draws on dataset with smaller Mahalanobis distance, also observed in Figure 5.
We next check the impact of the threshold a in the Step-1 OOD diagnostic by varying it from 0.01 to 0.5, specifically over the set [0.01, 0.05, 0.1,0.15, 0.2, 0.3, 0.4, 0.5], as shown in Figure 11b. As a increases, more test datasets are rejected, and the overall posterior quality of accepted amortized posterior draws generally improves as measured by the median and IQR of the posterior metrics (lower posterior metric values indicate higher quality). The quality of rejected amortized posterior draws also improves as a increases, while remaining consistently worse than that of the accepted amortized draws. 11 Overall, the posterior quality of the accepted amortized draws, in terms of median and IQR of W1 and MMTV, is not very sensitive to the threshold a, and 0.05 appears to be reasonable default choice.
These results support the use of the hihaalanobi-ddittncc--bbe in Step 1: it lends to flag the most. problematic dalasets and thereby improves the quality of accepted amortized posterior draws at negligible additional cost. At the same time, the residual low-quality posteriors at
OOD test as a lightweight first-linc diagnostic
For Bernoulli GLM, where test datascts and training datascts come from the same distribution, we observe slightly reverse trend as increases. This consistent with the corresponding result shown in Figure 1la and further suggests that the Mahalanobis distance good measure posterior quality for in-distribution datasets this case.
31

---

Published in Transactions on Machine Learning Research (02/2026)
GEV
Bernoulli GLM
Psychometric curve
Decision model
5.0
0.47
7 == -0.46
r 0.37
r 0.42
W1
2.5
0.0
0
10
20
4
10
0
10
1.0
0.64
%
r = -0.13
T' = 0.52
F = 0.35
MMTV
0.8
0.5
0.2
0.0
10
0
Mahalanobis distance
(a) Posterior quality metrics (W1 and MMTV) versus Mahalanobis distance for the four benchmark tasks.
GEV
Bernoulli GLM
Psychometric curve
Decision model
1.0
W7
0.5
0.0
1.0
0.8
MMTV
0.5
0.2
0.0
0.01 0.2 0.3
0.5 0.010.1 0.2 0.3 0.4
0.01 0.1 0.2
0.3
0.4
0.5
0.01 0.1
0.3 0.4 0.5
Threshold
Accepted
Rejected
(b) Sensitivity of the rejection threshold a in the OOD test. The median + IQR (shaded area) of the posterior quality metrics is shown separately for accepted and rejected datasets.
Figure 11: Relationship between amortized posterior quality metrics, Mahalanobis distance, and the OOD rejection threshold in Step 1. (a) Scatter plots of W1 and MMTV versus Mahalanobis distance, with Pearson correlation coefficient r reported in cach panel. (b) Effeec of varying the threshold a on the postcrior quality of accepted and rejected amortized posterior draws. See text in Appendix D for details.
small Mahalanobis distances underscore that this diagnostic cannot guarantee accuracy. For applications that require tighter accuracy guarantecs, it is therefore natural to enforce escalation to Step 2 (PSIS) irrespective of the OOD outcome, trading additional computation for more robust posterior approximation.
32

---

Published in Transactions on Machine Learning Research (02/2026)
E Amortized initialization for NUTS
In addition to ChEES-HMC, we evaluate the effectiveness of amortized posterior draws as initializations for the NUTS sampler. The experimental settings mirror those used for ChEES-HMC (Section 3.4), except that we launch only four chains, which is the typical configuration for NUTS. As shown in Figure 12, amortized initializations reduce the number of required warm-up iterations for both the GEV problem and the decision model. For the psychometric curve and Bernoulli GLM problems, all three initialization methods (amortized, PSIS-refined, and random) yield similar convergence behavior according to the R diagnostic (Vehtari et. al., 2021).
Notably, NUTS generally requires fewer warm-up iterations than ChEES- HMC across the evaluated problems, suggesting that while amortized initializations are still beneficial, the relative gain is more pronounced for ChEES-HMC, which runs many short chains in parallel.
GEV
Bernoulli GLM
10-1
10-1
10-2
10-2
10-3
10-3
Amortized
Psychometric curve
Decision model
Amortized PSTS
10
10-1
Random Init.
10
10-2
(02)
10
10-3
1050 100
200
300
500
1050
100
200
300
500
Warmup iterations
Warmup iterations
The effect of initialization for NUTS. The figure shows median+IIR across 20 test datasets. Using amortized posterior draws initializations for NUTS reduces the required warmup in the GEV and decision model tasks.
Figure 12:
33


Reviews:
[{"review_id": "nodrZebI2E", "reviewer": "Reviewer_KzLH", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces a three-stage amortized Bayesian workflow designed to make posterior inference scalable, while maintaining reliability. It transitions between pure amortized inference (fast but potentially unreliable) and MCMC (accurate but slow) trying to stay on a speed/accuracy \"pareto front\" of inference methods.\n\nThe workflow first tries pure amortized variational inference by training an affine coupling normalizing flow and doing posterior sampling. If that fails diagnostics, inference is triaged to an importance sampling step. Finally, if the IS step fails diagnostics inference is triaged to a full MCMC, initializing to the amortized inference samples.\n\nStrength:\n* Principled, diagnostic-driven pipeline\n* Reusing computations at each step\n* Modularity of pipelines, e.g. replacing type of neural density estimator in step 1 or alternative MCMC method in step 3\n\nWeaknesses (further discussed in requested changes below):\n* Hyperparameter choice for neural inference training is glossed over. How can we reliably set these?\n* Unclear what is claimed as contributions\n* Unclear how this work relates to other work", "claims_and_evidence_yn": "No", "claims_explanation": "For the models tested, the amortized Bayesian workflow works in practice in that it has significant time speedup with most instances using only amortized neural inference, fewer triaging to PSIS, and fewer still triaging to MCMC. They also provide clear evidence that initializing from the amortized inference samples speeds up MCMC convergence. \n\nOn the other hand, as the authors acknowledge: \"while our diagnostic for the amortized posterior draws in Step 1 is effective and highly efficient in practice, it remains an imperfect proxy for the true posterior approximation error and can occasionally result in the acceptance of poor-quality amortized draws\"\n\nAlso, it seems like the problems/models are chosen to be ones that it's expected amortized inference and PSIS work well. There's limited insight into why the workflow fails on a given model/dataset, etc. Overall, the experimental setup seems geared to artificially boost claims of a general purpose black-box inference workflow.", "audience_yn": "Yes", "audience_explanation": "Yes, many TMLR readers would be interested in methods / experiments surrounding more flexible, general purpose, robust, automated Bayesian modeling and inference.", "requested_changes": "Some changes that would be critical to securing my recommendation for acceptance:\n\n* Make clear in the first section what the claimed contributions are. Is there technical novelty in the triaging with diagnostics process, is it the empirical evidence for the inference speedup, etc.?\n\n* Be clearer on limitations of the workflow and what's shown in the paper. Does is apply to a narrow class of Bayesian models? What happens for high-dimensional models, discrete latent variables, etc.? Can you provide any deeper insight into the cases where the workflow failed in the experiments section and model classes where the workflow fails? Experiments on a wider class of models would be beneficial.\n\n* Make clear why the problem you're solving is an important one. Is computational speed still a bottleneck given modern compute? And is this really the case on the size of the models given in the experiments? Could reliable blackbox inference be a bigger bottleneck, as both neural VI and MCMC can fail on different models/datasets.\n\n* More depth in \"Related Work\" section. Many methods and areas are mentioned like a bucket list but it isn't clear how they relate to the workflow presented in the paper. I think it would strengthen the paper to discuss more related work on automated Bayesian inference systems (e.g. probabilistic programming libraries) and how the approach they take differs from the one given here. This obviously isn't the first work to propose a general purpose Bayesian inference workflow combining amortized VI with MCMC.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "eQAnxOxf08", "reviewer": "Reviewer_EoXK", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a principled workflow for Bayesian inference in the context of several data sets of observations. The workflow prioritizes inference speed, starting by training an amortized variational approximation and sequentially moving towards slower but more theoretically-grounded methods of inference (PSIS-SIR and ChEES-HMC, respectively). To this end, this work introduces multiple practical diagnostics to decide, for each observed data set, whether to stop the inferencial procedure and retain the generated samples or to continue to refine the samples via a slower inference algorithm.", "claims_and_evidence_yn": "Yes", "claims_explanation": "Yes. The experiments section shows clear speedups in inference without sacrificing (much) sampling quality, across several different experiments and with two different metrics. \n\nI really appreciate the inclusion of Section 3.4 and Figure 5. It provides solid evidence that reusing computation can really help (and at worst it will most likely not hurt).", "audience_yn": "Yes", "audience_explanation": "Yes. Practitioners with access to several data sets but limited computation specifically can benefit from potential orders of magnitude in speed up.", "requested_changes": "My main concerns are twofold:\n\n- the workflow could benefit from more clarity;\n- and it is unclear how some of the proposed diagnostics affect downstream sampling quality.\n\nFor the first point, see the additional comments. I highlighted multiple points on which I would like an expanded discussion or clarification thereof. In summary, while I understand that the paper might want to present the proposed workflow as succinctly as possible, I believe that including Appendices A and B in the main text (each paragraph in its corresponding subsection) would actually clarify the details of the workflow.\n\nFor the second point, I believe that the paper would greatly benefit from an ablation study of the user-specified parameters in the diagnostics to go from training phase to Step 1, and from Step 1 to Step 2. I expand on some of these details in the additional comments, but succinctly it’s unclear to me how the choices of the test statistic (in the training phase), the summary statistic for Mahalanobis distances (in Step 1), and the predefined acceptance threshold $1-\\alpha$ (also in Step 1) affect the quality of the resulting samples. I understand that some of these diagnostics might have been studied in the literature already (like the Pareto-$\\hat k$ of Step 2). In that case, instead of an ablation study, a discussion might be suitable (e.g., something like the discussion of using alternative diagnostics in Appendix B.2).", "broader_impact_concerns": "None.", "additional_comments": "The OOD check in Step 1 aims to detect atypical observed data sets. I have a few concerns with this approach.\n\n- Atypicality is measured only through a summary statistic, but it’s unclear how different choices of summary statistic affect the flagging of a data set as atypical. In other words, how robust are these summary statistics to atypicality detection?\n- This is the only check used to determine whether samples from the amortized approximation are accepted or sent to PSIS-SIR. Yet it’s unclear to me, from reading the paper, whether this is close to a sufficient (with high probability) condition for a data set to warrant its amortized draws to be rejected. That is, can a data set be considered typical (thus its draws from $q_\\phi$ be accepted) yet still lead to poor posterior draws? Figure 5 seems to suggest so: there are many outliers in the accepted data sets that compare poorly to the ground truth, particularly in the GEV experiment. Are there any comments on this on top of the discussion in Appendix B.2?\n- Figure 3 suggests that the training samples are much more concentrated around small values of Mahalanobis distances when compared to the test samples. Can you say anything about this in the context of the data sets in the experiments? E.g., from Table 1 I see that the GEV experiment resulted in almost half of the data sets getting flagged as atypical, which suggests to me that the test-based Mahalanobis distances where much larger than the training ones. Any intuition for this? Where I’m coming from: as a practitioner, I would like to understand how this part of the workflow is related to, e.g., the complexity of my posterior so that I can maybe adjust my predefined percentile threshold $1-\\alpha$ accordingly\n- I believe that expanding the discussion in Section 2.2.1 would make it much more useful for potential users to understand the OOD diagnostic as well as the trade-offs of using different diagnostics. Specifically, I recommend including Appendix B in the main text. Appendix B.1 provides much needed details about the OOD diagnostic, while the discussion in Appendix B.2 regarding alternative diagnostics is really important in my opinion.\n\nFrom Appendix A, it seems like the assessment of uniformity of the ranks of the test statistic is done visually? (E.g., Fig. 7 in Appendix C.) Is this visualization just for presentation purposes? Ideally you would want to automate checking uniformity (e.g., via a Kolmogorov-Smirnov test) rather than relying on a visual check. Please let me know if I’m misunderstanding.\n\nI think the paper would benefit from having a pseudo-code algorithm version of PSIS + SIR. Something like:\n- In: $(\\theta_1,\\dots,\\theta_S)$\n- [Pseudo-code steps] \n- Out: $\\hat k$ and $(\\theta_1’,\\dots, \\theta_S’)$\n- And then if $\\hat k >$ some threshold then I accept $\\theta’$, otherwise I move on to Step 3", "extra_scores": {}, "extra_text": {}}, {"review_id": "xj4AGxRQej", "reviewer": "Reviewer_8xGs", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a Bayesian inference workflow starting from a trained amortized inference artifact. Each stage of the workflow is coupled with a diagnostic indicating sufficient inference quality, or the need for continued refinement of the inference result. First the artifact training is verified with simulation-based calibration checking. In the first stage, at inference time, samples from the artifact are checked against an out-of-distribution (OOD) diagnostic. If the diagnostic is not satisfied the second stage corrects the posterior samples with importance sampling. The quality of the importance samples is verified with a Pareto distribution fit for the importance weights. If the diagnostic indicates poor sampling quality the posterior samples are finally corrected for in the third stage through a more expensive, parallel MCMC sampling scheme, ChEES-HMC. Using the amortized samples as initializations for the MCMC chains is one contribution of the paper. The final diagnostic is the $\\hat{R}$ statistic, which on failure triggers a fallback to the Standard Bayesian Workflow (Gelman et al. 2020). \n\nWhile each individual piece of this pipeline is well-established, the paper's contribution is the integration into a pipeline with solid default choices for each diagnostic. Whenever a diagnostic indicates sufficient posterior sample quality it can skip the more expensive follow up phases, aiming at an optimal trade-off between computation and sufficient posterior sample quality. Its choices are intuitive and easy to understand. While the combination of amortized inference with importance sampling or MCMC is well-established,  the paper provides a sound strategy for practitioners to follow. The paper experimentally validates the workflow on four problems from the Bayesian inference literature and demonstrates that each stage in the pipeline can improve inference results.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The stages of the workflow pipeline are well-established and convincingly composed. Each diagnostic is chosen from a well-defined statistic and pointers are given to alternative diagnostics where applicable. The experiments are satisfying, but could be more comprehensive both in size and number. Nonetheless, they are sufficient to demonstrate the workflow's applicability.", "audience_yn": "Yes", "audience_explanation": "Amortized inference is an emerging trend in Bayesian inference, but not yet widely integrated in many more traditional applications and probabilistic inference systems. Approaching the integration on the level of a workflow with guarantees expected in established Bayesian inference applications delivers the speedups achievable through amortized inference to the whole community. While many machine learning papers exist that describe how to train and check amortized inference artifacts in specific settings, few provide such general and practical workflow suggestions.", "requested_changes": "I think the paper is good as it is, but I find a few parts could be improved:\n\n1. The paper makes claims to move along a Pareto front between minimal compute and optimal inference quality. This is only touched upon in the beginning of the paper, but I think characterizing the trade-offs involved more carefully would be beneficial, since future refinements of the workflow will very likely be necessary. For instance there is a difference between inference speed and cost, and in many cases inference cost is more important than speed. Similarly, sometimes different types of hardware, e.g. large GPUs, might or might not be available. Considering only the speed dimension used in the paper, other variants of the workflow (or other workflows) could very well be better. For instance one variant that would be interesting to explore would be to adaptively increase sample size from the artifact if amortization + importance sampling are not sufficient, instead of using MCMC. Following Glöckler et al. 2022 [1], one could also deploy iterative refinement of the amortized artifact and skip MCMC completely, which is likely going to be faster than ChEES-HMC in many cases. Such trade-offs are touched upon a little bit in \"Modularity and practical flexibility\" in the Discussion Section 4, but could be significantly expanded.\n\n2. I think making each diagnostic calculation very explicit in the paper, such that it could be right away implemented in any context, could simplify the workflow suggestions. Some of this information is provided in Appendix A, but since I am not familiar with the details of each diagnostic, I would have to search the literature now to implement the workflow myself, which might slow the adoption of the workflow. In general, I think recipe-like papers are a great way to accelerate adoption.\n\n3. It would be interesting to account for training time of the amortization artifact separately in Table 1 as well.\n\n4. The paper more or less assumes that an amortized inference artifact is already available and can be trained. Providing a bit more context into how this can be done might prove very helpful for implementers not versed in training deep models for amortized inference. In fact this could prove to be the most critical part of the pipeline. Section 2.1 only touches on the inference objective and creation of synthetic data, but says little about the deep learning aspects of the procedure.\n\n[1] https://arxiv.org/abs/2203.04176", "broader_impact_concerns": "-", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)