Compositional simulation-based inference for time series
Authors: Manuel Gloeckler
Shoji Toyota
Kenji Fukumizu
Jakob H. Macke
Venue: ICLR 2025 Poster (2025)
Keywords: Simulation-based inference
Bayesian inference
time series
markovian simulators
Amortized Bayesian inference
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: Simulation-based inference for Markovian simulators leveraging the factorization

Abstract: Amortized simulation-based inference (SBI) methods train neural networks on simulated data to perform Bayesian inference. While this strategy avoids the need for tractable likelihoods, it often requires a large number of simulations and has been challenging to scale to time series data. Scientific simulators frequently emulate real-world dynamics through thousands of single-state transitions over time. We propose an SBI approach that can exploit such Markovian simulators by locally identifying parameters consistent with individual state transitions. We then compose these local results to obtain a posterior over parameters that align with the entire time series observation. We focus on applying this approach to neural posterior score estimation but also show how it can be applied, e.g., to neural likelihood (ratio) estimation. We demonstrate that our approach is more simulation-efficient than directly estimating the global posterior on several synthetic benchmark tasks and simulators used in ecology and epidemiology. Finally, we validate scalability and simulation efficiency of our approach by applying it to a high-dimensional Kolmogorov flow simulator with around one million data dimensions.

Published as a conference paper at ICLR 2025
COMPOSITIONAL FOR TIME SERIES
SIMULATION-BASED
INFERENCE
Manuel Gloeckler
Shoji Toyota"
University of Tübingen Tübingen, Germany manuel.
Kyushu University
Fukuoka, Japan
de
boyota@ait..yushuuaoao ip
Kenji Fukumizu
Jakob H. Macke
The Institute of Statistical Mathematics Tokyo, Japan
University of Tübingen & MPI-IS Tübingen
Tübingen, Germany jakob. macke.ic -treii ngen..de
fukumi zuê sm. a ip
ABSTRACT
Amortized simulation- baauee inference (SBI) methods train neural networks on simulated data to perform Bayesian inference. While this strategy avoids the need for tractable likelihoods, it often requires large number of simulations and has been challenging to scale time series data. Scientific simulators frequently em- ulate real-world dynamics through thousands of single-state transitions over time. We propose an SBI approach that exploit such Markovian simulators by lo- cally identifying parameters consistent with individual state transitions. We then compose these local results to obtain posterior over parameters that align with the entire time series observation. We focus on applying this approach to neural poste- rior score estimation but also show how it can be applied, e.g., to neural likelihood (ratio) estimation. We demonstrate that our approach is more cntulltion-efficinn than directly estimating the global posterior on several synthetic benchmark tasks and simulators used in ecology and epidemiology. Finally, we validate scalability and simulation efficiency of our approach by applying it to nigh-dimensional Kolmogorov flow simulator with around one million data dimensions.
INTRODUCTION
Numerical simulations are central approach for tackling problems in wide range of scientific and engineering disciplines, including physics (Brehmer Cranmer 2022; Dax et al. 202T), molecular dynamics (Hollingsworth & Dror 2018), neuroscience (Gonçalves ves et al. 2020 and climate sci- ence (Watson- Parris et al., (021)) Simulators often include at least some parameters that cannot be measured experimentally. Inferring such parameters from observed data a fundamental challenge. Bayesian inference provides a principled approach to identifying parameters that align with empiri- cal observations (Gelman et al. 20200.. Standard algorithms for Bayesian inference, such as Markov Chain Monte Carlo (MCMC) GGiks et al. 1995) and variational inference (Beal 2203). generally require access the likelihoods 0)x= However, for many simulators, directly evaluating the likelihood remains intractable, rendering conventional Bayesian approaches inapplicable. Yet, gen- erating synthetic data x p(z[0) feasible for numerical simulators. Simulation-base based inference (SBI) methods offer powerful alternative to perform Bayesian inference for such simulator models with intractable likelihoods (Cranmer et al. 2020).
Classical SBI methods, like Approximate Bayesian Computation (ABC) (Beaumont et al. 2002) and synthetic likelihoods (Wood) 20010), struggle to effectively scale to high- high-dimensionall simula- tions. To address this, SBI methods have been developed, which train neural networks to represent likelihoods (Papamakarios et al. 2019 Glöckler et al.] 2022 Boelts et al.) 2022). likelihood ra-
tios (Durkan al. [2020] Hermans et al.) 2020; 2022 Miller et al. (022). posteriors (Papamakarios & Murray 2016 Lueckmann et al. 2017
Greenberg et al. 2019 Deistler et al. 2022 Geffner et al.
"Equal Contributions

---

Published as a conference paper at ICLR 2025
10
true
approx.
Score accumulation
x11
10'0"
Inference:
Full time series
Training: Single transiioos
15
T=
15
1x. B]
(elo//xttiig Va, ((9a(( (1e)
Vo, log (0
Prior
Proposal and Transition
Local approximation
Global score spproximation
Fig 1: Illustration of Factorized Neural Score Estimation (FNSE). The goal is to perform param- eter inference on full time series model. The training process uses smaller subsets of single-stan state transitions initialized arbitrary proposal p(x'), with parameters sampled from prior distribution. During inference, the time series is divided into single state transitions, and each state transition is evaluated by the neural network to estimate local posterior scores. These local estimates are then aggregated to form global approximation, which is subsequently used to sample from the overall posterior distribution. Here, a denotes the diffusion time, and 0g is the associated noisy parameter.
Wildberger et al. 2023 Sharrock et al. 2024) or target several properties at once (Radev al.) 2023 Gloeckler et al.) 2024b). These methods allow for parameter inference without requiring ad- ditional simulations after training, making them more efficient than traditional approaches (Lueck- mann et al.) [202]]: they effectively amortize the cost of the simulation, and/or the full inference approach (Gershman & Goodman 2014 Le et al. 2017).
2023
However, applying these neural amortized methods to time series simulations can be challenging due to the high computational cost of repeated simulator calls. Running numerous time series posssill with varying sequence lengths- least wasteful. Instead, it seems advantageous to exploit the temporal structure these simulators: Many-if not most-scientific simulators for time series data are based on (stochastic) differential equations that model dynamics a system through state transitions over time (..g.. Se..SSedmeier et al.) (2016); Goswami 20swaaGo et al. (2006); SSrauss Effenberger (2017); [Siur et al. (2016)). directly model processes which are iteratively updated at each time-step. As consequence, such simulators have an inherently Markovian structure, which can be leveraged for efficient inference! particular, we can approach the problem on the level of single-state transitions, a simpler task that can be tackled with fewer simulations. For instance, in the Kolmogorov Flow simulator (Sec. [4.4). the simulation outputs reach nearly one million dimensions for time series of length hundred, which is out of range for most direct SBI approaches. In contrast, the dimensionality of single state transitions is only few thousand.
be computationally prohibitive at
In this work, we propose a method for efficient simulation-bel based inference for Markovian time series simulators. Unlike other neural SBI methods that require long sequences for training, approach locally estimates parameters consistent with a single state transition. By aggregating these local solutions, recover global posterior approximation for time series of arbitrary length, significantly reducing simulation costs (Fig. []. We focus on applying our strategy to posterior score-based approaches (FNSE, Factorized Neural Score Estimation) and additionally apply it to likelihood --ratoo) estimation methods (FNLE and FNRE). A related challenge has been addressed in prior studies (Geffner et al. 2023 Linhart al. 2024 Boelts al. [022] which had focused on efficient inference of i.i.d. (independently and identically distributed) observations. We here show that these approaches can be extended to Markovian simulators simulators-arguablly the dominant model- class for mechanistic simulators for time-series and empirically evaluate our approach on series of synthetic benchmark tasks, as well as simulators from ecology and epidemiology. Finally, we demonstrate the scalability and efficiency of this approach the high-dimensional Kolmogorov flow problem.
2 PRELIMINARIES
2.1 PROBLEM SETTING
We target simulators p p( () that output time series x (......<>) E Txd, Many simu- lators in scientific fields are implemented using (stochastic) differential equations, which inherently satisfy the Markov property. We focus on this class of simulators, where the likelihood takes the
2

---

Published as a conference paper at ICLR 2025
following form
p(a 10) p(x000)
II pl plaxt+1
++ lar 0).
To simplify the exposition, we constrain the initial condition to be independent of the parameter, p(²²10) = =p(x=), and assume that the simulation transition p(x"+1|=* ,0) does not depend on time 6. We relax both of these assumptions in Appendix [A.2] and also address the general case of higher-order Markov chains in Appendix A.22.3 We assume one can sample from the state-transition
function x1 [xx"+1 | 0), but not evaluate its probability density.
2.2
AMORTIZED METHODS FOR SIMULATION- -AASSE INFERENCE
Neural Posterior Estimation (NPE): Neural Posterior Estimation (NPE: Papamakarios & Murray
Lueckmann et al. (2017) Greenberg et al. (019);; Wildberger et al. (2023)) methods typ- ically use a conditional neural density estimator, e.g.. normalizing flow (Rezende & Mohamed,
(2016):
Papamakarios et al. 2017: Durkan et al.) 2019) to approximate the target posterior distribu- is trained via maximum likelihood to estimate p(0) (x:01 from a dataset .r: where (0., x:::) ) p(0 p(xx::T()) For observations of inhomogeneous size, NPE often incorporates an embedding network. Embedding networks can be designed as recurrent networks (Lueckmann et al. 2017), permutation-innarian architectures (Chan) ct al. 2018 Radev et al.) 2020), or attention- models (Schmitt et al., 2023) to respect known invariances in the data. At inference time, NPE yields amortized approximation of the posterior for any observation x 0:T therefore offering rapid evaluations or samples of the posterior.
2015
tion. The model 4o [0(z) xz
of parameter-data pairs
Neural likelihood(-ratio) estimation (NLE, NRE): Neural Likelihood( Ratio) Estimation
(NLE (Papamakarios et al.) 2019), NRE (Hermans et al. 2020 2022 Miller et al., 2022) meth-
ods do not directly approximate the posterior distribution but instead train a surrogate Io (a: (0)
for the likelihood (T () or likelihood ratio p( ax..:TTpp ( Once the likelihood is accessi- ble via the surrogate, standard techniques such as MCMC or variational inference can be applied to perform inference. These approaches circumvent the problem of not having access to the likelihood, but still leaves us with the drawbacks of standard inference methods, i.e., slow sampling and po- tential failure modes such as robustly handling multimodality (Glöckler al., 2022). Nonetheless, they also inherit the flexibility of standard inference techniques, such as the ability to handle mul- tiple independent (or conditionally independent) observations by simply multiplying the likelihood terms (Boelts et al. 2022).
Neural score estimation (NSE): Score-based diffusion models are a powerful tool for generating samples from a target distribution. Common diffusion models (Song & Ermon 2019 Song et al. 2021 Ho et al. 2020) are based on stochastic differential equations (SDEs) that can be expressed as dl. f(a) 22.00.dd.. da + g(a)dul) with being standard Wiener process. The drift and diffu- sion coefficients ff and g are chosen such that the solution to this SDE defines a diffusion process that transforms any initial distribution into simple noise distribution PA (0x) Samples from any target, such as the posterior (:T), can then be obtained by simulating the reverse diffusion process (Anderson) 1982)
(O M.A, o
de, = {f(a) - (a) Ve, logp((..). x )T) da + g(a)dk
where v is backward-in -inttr---i--eeeeeee Wiener process. In practice, we not access the analytic form of the score s (x:T) Vo, logp(0 le 6:1 but we can estimate it from samples using conditional
denoising score-matching (Hyvärinen & Dayan 2005 [ong et al. 202T)
L(Φ) E,, [>(a)]
(O)))(s((00())))))) (0.|æ"
Vo, log pigp(O... 0)
where À denotes a positive weighting function and p(0,10) = NN(0,; ((a))
This recently proposed approach has been highly successful across various tasks (Geffner et al. 2023 (Shrrock et al., [2024] Gloeekker et al. (024b).. offers a trade-off between the efficiency of the static NPE method and the flexibility of slower but more flexible NL(R)E method. By enabling
feasible post-hoc modifications and composability through appropriate adjustments to the backward
diffusion process, it bridges the gap between these two approaches (Geffner et al. 2023 Gloeckler
let al. 2024b).
'The diffusion time is denoted by a, against convention, distinguish from & in the time series.
3

---

Published as a conference paper at ICLR 2025
Algorithm Training
Algorithm 2 Evaluation
1: Input: prior p(0), proposal p(x²), tran- sition function T(x'* < 0), score net
1: Input: prior p(0), observation x??
compose method (see Sec. 3.2.2
$6.
2:
2: D
// Generate training dataset
3: def Sglob (0,, x ):
3: for i = 1 to N do
4: Slocal []
4: 0: ~ p(); x. ~ -(x"))
5: for t = 0 to T do
5: a: ~ T(x c 0.)
6: Stocal += [s x=+++11]]
6: DU {10,, o x )
7: Sglob =compose p p())
7: end for
8: return Syloh
8: Train Sy by minimizing Eq. using D
9:
Sample p(0)x..T) via Sylob
METHODS
3.1 GLOBAL INFERENCE FROM SINGLE-STEP TRANSITIONS
Direct estimation of the global target distribution (i.e., the posterior p(0 T.) in NPE or the like- lihood p(x0 () in NLE; Sec. 2.2 often requires a large number of simulations of the entire time series, leading to intractable, or at least very expensive, estimation problems. We mitigate the prob- lem of expensive simulation calls by leveraging the Markov factorization of the forward model: Since the simulator is completely specified by its state-transition probabilities, p(x*+' 0), these also contain all relevant information about the parameters. global target distribution (e.g.. the posterior) for a time series of variable size given a dataset of single-step transitions. Building on this observation, we aim to estimate local target scores 8tocal using single-step transition simulation data D {(0,, x; )] (Alg. | Sec. 3.2.1). Afterwards, we recover the global target distribution by aggregating these local solutions using a composition rule (compose, Alg. 2 Sec. 3.2.2). In the following, we apply this approach to NSE, NLE, and
therefore possible recover the
NRE (Sec. 3.2 3.3) and discuss the design of the required proposal p(xx)) (Sec. 3.4 for sampling single-state transitions.
3.2 FACTORIZED NEURAL SCORE ESTIMATION (FNSE)
3.2.1 LOCAL SCORE ESTIMATION
Assuming that p(x ()) satisfies Eq. in NSE can be factorized as
and p(*" 0) = =((x)), the global target Vo log p(0 (00)pp00)) xe ":T)
T-1 Σ 1=U
Ve logp(0)= 0.T
Vo log p( x x - (T - 1). Vo log (().
(3)
Here, p(0)2² 1="", x+1) is a posterior with x' following any proposal distribution x² ~ f(xx) and "+1
following the state transition x++1 T(x*+1/0' 0) (Appendix Sec. A. The factorization (3) implies that the global posterior is fully characterized by s(0) [ x²
Vo log (00 a'+
We can estimate this quantity using only single-state transitions by minimizing the loss
(()) = E.,,. 0., , [A(a)]lso(O. |æ", x(+1)-- Ve. log (ogp(0,(0)11)
(4)
given a proposal distribution x' p(x'). As a result, we can learn local score estimator by empirically minimizing this loss given a dataset of single-step simulations (Alg. D For globally amortized posterior approximation, the proposal distribution must at least satisfy two properties: (i) must have support at any point x' the simulator can attain, (ii) the proposed point x' must be independent of the parameter 0 involved in the transition to (+11 to ensure Eq. 3 remains valid (Appendix Sec. A. These constraints leave quite some flexibility in the choice of the proposal. Its choice will have significant impact, given finite simulation budget (Sec. (3.4) Our approach can readily be extended to higher-order Markov chains, which require proposing multiple past ob- servations (Appendix Sec. A.2.3. Between the two extreme cases of learning from single-state transitions, or full time series, there is a continuum of approaches, similarly to what was shown in Gettner et al. (2023) for independent observations. We can extend results to amortizing over any finite number of transitions, leading to a partially factorized approach (Appendix Sec. A.2.4).
4

---

Published as a conference paper at ICLR 2025
3.2.2 LOCAL SCORE COMPOSITION
At inference time, need to recover the global score Sylob 1x:0T) Vo, log ( (=) by composing all local scores Sf (0 1=)+++)
Vo, log | + estimated by the network, thereby defining the compose function in Alg. 2 Finally, the composed global score Sglab will
be passed to a diffusion sampler to obtain samples from the target posterior (())=!T).
While the composition introduced in Eq. Bis valid for 0, it becomes invalid for any 0 (i.e., each noisy posterior). The reason for this is the following: The likelihood for the 'noisy' parameter p(x::T')0, (0.) longer satisfies the Markov property, even if p( az.:71pp (0) does (Weilbach et al. 2023 Gloeckler et al., 2024b, Geffner et al. 2023; Rozet & Louppe, 2023). In the i.i.d. setting, this issue has been tackled by Linhart et al. (2024), and we will adapt their approach to the Markov setting. The correct global score is intractable in most cases, but they proposed an approximation that can be directly adapted for Markovian time series:
Vo, log (0gp(00/x= = A(Oa)-
Σ
'V-'Vo E 1 s (0,/x', x²+1)=++++))) + (1-11)))) T) 'Vo. log log
where AO.)) - Σt 2a1,141 + (1 (1-7) Here E, denotes the covariance matrix of 0,) (the usually tractable denoising prior) and In,.,t1 denotes the covariance matrix of p(0) O, 1+) (the denoising posterior), which we need estimate. This can be done by estimating the posterior covariance from samples of the local posteriors obtained via the diffusion model, referred to as GAUSS. Unless otherwise specified, we use this approximation as the default composition rule for FNSE. Alternatively, we can estimate it via Tweedie's 's moment projection using the Jacobian of the score estimator. referred to JAC (Linhart et al./ (2024), see Appendix B for details). In contrast, Geffner et al. (2023) addressed this issue (in the i.i.d. setting) through post-hoc sampling corrections (we referred to this uncorrected variant as FNPE, Appendix Sec. B.2 for details)
3.3 FACTORIZED LIKELIHOOD((RATIO) ESTIMATION (FNLE, FNRE)
For likelihood(-ratio) estimation, applying our approach is straightforward: by assumption, the like- lihood factorizes as shown in Eq. As consequence, we can directly learn the transition density p(a*+1 (e',)), similarly to the method described in Alg. 1. by replacing the score matching loss with the appropriate likelihood (or likelihood-ratio) loss (Papamakarios et al. [019 Hermans al. 2020 20022 Miller et af.) 2022). Once the transition density is obtained, the global log-likelihood approximation Eylab be computed directly from Eq. [T] (i.e., by simply summing up the local approximations). Classical MCMC techniques can be employed for sampling. We use reference im- plementations of NLE and NRE implemented in the shi package (Tejero-Cantero al. 2020b
Boelts et al. 2024), adapted to the Markovian setting.
3.4 CONSTRUCTION OF THE PROPOSAL DISTRIBUTION
A proposal p(x*) for single state transitions needs to satisfy conditions (i) and (ii) in Sec. 3.2.1
but solely relying on these properties does not necessarily make up a "good" proposal. Essentially, it specifies which regions of the data domain are represented in the training dataset. Therefore, better local posterior approximation the proposal compared to those that are less likely. This provides an opportunity and challenge to design appropriate proposals for given simulator.
expected for states that are more likely to be generated by
For an amortized posterior approximation, the proposal distribution design should be guided by the prior predictive distribution. The resulting trajectories o: ~ pT= encompass likely states x' and we can use this to construct reasonable proposal p(x') (e.g. by simply proposing these at random)-bowever however, this approach comes at the cost of upfront simulations. Alternatively, we can construct the proposal p(x') based on observations available prior to training. This will focus sim- ulation resources on the intended application but might not yield a good amortized approximation. For further details on proposal design, Appendix Sec.
5

---

Published as a conference paper at ICLR 2025
4 EMPIRICAL RESULTS
4.1 EVALUATION APPROACH
We empirically assess the accuracy and computational efficiency of the proposed approach, com- paring it against non-factorized NPE with appropriately chosen embedding net baseline. We note that the simulation budget is not determined by the number of simulator calls N, but rather by the number of calls to the state-transition function, i.e., T N. We, therefore, configure the NPE baseline with T max = 10 steps, i.c., if the simulation budget is 10k, the NPE baseline used 1000 simulations, each of which simulated for 10 steps (experiments with longer segments in Appendix E). We additionally use data augmentation by duplicating shortened variants (T < T of these simulations. This aims to help the RNN emhedding net to generalize to different sequence lengths. To also amortize over different initial conditions, the initial condition for each ten-step simulation was drawn from the proposal also used in the factorized methods. As main metrics for comparison, we use sliced Wasserstein distance (sW,) and Classifier two sample test accuracy (C2ST) on refer- ence posterior samples Lopez-Pazz & Oquab 2016 Kolouri et al., 2019 Bischoff et al. 2024). We average the estimated value over total of 10 randomly drawn observations. The whole process, i.e., training, sampling, and evaluation, was repeated five times, starting from different random seeds.
We begin by evaluating the methods simulators (Sec. (.2).. Next, we apply the approach to classical models from ecology and epidemi- ology, including the stochastic Lotka-Volters and SIR models (Sec. 4.3). Finally, we demonstrate the scalability of the method on a large-scale Kolmogorov flow task, where we perform inference on very high-dimennn dimensional data using only 200k simulator steps (Sec. 4.4
of newly designed benchmark tasks for Markovian
4.2 BENCHMARKS
To investigate several properties the proposed approach, we develop several synthetic tasks of first-order Markovian time series with associated reference posterior samplers (details in Ap- pendix[.II..
Gaussian RW: A Gaussian Random Walk of form x = a - 0 + for (0,, < with x² E gd and E Rdd This simple task offers an analytic Gaussian posterior.
Mixture RW: A Mixture of Gaussian Random Walk of form a'+1 x' + uu 0 + for (0, I), u ~ Unif({ 1., 1.})) and x', E Rd. By design, this task has a mixture of Gaussian transition density and non- Gaussinn bimodal posterior.
Periodic/Linear SDE: Both tasks are governed by the linear SDE dx' = A(0)x' + B(0)dw" The periodic SDE is oscillatory with 0, xx E RR² and posterior with four modes. The linear SDE involves A and B with x² E R and 18
Double well: A nonlinear SDE dæ' potential with modes position depending on 0,, O2 (Singer 2002 Cai et al. 2023).
0, + fg (x adw'" which samples from double-well
We first examine the scalability of all evaluated methods ality of the time series while maintaining a fixed simulation budget. Specifically, we assess each method on the Gaussian random walk (RW) example in 1, 2, and 10 dimensions, using total of 10k simulations (Fig. 2 Both FNLE and FNSE show better simulation efficiency than the NPE baselines. FNRE is more sensitive local errors, leading a decline in performance over a long time series. A similar observation was made by Geffner et al. (202)) for i.i.d. data. Notably, this issue is reduced in the FNLE approach, which was not analyzed in their work but successfully ap- plied by Boelts et al. (2022). On the other hand, the computational cost escalates significantly on long time series, as evaluating the global likelihood requires T forward passes for each iteration of MCMC sampling, making FNLE and FNRE relatively slow (especially NLE, as a forward pass through the normalizing flow in FNLE is more costly than the classifier used in FNRE, Fig. 2 In contrast, the sampling cost for FSNE remains moderate due to its more efficient sampling method, while performance remains similar to or better than FNLE.
function of the length and dimension-
In contrast to FNLE and FNRE. the approximative the sampling quality for FNSE. Therefore, we investigate the scaling behavior of different score composition techniques for long time series (Fig. 2p. Appendix Fig. (8)). Previous work within
composition (for 0) might impact

---

Published as a conference paper at ICLR 2025
a
Gaussian RW (td)
Gaussian RW (2d)
Gaussian RW (10d)
Gauts RW [2d]
2-1
2
Sampling time (9)
10a
MS
2
2-3
2
104
2
[5d]
100
sWi
10" 10 10
10
10" 10 10 10
10"
10
10
10
10
Transitions
Transitions
Transitions
50
100
500
NPE
FNLE
FNRE
FNSE
Transitions GAUSS
FNPE
JAC
Periodic SDE
Linear SDE
Mixture RW (5d)
Double well
SW: (100k sim.) sW. (10k Sim.)
(10d)
2
2
2
T
10
100
sW.
2
2
Mixture (5d)
2
2-4
20
21
SW1
23
2
2
0
100
10c
102
10a
10a
102
100
10: Transilions
102
100
10a
10a
10-1
10
10
10a
Transilions
Transitions
Transitions
Proposal std.
NPE
FNLE
FNRE
FNSE
Fig 2: Benchmarks: We validate our method on Gaussian Random Walk with different dimen- sions for different lengths (i.e. Transitions), also tracking sampling times (a). We assess FNSE score accumulation over Gaussian RW and Periodic SDE tasks using fixed Euler-Maruyama sam- pler (b). We compare methods across tasks and transition steps (c). Finally, we examine the effect of the proposal on NFSE trained with 10k simulations from normal proposal of varying standard deviation (d).
the i.i.d. case did not investigate performance beyond thirty samples (Geffner et al. 2023 Linhart al., 2024). Specifically, for time series exceeding length of T > 100, the FNPE method starts to become numerically unstable due to diverging trajectories. Although the assumptions of the GAUSS/IAC approximations are violated (as the Mixture RW has non-Gaussian posterior), these methods still tend to perform well in practice in our implementation. Even when Langevin correc- tions (Song al. 2021 Geffner et al.) 2023) are introduced, the GAUSS and JAC method remain superior, especially for large 7 (Fig. AT0 Appendix Sec. B.2
Next, we conducted benchmarks across several tasks. Overall, FNSE and FNLE show higher ac- curacy compared to FNRE and NPE in most tasks (Fig. 2r, Appendix Fig. (c). An exception is the double-well task, where NPE performs relatively well due the SDE's rapid convergence to stationary distribution, enabling the RNN to learn generalizable summary statistic for infer- ence. FNSE strongly outperforms other methods in cases of multimodal transition densities (c.g., Mixture RW) and performs comparably to likelihood-hasee approaches in simpler cases (e.g., Peri- odic/L.nnear SDE) even the posterior is strongly multimodal (Periodic SDE). Learning from single- step transitions is the extreme end of spectrum of possibilities. We, therefore, also evaluate this benchmark on partially factorized variants using five transitions each (Appendix Fig. 9. Overall, this could improve (e.g., Double Well) hurt (e.g., Mixture RW) performance relative the NPE baseline. We also compared method to broader spectrum of applicable baselines (Chen et al. (22T 2023); Sharrock al. (2024). Appendix Table AB and plot our results also against training budget (Appendix Fig. A6).
Finally, we investigate the impact of the proposal distribution, specifically for the FNSE method trained over 10k simulation steps (Fig. 2i. extended results in Fig. ALT.. We use a Normal dis- tribution centered at zero, adjusting its standard deviation ranging from a 0.1 (too narrow) to 1000 (too wide). Our results indicate relatively broad range of standard deviations of the pro- posal that result in good performance. If the distribution is too narrow, the observation may reach values outside the training domain, while if is too wide, the model is trained on values not observed during evaluation. In addition, we compared our default choice to proposals constructed using prior predictive samples (using a subset of the simulation budget, Sec. 3.4 Appendix Sec. E.4). Overall,

---

Published as a conference paper at ICLR 2025
a
Latka Volterra
b
SIR
posterior predictive
posterior predictive
posterior
5
1.5
2
60
0-
3
0
-1.0
0-
Time posterior
100
3
0
50
100
-3
60
3
-3
0:
3
-3
82
83
Time
Xg/9c
Xg/8.
T-10
T-100
C2ST
sW;
d
C2ST
sW1
Lotka Volterra
1.0
10
1.0
10a
10"1
10
10
100
100
10
100
10
100
Transilions
Transitions
Transitions
Transitions
NPE
FNLE
FNRE
FNSE
Fig 3: Lotka Volterra and SIR experiments The FNSE approximate posterior (predictive) of the best performing FNSE model using 100k transition steps to train, visualized on subsequences from fixed observation and associated true parameter. (a, b). We then show the quantitative performance in terms of C2ST and sW1 for each task on ten randomly selected observations; each run is repeated five times (c. d).
this leads to similar results even if the simulation budget used to construct the proposal is subtracted from the training budget (Tab. [
4.3 LOTKA VOLTERRA AND SIR
To further evaluate our method, we tested it on two famous models from ecology and epidemiology: the Lotka- Volterra and
(SIR) simulators. The Lotka Volterra simu- lator models predator-prey dynamics through four parameters that govern prey and predator growth, hunting rates, and mortality. The SIR simulator is a fundamental model for understanding infectious disease spread. Although commonly used for benchmarking SBI fixed-size observations Lueck- mann et al. 2021), we adapted these models for our analysis (details in Appendix Sec. D.I.
We begin by demonstrating how the posteriors evolve more time points are observed using the NFSE method (Fig. bb.... In both tasks, we find that the posteriors converge to the true parameter values, and the posterior predictive simulations increasingly align with the time series observation. Notably, unlike in previous i.i.d. scenarios (Geffner et al. 2023 Linhart et al.] [202)), the poste- rior uncertainty about the parameters does not necessarily decrease as additional time points are included. This is particularly evident in the SIR task, where the posterior remains largely unchanged between T 10 and T 100 (Fig. BD). This can be explained by the fact that the initial dy- namics contain significant information about the parameters, while the later dynamics do not, i.e., the infected population consistently declines to zero, and the susceptible and recovered populations converge to steady-state values that are independent of the parameterization.
We trained all methods on this task using 100k transition simulations. We then evaluate the accuracy of all models for performing inference sequences of length T = 1, 10, 100. For both tasks, we see that the factorized methods (except FNRE in Lotka Volterra) generally outperform the NPE baseline (Fig. bdd, also Appendix Fig. 6 per training budget, Appendix Tab. 3 for additional baselines). NESE, in particular, significantly outperforms all other approaches on the SIR task. We additionally performed iimulation-basee calibration analysis (Talts et al. 2018) for FNSE and NPE, which, in line with these results, showed that FNSE better calibrated than NPE (Fig. A7).. Notably, on the SIR task, found NPE-in contrast FNSE-to be unable generalize to longer time
8

---

Published as a conference paper at ICLR 2025
X1,
Xzo
Xsp
Xap
Xso
Xao
Xzp
Xao
Xao
3
8,
12
T=10
0.8
P
1.2
pror
1
post.
[
%
10
T=100
prior
post.
Xol 2-10*
Fig 4: Kolmogorov flow experiment: Single example observation of length hundred (top row). We visualize the posterior distribution condition the observation up / 10 (a, left top), along with a quantitative comparison of the mean absolute error (MAE) between posterior predictive samples and prior predictive samples on fifty different observations (left, bottom), including the one shown above. The vorticity of two selected predictive samples is visualized on the right. This analysis is then repeated for the entire observation (b).
series. In addition, compared default proposal to constructed using simulations. For the Lotka Volterra task, this constructed proposal performed significantly better than our simple default (Appendix Tab. 2 Fig.
4.4 KOLMOGOROV FLOW
Finally, we consider task with very high dimensional observations, previously considered in the context of data-assimilation (Rozet & Louppe, 2003).. The simulator models incompressible fluid dynamics governed by the Navier-Stokes equations:
du Øt
- Vu+
Vptf P
V u
Re
where the velocity field, Re the Reynolds number, p the fluid density, p the pressure field, and f is the external forcing. We added Gaussian noise with standard deviation a 5- 10-3 after each transition. Following Kochkov et al. (202T); Rozet Louppe (2023), we consider two-dimensional domain [0, 2T] with periodic boundary conditions and an external forcing f that corresponds to Kolmogorov forcing with linear damping Chandler & Kerswell (2013); Boffetta &
Ecke 2012). We take the Reynolds number and density free parameters 0 (Re,
In contrast to Rozet Louppe (202)), we consider the problem of parameter inference on 6. We apply FNSE to this problem, using only 200k transition evaluations, which successfully recovers amortized posterior estimator that generalizes to long-term observations (Fig. 4p. Notably, in such high-dimensional state space, simple proposal distribution would be ineffective. But, we can both the initial distribution and simulator to propose a variety of feasible states (details in Appendix 100 are well-
Sec. E.5 The posterior distributions for observations with both T 10 and concentrated around the true parameter values that generated the specific observations (Fig. b).. Posterior predictive samples are, on average, significantly closer the observed data compared to prior predictives (in mean absolute enror calculated over 1000 predictive samples from 50 differ- ent observations). Although the difference less pronounced for to the divergence for slight parameter modification on long simulations. To further investigate the performance, we perform simulation- based calibration analysis (Talts et al. (0200 (Fig. 12). Over- all, given the constraints the simulation budget, the posterior calibration okay but deteriorates for larger time series, indicating more simulations might be necessary to improve performance on T = 100.
100, this primarily due

---

Published as a conference paper at ICLR 2025
5 DISCUSSION
5.1 RELATED WORK
We build upon previous work in the ii.d. observation setting (Geffner et al., 2023 Linhart et al, 2024 Boelts et al. 2022), extending these methods to Markovian time series. In the context of data assimilation, Rozet & Louppe (2023) introduced method to learn a local score network to estimate p(x"-T) using a neighborhood o' Kit+ In contrast, we focus on solving a local inverse problem that leads to a global solution. Other work aims incorporate the task structure (approximately) into the estimating neural network architecture (Weilbach et al.
2020 2023; Gloeckler et al., 2024a). A
different line of work aims tackle high-dimensional SBI using summary statistics. These can be handcrafted (Alsing & Wandelt 2018). learned (Chen et al., 202 (023) or based on path signatures (Dyer et al., 202 2022). However, estimating statistics that generalize across sequence lengths from sparse training can be challenging. We enhance the simulation efficiency of amortized SBI by leveraging the Markovian structure, which is different from sequential training schemes (Glöckler et al. 2022| Sharrock et al. 2024) Greenberg et al. 2019) Durkan et al., 2020; Deistler et al. 2022 Gutmann et al. 2016] Warne al. 2022). which adaptively make simulations more informative for specific observations (but are thus not amortized). Similarly, other lines of work aim to distribute simulations based on cost function (Bharti et al. (2024) and increase simulation diversity (Niu et al. 2023). The simulation efficiency of our approach could be further enhanced by incorporating such techniques. Recently, a simulation -bassd conditional kernel density approximation for inference in SDEs was proposed (Cai al. 2023), which can be viewed
special of FNLE, substituting normalizing flows with kernel density estimation. Efficient inference in time series models was also previously approached in the ABC setting (Warne et al. 2022
Prangle 2016 Bernton et al. 2019
Dyer et al. 2021b).
2000 Arasaratnam & Haykin, 2009; Julier & Uhlmann 1997) is inherently related to FNLE, as these methods can an online (non-amortized) manner, but do require the availability of likelihoods (and are used to perform inference e.g. (Brouste et al. 2014)). The field of amortized inference with tractable likelihoods is inherently related (Le al. 2017 Choi et al.)
Moreover, any Kalman filtering method (Kalman 1960 Wan & Van Der Merwe
approximate the marginal log- ligeehoodd
2019 Ganguly et al. 2023 Margossian & Blei 2024).
5.2 LIMITATIONS AND FUTURE WORK
Our primary focus here was on inference for Markovian stochastic processes that share the same transition distribution across time. We also sketched extensions to cases where it varies with time, well as to processes with parameterized initial distributions. We note that Hidden Markov Mod- els, although widely applied practical scenarios, do not fall within the scope of our method. In combination with data assimilation techniques (Rozet & Louppe 2023) the Markovian hidden state can be recovered to which our technique can be directly applied. Extending our approach to more general probabilistic models remains direction for future research.
The score corrections caused by score aggregation described in Subsection 3.2.2 could be further improved. In their methods, [inhart al. (2024) approximates the data distribution by single Gaussian distribution, which can be violated in real-world scenarios. Our results demonstrate that this approximation empirically still performs well even if these assumptions are violated.
5.3 CONCLUSIONS
We introduced wawpulllatioo-fffiiieeeee approach for amortized inference in Markovian simulators. Al- though flexible embedding networks are commonly used for high-dimensional time series (Lueck-) mann et al. 2017 Radev et al. 2020 Schmitt et al. 2023), they often demand extensive amounts of simulations and can be fragile when faced with data perturbations (Cannon et al. 2022 Gloeck- ler et al. 2023). Success in these methods hinges on the embedding network's 's ability to capture robust and generalizable representations. In contrast, our approach decomposes the inference task into smaller, locally solvable problems, reducing computational costs and enhancing scalability for large-scale, complex simulations.
10

---

Published as a conference paper at ICLR 2025
SOFTWARE AND DATA
We use JAX (Bradbury et al. 2018) as the computational backbone and hydra Yadan 2019)
track configurations. For reference implementation of baselines, sbi Tejero-Cantert et al., 2020a; [Boett et al. 2024). Code to reproduce the experiments can be found in https:
//github. com/ www.ccooommmmccec.cc http//macelab/maaa. markasbbamakksbil//
AUTHOR CONTRIBUTIONS
Conceptualization: MG, ST, KF, JHM. Methodology: MG, ST. Software: MG. Investigation and Analysis (Benchmarks)) MG, ST. Investigation and Analysis (Lotka Volterra, SIR and Kolmogorov Flow): MG. Visualization: MG. Writing: MG, ST. Writing (Review & Editing): KF, JHM. Funding acquisition: ST, KF, JHM. Supervision: KF, JHM.
ACKNOWLEDGMENTS
We thank all members of the Mackelab, as well Keisuke Yano and Yoshinobu Kawahara for discussions. MG and JHM are supported by the German Research Foundation (DFG) through Ger- many's Excellence Strategy (EXC-Number 2064/1, Project number 390727645), the German Fed- eral Ministry of Education and Research (Tubingen AI Center). the "Certification and Foundations of Safe Machine Learning Systems Healthcare" project funded by the Carl Zeiss Foundation, and the European Union (ERC, DeepCoMechTons 101089288). MG is member of the International Max Planck Research School for Intelligent Systems (IMPRS-IS), ST and KF are supported by JST CREST JPMJCR2015. ST is supported by Grant-in-Aid for Research Activity Start-up 23K19966 and Grant-in-Ar Aid for Early-Cared Scientists 24K20750. KF is partially supported by JSPS Grant- in-Aid for Transformative Research Areas (A) 22H05106.
References
Justin Alsing and Benjamin Wandelt. Generalized massive optimal data compression. Monthly Notices of the Royal Astronomical Society: Letters, 476(1):L600-66)) February 2018. ISSN 1745- 3933.
Brian DO Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Ap-
plications, 12(3):313-326, 1982.
Christophe Andrieu and Gareth O Roberts. The pseudo-marginal approach for efficient monte carlo computations. The Annals of Statises, 2009.
Ienkaran Arasaratnam and Simon Haykin. Cubature kalman filters. IEEE Transactions automatic control, 1254-129,,, 2009.
Matthew Beal. Variational Algorithms for Approximate Bayesian Inference. PhD thesis, Gatsby Computational Neuroscience Unit, University College London, 2003.
Mark Beaumont, Wenyang Zhang, and David Balding. Approximate bayesian computation in population genetics. Genetics, 162(4):2025-2035 2002.
Espen Bernton, Pierre E. Jacob, Mathieu Gerber, and Christian Robert. Approximate bayesian computation with the wasserstein distance. Journal of the Royal Statistical Society Series B: Statistical Methodology, 81(2):235-269, 02 2019. ISSN 1369-7412. doi: 10.1111/rsb.12311.. 12312
URL htt.httts///dttt http://dooi.ooo /ogp:::///oo...//...11111
Michael Betancourt. Hamiltonian monte carlo for bayesian inference.
arXiv preprint
1711v01701.024344 2017.
Ayush Bharti, Daolang Huang, Samuel Kaski, and François-Xave Xavier Briol. Cost-aware simulation- based inference. arXiv preprint arXiv:2410.07930 2024.
Sebastian Bischoff, Alana Darcher, Michael Deistler, Richard Gao, Franziska Gerken, Manuel Gloeckler, Lisa Haxel, Jaivardhan Kapoor, Janne K Lappalainen, Jakob H Macke, et al. A practi- cal guide to sample- sample-based statistical distances for evaluating generative models in science. Trans- actions Machine Learning Research, 2024.
11

---

Published as a conference paper at ICLR 2025
Jan Boelts, Jan-Matthis Lueckmann, Richard Gao, and Jakob H Macke. Flexible and efficient simulation-bast -bastd inference for models of decision- decision-making. Elife, 11:e77220, 2022.
Jan Boelts, Michael Deistler, Manuel Gloeckler, Álvaro Tejero-Cantero, Jan-Mathis Lueckmann, Guy Moss, Peter Steinbach, Thomas Moreau, Fabio Muratore, Julia Linhart, Conor Durkan, Julius Vetter, Benjamin Kurt Miller, Maternus Herold, Abolfazl Ziaeemehr, Matthijs Pals, Theo Gruner, Sebastian Bischoff, Nastya Krouglova, Richard Gao, Janne K. Lappalainen. Bálint Muesányi, Felix Pei, Auguste Schulz, Zinovia Stefanidi, Pedro Rodrigues, Comelius Schröder. Faried Abu Zaid, Jonas Beck, Jaivardhan Kapoor, David S. Greenberg, Pedro J. Gonçalves, and Jakob H. Macke. shi reloaded: toolkit for simulation-base inference workflows, 2024.
Guido Boffetta and Robert Ecke. Two-dimensional turbulence. Annual review of fluid mechanics, 44(1):427-4511 2012.
Benjamin Boys, Mark Girolami, Jakiw Pidstrigach, Sebastian Reich, Alan Mosca, and O. Deniz Akyildiz. Tweedie moment projected diffusions for inverse problems, 2024.
James Bradbury, Roy Frostig. Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne and Qiao Zhang. JAX: composable transformations of PythonNummPy programs, 2018.
Johann Brehmer and Kyle Cranmer. Simulation- inference methods for particle physics. In Artificial Intelligence for High Energy Physics, pp. 579-611. World Scientific, 2022.
Alexandre Brouste, Masaaki Fukasawa, Hideitsu Hino, Stefano Iacus, Kengo Kamatani, Yuta Koike, Hiroki Masuda, Ryosuke Nomura, Teppei Ogihara, Yasutaka Shimuzu, Masayuki Uchida, and Nakahiro Yoshida. The yuima project: A computational framework for simulation and inference of stochastic differential equations. Journal of Statistical Software, 57(4):1-51, 2014.
Xin Cai, Jingyu Yang, Zhibao Li, and Hongqiao Wang. Simulation-based transition density approx- imation for the inference of sde models. arXiv preprint arXiv:2401.02529 2023.
Patrick Cannon, Daniel Ward, and Sebastian M. Schmon. Investigating the impact of model mis- specification in neural simulation- inference, 2022.
Jeffrey Chan, Valerio Perrone, Jeffrey P Spence, Paul A Jenkins, Sara Mathieson, and Yun S Song. A Likelihood-Free inference framework for population genetic data using exchangeable neural networks. Advances in neural information processing systems, 31 31:8594-8605, 2018.
Gary J. Chandler and Rich R. Kerswell. Invariant recurrent solutions embedded in turbulent two- dimensional kolmogorov flow. Journal of Fluid Mechanics, 722:554-595, 2013.
Yanzhi Chen, Dinghuai Zhang, Michael Gutmann, Aaron Courville, and Zhanxing Zhu. Neural approximate sufficient statistics for implicit models, 2021.
Yanzhi Chen, Michael U Gutmann, and Adrian Weller. learning summary statistics necessary for likelihood-free inference? In International Conference on Machine Learning, pp. 4529-4544, PMLR, 2023.
Kyunghyun Cho, Bart van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Hol- ger Schwenk, and Yoshua Bengio. Learning phrase representations using RNN encoder-decoder for statistical machine translation. In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), 2014.
Kristy Choi, Mike Wu, Noah Goodman, and Stefano Ermon. Meta-amortized variational inference and learning. In International Conference on Learning Representation, 2019.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation namula inference. Proceedings of the National Academy of Sciences, 1 117(49)33055-30006 2020.
Maximilian Dax. Stephen R. Green, Jonathan Gair, Jakob H. Macke, Alessandra Buonanno. and Bernhard Schölkopf. Real-time gravitational wave science with neural posterior estimation. Phys. Rev. Lett., 127:241 103, Dec 2021.
12

---

Published as a conference paper at ICLR 2025
Michael Deistler, Pedro J. Goncalves, and Jakob H. Macke. Truncated proposals for scalable and hassle-free simulation-base inference. In Advances in Neural Information Processing Systems, 2022.
Arnaud Doucet, Adam M Johansen, et al. A tutorial on particle filtering and smoothing: Fifteen years later. Handbook of nonlinear filtering, 12(656-704):3, 2009.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Ad- vances in neural information processing systems, 32, 2019.
Conor Durkan, lain Murray, and George Papamakarios. On contrastive learning for likelihood inference. In International Conference on Machine Learning. pp- 2771-2781. PMLR, 2020.
Joel Dyer, Patrick Cannon, and Sebastian Schmon. Approximate bayesian computation with path signatures. arXiv preprint arXiv:2106. 12555, 2021a.
Joel Dyer, Patrick Cannon, and Sebastian M Schmon. Approximate bayesian computation with path signatures. arXiv preprint
2021b.
Joel Dyer, Patrick Cannon, and Sebastian M Schmon. Amortised likelihood inference for expensive time-series simulators with signatured ratio estimation. In International Conference Artificial Intelligence and Statistics, pp. 11131-11144. PMLR, 2022.
Ankush Ganguly, Sanjana Jain, and Ukrit Watchareeruetai. Amortized variational inference: A systematic review. Journal of Artificial Intelligence Research, 78:167-215, October 2023. ISSN 1076-9757.
Tomas Geffner, George Papamakarios, and Andriy Mnih. Compositional score modeling for simulation-based inference. In International Conference Machine Learning, pp.
PMLR, 2023.
Andrew Gelman, Aki Vehtari, Daniel Simpson, Charles C Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner, and Martin Modrák. Bayesian workflow. arXiv preprint 2rXiv:201..0/000 2020.
Samuel Gershman and Noah Goodman. Amortized inference in probabilistic reasoning. In Pro- ceedings of the annual meeting of the cognitive science society, volume 36, 2014.
W.R. Gilks, S. Richardson, and D. Spiegelhalter. Markov Chain Monte Carlo in Practice. Chapman & Hall/CRC Interdisciplinary Statistics. Taylor & Francis, 1995. ISBN 9780412055515.
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, and Jakob H. Macke. Adversarial robustness of amortized bayesian inference, 2023.
Manuel Gloeckler, Michael Deistler, Christian Weilbach, Frank Wood, and Jakob H. Macke, All-in- one simulation-basee inference, 2024a.
Manuel Gloeckler, Michael Deistler, Christian Dietrich Weilbach, Frank Wood, and Jakob H. Macke. All-in-one simulation-basse inference. In Proceedings of the 41 International Con- ference Machine Learning, 2024b.
Pedro Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William Podlaski, Sara Haddad, Tim Vogels, et al. Training deep neural density estimators to identify mechanistic models of neural dynamics. Elife, 9:e56261, 2020.
A Goswami, BV Rao, A Goswami, and BV Rao. Markov models in genetics. A Course in Applied Stochastic Processes, pp. 179-170, 2006.
David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood-free inference. International Conference on Machine Learning, pp. 2404-2414, PMLR, 2019.
13

---

Published as a conference paper at ICLR 2025
Michael U Gutmann, Jukka Cor, al. Bayesian optimization for likelihood- free inference of simulator- baase statistical models. Journal of Machine Learning Research, 17(125):1-47, 2016.
Yoshihiro Hayashi, Junichiro Shiomi, Junko Morikawa, and Ryo Yoshida. RadonPy: auto- mated physical property calculation using all-atom classical molecular dynamics simulations for polymer informatics. npj Computational Materials, 8(1):222, November 2022. ISSN 2057-3960. doi: 10. 10.103//411244022-000000 841524-022-00908
URL https:// org/10 1038/
Joeri Hermans, Volodimir Begy, and Gilles Louppe. Likelihood-fres meme with amortized approxi- mate ratio estimators. In International Conference on Machine Learning, pp. 2020.
PMLR,
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, 2022. ISSN 2835-8856.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, volume 33, 2020.
Scott A. Hollingsworth and Ron O. Dror. Molecular dynamics simulation for all. Neuron (Cam- bridge, Mass.), ). 99(6): 129-1 2018. ISSN 0896-6273.
Aapo Hyvärinen and Peter Dayan, Estimation of non-normalized statistical models by score match- ing. Journal of Machine Learning Research, 6(4). 2005.
Simon Julier and Jeffrey K Uhlmann. New extension of the kalman filter to nonlinear systems. In Signal processing, sensor fusion, and target recognition VI, volume 3068, pp. 182- 193. Spie, 1997.
Rudolph Emil Kalman. A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82, 1960.
Tero Karras, Miik Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion- based generative models, 2022.
Samar Khatiwala. Efficient spin-up of earth system models using sequence acceleration. Sci- Advances, 10(18):cadn2839, 2024. doi: 10.1126//ciadvaadn2289, URL pattps://ve
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan computational fluid dynamics. Proceedings of the National Academy of Sciences, (18(21):e10118811 18, 2021.
Hoyer. Machine
Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, and Gustavo Rohde. Generalized sliced wasserstein distances. Advances in neural information processing systems, 32, 2019.
Tuan Anh Le, Atilim Gunes Baydin, and Frank Wood. Inference compilation and universal proba- bilistic programming. In Artificial Intelligence and Statistics, pp. 1338-1348. PMLR, 2017.
Julia Linhart, Gabriel Victorino Cardoso, Alexandre Gramfort, Sylvain Le Corff, and Pedro L. C. Rodrigues. Diffuusio posterior sampling for simulation- inference in tall data settings, 2024.
David Lopez-Paz and Maxime Oquab. Revisiting classifier two-sample tests. arXiv preprint arXiv: 06540, 2016.
Jan-Mauhis Lueckmann, Pedro J Goncalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher. and Jakob H Macke. Flexible statistical inference for mechanistic models of neural dynamics. Advances in neural information processing systems, 30, 2017.
Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke. Bench- marking simulation-base inference. In International Conference on Artificial Intelligence and Statisties, pp. 343-351. PMLR, 2021.
14

---

Published as a conference paper at ICLR 2025
Charles C. Margossian and David M. Blei. Amortized variational inference: When and why?, 2024.
Benjamin K Miller, Christoph Weniger, and Patrick Forré. Contrastive neural ratio estimation. Ad- vances Neural Information Processing Systems, 35:3262- 3278, 2022.
Ziang Niu, Johanna Meier, and Srançois-Xaviee Briol. Discrepancy-bassee inference for intractable generative models using Quasi-Monte Carlo. Electronic Journal of Statistics, 17(1):1411 1456, 2023. doi: 10. 00.121/23-EJS21133
URL bttps:://///oo
Franco Ormeño and Ignacio J. General. Convergence and equilibrium in molecular dynamics sim- ulations. Communications Chemistry, 7(1):26, February 2024. ISSN 2399-3669. doi: 10. 542004-02--01114-3 URL arg.ioog//iii 100 1t2/:1100000000 1038/ 112000000-244--1111- 024-01114-5
George Papamakarios and lain Murray. Fast e-free inference of simulation models with bayesian conditional density estimation. Advances in neural information processing systems, 29, 2016.
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. Advances in neural information processing systems, 30, 2017.
George Papamakarios, David Sterratt, and lain Murray. Sequential neural likelihood: Fast likelihood-fie free inference with autoregressive flows. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 837-848. PMLR, 2019.
Januara Dennis Prangle. Lazy abc. Statistics and Computing, 26(1-2):111-185, January 2016. ISSN 0960- 3174.
Stefan Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone, and Ullrich Köthe. Bayesflow: Learning complex stochastic models with invertible neural networks. IEFE transactions on neural networks and learning systems, 33(())1452-144 2020.
Stefan T. Radev. Marvin Schmitt, Valentin Pratz, Umberto Picchini, Ullrich Köthe, and Paul- Christian Bürkner. Jana: Jointly amortized neural approximation of complex Bayesian models. In Robin Evans and Ily Shpitser (eds.), Proceedings of the Thirty-Ninth Conference on Uncer- tainty in Artificial Intelligence, volume 216 of Proceedings of Machine Learning Research, pp- 1695-1706. PMLR, 31 Jul-04 Aug 2023.
Jishnu Ray Chowdhury and Cornelia Caragea. Recursion in recursion: Two-level nested recursion for length generalization with scalability. Advances in Neural Information Processing Systems, 36, 2024.
Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows In Proceedings of the 32nd International Conference on International Conference on Machine Learning-Vone Volume 37, pp. 1530-1538. JMLR. org, 2015.
François Rozet and Gilles Louppe. Score-based data assimilation. In Thirty-seventh Conference on Neurul Information Processing Systems, 2023.
Marvin Schmitt, Stefan T. Radev, and Paul-Christian Bürkner. Fuse it or lose it: Deep fusion for multimodal iimulation-based inference, 2023.
Katrin Sedlmeier, Sehastian Mieruch, Gerd Schädler, and Christoph Kottmeier. Compound extremes in a changing climate-a markov chain approach. Nonlinear processes in geophysics, (3(6:375 390, 2016.
Louis Sharrock, Jack Simons, Song Liu, and Mark Beaumont. Sequential neural score estimation: Likelihood--free inference with conditional score based diffusion models. In Proceedings of the s International Conference on Machine Learning. 2024.
Hermann Singer. Parameter estimation of nonlinear stochastic differential equations: simulated maximum likelihood versus extended kalman filter and itô-taylor expansion. Journal of Compu- tational and Graphical Statistics, 11(4):972-995, 2002.
Anshul Sirur, David De Sancho, and Robert B Best. Markov state models of protein misfolding.
The Journal of Chemical Physics, 144(7), 2016.
15

---

Published as a conference paper at ICLR 2025
Anders Sjöberg. Jakob Lindqvist, Magnus Önnheim, Mats Jirstrand, and Lennart Svensson. Meme- correction of score-based diffusion models for model composition, 2024.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Bus E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Yang Song, Jascha Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In Interna- tional Conference on Learning Representations, 2021.
R. Du Toit Strauss and Frederic Effenberger. A hitch-hiker's guide to stochastic differential equa- tions: Solution methods for energetic particle transport in space physics and astrophysics. Space Science Reviews, 21 212(1-2):151-192, March 2017. ISSN 1572-9672.
Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman. Val- idating bayesian inference algorithms with iauulttion-bsee calibration. arXiv preprint arXiv: 1804.06888,,8 2018.
Scan Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari. and Andrew Gelman. Validating bayesian inference algorithms with simulation-basee calibration, 2020.
Alvaro Pejero-Cantero, Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pe- dro J. Gonçalves, David S. Greenberg, and Jakob H. Macke. sbi: toolkit for simulation-bed inference. Journal of Open Source Software, 5(52):2505, 2020a.
Alvaro
Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pe- dro J. Gonçalves, David S. Greenberg, and Jakob H. Macke. sbi: A toolkit for simulation-basee inference. Journal of Open Source Software, 5(52):2505, 2020b.
Eric A Wan and Rudolph Van Der Merwe. The unscented kalman filter for nonlinear estimation. In Proceedings of the IEEE 2000 adaptive systems for signal processing, communications, and control symposium (Cat. No. 00EX373), pp. 153- leee, 2000.
David J. Warne, Thomas Prescott, Ruth E. Baker, and Matthew J. Simpson. Multifidelity multi- level monte carlo to accelerate approximate bayesian parameter inference for partially observed stochastic processes. Journal of Computational Physics, 469:111543, 2022. ISSN 0021-9991.
doi: https//dbb.oor1110.01///...g2p2p 2022.1111541. URL https: www. sciencedirect.
com/
D. Watson-Parris, A. Williams, L. Deaconu, and P. Stier. Model calibration using esem v1.1.0 an open, scalable earth system emulator. Geoscientific Model Development, 14(12):7659-7672 2021.
Christian Weilbach, Boyan Beronov, Frank Wood, and William Harvey. Structured conditional con- tinuous normalizing flows for efficient amortized inference in graphical models. In International Conference on Artificial Intelligence and Statistics, pp. 4441-4451. PMLR, 2020.
Christian Dietrich Weilbach, William Harvey, and Frank Wood. Graphically structured diffusion models. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), Proceedings of the 40th International Conference Ma- chine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 36887-36909. PMLR, 23-29 Jul 2023.
Jonas Bernhard Wildherger, Maximilian Dax, Simon Buchholz, Stephen R Green, Jakob Macke, and Bernhard Schölkopf. Flow matching for scalable Aimulation-basee inference. In Thirty- seventh Conference on Neural Information Processing Systems, 2023,
Simon N Wood. Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466 102-00110 2010.
Omry Yadan. Hydra- framework for elegantly configuring complex applications. Github, 2019.
16

---

Published as a conference paper at ICLR 2025
Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov, and Alexander Smola. Deep sets, 2018.
Yi Zhang, Arturs Backurs, Sébastien Bubeck, Ronen Eldan, Suriya Gunasekar, and Tal Wag- ner. Unveiling transformers with lego: synthetic reasoning task, 2022. URL https://arxiv, org/abx/226604430,, 2022.
17

---

Published as a conference paper at ICLR 2025
A METHODS
A.1 DERIVATION OF SCORE FACTORIZATION
3
For 0: following the simulator p(x" (0), the global score Ve logp(0 x0:T) is given by
Vo log ((0(x = Vo Volog {p(0(=u1) p( }} = Vo log {p(0) p(x (0)}
lolog p(0) p(x")) p( lao", 00
1()= + V log Volog {p(0) p(x" 00)
p(a 1xe",
|xe, 0)
r=0
Volog {p(0)=") p(²")} V log II p (x'
T-1 Σ
Vologp(0l=)) + V lo + lx'
(5)
Let xt follow any proposal distribution p(x'), and x+ follow the simulation transition p(a xt 0) given æ'. Noting that the replacement does affect the conditional distribution p(æ ,²,0), 0). we may assume x. ~ p(x') in (
By Chain Rule, we have
p(x ++1 01x²)= = p(a - |x'" 0) p(0)x²))
(6)
1 x++11 ,,),,,x)),,))))) 5(x |x"), (("),,
where P means x' follows the proposal p(x'). This leads us to the equality
p( x. 0) p(0 x²) 7(x'+1
(x0) p(0 ++1 xx) 5( x²)
p 0 x
p(0)
Here, the second equality is derived from the equality p(0]x²) p(0). to which the fact that p(x')
does not involve a parameter leads. Substituting this equality into Eq. 5 we have
p(0 x x²) p(x lx))
Vo Polgp(0(
+
Ve log
p(0)
T-1
Va log (((0))))))))
+
Vo log [+ ;,++,,)')))))))))),))) - T Vo logr(0)
(7)
0=0
Ve log 5(0 x0 +,) -T Vo log
Which concludes the proof. Here, the final equality is derived from the assumption p(a:"10) p(x").
A.2 EXTENSIONS
A.2.1 TIME INHOMOGENEOUS SIMULATION
The assumption that p(x'+1 +1 [= 0) does not depend on t, referred to as time-homogeneous transi- tion simulations, in the main manuscript be relaxed. Regarding FSNE, the score decomposition (3 also holds for a time- time-inhoman iame--nhooggeouut simulators where p( x O) varies among the time;
the score Vo log p (0)m', +.+1), however, then depends time t. To address the issue, train
18

---

The condition on initial state p(=" (0) can also be relaxed. Regarding FNSE, the global score Published as a conference paper at ICLR 2025
time- timeiinooogeeeeus score (0(0 x++1 .t) using training data D {O,, xz 2 ,t:],, where is drawn from random number generator, æ is drawn from a proposal distribution p(æ)), and aff+1 T(x |- 0,). FNPE and FNRE can be handled in the same way as FNLE by adding
as variable. Experimental evaluations are presented in Appendix E.6
A.2.2 INITIAL PARAMETER ESTIMATION
Ve log ( " can be factorized by
V 10 [000000 Ve logr(0)<>))
Va Foloog(00> + T Vo Vologp(0)
as proven Eq. This implies that the global score is obtained by estimating the two scores Vo log p(0ul=0))) and Ve log p(0 ( af separately. The estimation of these two scores can be combined by modeling a score function that outputs different scores depending on the number of input variables. FNPE and FNRE can be extended by estimating the two likelihoods p(æī(0) and p(ar) [1 a f), similar to FNSE.
A.2.3 HIGHER ORDER MARKOV CHAINS
Our proposal can be generalized to simulators with high order markov chain
p 0)) p(x 0mm-1m 00)
II
p( [+1
x m+ +++++1 ,0)
(8)
1=m-1
with its degree m. Under the setting, the score Vo log p (0 (:T) has the factorization
Ve logp(0)az = Vo log
x x.-m-1)
00m-1)
T-m-1
Σ ビ=m-1
Vo log 100 (0000mm++ -T -mmmm 1). V log
which implies that the global score be recovered by the two scores V log 0m and Ve log =/))====+++++)) Here, x m==m+++ follows a proposal p(z'-mit+1 and x m- follows the simulation transition (g 0) in the score Va log (0 The estimation of these two scores can be merged same as the initial parameter estimation case in the last section.
A.2.4 PARTIAL FACTORIZATION
Instead of considering just a single transition, the methodology naturally extends multi- ple transitions, also explored by Geffner et al. (2023). In PFNSE with M, target Ve (ogg(())' .......x*+... f+M ), while in PFNLE or PFNRE, we target p(=',... M (0) (up to con- stants), using dataset of M-step simulations. A key advantage of partially factorized methods is that they require fewer network evaluations for inference on a fixed-size observation, which can ac- celerate sampling and reduce the accumulation of local errors over time. Moreover, longer training simulations potentially lead to superior performance in long-term predictions as fewer local errors accumulate. However, they may require more simulations for effective training in return for these desirable properties.
A.2.5 MISSING VALUES IN TIME SERIES OBSERVATION
Our method is capable of handling missing data. The proposed factorization leverages the Markov properties time series simulators at observation points, which remain intact even in the presence of missing data. To preserve the amortized property regarding missing points, include the time step size At as an additional imput variable to the model. This allows to skip over timepoints that are missing and generally allows for irregular time grids.
19

---

Published as a conference paper at ICLR 2025
B SCORE COMPOSITION RULES
Our score decomposition 13 relies on the Markov assumption of the simulator ((x::T(() as shown in its proof (Appendix A. On the other hand, p(> [0,] p(x" ())(010a)) no longer satisfies the Markov property for > O, which renders Eq. invalli for any > 0. In this section, review the two existing remedies for tackling this problem.
B.1 METHOD BY GEFENER ET AL.. 2023
Geffner et al. 2023) avoids the issue by constructing an alternative sequence of distributions qa (0 < S A) that satisfies qo p(0)=0:7)) and qA (N:0:AA HA, o I). Adapting their proposal to the Markovian case, we obtain
T-1 Σ 80(0, lo', 1 (1- T) (1-1)(A a) V log (()).
Ve, (0.g((00
While this approach is generally applicable, it renders the backward SDE (Eq. 2 invalid as the marginals no longer follow the corresponding forward diffusion process. This requires the use of MCMC-based sampling corrections (in the sense of predictor-correct diffusion samplers (Song et al.; 2021 ), such as unadjusted Langevin dynamics (Geffner et al.) (222)) or similar approaches (Sjöberg
al. 2024).
B.2 METHOD BY LINHART ET AL. (2024)
Linhart et al. (2024) have recently proposed an alternative approach by approximating the marginal scores derived from the forward diffusion process. This enables the utilization of the reverse SDE and consequently allows for the application of standard diffusion-based sampling techniques. They showed that the correct score could indeed be written as suggested in Eq. plus an additional but intractable correction term of the form
Vol(O, Vel(0,x00T) ) log p(0 1- II 0 x', x +d0
)d0
To estimate it analytically, they employed Gaussian approximations, p(0 x', +1+ N(0; Ma,t,t11 x², +1) Eu,i (0a, x², (., 2²+1)) 0.) = (0; na (0.). En),
making the additional terms tractable. Eventually, have
Vo, log 1og D:T) R A(O))
Σ 1=0
E De. e.,,441 (0.l x1) + (1-T)E, 'Ve, logp(0.), Vo,
where A(O.) = Et 2.1.141 (1- This approach requires specifying two "hyper- parameters": the denoising prior covariance Σ. in 0,) and denoising posterior covariance Ea,,1411 in p(0)0,, x', x²+1), ). The denoising prior covariance is typically known analytically, along with the marginal prior score. Specifically, given a Gaussian prior with covariance Eg we have Σa (Eg² s(a)2/a(a)) I by Gaussian marginalization rules with the perturbation kernel p(0.(0) = (e s(a)0, (a) I)) specific the diffusion model (i.e. Appendix Sec. D.2).
In contrast, the denoising posterior covariance must be estimated. To address this, GAUSS assumes Gaussian clean posterior (i.c., at = 0) and analytically computes the denoising posterior covari- ances (analogous to the prior, but with an estimate of " This estimate can be obtained through samples, given that we have diffusion model that indeed can sample from p(0]æ*, x(+1). Alternatively, JAC directly estimates the denoising covariance iteratively using Tweedie's Moment projection (Boys et al. 2024), leveraging the Jacobian of the score network:
m(a)²
E Ea
s(a) (I Vo,
S (0(0a . a
If the Gaussian assumption is satisfied, both approaches are theoretically equivalent (neglecting numerical approximation errors). Yet. if these are violated, these yield different approximations to the denoising covariance matrices and hence will behave differently.
20

---

Published as a conference paper at ICLR 2025
Lotka Valterra (Phase Space)
b
Proposal
Postarior epproximation
12
12
FNSE( 100k)
0.3
Predstor
Predstar
FNSE (95 (+ proposal)
0.1
0-
10
100
Prey
Proy
Transitions
Fig 5: Proposal: A set of trajectories within the phase space of the Lotka Volterra task for different parameterizations sampled from the prior (constructed using total of 5k simulation steps) (a). The proposal was constructed by randomly sampling noisy points from the state space trajectories (b). Posterior approximation in sliced Wasserstein distance constructed by improved proposal compared to our "naively" chosen baseline (c).
B.3 PRACTICAL IMPLEMENTATION DETAILS
The GAUSS and JAC approximation is only valid if A is positive definite, as the derivation of this form depends on it. Ignoring this can lead to numerical instabilities. Unlike Linhart et al.) (2024). we do not clip the diffusion process to specific prior range to avoid numerical issues; we adapt all estimated posterior precision matrices to ensure positive definiteness of A. This adaptation resolves the problem even for long time series (i.e., T > 100; see Appendix E.3 for details). Importantly, since the GAUSS method is derived under Gaussian assumptions, presupposes that the poste- rior covariance matrices (and their corresponding scores) are smaller than those of the prior due to observational constraints. This assumption may not hold true, particularly in multi-modal tasks. For instance, in the Periodic SDE task, four modes are symmetrically positioned around the origin. While the variances of each mode may decrease, the total variance often still increases, especially when the modes distant from the origin. This situation thus invalidating the approximation. Similarly, for JAC, the estimated posterior covariance might not even be positive and definite as the Jacobian of the scoring network might badly represent the Jacobian of the true score.
render A(0a) non non-positive definite,
To address these challenges, we make minimal adjustments the initially estimated Eg,.e+1 to ensure that A(O.)) remains positive definite. We achieve this by considering its eigenvalue decom- position A = VAV-1, allowing us to identify the minimal adjustment required to enforce positive A_ definiteness: A_ = -V min(A, A, 0)V-1 Consequently, we update Σー1 a.t. ET -t,, + T + cT for a small nugget e. The A computed with these adjusted covariance matrices is guaranteed to be positive definite. The JAC method encounters similar issues, and we apply an analogous adjustment to its initial estimates.
We utilize the GAUSS composition rule our default method: We estimate the posterior precision matrices E as recommended by Linhart et al. (2024). For each posterior p(0) [|" x²+)), we produce 500 dimo samples in parallel from the diffusion model. These samples then used to estimate the associated covariance matrix. Although this approach introduces a slight constant overhead prior to each sampling run, it does require the iterative estimation of Jacobians, as is necessary with the JAC method. On the other hand, within our settings, the dimension of the parameter is often much smaller than the data dimension, rendering computation of the Jacobian not too expensive.
C PRACTICAL GUIDE TO PROPOSAL CONSTRUCTION
In this section, we provided a practical guideline for constructing the proposal p(x'). In general, the proposal will be task-specific. For an amortized approximation, a valid proposal p(x') should satisfy the conditions: (i) It must have support any point x' the simulator can attain (at evaluation time). (ii) The proposed point x' must be independent of the parameter involved in the transition to a I1 ensure Eq. 3 remains valid (Appendix Sec. A. But solety relying these properties
21

---

Published as a conference paper at ICLR 2025
Periodic SDE
Lincar SDE
Mixture RW
Double Well
Task
T=1
T=11
T=101
T=1
T=11
T=101
T=1
T=11
T=101
T=1
T=11
T=101
FNSE
0.53
0.60
0.76
0.53
0.70
0.90
0.52
0.69
0.69
0,50
0.59
0.81
FNSE (pred. prop) 0.53
0.60
0.71
0.53
0.70
0.90
0.52
0.69
0.75
0.50
0.57
0.80
Table I: Predictive proposal: Posterior approximation for different proposals given as C2ST metric for each benchmark task. Here, we compare our default proposal (trained with 100k step simula- tions), chosen by domain knowledge, and the prior predictive construction (using 5k simulation steps) training on the remaining 95k step simulations.
Lotka Volterra
SIR
T=1
T=11
T=101
T=1
T=11
T=101
FNSE FNSE (pred. prop)
0.57
0.78
0.94
0.67
0.78
0.74
0.52
0.65
0.88
0.67
0.79
0.74
Table 2: Predictive proposal: Posterior approximation for different proposals given as C2ST metric for Lotka Voletrra and SIR. Here, we compare our default proposal (trained with 100k step simu- lations), chosen by domain knowledge, and the prior predictive construction (using 5k simulation) trained on 95k step simulations.
does necessarily make up "good" proposal. Essentially, the proposal specifies which regions of the data domain are represented in the training dataset. On finite simulation budget, we thus would expect better posterior approximation for states that are more likely to be generated by the proposal than those that are less likely. This provides an opportunity and challenge to improve the proposal distribution design for given simulator.
This interpretation makes it fairly intuitive to know what might be good proposal and domain knowledge can be used to design a proposal from scratch. applicants know that the dynamics are bounded (or at least likely contained) within certain region, then one can design proposal that roughly covers it. the applicant knows that the simulator has stationary distributions, then these would serve good proposals. Yet, in general, it might be hard suboptimal to construct proposal this naive way. The default proposals used in our analysis are indeed simple distributions (e.g., Gaussians) chosen based such knowledge. Yet we here also provide some general approaches to select reasonable proposals in Sec. C.1 and Sec.
C.1 PRIOR PREDICTIVE CONSTRUCTION
Condition (ii) does not generally exclude the use of the simulator to construct a good proposal and we hence can use some simulation outputs to construct good proposal. As an example, let's consider the Lotka- Volterra task. By sampling from the prior predictive (i.e., sampling parameters from the prior and simulating data), we can visualize the trajectories ":7 of simulation outputs within its phase space (Fig. 5p). These states are likely to encompass most states, and we can use this to construct an effective proposal for our factorized methods, e.g., by simply sampling noisy points phase space trajectories uniformly (Fig. 5p). For example, one can:
(A) Sample a few parameters [ from the prior p(0)-
(B) Perform time series simulations <.:TT p(x (0:) for each parameters until a certain time x!"
(C) Estimate the proposal ((x'), e.g., by randomly selecting a state of from N (T + 1) states with some noise (which corresponds to a kernel density estimate).
What condition (ii) prevents that we cannot directly leverage transition from prior predictive sim- ulations (0; ....T) each transition will depend the same parameter, and thus, each triplet (0:, 2z x/+ /+1, will not satisfy that æ is independent of 0, By constructing the proposal above, apriori, we avoid this problem. Many time series simulators are often designed to ensure that their
outputs change minimally for large 7 (Hayashi et al. 2022 Khatiwala 2024 Ormeño & General
22

---

Published as a conference paper at ICLR 2025
2024 (a.k.a. reach stable states stationarity), number of simulation trajectories for not too many transitions can, thus, be sufficient to explore the phase space.
Attentive readers might recognize that doing so will itself cost simulations. A quantity aimed to reduce. To investigate if better proposal is worth this cost, we tested this on the Lotka Volterra task, comparing it with our default proposal x' Unif(O, 10), heuristically chosen by domain knowledge that trajectories will likely oscillate between zero and ten. We equalized the simulation costs for both experiments: 100k simulation steps are used to train the model in the heuristic proposal, whereas, for the prior predictive approach, 5k are allocated to construct the proposal, and 95k for training. The prior predictive construction further improved the performance, even if this cost is subtracted from the training budget (Fig. 55) We further ran this experiment for all other tasks, which showed that this approach indeed leads to proposals that perform similarly, sometimes better, than the default proposals chosen by us [Tab. Tab. An alternative and related approach, also designed to explore the phase space (but not requiring additional simulations), is the routine applied in the Kolmogorov flow experiment (Appendix Sec. E..). Here, we simply use different parameter each time perform a transition hence avoiding violating condition (ii).
C.2 TARGET-BASED CONSTRUCTION
We can also create the proposal p(x') using any sort of (real/synthetic) set of observations x:: accessible prior to training. For amortized inference we can proceed analogously to Sec. C.1
we are only interested in inference for single observation as in a sequential SBI setting Pa pamakarios & Murray) 2016 Greenberg et al., 2019 [Sharrock et al.) [2024) Deistler et al.) 2022). we can construct the proposal using this single observation by simply proposing observed states ay
at random. The full amortization property is lost in this way (as the support is constrained), but the simulation budget can be targeted on this more specific inference problem, which can lead to improved performance.
D EXPERIMENT DETAIL.
D.1 TASKS
In this section, we elaborate on all tasks used in the analysis in detail. For all tasks, we used standard normal prior but eventually reparameterized them certain constraints if necessary for the task (see individual tasks for specification).
Note that the training aims to also amortize over the initial distribution p(="). We evaluate the performance always on specific and fixed initial distribution (see individual tasks for specification), which is also different from the proposal p(x").
This task serves as simple baseline. It serves as an extension of the Gaussian Linear task introduced by Lueckmann et al. (02I) to the Markovian case. is defined by a Gaussian Random Walk of form
Gaussian RW:
x a- x 0 + for c~ c~,000 I), a = 0.9,
with x 0 E Fdd and 0 E Rd. This task offers an analytic solution for the posterior, which is Gaussian. Dimension can be set as wanted, we choose d 2, 10. As proposal, we choose p(x') = Nix=; vIOI), which is motivated by the fact that it the variance of the corresponding stationary distribution. For evaluation, we fixed p(æ") = 6(x²) (i.e., a point mass at zero).
Mixture RW: A Mixture of Gaussian Random Walk of form
f+1 x f 0+c for I(, u Unif(((-1,1)
with x',0 0 E Rdd We choose = 2, As proposal, we choose p(x') = NNa=;; v10T). For evaluation, we fixed p(=")= 5(x²) (i.e., point mass at zero).
By design, this task has a tractable transition density, which is a mixture of Gaussians, hence allow- ing exact computation of the marginal -oellike. In contrast to the previous task, this offers
23

---

Published as a conference paper at ICLR 2025
non-Gaussian bimodal posterior. We generate reference posterior samples using Hamiltonian Monte Carlo (HMC, (Betancourt) 2017)) with integration steps. To avoid mode collapse, we run 100 par- allel chains, each initialized different samples obtained by importance resampling of 50 initial samples from the prior. We run the sampler for 600 burn-in iterations, during which the integration step size is adapted accept around 80% of the proposed steps. We additionally thin the accepted samples in the chain by a factor of six to avoid autocorrelation.
Periodic SDE:
The periodic SDE is a linear SDE parameterized as follows:
0 03 02
dx
x² dt + 0.11 d.1,,
0
which gives us a stochastic oscillator. Each transition corresponds to simulating this SDE for 0.1 ms, using 20 Euler Maruyama steps. As proposal, choose p(x') (=="; 0, I). For evaluation, we fixed the initial distribution to p(x") = 6(( (.55 0. to quickly start with oscillations.
The corresponding marginal log-likelihood is computed using Kalman filter Kalman (1960). Ref- erence posterior samples are obtained in combination with HMC, using procedure similar to the one discussed above. This task offers multimodal posterior (4 modes) with all modes being point symmetric around the origin.
Linear SDE:
The linear SDE task is given by
x² (A(0) - + 0.5T)
where every entry of A and is directly parameterized by 0 R8 and x' E R3. Reference posterior samples are obtained in the same way for the previous task. As proposal we choose p(xx' = (aa",, 0, I). For evaluation, we fixed the initial distribution to =
Double well:
This is nonlinear SDE
dx' 0 x + 02 (z')² + a dw',
which samples from double-well potential with modes position depending 0,, 02 (Singer 2002). As proposal we use p(x²) p(=o) 81(1, 11 1, ) We use combination of particle filter (Doucet et al. 2009) and geeudo-marginal Metropolis- Hastiins MCMC (Andrieu & Roberts 2009) sampler to generate reference posterior samples.
2.5, as intial distirbution we at evaluation we
Lotka- Volterra: We use stochastic Lotka- model, a classic system of differential equa- tions used to describe predator-prey interactions. We use a stochastic variant, where the amount of noise also depends on the population of each species. The model dynamics are given by the following equations:
dx
x¹ - -B. x/ y'dtt ox dw
dt
dy' = - 3/ 8 x' of dt ox'y dwž.
dt
where x represents the prey population and y represents the predator population. The parameters a, G, 7 and control the interaction between the species. The noise hyperparameter was set to 0.05. The dynamic are constrained to remain positive. As a proposal we use p(x*) ) Unif(O,, 10). As evaluation time the initial distribution fixed to p(x") 8(x (1.,
To sample from the posterior by traditional means, we use a Particle Filter to obtain a stochas- tic estimator of the marginal likelihood. We use this stochastic estimator in pseuio-marginal- like Satropolis-Hasting MCMC algorithm to obtain reference samples from the posterior distribu-
tion (Doucet al. 2009 Andrieu & Roberts 2009).
24

---

Published as a conference paper at ICLR 2025
SIR: We use stochastic SIR method defined through the following equations:
ds
di - - 0.00duu..
dl' di
BS - - 001 1d02
dR' = + 0.01du3.
di
here S, I, and R represent the susceptible, infected, and recovered populations. We use as a proposal f(x') Unif(0,55) 5). As initial distribution, sample 5 UUnif(2.5,5). 70 Onif 2.5) and R ~ 8(R)).
To sample from the posterior by traditional means, estimator of the marginal likelihood. We use this stochastic estimator in a pseudo-marginal-iik Metropolis-Hastings MCMC algorithm obtain reference samples from the posterior distribution.
Particle Filter to obtain stochastic
D.2 TRAINING AND EVALUATION
For the NPE baseline, we utilize a 5-layer neural spline flow (Durkan et al. 2019), with each layer parameterized by a 2-layer MLP having hidden dimension of 50, Additionally, employ Gated Recurrent Unit (GRU) network Cho et al. (2014) as the embedding network, also with hidden dimension of 50. To enhance the RNN embedding network's ability to generalize across different sequence lengths, we apply data augmentation. Specifically. for each subsequence T < Tmmx = 10, we duplicate 10% of randomly selected parameter-data pairs and shorten the corresponding data time series to length T. This does not add any more simulation calls but directly trains the neural net for sequences < T maxx
For FNLE and FNRE, use adapted reference implementation implemented in the shi pack- age (Tejero-Cantero et al. 2020b). In FNLE, we use a 5-layer Masked Autoregressive Flow (Pa- pamakarios et al., 2017), each parameterized by 2-layer MLP with a hidden dimension of 50. In FNRE, we use resnet classifier with two blocks each considering 2 layers with hidden dimen- sion of 50. As MCMC sampling algorithm, mode-collapse, we run 100 parallel chains.
per-axis slice sampling algorithm. To avoid
Both approaches trained with training batch size of 1000 until convergence, determined by the default early stopping routine.
For FNSE
custom implementation in JAX (Bradbury et al. 2018). In all experiments,
use the Variance Preserving SDE (Song et al. 2021) using
f(t) -05 (Bmin + t (Amax Smin )m
g v //Bin + + (Bmax Pmin)
We Amin 0.1, and 3max 10 for all experiments. Both for the time interval [10 1.]. The associated conditional means and variances can be derived from this SDE (Song et al. 2021).
For the score estimation network, use a 5-layer MLP with hidden dimension of 50 and GELU activations. The diffusion time is embedded using random Fourier embedding. Precondition to the scoring network by performing time-dependent z- scoring (Karras al. 2022). We the denoising score matching loss with weighting function in (Song et al. (2021).
We use an AdamW optimizer with a learning rate of 5 10-4 with a cosine schedule and a training batch size of 1000. Similar to the SBI routine, we use early stopping, but with a maximum number of 5000 epochs.
E ADDITIONAL EXPERIMENTS
E. EXTENDED EXPERIMENTAL. COMPARISON
In this section, we extended the empirical evaluation and collected the results Table and additionally added the following baseline comparisons:
We
25

---

Published as a conference paper at ICLR 2025
NPE' has an alternative embedding network (no RNN). Inspired by the fact that fac- torized methods work on all windows (x', x(+1), we embed all such windows into 200 dim feature vectors (using MLP). All feature vectors accumulated (using summa- tion) and mapped into a final 50-dim summary statistics (using another MLP). This is also similar to the embedding networks used in the i.i.d. case (Zaheer et. 2018).
NPE50 corresponds to our standard NPE with RNN embedding net. By default, we do train NPE on sequences of length ten; here we extended this to 50 (although this means that the number of joint samples decreases accordingly, i... for 10k, only 200 for 100k 2000 to ensure an overall equivalent training budget).
NPE+ with pre-trained neural sufficient summary statistic. We follow (Chen et al. 2021 infomax learning approach using the distance correlation objective to pre-train the embed- ding net.
NPE+ with pre-trained sliced neural sufficient summary statistic (Chen et al. 2023)
NLE+: NLE with neural sufficient summary statistic (Chen et al. 2021).
NLE++ NLE with neural sliced sufficient summary statistic (Chen et al. 2023
NSE Neural posterior score estimation (Sharrock et al. 2024 Geffner et
2023 (analo-
gous to standard NPE baseline, but with score estimation).
Periodic SDE
Linear SDE
Mixture RW
Double Well
Task
T=1
T=10
T=100
T=1
T=10
T=100
T=1
T=10
T=100
T=1
T=10
T=100
NPE
0.82
0.84
0.96
0.72
0.93
0.97
0.92
0.99
1.00
0.64
0.67
0.81
NPE'
0.78
0.91
0.99
0.70
0.93
0.99
0.91
0.99
1.00
0.55
0.72
0.99
NPE50
0.80
0.93
0.98
0.75
0.93
0.97
0.92
0.99
1.00
0.58
0.68
0.85
NPE+
0.80
0.89
0.99
0.75
0.92
0.97
0.91
0.99
1.00
0.67
0.65
0.91
NPE+
0.79
0.90
0.97
0.74
0.92
0.97
0.91
0.99
1.00
0.64
0.64
0.83
NLE+
0.80
0.88
0.97
0.76
0.92
0.97
0.96
0.99
1.00
0.69
0.64
0.78
NLE++
0.75
0.93
0.98
0.76
0.93
0.97
0.95
1.00
0.55
0.60
0.81
NSE
0.77
0.92
0.97
0.68
0,93
0.97
0.99
1.00
0.60
0.59
0.79
FNLE
0.57
0.66
0.83
0.57
0.78
0.98
0.71
0.89
0.98
0.52
0.64
0.87
FNRE
0.63
0.89
0.98
0.69
0.92
0.99
0.91
0.99
1.00
0.64
0.89
0.97
FNSE
0.56
0.65
0.85
0.56
0.82
0.97
0.57
0.78
0.92
0.52
0.68
0.89
NPE
0.68
0.63
0.95
0.70
0.84
0.96
0.91
0.99
1.00
0,62
0.64
0.76
NPE
0.62
0.62
0.99
0.67
0.74
1.00
0.91
0.99
1.00
0.62
0.62
0.99
NPE50
0.81
0.83
0.87
0.71
0.92
0.97
0.91
0.99
1.00
0.63
0.69
0.76
NPE+
0.72
0.69
0.99
0.66
0.89
0.96
0.91
0.99
1.00
0.60
0.60
0.80
NPE+
0.67
0.65
0.95
0.71
0.91
0.97
0.91
0.99
1.00
0.65
0.61
0.78
NLE+
0.67
0.82
0.97
0.71
0,91
0.97
0.99
1.00
0.55
0.58
0.81
NLE++
0.73
0.8
0.96
0.70
0.92
0.97
0.92
0.99
1.00
0.55
0.58
0.79
NSE
0.71
0.65
0.96
0.66
0.77
0.93
0.99
1.00
0.54
0.54
0.77
FNLE 0.54
0.60
0.72
0.53
0.69
0.89
0.61
0.76
0.89
0.50
0.54
0.75
PNRE 0.55
0.75
0.92
0.54
0.81
0.99
0.70
0.86
0.93
0.50
0.74
0,90
FNSE 0.53
0.60
0.76
0.53
0.70
0.90
0.52
0.69
0.69
0.50
0.59
0.81
Table 3: Extended baseline comparison: Performance per method given C2ST metric for each benchmark task. The top half of the table reports results trained on 10k step simulations, and the bottom half on 100k. Best performing methods are marked in bold.
Overall, all "global" methods that aim to reduce the time series into a statically sized summary statistic do behave similarly within our evaluation suite (Table 3 4. A fundamental challenge with these approaches is sequence length generalization, i.c.. finding statistics that generalize to larger lengths than observed training, which challenging task (Ray Chowdhury & Caragea 2024 /Zhng et al.) 2022). In contrast, factorized methods do not have this problem. Instead, their main source of error is due to the accumulation local approximation emrors. This can lead to deteriorating performance and is especially visible for NRE within our evaluation. Notably, an increase in C2ST is, to some degree, expected, given that posteriors with more observations likely contract, making slight deviations casier to detect by classifier. In addition, results often plotted against the training budget within SBI benchmark (Lueckmann et al., (017) and not by observation size. We thus visualize results also in this more similar fashion (see Fig. Ab). We
26

---

Published as a conference paper at ICLR 2025
Lotka Volterra
SIR
T=1
T=10
T=100
T=1
T=10
T=100
NPE
0.88
0.96
0.99
0.98
0.99
1.00
NPE
0.84
0.96
1.00
0.97
0.99
1,00
NPE+
0.98
0.99
1.00
NPE+
0.88
0.95
0.99
0.98
1.00
1.00
NLE+
0.92
0.96
1.00
0.98
0.99
1.00
NLE+
0.85
0.96
0.99
0.97
0.99
1,00
NSE
0.84
FNLE
0.89
0.98
0.93
0.96
0,96
FNRE
0.83
0.99
0.99
0.94
0.98
1,00
FNSE
0.66
0.90
0.99
0.84
0.89
0.85
NPE
0.85
0.93
0.99
0,96
0.95
1.00
NPE
0.83
0.87
0.99
0.98
0.99
1,00
NPE+
0.87
0.86
0.99
0.96
0.95
1,00
NPE+
0.96
0.92
1,00
NLE+
0.87
0.99
0,96
0.96
1,00
NLE++
0.84
0.95
0.99
0.96
0.97
1,00
NSE
0.84
0.93
0.99
0.96
0.95
1,00
FNLE
0.60
0.78
0.96
0.93
0.95
0.97
FNRE
0.81
0.89
FNSE
0.57
0.78
0.94
0.67
0.78
0.74
Table 4: Extended baseline comparison for Lotka Volterra and SIR: Performance per method
given as C2ST metric for Lotka Volterra and SIR task. The top half of the table reports results trained 10k step simulations, and the bottom half 100k. Best performing methods are marked in bold.
Periodic SDE
Linear SDE
Mixture RW (5d)
Double well
Lotka Volterra
SIR
2a
2
2
2
2
[T=1]
2-1
2-
-
2-1
-
2-
10 10 10
10* 10 10
10 10
10
10 10 10
10 10a
10
10a 10a 10
SW [T=10]
20
20
2
2-1
10
10
10c
10 10
106
10 10
106
10a 10 106
10 10a 108
10 10a
10
2a
2
2a
2a
(T=10)
20
20
2-1
10a 10a
10c
10 10a 10c
10a Training budget
10a 10
10
10a 10c
10a
10a 10c
10a 10a 10
Training budget
Training budget
Training budget
Training budget
Training budget
Fig 6: Posterior approximation plotted against training budget:
Each column corresponds
task. Each row shows the performance in sliced Wasserstein distance for test observations generated using 1,10, and 100 transitions. contrast to the main figures, we see the x-axis now specifies the training budget, i.c., the number of transition evaluations that can be used to simulate the training data.
performed this training budget (equaling the number of transition step evaluations) of 10 105 and 10 (which corresponds to 103, 104, and 10 "full" simulation for NPE) as in Lueckmann et al. 202T).
27

---

Published as a conference paper at ICLR 2025
Lotka Volterra
SIR
02
FNSE
Emp. CDF
Emp. CDF
1
0
1 0
0
Emp. CDF
Emp. CDF
NPE
0
Fractional
Fractional rank
Fractional rank
Fractional
Fractional
Fractional rank
Transitions
T=1
T=5
T=10
T=20
T=50
T=100
Fig 7: We perform simulation- based calibration analysis (SBC) Talts et al.) 2018) for the Lotka Volterra and SIR task trained using 100k step simulations using FNSE and NPE. We used 1000 test simulations for evaluation. The shaded Grey area indicates the expected results under the true posterior. Colors indicate the number of transitions of the test simulations.
E.2 CALIBRATION ANALYSIS LOTKA VOLTERRA AND SIR
We performed simulation- calibration analysis Talts et al. 2018) for the Lotka-Volterra and SIR tasks. Consistent with our evaluation using C2ST and sW, we observe that FNSE is better calibrated than NPE across different sequence lengths (Fig. 7ab). This difference is especially pro- nounced the SIR task, where NPE completely fails to generalize to sequence lengths larger than those observed during training (T' > 10). In contrast, the effect is less prominent in the Lotka- Volterra task. We hypothesize that, due to the periodic nature of the Lotka- Volterra dynamics, the embedding network more easily identifies representations that generalize well across sequence lengths.
E.3 COMPOSITION METHODS
We compared the results of the different implemented score composition methods. Overall this showed a clear advantage of our implementations of GAUSS/JAC over FNPE (Fig. [0]. Even Langevin corrections are introduced FNPE does tend to perform worse and diverge given a larger number of observations.
The performance of both the GAUSS and JAC methods is quite similar (Fig. 10 first column), espe- cially in the Gaussian RW task. This is expected as both are equivalent if Gaussian assumptions are not violated. Notably, for samplers with fewer steps on non- GanGGussian tasks, JAC tends to perform slightly better than GAUSS on long time series (Fig. [0b, columns two onwards). The imposed reg- ularity conditions on the posterior covariances appear to resolve the numerical challenges associated with this approach, as observed by Linhart et al. (204)).
E.4 PROPOSAL DISTRIBUTION
The choice of proposal distribution can have a significant impact on a finite simulation budget. In benchmark tasks, we primarily selected simple yet reasonable distributions based our prior knowledge of the dynamics involved. However, we examine the impact of this choice more closely in the Gaussian RW and Mixture RW tasks (see Fig. 11
We set the proposal distribution to a Gaussian with a zero mean and variable variance, ranging from very narrow to quite wide. Importantly, at evaluation, we initialize the distribution point mass at zero. For a very narrow proposal, we thus anticipate strong performance at T 1. However, as T increases, performance may decline because the time series dynamics could extend beyond the
28

---

Published as a conference paper at ICLR 2025
tigh-probability support of the proposal. Conversely, excessively wide proposals may also hinder performance the values at the observation time fall outside the training set. It is crucial to note that this performance evaluation is relative to our specific evaluation pipeline. In both tasks, evaluations are conducted on random walks that consistently start from zero, thus constraining the range of attainable values given maximum number of steps.
Overall, there is indeed an optimal proposal for all tasks, which is close to, but not identical to, our intuitively chosen options. As expected, the significance of the proposal also increases with dimensionality.
It important emphasize that the proposals we utilized generally not optimal. Other tasks, such as the Periodic/Linear SDE, Double Well, SIR, and Lotka- Volterra models, exhibit high correlation among different variables, which the proposal we chose does not capture. See Appendix Sec. S.4.3.4 for some examples. Nevertheless, to ensure fair comparison with the NPE baseline, we opted for relatively simple proposals.
E.5 KOLMOGOROV FLOW
For the Kolmogorov flow, modify the score architecture. We utilize convolutional embedding network composed of four blocks, each containing a Conv2I) layer, GroupNorm, and GeLU activa- tion. The output is then processed by two-layer MLP to produce 100-dimensional embedding. This embedding is generated from the two frames passed to the score network and subsequently concatenated. On this embedding, we employ the same five-layer MLP that we for all other tasks, but we increase the hidden dimension to 400.
is important note that in the NPE approach, we would need construct an embedding network not only for individual "images" but also for complete video of arbitrary length, which poses significant challenges.
As an initial distribution, we utilize a filtered velocity field provided by jax-cfd library, with a maximum velocity and peak wave number of three. Our approach closely follows the methodology outlined by Rozet & Louppe| (2023). The simulator performs 100 solve steps per transition, with step size of 10- seconds; each transition thus emulates 10 milliseconds. To introduce stochasticity into the PDE, we add small amounts of Gaussian noise with standard deviation of 5 10-3 after each transition. We used the jax-cfd library (Kochkov et al. 202I) to solve the Navier-Stokes equations on 64 x 64 grid.
Designing an effective proposal distribution in a high-dimensional setting is challenging. Simply using the initial distribution is unlikely to yield satisfactory performance, over time. Ideally, we would like to obtain samples from p( p(x*) for a range of different t; however, this approach would require a significant simulation budget to step to each I.
the dynamics evolve
To address this challenge while still generating diverse set of simulations, we employ the following scheme:
(i) Sample an initial value x ~ p)=))
(ii) For f = O,...,
Sample 0, ~ p(0)
Perform transition 0z+1 (2" 0:)
Notably, this procedure does not violate the requirement that p(x²") is independent 0, as each transition is performed with different, independent parameter. This approach is motivated by the
fact that, at best, we would like proper samples from the associated prior predictive distribution (Appendix Sec.(c)
For 1,000 initializations, we run this procedure for 100 steps, and for 10,000 initializations, we run it for 10 steps, totaling exactly 200,000 transition emulations.
To assess the average performance across different observations and observation lengths, we con- ducted iouuulation-baeee calibration analysis (Talts et al. 2018) (see Fig. 12). This analysis eval- uates whether the true parameter ranks across posterior samples are uniformly distributed, as would be expected under the true posterior distribution. In other words, examines whether the estimated
29

---

Published as a conference paper at ICLR 2025
posteriors adequately cover the true parameters. The results indicate generally good, though not perfect, calibration the estimated posteriors. However, the calibration tends to deteriorate with longer observation lengths. In this scenario, the posterior becomes very narrow, making the metric sensitive to such narrow distributions (which increases the likelihood missing the true posterior). Furthermore, this suggests that for certain observations, the approximate posterior may be biased. Given the limited amount of training data, such outcomes are to be expected.
E.6 TIME-INHOMOGENOUS GAUSSIAN RANDOM WORK
To evaluate the performance of the time-inhomogeneous extension described Appendix A.2.1 we consider the Time- Inhommeenoouu Gaussian Random Walk (TI Gaussian RW)
1 1+1) a: + for <~ I), 0.9, an 0 E Rd a +
E++1
with æ, multiplied by the time-inhomogeneous coefficient (a (+++ We demonstrate the scal- ability of our methods with dimension d ranging from 1,2, to 10, and the number of observation transitions ranging from 1, 10, to 100. As a comparison metric, we use C2ST. The simulation bud- is fixed at 10k and 100k, consistent with the problem settings in the main section. Note that the training data include the time variable I, and therefore, the training data for estimating the transi- tion distribution p( +1 x O) reduced to 100 compared to the time- meme-hoogeneouss setting in the main section. The maximum time length Tmax in training set to 100.
The experimental results presented in Fig. 13 perform the NPE baseline overall. In particular, it has been confirmed that FNSE scales better for longer observation times compared to other methods.
We can confirm that the proposed methods out-
30

---

Published as a conference paper at ICLR 2025
Gaussian RW (1d)
Gaussian RW (2d)
Gaussian (10d)
Gauss RW [2d]
10a
C2ST
C2ST
Sampling
Mixture RW [5d]
10a
C2ST
10
10
10a
10-1
0.5
Transitions
Transitions
Transitions
50
100
500
NPE
FNLE
FNRE
FNSE
Transitions GAUSS
FNPE
JAC
Periodic SDE
Linear SDE
Mixture RW (5d)
Double well
d
(10d)
1.0
10
100
10k sim. C2ST
C2ST
0.5
Mixture RW (5d)
1.0
1.0
100k sim.
C2ST
102
10a
10a
10
10a
10
10a
10
100
10a
Transilions
Transitions
Transitions
Transilions
Proposal std.
NPE
FNLE
FNRE
FNSE
Fig 8: Benchmark with C2ST metrie: We validate our method on Gaussian Random Walk with different dimensions for different lengths (i.e. Transitions), also tracking sampling times (a). We assess FNSE score accumulation over Gaussian RW and Periodic SDE tasks using fixed Eu- ler-Maruyama sampler (b). We compare methods across tasks and transition steps (c). Finally, we examine the effect of the proposal NFSE trained with 10k simulations from normal proposal of varying standard deviation (d).
31

---

Published as a conference paper at ICLR 2025
Periodic SDE
Linear SDE
Mixture RW (5d)
Double
50k sim.
C2ST
1.0
500k sim. C2ST
0.5
10
10
10a
10
10'
10a
10
10
10a
10"
10'
102
Transitions
Transitions
Transitions
Transitions
b
Perlodic
Linear SDE
Mixture (5d)
Double well
sW; (500k sim.) sW; (50k sim.)
2
2
2
2/3
2:1
2D
2
2/3
2:3
10?
10a
10a
102
10a
10
107
10a
Transitions
Transitions
Transitions
Transitions
NPE
PFNLE
PFNRE
PFNSE
Fig 9: Benchmark partially factorized methods: We show the benchmark performance using
partially factorized methods PFNLE, PENRE, and PFNSE. Each is trained with 5 steps. The perfor- mance shown with respect C2ST (a) and sliced Wasserstein distance (b).
32

---

Published as a conference paper at ICLR 2025
Gaussian RW
Periodic SDE
Mixture RW (2d)
Mixture RW (5d)
2-a
2
2-1
2
EMBRO SWr
2
2-3
2
Langevingg 22
2-1
2
2-
2
10
50
100
800
10
50
100
500
10
50
100
500
5
10
50
100
500
Gaussian (2d)
Perlodic
Mixture (2d)
Mixture RW (5d)
EM1100)
2
2
2-3
2²
Langeving
2-4
2:3
2-3
2
500
100 500
Transilions
Transitons
Transitions
Transitions
FNPE
GAUSS
Performance of score estimators trained with 10k simulations using different diffusion samplers (default and fewer steps) and score composition methods. Each column corresponds to a different task, and each row represents either a standard Euler-Maruyana diffusion sampler one equipped with Langevin corrector, both performing steps. The analysis is conducted for two time discretizations: (a) 500 steps and (b) 100 steps.
Fig 10: Benchmark score composition methods:
Gaussian RW (1d)
Gaussian RW (2d)
Gaussian RW (10d)
Mixture RW (2d)
Mixture RW (5d)
1.0
1.0
1.0
1.0
1.0
C2ST
0.5
0.5
0.5
0.5
2
2
2-1
2
2-1
2-1
2-3
2
2
2-3
2
10 10
10d
10a
103
10
1
10
10a
10a
10
1
100
10
102
10a
10"
10
10
10a
10a
100
10
10a
10a
Proposal
Proposal
Proposal
Proposal
Proposal
10 100
Fig 11: Benchmark proposal sensitivity: Each column represents a different task. The top row displays results using the C2ST metric, while the bottom row shows the sliced Wasserstein distance. FNSE models are trained with a fixed budget of 10k simulations, with each using different proposal distributions characterized by a mean of and varying standard deviations (x-axis).
33

---

Published as a conference paper at ICLR 2025
Density
Reynolds Re
1.0
T=1
empirical CDF
T=10
T=20
0.5
T=50
T=100
0.0
500
1000
500
1000
posterior rank
posterior rank
Fig 12: Simulation-based calibration Kolmogorov flow: The simulation-basee calibration results present the empirical cumulative density functions (CDF) of the ground-truth parameters, ranked according to the inferred posteriors derived from 100 different observations. well-calibrated pos- terior should exhibit uniformly distributed ranks, highlighted by the shaded gray area. Repeated for data emulated for T=1,10,20,50, steps).
Gaussian RW (1d)
TI Gaussian RW (2d)
TI Gaussian RW (10d)
1.0
10k sim.
C2ST
0.5
1.0
100k sim.
C2ST
0.5
100
101
102
100
10
102
10
10
102
Transitions
Transitions
Transitions
NPE
FNLE
FNRE
FNSE
Fig 13: Time inhomogenous Gaussian random walk: We train all methods on 10k (top row) and 100k (bottom row) simulation the time-inhomogenous Gaussian RW talk. All methods have also been adapted to input the current time (as given by integers ranging from zero to 100).
34


Reviews:
[{"review_id": "KgLgnv38ci", "reviewer": "Reviewer_bzn4", "summary": "This paper proposes a simulation-based inference (SBI) method for state-space models where the transition dynamics is Markovian. The core idea is to simulate many single-state transitions and then aggregate them instead of simulating entire trajectories of time-series, so as to reduce the total number of simulator calls.", "questions": "* Is there a reason why comparison with NSE is not done? \n* Can FNSE be also made sequential in a similar manner as other simulation-efficient SBI methods such as SNPE/SNLE/SNRE? \n* Which score composition method is used in the experiments when implementing FNSE? (I might have missed if this is mentioned in the appendix)", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The paper addresses the relevant problem of computational complexity in SBI which arises when the simulator is costly to sample from. \n\nThe experimental evaluation seems thorough, and the results are significant for FNSE.", "weaknesses": "My main concern is related to the novelty of the paper. The main technical contribution lies in deriving the factorized NSE method, which as the authors mention, is an extension of the setting of the Geffner et al. (2023) paper to the Markov setting. An in-depth analysis/discussion around selecting the proposal $\\tilde{p}(\\mathbf x^t)$ and providing recommendations on its choice would have strengthened the paper significantly and added to the technical contribution. The paragraph discussing the proposal is kind of confusing. The authors say that \"there is lots of flexibility in the choice of the proposal...\", but also that \"...designing a good proposal can be challenging\". I appreciate the experiments evaluating the sensitivity of the results to the choice of proposal, but it is not clear how the proposals are chosen in the first place.\n\nThe clarity of writing can be improved by providing examples, especially when it comes to motivating the setting, the core idea, and the assumptions. For instance, the authors say \"...given sufficiently many single-step transitions, it should be possible to infer the global target...\". A discussion around the cases in which this does (and does not) hold would help the reader understand the scope of this work. Another example is the requirement that \"...the proposed state must be independent of the parameters involved in the current state transition...\", where it is not clear how restrictive this requirement is, what kind of cases does it cover, etc. Similarly, some examples of expensive real-world scientific simulators in para 3 of the intro would be nice to have.\n\nOn the topic of clarity, I found Section 3.2.2 difficult to follow. In particular, it is not clear when the authors are talking about background of previous work (Geffner 2023 and Linhart 2024), and which part is their contribution. Perhaps some of the discussion regarding the implementation details can be moved elsewhere, as they are not easy to follow without reading those two works anyway.  \n\nThe paper is presented as a sample-efficient method for Markovian time-series simulators. However, the experiments measure the performance as a function of the number of transitions. To conclude that the proposed method is \"simulation-efficient\", I would have expected performance to be plotted as a function the number of calls to the simulator (or computational cost) for a fixed number of transitions.  \n\nWould be nice to have a small discussion on other sample-efficient SBI methods in the related works (e.g. based on sequential sampling, Bayesian optimization, etc.).\n\nMinor comment: Typo in line 125-126 (the word \"approaches\")", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wCQuR2qMmy", "reviewer": "Reviewer_g3YX", "summary": "The paper applies compositional score-matching, previously introduced for SBI on exchangeable models, to non-exchangeable, Markovian models. The appealing aspect of the methodology is that a score estimator can provide good posterior estimates without simulating the full process for each batch instance, but instead simulating single transitions per batch instance and learning to aggregate across time. The hope is that one can achieve the same posterior accuracy with less than $N * T$ evaluations, where $N$ is the simulation budget and $T$ is the maximum time horizon of the simulator.", "questions": "Why does sW seem to, counterintuitively, increase for increasing number of transitions (Figure 2)?\nWill the method be robust to missing or corrupted data?", "limitations": "", "rating": 5, "confidence": 5, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The paper is self-contained, easy to follow, and it attempts to address a challenging problem that is relevant to the field as a whole. \n- The method is applicable to different classes of SBI methods (e.g., likelihood approximation, likelihood ratio approximation, direct posterior estimation).\n- The idea is original and can stimulate further research into efficient SBI on dynamic models, especially when simulation budgets are scarce.", "weaknesses": "- As far as I understand it, the proposed method is a straightforward extension of FNPE with an additional input to the score estimator (Eq.6 in [1]). As such, the claim that a new “general SBI framework” is proposed requires some calibration. In contrast, the authors could highlight and extend the empirical aspects of the work.\n\n- While the basic idea is rather appealing, I find it a bit unconvincing that FNPE can robustly approximate the correct global posteriors with sparse training (the same goes for FNLE). There is an aspect of the methodology that strikes me as magical: suppose that most information about $\\theta$ in a time-varying signal is contained in later time segments (e.g., as in certain non-stationary signals) or subject to latent transitions (e.g., as in HMMs). It would be incredibly hard to get to that information if the simulation is not evolved for longer $T$, but the factorized approach is supposed to somehow get enough training signal from very sparse simulations in all relevant time windows. How is it possible for the score network to properly learn the correct composition without any assumptions on signal stationarity or smoothness? I assume the reasonable performance on the toy examples can be attributed to the rather short signal lengths and the use of simple, low-dimensional models (Fig.9 reveals striking miscalibration for longer time series on the only challenging model, confirming my fears). \n\n- Overall, the evaluation lacks robustness and I am concerned that it is subject to randomness given the extremely limited number of test simulations used (10!). Since the work is explicitly situated in an amortized inference context and pitched as such, the evaluation could have been much more comprehensive, featuring hundreds or even thousands of test simulations and at least some practically relevant metrics with an absolute interpretation, such as calibration error or other measures used for evaluating computational faithfulness in Bayesian analysis ([2]). In addition, there are no ablation studies (e.g., one such study could vary the simulation budget or use an adaptive solver for better-than-uniform sampling, another study can consider extremely long time horizons $T$), all models only have a few parameters (SIR and LV are toy models from the 60s, but presented as distinct from the toy model section), and it is hard to tell if the performance on the only non-trivial model (4.4) is acceptable given the lack of ground-truth and the poor calibration for not not even that long time horizons ($T=100$).\n\n- It may be helpful to clarify some points in Section 2.2 regarding related work on amortized inference. It might strengthen the section to highlight that amortized Bayesian inference is explicitly introduced in [3] and extensively discussed and expanded upon in inference compilation methods [4, 5, along with additional references therein]. Including this line of research could offer a more comprehensive overview of the field. Additionally, [6, 7] don't seem to introduce, discuss, or validate amortized methods; instead, they focus on sequential likelihood-free techniques, such as SNPE [7], without suggesting or claiming to use amortization. In fact, [6] cursorily note that learning over the prior predictive is possible but ultimately dismiss it as “grossly inefficient” (p.3). \n\n[1] Geffner, T., Papamakarios, G., & Mnih, A. (2023, July). Compositional score modeling for simulation-based inference. In International Conference on Machine Learning (pp. 11098-11116). PMLR.\n\n[2] Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., ... & Modrák, M. (2020). Bayesian workflow. arXiv preprint arXiv:2011.01808.\n\n[3] Gershman, S., & Goodman, N. (2014). Amortized inference in probabilistic reasoning. In Proceedings of the Annual Meeting of the Cognitive Science Society (Vol. 36, No. 36).\n\n[4] Le, T. A., Baydin, A. G., & Wood, F. (2017). Inference compilation and universal probabilistic programming. In Artificial Intelligence and Statistics (pp. 1338-1348). PMLR.\n\n[5] Wu, M., Choi, K., Goodman, N., & Ermon, S. (2020). Meta-amortized variational inference and learning. In Proceedings of the AAAI Conference on Artificial Intelligence.\n\n[6] Papamakarios, G., & Murray, I. (2016). Fast ε-free inference of simulation models with bayesian conditional density estimation. Advances in Neural Information Processing Systems, 29.\n\n[7] Lueckmann, J. M., Goncalves, P. J., Bassetto, G., Öcal, K., Nonnenmacher, M., & Macke, J. H. (2017). Flexible statistical inference for mechanistic models of neural dynamics. Advances in Neural Information Processing Systems, 30.\n\n\n*Minor*\nSome attention to typos and neologisms is needed (e.g., P2L107 -> one can sample, …, to its probability density?, P7L377, etc.).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "m6g7XAJqyI", "reviewer": "Reviewer_Z2NW", "summary": "The paper proposes a class of methods to perform neural simulation-based inference on time-series data when the forward model has Markovian structure, i.e. when the next step only depends on the current one. This is the case in many physical systems defined by PDEs, for example. The method exploits this structure of the simulator / forward model to learn locally (in timestep) amortized estimators for global parameters of interest that describe the system in question. The paper exploits many recent advances in amortized simulation-based inference and adapts them to the Markovian time-series setting.", "questions": "- How sensitive is the more complex Kolmogorov flow example to the choice of proposal distribution?\n- Beyond the empirical results demonstrated, do the authors have a deeper understanding of the robustness of score composition rules?\n- The generalization beyond the number of timesteps seen is one of the most useful properties of the method -- how do the generalization capabilities depend on the number of timesteps used during training (e.g., max 10 vs 100 timesteps during training, then 1000 during evaluations)?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- The paper is exceptionally well-written. Secs. 2 and 3 set up the problem well, and describe the contributions in the context of existing literature. I feel like I learned a lot about the surrounding field beyond the specific contributions of the paper by reading these sections. Limitations are discussed upfront, and some extension settings explored in an initial way.\n- The experiments in the paper are well-motivated. The paper looks at toy models, standard benchmarks, as well as more challenging simulators that highlight some of the specific advantages of the method (e.g., generalization beyond the number of steps used during training).\n- The paper explores a timely question of broad scientific as well practical real-world interest. The generalization capabilities of the method beyond the number of timesteps trained for seem especially practical (and maybe should be highlighted more) -- many applications are bottlenecked by the inability to practically train on the longer-horizon timescales needed during deployment.", "weaknesses": "While the core methods are described nicely, some practical aspects receive less thorough treatment. Specifically, the proposal distribution seems like a critical design choice (especially for complex problems like Kolmogorov flow), and the score composition rules prove surprisingly robust even when their assumptions are violated (e.g. GAUSS in non-Gaussian setting). There is limited discussion of why, and limited principled guidance on these implementation choices that practitioners might need to apply the method to new domains. A more thorough discussion here could significantly strengthen the paper in particular for practitioners across domains.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Hlycc2ioR0", "reviewer": "Reviewer_bjS2", "summary": "The authors introduce a novel approach for simulation-based inference (SBI) for time series models of which the joint probability density is Markovian. The work builds on the recently introduced framework of score-based generative models for SBI which aims to approximate the score of the posterior distribution when multiple data points (in this case time points) are conditioned on. They show that their approach is beneficial over a state-of-the-art method in time series benchmarks.", "questions": "- Is there any reason why NPE performs so badly in the C2ST statistic in Figs 3c and 3d? Previously literature, e.g., [1-2], seems to report different results.\n- Could you kindly explain the exact experimental setup in Fig2a and Fig2c? Does, e.g., $100$ transitions mean that NPE received $100$ data samples of length $100$ while for $10$ transitions it received $1000$ samples?\n\n#### References\n- [1] Flow Matching for Scalable Simulation-Based Inference, Figure 4, https://arxiv.org/pdf/2305.17161\n- [2] Benchmarking Simulation-Based Inference, Figure 2, https://arxiv.org/abs/2101.04653", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The paper proposes an intuitive approach for SBI for (very) high-dimensional time series.\n- The topic (SBI for high-dimensional time series data) is very timely and relevant given the (seemingly) increased interested of applied researchers.\n- The paper is well written and easy to follow.", "weaknesses": "- Despite the intuitive idea, the proposed method is very incremental. As far as I can see, instead of modelling the score $s(\\theta_a|x^t)$ the paper proposes to model $s(\\theta_a|x^t, x^{t-1})$, i.e., adding one variable to the score network input, in order to make the method amendable for time series? \n- The evaluations are in my opinion not fully convincing. \n    - The authors chose to compare their approach to exactly one other baseline which I find a bit insufficient given the wealth of SBI methods and the fact that the authors propose approaches for local FNLE, FNRE and FNSE.\n    -  The baseline NPE uses an RNN as an embedding network. Given the authors assume Markovianity, they could have, e.g., just computed summary statistics and used this as conditioning variables, instead of running a costly RNN (which also needs to be trained in addition to the flow layers). Alternatively approaches that reduce the dimensionality more efficiently like neural sufficient statistics could have been benchmarked for a more thorough evaluation [1,3].\n  - Relatedly, there is work on high-dimensional SBI that the authors could potentially have included in their evaluations: [1-4]\n\n#### References\n- [1] Generalized massive optimal data compression, 2017, https://arxiv.org/abs/1712.00012\n- [2] Neural Approximate Sufficient Statistics for Implicit Models, 2021, https://arxiv.org/abs/2010.10079\n- [3] Is Learning Summary Statistics Necessary for Likelihood-free Inference?, 2023, https://proceedings.mlr.press/v202/chen23h.html\n- [4] Simulation-based inference using surjective sequential neural likelihood estimation, 2023, https://arxiv.org/abs/2308.01054", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "2euoR9bvSw", "reviewer": "Reviewer_Ck63", "summary": "This work proposes a set of Simulation-Based Inference (SBI) methods to perform approxiamte Bayesian (parameter) inference given an\nobservation that takes the from of a specific class of timeseries (finite-length homogenous Markov chains).\nThese methods are adaptations of existing score, likelihood, and ratio-based estimation methods, whose original formulation\nassume i.i.d observations, and not observations arising from a Markov Chain.\n\nThe methods do so by estimating ratio/likelihood/scores \"locally\", e.g. to the individual factors of the Markov chain.\nThey then produce a \"global\" posterior estimate (e.g. an estimate of the posterior of the parameters given the entire Markov chain sample)\nby relying on the factorization structure of the Markov chain probability. The composition is straightforward in the case of ratio and likelihood methods,\nwhile when using scores, additional approximations are necessary due to the non-Markovian structure of the blurred posteriors estimated in the previous step.\n\nThese approximations are direct adaptations of the one already used in prior work:\n- one performs iterative Langevin-based sampling of an annealing path of distributions bridging from a standard Gaussian distribution to the true posterior, and whose intermediate distributions are precisely the SDE solutions whose scores were estimated.\n- one produces an approximation of the scores of the SDE initialized at the posterior of the parameter given the entire Markov trajectory, which can then be used to simulate trajectories from the time-reversed SDE, whose final iterate is marginally distributed according to the true posterior.\n\nThe performance of the method is investigated on a set of experiments, which include standard benchmark models and a Kolmogorov flow model with high dimensional observations. The simulators are adapted (using noise injection) to fit the formalism of stochastic inference methods.", "questions": "Is there a particular motivation behind using Sliced Wasserstein Distance for evaluation? If using it as the main criterion, it would be nice to provide an intuition as of what features of differences in distribution this distance is capturing. Otherwise, I'd suggest using the Sinkhorn Divergence instead -- with a small regularization parameter, it is essentially a direct approximation of the Wasserstein distance.\n\n**Minor**\n\n- The use of \"time-variant\" and \"time invariant\" should be replaced by the more standard \"(time) homogenous/unhomogenous\".\n- \"MCMC-based sampling corrections, such as unadjusted Langevin dynamics\". The ULA doesn't include a correction step (it is \"undajusted\"). What do the authors mean by that?\n- To produce symmetric quotation marks, consider using ``text'' instead of \"text\" in latex.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper proposes a conceptually simple, attractive adaptation of SBI methods to Markov chains. For a total simulation budget, the shift towards learning the transition factor allows to increase the number of available observations. Moreover, learning the transition probabilities may be a simpler problem than learning the parameter-to-whole-timeseries relationship, as the observation is lower dimensional in the former case.", "weaknesses": "**Regarding Presentation**\n\nThe presentation of the method needs some improvements. \n\nFirst, regarding the the sampling part, the presentation of FNPE and Gauss do not define precisely key details of the approximations.\n- In FNPE, where is $q_a$ used for instance?\n- In GAUSS/JAC, What is $\\Sigma_{a,t, t+1}$?\n\nI had to look at each of the separate papers that introduced such approximations in the iid case to understand what was going on.\nThe paper should be updated to define all terms properly -- most of will constitute background. Even though the adaptations made by the method are almost independent from the these terms, the paper should remain self contained as far as the core method is concerned. To save space, the FNPE approximation method could be placed in the appendix, as it is not used in the experiments.\n\nSecond, the local score estimation section could also be made clearer: for instance, this section\nstarts by breaking down the true (unblurred) posterior (Equation (3)), highlighting how knowing the score of each\nfactor is enough to get the score of the resulting posterior. However, what is actually needed for sampling is the blurred posterior,\nwhich cannot be factorized as Equation (3). Equation (3) thus ends up not being very relevant to the method.\n\n\n**Regarding experiments**: the presentation of the Kolmogorov flow experiments is missing details: for instance, it does not include any mention of stochasticity, which is confusing given that the method is about estimating probability distributions.\n\n\n**Regarding performance**\n- The main weakness of the method right now is that the number of timesteps investigated in the paper remains smaller (up to 100)\ncompared to the regime of some application of interest, like neuroscience.\n- In the Kolmogorov flow experiments, the MAE of the posterior predictive significantly increases for when using 100 instead of 10. It would be nice if the authors commented on this point in the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "G0nnCYpGZy", "reviewer": "Reviewer_MmYK", "summary": "The authors consider performing simulation-based inference for Markovian simulators. They propose two primary methods, based on score matching, and likelihood (ratio) estimation. In both cases, they utilize the Markovian structure to increase the efficiency of inference, by performing inference using data from single-state transitions, rather than requiring whole sequence simulation. In the likelihood case, the global likelihood is the product of the transition functions $p(x^{t+1}|x^t, \\theta)$ (eq. 1), which once estimated can be used with MCMC to infer the posterior. Similarly, for estimating the score, they again use the Markovian structure to factorize the score (eq. 3), and use a \"local\" score estimator to approximate the global score. In both cases, a proposal distribution is used, $\\tilde{p}(x^t)$, from which a sample is drawn, and a single transition is performed, which is used for learning.", "questions": "Is $x^t$ in the SIR example the three states $\\{S^t, I^t, R^t\\}$? Correct me if I am misunderstanding, but if this is the case, the proposal is  $\\tilde{p}(x^t) = \\text{Unif}(0, 5)$, presumably independent for for each state S, I and R. At any given time $t$, I would expect the states to be correlated (e.g. more recovered likely implies less infected and susceptible individuals). The use of an independent uniform does not capture this. Further, I presume the method is only valid  if we observe all three states?\n\nAs mentioned above: Should we aim for $\\tilde{p}(x^t)$ to resemble a prior predictive simulation for a randomly chosen $t$? Could the proposal $\\tilde{p}(x^t)$ be sequentially improved as the posterior is learned?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The method is to the best of my knowledge novel and provides a solution to utilizing the known dependency structure of Markovian simulators; knowledge which is often neglected in simulation-based inference.  The application area is an important area and handling sequence data in simulation based inference is known to be often challenging (in terms of performance and computational cost). The experiments are convincing and consider a set of familiar but interesting models.", "weaknesses": "My main criticism is that the paper should include a clearer description of what might constitute a \"good\" proposal distribution $\\tilde{p}(x^t)$, especially when it is first introduced. I further feel that the requirement for choosing this proposal is a major problem for the practical utility of the introduced methods, and the authors should better address this problem. E.g. should we aim for it to resemble a prior predictive simulation for a randomly chosen $t$? Could the proposal $\\tilde{p}(x^t)$ be sequentially improved as the posterior is learned? The use of hand-picked simple distributions is somewhat concerning for broader applicability.\n\nSome smaller issues:\n- The names of the contributed methods, at least FNSE, should be introduced earlier, for example, in the final paragraph of the introduction (FNSE is in Figure 1, but it's easy to miss the name). \n- \"except FNRE; in LV\": I don't believe LV abreviation was introduced, presumably Lotka-Volterra.\n- Text size in figures is often a little too small.\n- There is a reasonable degree of entangling of previous work in section 3.2.2, which made it hard to read. For example, the abbreviation FNPE is introduced without definition, initially implied to refer to the method by Geffner. It is then introduced and thereafter used as an adaptation of the Geffner method. It also refers to a score-based method, so the use of NPE is confusing as it implies a neural posterior estimation method, using e.g. flows.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)