Sample-Size Extrapolation in Simulation-Based Inference via Statistical Scaling

S AMPLE -S IZE E XTRAPOLATION IN S IMULATION BASED I NFERENCE VIA S TATISTICAL S CALING Shiyi Sun1

Lei Fang2

Geoff K. Nicholls1

arXiv:2610.16820v1 [stat.ML] 12 Oct 2026

A BSTRACT In simulation-based inference (SBI), neural posterior estimators trained on a bounded range of sample sizes can become inaccurate on larger datasets because their parameterizations do not specify how posterior location, uncertainty and shape must change as observations accumulate. We propose center-estimated NPE (CE-NPE), which combines a statistical center with an explicit contraction rate and learns a local center correction and residual densities in posterior-scale coordinates. The center is a statistical estimate, such as a closed-form statistic or a simulatorbased moment match, which does not require an evaluable likelihood. Our analysis explains how contraction amplifies center error and decomposes Gaussian posterior error into shape, mean and covariance contributions. On six regular benchmarks, CE-NPE gives the best reference-posterior agreement of all tested amortized estimators at every out-of-range size, up to five times the largest training size, without retraining; it also fits better within the training range on most of them. Matched controls and nonregular examples characterize the scope of these gains and identify the conditions on which they depend: an appropriate contraction rate, an accurate center, and a stable residual representation.

1

I NTRODUCTION

Simulation-based inference (SBI) targets models whose likelihood cannot be evaluated but which can be simulated (Cranmer et al., 2020). Neural posterior estimation (NPE), its most direct form, amortizes Bayesian inference by learning a conditional density from simulations and reusing it on new datasets (Papamakarios & Murray, 2016; Greenberg et al., 2019; Lueckmann et al., 2021). In deployment, however, the estimator may be asked to condition on n observations when n exceeds every sample size seen during training. Permutation-invariant encoders such as Deep Sets accommodate this change in cardinality (Zaheer et al., 2017), but the posterior must also change its location, width and shape as evidence accumulates, and providing the network with n as an input (Radev et al., 2022; Habermann et al., 2025; Mittal et al., 2025) leaves these relationships to be learned over a bounded training range. Accepting arbitrary dataset sizes is therefore an architectural property, whereas extrapolating posterior uncertainty is a statistical one. We ask how posterior estimators trained on small datasets can extend to larger sample sizes without retraining. Statistical asymptotics describe how a posterior should change with n, and therefore supply structure for this problem. Under regularity conditions, the Bernstein–von Mises theorem describes an approximately Gaussian posterior with width of order n−1/2 after appropriate centering (van der Vaart, 1998, Ch. 10). More generally, a justified contraction rate n−α , which nonregular models can change (Bochkina & Green, 2014), suggests learning in local coordinates nα (θ − cn ). Working in these coordinates removes the need for the network to infer the leading contraction rate from the training range, but it makes center accuracy critical: a location error δn becomes a shift nα δn in those coordinates, so a small fixed bias can become large relative to posterior uncertainty. Moreover, even a center that contracts at the correct rate leaves errors in covariance and residual shape, so scaling alone does not ensure calibrated inference. We propose center-estimated NPE (CE-NPE), which anchors the posterior at a statistical estimator cstat (x1:n ) and learns only a correction, multiplied by n−α . This design places the learned adjustment 1 Department of Statistics, University of Oxford, Oxford, United Kingdom. 2 Mathematical Institute, University of Oxford, Oxford, United Kingdom. Correspondence to: Shiyi Sun <shiyi.sun@spc.ox.ac.uk>.

1

(i) free conditioning: qϕ (θ | sn , log n) scale learned, center learned, shape learned

x1:n

sn = eϕ (x1:n )

(ii) structural scaling: qϕ (θ | sn , log n) = rϕ nα (θ − cn ) | sn , log n nαd scale imposed, center learned, shape learned

(iii) center-estimated (CE-NPE): cn = cstat (x1:n ) + n−α ∆ϕ (sn , log n) scale imposed, center anchored at a statistic plus a learned local correction, shape learned

Figure 1: Progressively removing what the network must extrapolate. All three constructions share a permutation-invariant encoder and differ in which ingredients of the local coordinates zn = nα (θ − cn ) are learned from the bounded training range and which are imposed: (i) free conditioning learns scale, center, and shape; (ii) structural scaling imposes the scale; (iii) CE-NPE additionally anchors the center at a statistical estimator.

on the posterior’s contraction scale and thereby reduces what must be extrapolated in parameter space; because it supplies structure that would otherwise have to be learned from simulations, it can also improve fitting within the training range. The anchor can be a closed-form statistic or a simulator-based moment estimate, so it does not require an evaluable likelihood, and the network still learns uncertainty, dependence and residual shape. Figure 1 relates CE-NPE to free conditioning and to structural scaling with a learned center, which serves as an intermediate ablation. Our experiments examine both the benefit of this parameterization and its source: six adapted regular benchmarks test reference-posterior agreement beyond the training range, matched-center controls assess what scaling adds beyond the statistical estimator, best-effort baselines assess sensitivity to training effort, and nonregular examples test the assumptions behind the local representation. Throughout, we report in-range performance alongside extrapolation results, separating errors already present within the training range from those that arise beyond it. Contributions. • A structured parameterization. We introduce CE-NPE, which combines a statistical center with an explicit contraction rate and learns the center correction and residual density in local coordinates. The center need not evaluate the likelihood and may be constructed from simple statistics or simulator access, so the estimator remains likelihood-free, and the parameterization reduces what the network must extrapolate beyond the training sample sizes (Section 3). The ingredients are classical; the contribution is to place them inside an amortized estimator and to show that its out-of-range gains are attributable to this structure rather than to network capacity. • Explaining the errors that remain after scaling. We show how posterior contraction amplifies center error and separate Gaussian posterior error into shape, mean and covariance contributions. These results explain why correct scaling and a contracting center alone do not guarantee accurate, calibrated inference (Section 4). • Testing the benefits and limits of structural extrapolation. We test whether explicit statistical structure improves amortized posterior estimation beyond the training sample-size range, under matched architectures and simulation budgets. On six regular benchmarks (Lueckmann et al., 2021), center anchoring gives the best reference agreement at every out-of-range size against free conditioning, structural scaling with a learned center, and a compositional score-based estimator (Geffner et al., 2023). Fixed-center controls isolate the contribution of scaling, while best-effort baselines assess sensitivity to training effort. Four nonregular prototypes with a faster rate, a boundary, a ridge and persistent multimodality show where the imposed structure fails, and a real-data application to 1000 Genomes windows (The 1000 Genomes Project Consortium, 2015) assesses consistency and held-out prediction without a reference posterior (Sections 5 and 6; Appendices D–G). 2

2

R ELATED WORK

Neural simulation-based inference. Neural SBI learns posterior densities, likelihoods, or likelihood ratios from simulations (Papamakarios & Murray, 2016; Greenberg et al., 2019; Papamakarios et al., 2019; Hermans et al., 2020); see Cranmer et al. (2020) for a survey. Neural posterior estimation fits a conditional density estimator to simulated parameter–data pairs, either sequentially around one observation (Greenberg et al., 2019; Deistler et al., 2022) or amortized over the prior, and its density family has moved from normalizing flows (Durkan et al., 2019) to diffusion heads (Chen et al., 2025) and to transformers that serve many conditional queries with one model (Gloeckler et al., 2024). Amortization has also been extended to misspecified simulators (Ward et al., 2022) and beyond the standard Bayesian target, to generalized Bayesian posteriors indexed by a loss temperature (Sun et al., 2026). Benchmarks (Lueckmann et al., 2021) and calibration studies (Hermans et al., 2022) evaluate these estimators at fixed sample sizes. We study a different axis: how an amortized estimator trained on a bounded range of sample sizes transfers to larger datasets, and we introduce explicit statistical structure into its dependence on sample size. Inference with variable-size datasets. Deep Sets and Set Transformer provide permutationinvariant representations of variable-size inputs (Zaheer et al., 2017; Lee et al., 2019), and Chan et al. (2018) apply exchangeable networks to likelihood-free inference from population genetic data. BayesFlow combines learned summaries with amortized posterior estimation, trains across variable numbers of observations, and shows posterior contraction empirically as observations accumulate (Radev et al., 2022). More recent amortized estimators condition transformers on sets of observations of varying size, within multilevel models (Habermann et al., 2025) and as in-context posterior estimation (Mittal et al., 2025). These architectures accommodate changing dataset size, but do not themselves prescribe a contraction rate. CE-NPE retains a set encoder and adds a specified rate and a statistical center to the posterior parameterization. Compositional inference. For conditionally independent observations, learned likelihoods or likelihood ratios can be combined across observations to perform inference at different sample sizes (Papamakarios et al., 2019; Hermans et al., 2020). Compositional score methods similarly construct posterior samplers from networks trained on individual observations or small subsets (Geffner et al., 2023; Linhart et al., 2026). These methods exploit factorization to accommodate additional observations, but require the likelihood to factorize across subsets. CE-NPE instead models the full-data posterior in centered, rescaled coordinates, relying on an accurate statistical anchor and a stable residual representation rather than on a factorization (comparison with PF-NPSE in Appendix F.6). Regression adjustments and learned summaries. ABC regression adjustment modifies simulated parameter draws to account for discrepancies between simulated and observed summaries (Beaumont et al., 2002), and semi-automatic ABC uses regression estimates of posterior means as summary statistics (Fearnhead & Prangle, 2012). CE-NPE uses an estimator for a different purpose, anchoring the coordinates in which an amortized posterior is learned, with a learned correction expressed in local coordinates and mapped to parameter space through n−α . Asymptotic structure in posterior estimation. Under regularity √ conditions, the Bernstein–von Mises theorem motivates centering and rescaling a posterior by n (van der Vaart, 1998, Ch. 10). Asymptotic reasoning also informs parametric choices for amortized posterior approximation (Maceda et al., 2026), while Frazier et al. (2024) study the statistical accuracy of neural posterior and likelihood estimation. Our focus is transfer from a bounded training range in sample size: a justified contraction rate defines local coordinates, and we analyze the remaining errors from centering, covariance and residual shape.

3

M ETHODS i.i.d.

Let θ ∈ Θ ⊆ Rd have prior p(θ), and let x1:n ∼ p(· | θ), where simulation from the generative model is available although the likelihood need not be evaluable. We train amortized estimators on simulated datasets with n ≤ nmax and evaluate them at sample sizes n > nmax . Every estimator is 3

trained with the standard NPE objective, the log density assigned to the simulated parameter given the simulated dataset and its size (Papamakarios & Murray, 2016; Greenberg et al., 2019). Because the observations are exchangeable, all estimators share the same permutation-invariant Deep Sets encoder sn = eϕ (x1:n ) (Zaheer et al., 2017; Radev et al., 2022); the comparison focuses on how the posterior parameterization depends on n, given a representation that already accepts variable-size datasets. The three constructions of Figure 1 differ only in how much of that dependence is learned and how much is imposed. 3.1

F REE CONDITIONING

The most direct approach is to provide sample size as an additional conditioning variable: qϕfree (θ | x1:n , n) = qϕ (θ | sn , log n).

(1)

The logarithm only places the sample-size input on a convenient numerical scale. Apart from this explicit input, the construction is the standard one when a set encoder is combined with amortized posterior estimation (Zaheer et al., 2017; Radev et al., 2022). It imposes no constraint on how posterior location, scale, or shape evolve with n: any contraction must be inferred from the values of n seen in training, allowing flexible interpolation within the training range but requiring the learned dependence to be extended beyond it. 3.2

S TRUCTURAL SCALING

In regular parametric problems the Bernstein–von Mises theorem implies that posterior uncertainty shrinks at the leading rate n−1/2 (van der Vaart, 1998, Ch. 10). More generally, suppose the relevant contraction rate is n−α , with α = 1/2 in the regular case, and define the rescaled variable zn = nα (θ − cn ),

(2)

where cn is a data-dependent center. In the basic structural estimator the center is learned from the encoded dataset, cn = cϕ (sn , log n). Rather than modelling θ directly, we learn a conditional density rϕ for zn : zn | x1:n ∼ rϕ ( · | sn , log n). (3) The corresponding density in the original parameterization is qϕscale (θ | x1:n , n) = rϕ (nα (θ − cn ) | sn , log n) nαd .

(4)

Sample size therefore enters through the coordinate transformation and its Jacobian, as well as through the log n input to the network. The transformation imposes a leading contraction factor, while the network learns the center and residual distribution. Because a size-dependent residual can still change its own width, the explicit factor alone does not fix the fitted contraction rate. 3.3

C ENTER - ESTIMATED NPE

With the rate imposed, the center cn must still extrapolate, and in the rescaled coordinates a center error δn becomes nα δn (Section 4). CE-NPE therefore anchors the center at a problem-specific estimator cstat (x1:n ) and lets the network ∆ϕ (sn , log n) learn only a correction: cn = cstat (x1:n ) + n−α ∆ϕ (sn , log n).

(5)

The factor n−α expresses the learned correction in the same local scale as the posterior uncertainty; Proposition 1 (Appendix B) shows that with this factor a bounded correction preserves the anchor’s error order, which a correction in parameter units does not. Substituting this center into zn = nα (θ − cn ) gives zn = nα [θ − cstat (x1:n )] − ∆ϕ (sn , log n). (6) The center is only a coordinate anchor: it need not be a MAP estimate, require an evaluable likelihood, or approximate posterior uncertainty or shape. Its key requirement is to track posterior location at the contraction rate, with error Op (n−α ), so that the error remains O(1) in rescaled coordinates (Section 4). In our benchmarks, cstat ranges from simple closed-form statistics such as x̄ to EM or simulator-based moment matching (Appendix C). 4

3.4

R ESIDUAL DENSITY FAMILIES

The structured parameterization does not require the rescaled posterior to be Gaussian, and we consider two residual families. The first is a Gaussian conditional density, rϕ (zn | sn , log n) = N (zn ; µϕ (sn , log n), Σϕ (sn , log n)) ,

(7)

which provides a simple residual family and lets us assess how far Gaussianity suffices after centering and scaling. This is not Gaussian variational inference: mean and covariance are amortized functions of (sn , log n), and in parameter space the density has center cn + n−α µϕ and covariance n−2α Σϕ . The second family uses a neural spline flow (NSF) (Durkan et al., 2019), which can represent non-Gaussian residual shapes. The free heads and the CE residual receive log n as a size input, so the conditional residual can account for finite-sample changes in shape and covariance, while the deterministic factor n−α sets the contraction scale; the size inputs of the saved runs are documented in Appendix A.2.

4

W HY SCALING ALONE IS NOT ENOUGH

The factor n−α fixes only the reference contraction scale; accurate extrapolation also requires that posterior location, uncertainty and shape remain correct as n grows. 4.1

K EEPING CENTER ERROR ON THE POSTERIOR SCALE

Let mn be the true posterior mean and cn the center of the residual coordinates. If posterior width decreases as n−α , a location error cn − mn becomes a shift nα (cn − mn ) on the posterior scale: the same absolute error becomes five times larger in posterior-standard-deviation units when the posterior becomes five times narrower (Figure 2). The residual decomposes as nα (θ − cn ) = | {z }

residual to be learned

nα (θ − mn ) | {z }

variation around the true mean

+

nα (mn − cn ) | {z }

,

(8)

shift caused by center error

Keeping this shift bounded requires the center error to decrease at least as fast as the posterior width. Otherwise, the residual network must learn a location shift that grows with sample size, even after the contraction scale has been removed. true mean mn

true mean mn

error ≈ 1.1 SD

error ≈ 5.6 SD

n=N

n = 25N

Figure 2: The same location error becomes more serious as uncertainty shrinks. CE-NPE uses cn = cstat + n−α ∆n to express the learned correction on this same scale (Equation 5). Proposition 1 (Appendix B) shows that, for a fixed trained network, if cstat − mn = Op (n−α ) and ∆n = Op (1) over repeated datasets, then cn − mn = Op (n−α ). A bounded local correction therefore preserves the anchor’s error order; a correction added directly in parameter units must itself shrink at rate n−α or faster to provide the same guarantee. Neither condition is enforced by training, and a bounded local correction cannot remove a persistent anchor bias; the anchor itself must therefore contract. 4.2

S EPARATING LOCATION , COVARIANCE , AND SHAPE ERRORS

Keeping the shift bounded is necessary but not sufficient: an error that remains bounded relative to posterior width need not be small. For example, two one-dimensional Gaussians with equal variances and means separated by one standard deviation have forward KL divergence 1/2, even as both variances shrink. 5

Let Pn denote the true posterior, Qn a fitted Gaussian posterior, and Gn the Gaussian with the same mean and covariance as Pn . Under the density, entropy, and covariance conditions stated in Proposition 3 (Appendix B.4), KL(Pn ∥Qn ) = KL(Pn ∥Gn ) + | {z } non-Gaussian shape

KL(Gn ∥Qn ) | {z }

.

(9)

mean and covariance mismatch

The first term measures the error that remains even for the best Gaussian approximation (Proposition 3). The second measures errors in the fitted mean, marginal widths, and correlations. Because the learned residual can have a nonzero mean, the fitted mean need not equal cn ; controlling the center addresses only part of posterior accuracy. For a fixed dataset, free conditioning, structural scaling, and statistical anchoring do not change the best achievable approximation within an unrestricted Gaussian family; they change how its parameters are learned and extrapolated. A spline residual can reduce the Gaussian shape error, but greater flexibility does not guarantee that training achieves this. The experiments therefore assess posterior-mean error, marginal width, and joint posterior agreement separately. Both analyses hold at any fixed sample size; whether the fitted correction, covariance and residual shape remain accurate beyond the training range is an empirical question, which Section 5 addresses by reporting fit at the training boundary and beyond it separately.

5

E XPERIMENTS

5.1

E XPERIMENTAL SETUP

Estimators and protocol. We compare five estimators under matched encoder widths, optimization schedules, and training simulation budgets: Free-G and Free-NSF, which learn the posterior through free conditioning; Scale-G, which imposes the contraction scale but learns the center; and CE-G and CE-NSF, which combine the imposed scale with a statistical anchor and a learned local correction. The suffixes G and NSF denote Gaussian and spline density heads, respectively. Training sizes are sampled uniformly from n ∈ {5, . . . , 20}, and evaluation uses n ∈ {10, 20, 50, 100}, with datasets paired across estimators and seeds. Out-of-range sizes are used only for evaluation. All regular benchmarks impose the regular rate α = 1/2. When no asymptotic argument supplies the rate, it can be estimated from how the anchor contracts across the training sizes; on our tasks this estimate serves as a check on the assumed rate rather than a substitute for it (Appendix F.4). Each training step draws a fresh batch of 128 simulated datasets for up to 4000 steps, evaluation uses 400 paired datasets per size, and reported values are means over three seeds (five on the mixture). Architectural and implementation details are given in Appendix A. Tasks. We use six adapted regular benchmarks from the sbibm suite (Lueckmann et al., 2021): Gaussian linear with Gaussian and uniform priors, Gaussian mixture, Bernoulli GLM, SIR, and Lotka– Volterra. Each dataset contains n independent simulator replicates. The SIR and Lotka–Volterra priors are truncated to a regular regime in which the data remain informative about every parameter. Benchmark definitions, centers, and reference computations are given in Appendix C. Symmetric multimodal benchmarks and nonregular stress tests are reported in Appendix E. Evaluation metrics. We assess posterior location, marginal width, and joint agreement with exact or numerically computed reference posteriors. Center error measures the fitted posterior-mean error in reference-standard-deviation units, averaged over coordinates; the SD ratio targets 1. Joint 95% HPD coverage (Hermans et al., 2022) targets 0.95, and C2ST (Lopez-Paz & Oquab, 2017) targets 0.5. These complementary diagnostics distinguish improvements in location and contraction from remaining errors in joint posterior structure (Section 4.2). KL divergence, MMD (Gretton et al., 2012) and simulation-based calibration (Talts et al., 2018) appear in the appendix studies that use them. 5.2

D OES STATISTICAL STRUCTURE IMPROVE SAMPLE - SIZE EXTRAPOLATION ?

Overall performance. CE-G and CE-NSF achieve the lowest C2ST at both out-of-range sizes on all six regular benchmarks (Figure 3; full results at every size in Appendix D), and they also improve 6

Gaussian Gaussian linear Gaussian linear mixture (uniform)

1.5

10 20

2 10 20

Bernoulli GLM SIR

50 100

5 2 1 10 20

20 10 5 2 1

10 20

50 100

50 100

5 2 1

Lotka Volterra

50 100

5

1

50 20 10 5 2 1

1.0

0.5 0.2 0.1 0.05

1.2 1

center error (ref. SD)

SD ratio

10 20

10 20

50 100

50 100

sample size n Free-G

5 2 1 0.5 0.2 0.1 0.05 0.02 20 10 5 2 1 0.5 0.2 0.1 20 10 5 2 1 0.5 0.2 50 20 10 5 2 1 0.5 0.2 0.1

50 100

0.0

10 20

50 100

50 100

50 100

10 20

50 100

50 100

10 20

50 100

50 100

50 100

10 20

50 100

10 20

50 100

0.6 10 20

50 100 1.0 0.8 0.6

10 20

50 100 1.0 0.8

0.5 10 20

50 100

0.8

1.0

0.0

10 20

1.0

0.5 10 20

50 100

0.6

1.0

0.0

10 20

0.8

0.5 10 20

50 100

1.0

1.0

0.0

10 20

0.6

0.5 10 20

50 100

0.8

1.0

0.0

10 20 1.0

0.5 10 20

C2ST

0.6

1.0

0.0

1.0 0.8

0.5 10 20

2 1 0.5 0.2 0.1

95% joint coverage

0.6 10 20

50 100

sample size n sample size n Scale-G CE-G Free-NSF

sample size n CE-NSF

Figure 3: Results on the six regular benchmarks. Seed means with ±1 seed-SD bands. The shaded band is the training range; dotted lines mark the targets (1, 0.95, and 0.5). SD ratio and center error are on log axes.

in-range agreement on most tasks. The in-range gain is expected: at n between 5 and 20 the posterior is already far narrower than the prior on these benchmarks, so the anchor and the imposed scale supply the location and width that a free network must otherwise learn. The advantage does not, however, extend uniformly to calibration: at n = 100, both CE estimators remain over-dispersed on SIR, where only CE-G maintains near-nominal joint coverage, and on Lotka–Volterra their marginal widths are close to the reference yet both under-cover. Roles of scaling and centering. The intermediate ablation separates the two ingredients. Scale-G keeps its relative width approximately constant on five benchmarks and slows its growth on SIR, but because its center is still learned, its posterior-mean error in reference SD units increases with sample size on every benchmark. Anchoring the center removes this growth: CE-G improves both center error and C2ST over the Gaussian baselines at every out-of-range size. A fixed-center mixture control isolates the benefit of scaling more directly: with the same EM posterior mean, imposing the contraction factor substantially reduces C2ST and KL at n = 100, although coverage moves farther from nominal (Appendix F.3). This control also outperforms full CE-G in C2ST and KL, so on the mixture the gain comes from explicit scaling around an accurate center, and the learned location adjustment adds no further improvement. The same controls on SIR and Lotka–Volterra lead to the 7

same conclusion: the fixed-center, imposed-scale control matches or slightly improves on CE-G at every size, and the in-range advantage of CE-NSF narrows or reverses at 5× (Appendix F.3). Residual flexibility. The spline head adds a further, task-dependent benefit: CE-NSF improves on CE-G in width and C2ST on the mixture and matches it on the other benchmarks, whereas Free-NSF, which lacks the centered coordinates, remains inaccurate. Flexibility therefore helps once the coordinates are right, but does not substitute for them. 5.3

R EAL DATA : MUTATION RATE AND GROWTH FROM 1000 G ENOMES

The regular benchmarks all have reference posteriors; population genetics provides a real-data application without one. Under a neutral coalescent model with growth, we infer the per-kb mutation rate and growth rate from genomic-window summaries, using a simulator-based moment estimator as the statistical anchor. The data are summaries of non-overlapping 10 kb windows on chromosome 22 of the 1000 Genomes phase-3 panel, and a dataset of size n is a subsample of n windows. The five estimators follow the main training-size protocol. With no reference posterior, we use two complementary criteria (Figure 4): consistency with the moment fit to all windows, which shares the center’s moment functional and therefore favors CE-NPE, and posterior-predictive energy distance to held-out windows, which is independent of that target. As n grows to five times the training maximum, the posteriors of CE-G and CE-NSF contract onto the full-data moment fit, whereas Scale-G and Free-NSF contract onto a lower growth rate and Free-G barely contracts (Figure 4a); accordingly, only the CE estimators keep their consistency with the fit while the others drift away from it (Figure 4b). This stability costs nothing in predictive accuracy: on the anchor-independent criterion, their held-out predictions match the best free baseline at every size. In a separate in-model evaluation, CE has the lowest SBC KS statistics, although its joint coverage remains below nominal (Appendix G).

growth rate

(a)

Free-G

Free-NSF

Scale-G

CE-G

CE-NSF

1.0 1.5 2.0

1.0 1.5 2.0

1.0 1.5 2.0

1.0 1.5 2.0

1.0 1.5 2.0

1.5 1.0 0.5 0.0

(per kb)

(per kb)

(per kb)

(per kb)

(per kb)

n = 10 n = 20 n = 50 n = 100 full-data fit (b) consistency with full-data fit distance to fit (in posterior SD) held-out energy distance 0.75 0.50 0.25

0.15

2

0.10

1 10

20

50

sample size n

Free-G

100

10

20

50

sample size n

Scale-G

CE-G

100

10

Free-NSF

20

50

sample size n

100

CE-NSF

Figure 4: Posterior contraction on real 1000 Genomes windows. (a) The typical posterior of each estimator as n grows from 10 to 100, with the moment fit to all windows as the cross. (b) Against n: consistency, the fraction of datasets whose 95% intervals contain the full-data fit (higher is better); distance, how far the fit lies from the posterior mean in posterior SDs (lower is better); energy distance between posterior-predictive and held-out real windows (lower is better; independent of the anchor). Tables in Appendix G. 8

6

D ISCUSSION AND LIMITATIONS

Discussion. Taken together, the experiments support combining a statistical anchor with an explicit contraction scale: this structure reduces the sample-size dependence that a neural posterior estimator must learn, and it changes how posterior parameters are learned and extrapolated rather than enlarging the Gaussian approximation family (Section 4.2). The fixed-center controls clarify the source of these gains: with covariance still learned, anchoring at the statistic and imposing the scale match or improve on full CE-G in out-of-range C2ST on the mixture, SIR, and Lotka–Volterra (Appendix F.3). CE-NPE is accordingly intended for datasets comprising independent simulator replicates conditional on the parameter, when a contraction rate is justified and an affordable statistical center tracks posterior location on that scale. Alternative explanations. Because several free baselines already fit poorly within the training range, we ran three controls to check that the out-of-range gap is not an artifact of under-training, of the particular training maximum, or of the choice of comparator (Appendix F). First, best-effort free baselines, retrained on every benchmark with a tuned learning rate, three times the step budget and longer patience, narrow the in-range gap on four benchmarks but leave the out-of-range ordering unchanged on all six. Second, retraining every estimator on SIR with a four times larger training maximum preserves CE’s C2ST advantage at 5N , although absolute accuracy and coverage change. Third, the compositional comparator PF-NPSE, trained on the same number of simulated datasets, approaches CE-NSF on Gaussian linear but remains less accurate on the other regular benchmarks, with several failures already present in range. Neither training effort, training range nor composition therefore accounts for the gap. Scope and failure modes. The nonregular prototypes probe the rate, centering, and residual assumptions (Appendix E). On T1, whose posterior contracts at rate n−1 to a Gamma limit (Bochkina & Green, 2014), imposing the regular exponent produces increasing width mismatch; correcting the exponent substantially improves agreement (Table 11). T2 instead exposes a learning failure at a boundary: neither fitted residual recovers the shape, although a Gaussian renormalized on the parameter support can represent the exact posterior. A persistent ridge or separated modes, as in T3 and T4, do not become a stable local posterior under contraction about a single center. The two symmetric sbibm benchmarks, SLCP and two moons, behave in the same way: every single-center estimator fails at every size, including within the training range (Appendix E); a flexible flow can represent such distributions but must learn their changing geometry in the rescaled coordinates. Practical and empirical limitations. CE-NPE is only partially amortized: simulator-based anchors can require numerical optimization at each query, adding task-specific computation beyond the neural estimator (Appendix C). The experiments cover finite simulation budgets and extrapolation up to five times the training maximum; additional training on every benchmark and a second training maximum on SIR test this sensitivity (Appendix F). Better agreement than the tested baselines does not by itself ensure calibrated credible regions, as the regular-benchmark and in-model population-genetic results show (Appendices D and G); calibration remains a separate concern in neural SBI (Hermans et al., 2022).

7

C ONCLUSION

Accepting arbitrary dataset sizes is architectural; extrapolating posterior uncertainty is statistical. Explicitly encoding a justified posterior contraction rate removes the leading scale-extrapolation burden, a contracting statistical center controls the amplified location error, and once these first-order effects are stabilized, residual density learning targets the remaining covariance and shape errors. On six regular benchmarks CE-NPE gives the best reference agreement of the tested amortized estimators at every out-of-range size, up to five times the largest training size, without retraining; on real genomic windows CE-NPE stays consistent with the full-data fit at the same relative sizes at no cost in held-out predictive accuracy. Where the rate or the center is not justified, as in the nonregular prototypes, the same analysis identifies which ingredient fails and points toward representations adapted to the posterior geometry. 9

AI USE STATEMENT LLM-based assistants helped polish the language and assisted with experiment, analysis, and figure code. The authors verified all content and take responsibility for the manuscript, code, and results. R EPRODUCIBILITY STATEMENT Task-specific protocols, reference calculations and per-case results are provided in the appendix and accompanying code (supplementary material). Every table can be regenerated from the included per-case files. Where an evaluation uses a subset of the paired cases, the appendix states the subset.

R EFERENCES Mark A. Beaumont, Wenyang Zhang, and David J. Balding. Approximate bayesian computation in population genetics. Genetics, 162(4):2025–2035, 2002. Mark A. Beaumont, Jean-Marie Cornuet, Jean-Michel Marin, and Christian P. Robert. Adaptive approximate Bayesian computation. Biometrika, 96(4):983–990, 2009. Natalia A. Bochkina and Peter J. Green. The Bernstein–von Mises theorem and nonregular models. The Annals of Statistics, 42(5):1850 – 1878, 2014. doi: 10.1214/14-AOS1239. URL https: //doi.org/10.1214/14-AOS1239. 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. In Advances in Neural Information Processing Systems, volume 31, 2018. Tianyu Chen, Vansh Bansal, and James G. Scott. Conditional diffusions for amortized neural posterior estimation. In Proceedings of the 28th International Conference on Artificial Intelligence and Statistics, volume 258 of PMLR, pp. 2377–2385, 2025. Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117(48):30055–30062, 2020. Michael Deistler, Pedro J. Goncalves, and Jakob H. Macke. Truncated proposals for scalable and hassle-free simulation-based inference. In Advances in Neural Information Processing Systems, volume 35, 2022. Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios. Neural spline flows. In Advances in Neural Information Processing Systems, volume 32, 2019. Conor Durkan, Iain Murray, and George Papamakarios. On contrastive learning for likelihood-free inference. In International Conference on Machine Learning, pp. 2771–2781. PMLR, 2020. Stefan Elfwing, Eiji Uchibe, and Kenji Doya. Sigmoid-weighted linear units for neural network function approximation in reinforcement learning. arXiv preprint arXiv:1702.03118, 2017. URL https://arxiv.org/abs/1702.03118. Paul Fearnhead and Dennis Prangle. Constructing summary statistics for approximate bayesian computation: Semi-automatic approximate bayesian computation. Journal of the Royal Statistical Society: Series B, 74(3):419–474, 2012. David T. Frazier, Ryan Kelly, Christopher Drovandi, and David J. Warne. The statistical accuracy of neural posterior and likelihood estimation. arXiv preprint arXiv:2411.12068, 2024. URL https://arxiv.org/abs/2411.12068. Tomas Geffner, George Papamakarios, and Andriy Mnih. Compositional score modeling for simulation-based inference. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of PMLR, pp. 11098–11116, 2023. Manuel Gloeckler, Michael Deistler, Christian Weilbach, Frank Wood, and Jakob H. Macke. All-inone simulation-based inference. In Proceedings of the 41st International Conference on Machine Learning, 2024. 10

Pedro J. Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F. Podlaski, Sara A. Haddad, Tim P. Vogels, David S. Greenberg, and Jakob H. Macke. Training deep neural density estimators to identify mechanistic models of neural dynamics. eLife, 9:e56261, 2020. David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood-free inference. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of PMLR, pp. 2404–2414, 2019. Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. Journal of Machine Learning Research, 13:723–773, 2012. Daniel Habermann, Marvin Schmitt, Lars Kühmichel, Andreas Bulling, Stefan T. Radev, and PaulChristian Bürkner. Amortized bayesian multilevel models. Bayesian Analysis, 2025. Advance publication, doi:10.1214/25-BA1570. Joeri Hermans, Volodimir Begy, and Gilles Louppe. Likelihood-free mcmc with amortized approximate ratio estimators. In Proceedings of the 37th International Conference on Machine Learning, 2020. 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. William Ogilvy Kermack and Anderson G. McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London. Series A, 115(772):700–721, 1927. Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. John Frank Charles Kingman. The coalescent. Stochastic Processes and their Applications, 13(3): 235–248, 1982. Juho Lee, Yoonho Lee, Jungtaek Kim, Adam Kosiorek, Seungjin Choi, and Yee Whye Teh. Set transformer: A framework for attention-based permutation-invariant neural networks. In Proceedings of the 36th International Conference on Machine Learning, 2019. Julia Linhart, Gabriel Victorino Cardoso, Alexandre Gramfort, Sylvain Le Corff, and Pedro L. C. Rodrigues. Diffusion posterior sampling for simulation-based inference in tall data settings. Transactions on Machine Learning Research, 2026. URL https://arxiv.org/abs/2404. 07593. David Lopez-Paz and Maxime Oquab. Revisiting classifier two-sample tests. In International Conference on Learning Representations, 2017. Ilya Loshchilov and Frank Hutter. SGDR: Stochastic gradient descent with warm restarts. In International Conference on Learning Representations, 2017. URL https://arxiv.org/ abs/1608.03983. Alfred J. Lotka. Analytical note on certain rhythmic relations in organic systems. Proceedings of the National Academy of Sciences, 6(7):410–415, 1920. Jan-Matthis Lueckmann, Pedro J. Goncalves, Giacomo Bassetto, Kaan Ö.cal, Marcel Nonnenmacher, and Jakob H. Macke. Flexible statistical inference for mechanistic models of neural dynamics. In Advances in Neural Information Processing Systems, volume 30, 2017. Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke. Benchmarking simulation-based inference. In Proceedings of the 24th International Conference on Artificial Intelligence and Statistics, volume 130 of PMLR, pp. 343–351, 2021. Elliot Maceda, Emily C. Hector, Amanda Lenzi, and Brian J. Reich. Demonstrating the power and flexibility of variational assumptions for amortized neural posterior estimation in environmental applications. Environmetrics, 37(1):e70074, 2026. doi: 10.1002/env.70074. 11

Sarthak Mittal, Niels Leif Bracher, Guillaume Lajoie, Priyank Jaini, and Marcus Brubaker. Amortized in-context bayesian posterior estimation. arXiv preprint arXiv:2502.06601, 2025. George Papamakarios and Iain Murray. Fast ϵ-free inference of simulation models with bayesian conditional density estimation. In Advances in Neural Information Processing Systems, volume 29, 2016. George Papamakarios, David Sterratt, and Iain Murray. Sequential neural likelihood: Fast likelihoodfree inference with autoregressive flows. In Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, 2019. Stefan T. Radev, Ulf K. Mertens, Andreas Voss, Lynton Ardizzone, and Ullrich Köthe. Bayesflow: Learning complex stochastic models with invertible neural networks. IEEE Transactions on Neural Networks and Learning Systems, 33(4):1452–1466, 2022. Scott A. Sisson, Yanan Fan, and Mark M. Tanaka. Sequential Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences, 104(6):1760–1765, 2007. Shiyi Sun, Geoff K. Nicholls, and Jeong Eun Lee. Amortized simulation-based inference in generalized bayes via neural posterior estimation. In Proceedings of the 43rd International Conference on Machine Learning, 2026. Gábor J. Székely and Maria L. Rizzo. Energy statistics: A class of statistics based on distances. Journal of Statistical Planning and Inference, 143(8):1249–1272, 2013. doi: 10.1016/j.jspi.2013. 03.018. Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman. Validating bayesian inference algorithms with simulation-based calibration. arXiv preprint arXiv:1804.06788, 2018. The 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature, 526 (7571):68–74, 2015. Aad W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 1998. Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian M. Schmon. Robust neural posterior estimation and statistical model criticism. In Advances in Neural Information Processing Systems, volume 35, 2022. G. A. Watterson. On the number of segregating sites in genetical models without recombination. Theoretical Population Biology, 7(2):256–276, 1975. Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov, and Alexander Smola. Deep sets. In Advances in Neural Information Processing Systems, volume 30, 2017.

12

Appendix Roadmap. Appendix A gives the architecture and training protocol; Appendix B the statements and proofs behind Section 4; Appendix C the benchmark and reference definitions; Appendix D the six regular benchmarks at every size; Appendix E the nonregular stress tests; Appendix F the controls (best-effort baselines, a second training maximum, fixed-center, exponent estimation, size context, a compositional comparator); Appendix G the real-data application.

A

I MPLEMENTATION AND TRAINING PROTOCOL

The five estimators share the encoder, the optimization protocol and the training budget; they differ only in the density head and in how the center is obtained (Section 3). A.1

A RCHITECTURE

All five estimators use the same Deep Sets encoder (Zaheer et al., 2017): • Per-observation MLP: xi ∈ Rdx 7→ hi ∈ R64 , with widths dx → 64 → 64 → 64. Pn • Mean pooling: h̄ = n−1 i=1 hi . • Aggregation MLP: [h̄, tn ] 7→ sn ∈ R64 , with widths (64 + dt ) → 64 → 64. Here tn ∈ Rdt contains the sample mean, log SD, skewness and excess kurtosis of each coordinate, plus task-specific summaries for SIR and Lotka–Volterra (Appendix C). The density heads and center networks are: • Gaussian head: MLP with hidden widths (64, 64); outputs µ ∈ Rd and a lower-triangular L, with Σ = LL⊤ and Ljj = exp(clip(aj , −7, 5)). • NSF head (Durkan et al., 2019): 5 transforms, 12 bins, conditioner hidden widths (64, 64). • Center/correction: MLP with hidden widths (64, 64) and d outputs. All MLPs use SiLU (Elfwing et al., 2017) between linear layers and no output activation. Parameter coordinates. Free-G uses θ; Free-NSF uses (θ − µprior )/σprior , with coordinate-wise moments estimated from 20,000 prior draws. Scale-G, CE-G and CE-NSF use zn = nα (θ − cn ) (Equation 2) without further standardization. A.2

S IZE INPUTS

We use log n to express neural size conditioning in the main-text parameterization. Sample size also enters through the explicit contraction factor n−α and the statistical anchor computed from the observed dataset. The regular benchmarks use α = 1/2. Appendix F.5 compares log n, 1/n, and no size input to the residual network while holding the structured center and contraction factor fixed. On this Gamma diagnostic, the three choices give nearly identical contraction exponents and similar width and calibration diagnostics, with small differences in KL. This comparison concerns residual context; it does not establish equivalence of all size encodings or network components. A.3

T RAINING AND EVALUATION • Training: up to 4000 steps, each on a fresh batch of 128 simulated datasets, n drawn uniformly from {5, . . . , 20}; Adam (Kingma & Ba, 2015) with learning rate 10−3 , cosineannealed to zero over the run without restarts (Loshchilov & Hutter, 2017), and gradient norm clipped at 5. The objective is the negative log density of the simulated parameter under the estimator (Section 3).

13

• Validation: 192 fixed cases at each n ∈ {5, 10, 15, 20} every 200 steps. Validation NLL is averaged across these sizes; the best checkpoint is retained. An improvement must exceed 10−4 , and training stops after 20 validation checks without such an improvement. • Evaluation: 400 paired cases at n ∈ {10, 20, 50, 100}, identical across methods and seeds; 256 posterior draws for calibration metrics, 1024 for reference metrics (Section 5.1). • Out-of-range sizes are never used for early stopping, hyperparameter selection or post-hoc calibration. • Seeds: three on the regular benchmarks, except the mixture, which uses five for C2ST and coverage and two for the saved width and center-error summaries; the prototypes use ten. Values are means across seeds, with sample SDs where shown; table captions state exceptions. The prototypes of Appendix E use the same training range and evaluation grid with a smaller network and their own budget, stated there. Two protocol variations are reported separately: best-effort free baselines (learning-rate sweep, three times the budget; Appendix F.1) and a second training maximum N = 80 on SIR (training n ∈ {20, . . . , 80}, evaluation n ∈ {40, 80, 200, 400}; Appendix F.2).

14

B

T HEORETICAL SUPPORT

This appendix collects the statements and proofs behind Section 4. The first part concerns the center: what a contracting anchor means, why the learned correction is multiplied by the contraction factor, and how center and rate errors change with sample size. The second part concerns the fitted posterior as a whole: an exact decomposition of its error into shape, mean and covariance terms. Throughout, the network is fixed after training. B.1

N OTATION AND THE MEANING OF A CONTRACTING CENTER

Let Pn be the true posterior for a dataset of size n, with mean mn and covariance n−2α Σn , where Σn → Σ ≻ 0. The exponent α > 0 specifies the contraction rate; Σn specifies the remaining covariance. All probability statements concern repeated datasets at a fixed generating parameter θ0 after the network has been trained. We distinguish three locations. The statistical anchor is cstat (x1:n ). The residual coordinate is centered at cn = cstat + n−α ∆n , where ∆n is the network’s local correction. The mean of the final fitted posterior is an . If its residual has mean µn , then an = cstat + n−α (∆n + µn ).

(10)

Thus, the anchor and the fitted posterior mean need not coincide. Only the sum of the two learned location terms affects an ; an implementation may combine them into a single output. The condition nα (cstat − mn ) = Op (1) means that the anchor’s error remains bounded in probability after multiplying by nα . Because Σn has a positive-definite limit, this is equivalent to bounded error in posterior-standard-deviation units. It does not require that the standardized error tend to zero. B.2

W HY MULTIPLY THE CORRECTION BY THE CONTRACTION FACTOR ?

Proposition 1 (Scaling of center errors). Suppose nα (cstat − mn ) = Op (1) and the fitted local correction satisfies ∆n = Op (1). Then cn = cstat + n−α ∆n

=⇒

nα (cn − mn ) = Op (1).

If instead the correction is applied in parameter units, then cn = cstat + ∆n

=⇒

nα (cn − mn ) = nα ∆n + Op (1).

In particular, a correction converging to a nonzero constant makes the standardized center error diverge. A parameter-unit correction of order Op (n−α ) retains bounded error. Proof. Subtract mn and multiply each expression by nα . For the local correction the result is nα (cstat − mn ) + ∆n , a sum of two tight terms. For the parameter-unit correction the second term is nα ∆n , giving the stated alternatives. This is a sufficient-condition statement about the correction, not proof that a particular training procedure meets it. The result also requires a contracting anchor. If the anchor retains a fixed bias, a bounded local correction becomes too small in parameter units to cancel it. For the final mean an , the same argument applies to ∆n + µn in Equation 10. B.3

H OW FAR CAN CENTER AND RATE ERRORS GROW ?

Set N = nmax and write a deployment size as n = f N , with f ≥ 1. The following calculation assumes specified power laws over the range of interest. The width σn is a scalar reference scale; it may describe one parameter direction. It is not a substitute for checking the full covariance in a multidimensional problem. Proposition 2 (Center and rate sensitivity). Suppose ∥af N − mf N ∥ ≤ eN f −γ and σf N = σN f −α . Then ∥af N − mf N ∥ eN α−γ ≤ f . (11) σf N σN 15

If the center-error law is an equality with eN > 0 and α > γ, a tolerance τ ≥ eN /σN is reached at fcenter = (τ σN /eN )1/(α−γ) . α Separately, suppose the fitted width is σ bf N = σN f −b , so that it matches the reference at N . With η=α b − α, σ bf N σ bf N log ≤ κ ⇐⇒ |η| log f ≤ κ. (12) = f −η , σf N σf N

Proof. Divide the center-error bound by the reference width. Under equality, solve for f at the specified tolerance. For the width statement, divide the two power laws and take logarithms. Equation 11 has three cases: the bound on standardized error decreases if γ > α, stays constant if γ = α, and grows if γ < α. These describe the error itself when the assumed center-error law is an equality; a growing upper bound alone √ does not establish that the error grows. With a fixed offset and root-n contraction, the error grows as f under equality. For example, increasing sample size from N to 25N multiplies the standardized error by five, even though the absolute offset has not changed. The threshold fcenter is therefore conditional on the assumed power laws, not an extrapolation range guaranteed by training. Here the error concerns the final fitted mean an , including any residual mean, rather than the coordinate anchor alone. Equation 12 says that an exponent that is too large produces a posterior that becomes too narrow, and conversely for an exponent that is too small. The T1 prototype (Appendix E) is a direct test: its posterior contracts√at rate n−1 while every estimator imposes n−1/2 , so η = −1/2 and the predicted over-dispersion is f . The measured factors on T1 are 1.56 and 2.21 at f = 2.5 and f = 5, within 2% of the prediction for both residual families. The fitted-width assumption is essential. A size-dependent residual can change its own width and partly undo an imposed exponent; the coordinate transformation alone does not imply that the fitted width follows the stated power law. These calculations also omit other sources of posterior error. They explain sensitivity under specified conditions, rather than establishing a deployment guarantee. B.4

A N EXACT DECOMPOSITION FOR A G AUSSIAN APPROXIMATION

The true posterior Pn has mean mn and covariance n−2α Σn . The fitted Gaussian Qn has mean an and covariance n−2α An . The rate factor is shared in this notation, but An is allowed to depend on n, so the notation does not assume correct fitted contraction. Proposition 3 (Exact error decomposition). Suppose Pn has a density with finite differential entropy and positive-definite covariance, and An ≻ 0. Define Gn = N (mn , n−2α Σn ), −1/2

and set hn = nα (an − mn ) and En = An KL(Pn ∥Qn ) =

Jn = KL(Pn ∥Gn ), −1/2

Σ n An

. Then

Jn |{z}

−1 + 21 h⊤ n An hn | {z }

Gaussian shape error

mean error 1 + 2 {tr(En ) − d − log det En } .

|

{z

covariance error

(13)

}

All three terms are nonnegative. If the eigenvalues of An and En stay in fixed compact subsets of (0, ∞), there are constants 0 < c < C < ∞ such that c{∥hn ∥2 + ∥En − I∥2F } ≤ KL(Pn ∥Qn ) − Jn ≤ C{∥hn ∥2 + ∥En − I∥2F }. Proof. Add and subtract log gn inside EPn log(pn /qn ). The first term is Jn . The remaining expectation involves the log ratio of two Gaussian densities, which is quadratic in θ. Since Pn and Gn have the same first two moments, that expectation equals KL(Gn ∥Qn ). The Gaussian KL formula gives Equation 13. For the quadratic bounds, diagonalize En . Its covariance penalty is the sum of g(λ) = λ − 1 − log λ over its eigenvalues. Here g(1) = g ′ (1) = 0 and g ′′ (λ) = λ−2 . Taylor’s theorem on the stated compact interval bounds this sum above and below by constant multiples of ∥En − I∥2F . The eigenvalue bounds for An give the corresponding bounds for the mean term. 16

The first term is the error that remains even after choosing the best Gaussian mean and covariance. The second measures displacement in units of the fitted uncertainty. The third includes errors in widths, correlations and covariance orientation. None of these is removed merely by writing an explicit factor n−2α . For CE-G, an includes both learned location terms in Equation 10; substituting the statistical anchor alone would omit part of the error. If Pn is itself Gaussian, Jn = 0. For a flow approximation, KL remains invariant under an invertible coordinate change, but the Gaussian three-term formula does not apply to its general density. √ Bounded error can still be substantial. For Pn = N (0, 1/n) and Qn = N (a/ n, 1/n), the formula gives KL(Pn ∥Qn ) = a2 /2 for every n. Both distributions contract at the same rate, yet their relative displacement stays fixed. Bounded standardized mean error therefore does not imply vanishing posterior error or nominal coverage. The moment-matched Gaussian is the best Gaussian. Taking an = mn and An = Σn in Equation 13 makes the mean and covariance terms vanish, and both are strictly positive otherwise, so Gn is the unique minimizer of KL(Pn ∥ · ) over Gaussians and Jn is the Gaussian family’s irreducible error. Because an unrestricted Gaussian remains an unrestricted Gaussian after an invertible centering and scaling, Free-G, Scale-G and CE-G share this best achievable approximation at a fixed dataset; their fitted errors differ because they learn and extrapolate the mean and covariance differently. Only a more flexible residual can reduce Jn , and whether a trained spline does so is an empirical question (Section 5.2).

17

C

B ENCHMARK AND REFERENCE DEFINITIONS

C.1

G AUSSIAN LINEAR

Prior

θ ∼ N (0, 0.1 I10 )

Simulator

xi | θ ∼ N (θ, 0.1 I10 )

Dimensions

θ ∈ R10 , xi ∈ R10

Center

sample mean x̄

Reference posterior

conjugate Gaussian posterior

References

Lueckmann et al. (2021)

C.2

G AUSSIAN LINEAR ( UNIFORM )

Prior

θ ∼ U [−1, 1]10

Simulator

xi | θ ∼ N (θ, 0.1 I10 )

Dimensions

θ ∈ R10 , xi ∈ R10

Center

sample mean x̄

Reference posterior

product of truncated normals

References

Lueckmann et al. (2021)

C.3

G AUSSIAN MIXTURE

Prior

θ ∼ U [−10, 10]2

Simulator

xi | θ ∼ 21 N (θ, I) + 12 N (θ, 0.01 I)

Dimensions

θ ∈ R2 , xi ∈ R2

Center

ten EM steps from the coordinate-wise median

Reference posterior

grid integration, 201 × 201 points around the MAP

References

Sisson et al. (2007); Beaumont et al. (2009)

C.4

B ERNOULLI GLM

Inference of a ten-parameter generalized linear model with Bernoulli observations and a Gaussian prior whose covariance encourages smoothness by penalizing second-order differences in the filter; the observation is the sufficient statistic of the GLM.
Prior θ = (β, f ): β ∼ N (0, 2), f ∼ N 0, (F ⊤ F )−1 , with Fi,i−2 = 1, p Fi,i−1 = −2, Fi,i = 1+ (i − 1)/9, Fi,j = 0 otherwise, 1 ≤ i, j ≤ 9 PT Simulator xi = (xi,1 , . . . , xi,10 ), xi,1 = t=1 zt ,
1 xi,2:10 = xi,1 V z, zt ∼ Bern η(vt⊤ f + β) , η(·) = exp(·)/(1 + exp(·)), frozen input between time bins t − 8 and t: Dimensions

vt ∼ N (0, I), V = [v1 , . . . , vT ] θ ∈ R10 , xi ∈ R10

Fixed parameters

duration T = 100 time bins

Center

simulator method of moments: damped Gauss–Newton steps matching X ⊤ σ(Xθ) to the replicate mean

18

Reference posterior

adaptive importance sampling from the Newton MAP

References

Lueckmann et al. (2017); Gonçalves et al. (2020)

C.5

SLCP

A task designed to have a simple likelihood and a complex posterior: one replicate is four twodimensional points from a Gaussian whose mean and covariance are nonlinear functions of θ, so the posterior is a union of four mirror images. Prior θ ∼ U [−3, 3]5 Simulator

xi = (xi,1 , . . . , xi,4 ), xi,j ∼ N (mθ , Sθ ), where s2 ρs s
mθ = (θ1 , θ2 )⊤ , Sθ = ρs 1s s12 2 , 1 2

2

Dimensions

s1 = θ32 , s2 = θ42 , ρ = tanh θ5 θ ∈ R5 , xi ∈ R8

Center

moment estimator on the branch θ3 , θ4 > 0

Reference posterior

importance sampling on a fundamental domain

References

Papamakarios et al. (2019); Greenberg et al. (2019); Hermans et al. (2020); Durkan et al. (2020)

C.6

T WO MOONS

Simulator

θ ∼ U [−1, 1]2    √  r cos α + 0.25 −|θ1 + θ2 |/ √2 xi | θ = + , where α ∼ r sin α (−θ1 + θ2 )/ 2 U (−π/2, π/2), r ∼ N (0.1, 0.012 )

Dimensions

θ ∈ R2 , xi ∈ R2

Center

moment estimator on the branch θ1 + θ2 ≥ 0

Reference posterior

grid integration on a fundamental domain, expanded by the symmetry group

References

Greenberg et al. (2019)

Prior

C.7

SIR

An epidemiological model of susceptible S, infectious I and recovered R individuals; the task infers the contact rate β and the recovery rate γ from binomial samples of I at ten evenly spaced times. Prior β ∼ LogNormal(log 0.4, 0.5), γ ∼ LogNormal(log 0.125, 0.2),

Simulator

truncated to R0 = β/γ ∈ [1.5, 8] (about 12% of prior mass removed; below 1.5 no epidemic occurs within the horizon, above 8 the peak precedes the first observation and β is weakly identified). Inference is in (log β, log γ).
xi = (xi,1 , . . . , xi,10 ), xi,t ∼ Binomial 1000, I(t) N stored as fractions, where I is simulated from

Dimensions

dS SI dt = −β N dI SI dt = β N − γI dR dt = γI θ ∈ R2 , xi ∈ R10

19

Fixed parameters

population N = 106 , duration T (S(0), I(0), R(0)) = (N − 1, 1, 0)

Center

simulator method of moments: coarse grid, then three damped Gauss– Newton steps matching the simulator’s mean output to the replicate mean

Reference posterior

grid integration in whitened Laplace coordinates; Gaussian to high accuracy, fitted contraction exponent 0.51

References

Kermack & McKendrick (1927)

C.8

= 160, initial conditions

L OTKA –VOLTERRA

The classic model of two interacting species; the task infers four interaction parameters from both populations observed at ten evenly spaced times. Prior α ∼ LogNormal(−0.125, 0.5), β ∼ LogNormal(−3, 0.5), γ ∼ LogNormal(−0.125, 0.5), δ ∼ LogNormal(−3, 0.5),

Simulator

truncated to parameters for which neither species falls below e−5 individuals at any observation time (about 8% of prior mass removed; an effectively extinct species carries no information about its rates). Inference is in log-rate coordinates. xi,1,t ∼ LogNormal(log X(t), 0.1), xi,2,t ∼ LogNormal(log Y (t), 0.1) stored as log-counts, where X, Y are simulated from

Dimensions

dX dt = αX − βXY dY dt = −γY + δXY θ ∈ R4 , xi ∈ R20

Fixed parameters

duration T = 20 time units, initial conditions (X(0), Y (0)) = (30, 1)

Center

simulator method of moments: 134 grid, then four damped Gauss– Newton steps with multi-start

Reference posterior

adaptive importance sampling with a Student-t proposal at the Gauss– Newton MAP; Gaussian to high accuracy, fitted contraction exponent 0.50

References

Lotka (1920)

Simulator-based centers (SIR and Lotka–Volterra). Both centers are likelihood-free methodof-moments estimates computed per dataset from simulator runs only. Objective: the simulator’s mean output (the ODE solution with observation noise off) is fitted to the replicate mean by weighted least squares, with weights n over the replicates’ sample variance per output coordinate (floored at 10−3 of its mean), plus the Gaussian prior on θ as a quadratic penalty; the likelihood is never evaluated. Initialization: a coarse grid over the prior box, one batched ODE solve for the whole grid (41 × 41 points within ±3.5 prior SDs for SIR; 134 points for Lotka–Volterra), scored by the same objective. Refinement: damped Gauss–Newton with a central-difference Jacobian and step halving, three steps from the best grid point for SIR and four steps from each of the eight best grid points for Lotka–Volterra, keeping the best fit (the oscillatory mean map has local optima at the wrong phase). The same routine is applied at every n; no out-of-range data, reference posterior or tuning enters the center. The Gamma model is defined with the ablation that uses it (Appendix F).

20

D

F ULL REGULAR - BENCHMARK RESULTS

This appendix reports the six regular sbibm benchmarks at every evaluation size, in-range sizes next to the out-of-range sizes, under the unified protocol of Appendix A. The two symmetric multimodal benchmarks are in Appendix E. Table 9: Results on the six regular benchmarks at every evaluation size (three seeds, except the mixture: five for C2ST and coverage, two for SD ratio; training n ∈ {5, . . . , 20}). The two symmetric multimodal benchmarks are in Table 13. Bold marks the best entry in each column and block (closest to the target where the metric has one). C2ST estimator

95% joint coverage

SD ratio

n = 10 n = 20 n = 50 n = 100 n = 10 n = 20 n = 50 n = 100 n = 10 n = 20 n = 50 n = 100

Gaussian linear Free-G 0.57 Scale-G 0.55 CE-G 0.51 Free-NSF 0.55 CE-NSF 0.51

0.61 0.57 0.51 0.58 0.51

0.75 0.65 0.55 0.72 0.54

0.89 0.71 0.57 0.87 0.56

1.02 1.01 1.00 1.01 1.02

1.09 1.00 1.00 1.06 1.01

1.24 1.00 1.03 1.26 0.99

1.41 1.00 1.06 1.53 0.98

0.96 0.96 0.96 0.96 0.97

0.97 0.95 0.96 0.98 0.96

0.97 0.91 0.96 0.99 0.93

0.86 0.87 0.97 1.00 0.93

Gaussian linear (uniform) Free-G 0.84 0.88 Scale-G 0.66 0.67 CE-G 0.61 0.61 Free-NSF 0.58 0.62 CE-NSF 0.59 0.59

0.94 0.68 0.58 0.78 0.56

0.98 0.74 0.57 0.92 0.56

1.62 1.06 1.03 1.03 1.04

1.97 1.03 1.01 1.12 1.02

2.72 1.02 1.02 1.35 1.00

3.62 1.02 1.04 1.62 0.98

0.96 0.96 0.96 0.96 0.96

0.98 0.94 0.94 0.98 0.94

0.90 0.92 0.95 0.99 0.94

0.55 0.90 0.97 0.99 0.94

Gaussian mixture Free-G 0.87 Scale-G 0.87 CE-G 0.64 Free-NSF 0.87 CE-NSF 0.56

0.90 0.89 0.55 0.89 0.52

0.93 0.90 0.53 0.93 0.51

0.95 0.92 0.54 0.95 0.51

2.73 2.81 1.22 3.14 1.33

3.34 2.83 1.03 3.44 1.15

4.54 2.74 0.90 4.08 1.05

5.86 2.70 0.83 4.76 1.01

0.93 0.94 0.94 0.94 0.92

0.93 0.89 0.94 0.94 0.94

0.97 0.81 0.92 0.93 0.94

0.98 0.75 0.88 0.86 0.93

Bernoulli GLM Free-G 1.00 Scale-G 0.99 CE-G 0.70 Free-NSF 0.98 CE-NSF 0.67

1.00 1.00 0.67 0.99 0.64

1.00 1.00 0.66 1.00 0.63

1.00 1.00 0.68 1.00 0.65

6.56 3.80 1.05 2.24 1.05

8.88 4.04 1.04 2.97 1.04

13.98 4.10 1.07 4.42 1.05

19.49 4.07 1.10 6.02 1.08

0.93 0.94 0.92 0.94 0.95

0.94 0.82 0.94 0.95 0.95

0.93 0.24 0.95 0.96 0.95

0.88 0.06 0.97 0.91 0.95

SIR Free-G Scale-G CE-G Free-NSF CE-NSF

0.88 0.88 0.64 0.74 0.60

0.92 0.90 0.58 0.79 0.56

0.97 0.96 0.66 0.92 0.61

0.99 0.99 0.77 0.96 0.74

2.27 2.30 1.13 1.46 1.58

3.07 2.60 1.08 1.80 1.30

4.37 3.09 1.27 2.44 1.50

7.43 4.46 2.07 3.78 2.31

0.94 0.95 0.94 0.95 0.95

0.96 0.92 0.93 0.97 0.93

0.61 0.46 0.95 0.63 0.90

0.36 0.15 0.97 0.38 0.76

Lotka–Volterra Free-G 0.99 Scale-G 0.99 CE-G 0.74 Free-NSF 0.97 CE-NSF 0.71

0.99 0.99 0.71 0.98 0.66

1.00 1.00 0.67 0.99 0.66

1.00 1.00 0.74 1.00 0.76

23.60 24.42 1.17 2.93 1.57

34.30 25.25 1.11 3.84 1.43

52.05 24.07 0.99 5.43 1.19

75.06 24.11 0.90 7.45 1.03

0.94 0.95 0.94 0.95 0.97

0.96 0.87 0.93 0.97 0.94

0.70 0.63 0.89 0.90 0.84

0.26 0.33 0.79 0.56 0.66

21

E

N ONREGULAR STRESS TESTS

The four prototypes. Each prototype is a one- or two-dimensional model with an exact reference posterior, each chosen to probe one assumption of the structured representation: • T1, parameter-dependent support. xi ∼ U [0, θ], θ ∼ U [1, 3]. The posterior is ∝ θ−n on [max(Mn , 1), 3] with Mn = maxi xi and contracts at rate n−1 to a Gamma limit, the nonregular Bernstein–von Mises case of Bochkina & Green (2014). Violated: the imposed rate 1/2. • T2, boundary parameter. xi ∼ N (θ, 1) with θ ∈ [0, 3] and a truncated-exponential prior of rate 2. The posterior is TN(x̄ − 2/n, 1/n, [0, 3]): the rate 1/2 holds, but the shape is non-Gaussian whenever the mass sits against the boundary at 0. Probed: residual shape at a boundary, which the Gaussian head can represent after renormalization. • T3, singular parameterization. xi ∼ N (θ1 θ2 , 1), θ ∼ U [−2, 2]2 . The posterior concentrates on the hyperbola θ1 θ2 = x̄, a ridge that does not contract. Violated: a single local center. 2 • T4, √ persistent bimodality. xi ∼ N (θ , 1), θ ∼ U [−2, 2]. The posterior has two modes at ± x̄ whose separation does not shrink with n. Violated: a single local center. Reference posteriors are computed on composite grids covering the prior support, and every estimator’s density is renormalized on that support before scoring. The two symmetric sbibm benchmarks belong to the same family: SLCP has four mirror-image modes (as T4), and the support boundary of two moons contracts at rate n−1 (as T1). Protocol. The protocol is that of the main text (Appendix A) with three differences: ten seeds instead of three, a smaller network (two hidden layers of 64 units, 3000 steps), and a nominal coverage level of 90% with C2ST evaluated in closed form on the reference grid. The structured estimators impose the same rate α = 1/2 as on the regular benchmarks, and the center-estimated estimators use the closed-form center of each prototype (Mn for T1, max(x̄, 0) for T2, the nearest point on the ridge for T3, the positive mode for T4); for T1 the encoder additionally receives Mn as a fixed summary, since mean pooling cannot represent a maximum. The true rate (1 on T1) is used only in explicitly labeled oracle ablations. One prototype-specific quantity is added to the main-text metrics: for T3 and T4, the distance from the posterior center to the nearest ridge or mode. Table 10: Nonregular prototypes at n = 100 (5× the largest training size; ten seeds). Each block violates a different assumption of the structured representation. T1: parameter-dependent support

T2: boundary parameter

method

KL

SD ratio

center err.

C2ST

KL

SD ratio

center err.

C2ST

Free-G Scale-G CE-G Free-NSF CE-NSF

12.27 1.31 1.11 0.65 0.49

1.09 2.69 2.71 2.02 2.82

4.72 2.02 0.54 0.87 1.81

0.92 0.82 0.78 0.70 0.68

0.12 0.24 0.10 0.17 0.18

1.00 0.83 0.81 1.02 0.84

0.59 0.39 0.18 0.52 0.24

0.61 0.61 0.58 0.61 0.59

T3: singular parameterization

T4: persistent bimodality

method

KL

SD ratio

center err.

C2ST

KL

SD ratio

center err.

C2ST

Free-G Scale-G CE-G Free-NSF CE-NSF

2.34 7.65 8.40 0.67 26.81

0.84 0.36 0.36 0.95 0.36

0.01 0.01 0.40 0.11 0.50

0.90 0.95 0.92 0.70 0.89

2.28 5.01 5.67 0.14 23.48

0.85 0.39 0.39 1.03 0.39

0.01 0.00 0.43 0.01 0.60

0.86 0.93 0.90 0.61 0.79

Rate violation: T1. T1 contracts at rate n−1 while every estimator imposes n−1/2 . The estimated posterior is therefore too wide, by a factor that grows with n, consistent with the rate sensitivity of Proposition 2. Table 11 then changes only the imposed exponent: the default 1/2, the value 0.8 that the training-range estimate of Appendix F.4 supplies, and the oracle 1. Width, center error and KL 22

improve smoothly as the exponent approaches 1, while coverage stays close to nominal for the default and estimated exponents, because a too wide posterior fails gracefully, and falls slightly below it at the oracle exponent. An estimated rate therefore repairs most of the failure without being exact. Table 11: Result on T1 at n = 100. CE-G imposed α 1/2 (default) 0.8 (estimated) 1 (oracle ablation)

CE-NSF

SD ratio

center err.

cov. 90%

C2ST

KL

SD ratio

center err.

cov. 90%

C2ST

KL

2.71 1.40 0.92

0.54 0.25 0.15

0.97 0.91 0.86

0.78 0.72 0.69

1.10 0.67 0.48

2.82 1.42 0.92

1.81 0.59 0.23

0.97 0.92 0.85

0.68 0.59 0.56

0.49 0.14 0.03

Boundary shape: T2. The fixed-θ0 strata of Table 12 isolate the boundary; away from it (θ0 = 0.5) every estimator is accurate, and at it neither residual learns the truncated shape. This is a learning failure, not a representational one: every estimator’s density is renormalized on [0, 3] before scoring, and a Gaussian with mean x̄ − 2/n and variance 1/n so renormalized equals the true posterior exactly. The exact posterior’s own 90% HPD coverage of θ0 = 0 is 0.96 at n = 20 and 0.94 at n = 100, not the nominal 0.90, because the truth sits on the boundary. Table 12: Fixed-θ0 strata, n = 20 → 100 (ten seeds). T2 at the boundary (θ0 = 0), with the exact posterior as reference; T3/T4, distance from each estimator’s posterior mean to the nearest ridge or mode in units of its width, which grows with n for any single-center summary, the exact posterior mean included. Accuracy on T3 and T4 is read from Table 10. T2, θ0 = 0 stratum estimator Exact posterior Free-G Scale-G CE-G Free-NSF CE-NSF

T3

T4

SD ratio

center err.

cov. 90%

KL

ridge dist.

mode dist.

1.00→1.00 0.85→0.74 0.83→0.86 0.83→0.79 1.01→1.14 0.97→0.90

0.00→0.00 0.19→0.29 0.19→0.61 0.19→0.22 0.03→0.58 0.08→0.18

0.96→0.94 0.70→0.93 0.53→0.03 0.60→0.58 0.77→0.26 0.68→0.68

0.00→0.00 0.11→0.11 0.13→0.43 0.11→0.16 0.01→0.23 0.02→0.07

– 3.4→7.9 3.4→7.9 3.4→6.3 3.4→7.8 3.3→5.3

– 9.5→21.2 9.5→21.2 9.2→12.7 9.5→21.0 7.7→8.1

Table 13: The two symmetric multimodal sbibm benchmarks. C2ST

95% joint coverage

MMD

estimator n = 10 n = 20 n = 50 n = 100 n = 10 n = 20 n = 50 n = 100 n = 10 n = 20 n = 50 n = 100 SLCP Free-G Scale-G CE-G Free-NSF CE-NSF

0.98 0.98 0.98 0.75 0.91

0.98 0.98 0.98 0.79 0.97

0.99 1.00 0.98 0.83 0.98

0.99 1.00 0.99 0.91 0.98

0.064 0.058 0.040 0.004 0.147

0.073 0.142 0.072 0.004 0.290

0.084 0.363 0.144 0.007 0.427

0.093 0.532 0.208 0.013 0.514

0.96 0.97 0.95 0.92 0.27

0.99 0.92 0.94 0.98 0.15

0.94 0.20 0.87 0.96 0.13

0.87 0.00 0.69 0.82 0.09

Two moons Free-G Scale-G CE-G Free-NSF CE-NSF

0.98 0.98 0.97 0.80 0.80

0.98 0.98 0.98 0.81 0.81

0.99 0.99 0.99 0.90 0.94

0.99 0.99 0.99 0.97 0.98

0.129 0.109 0.112 0.014 0.016

0.141 0.223 0.195 0.012 0.011

0.165 0.475 0.347 0.018 0.119

0.197 0.660 0.432 0.050 0.221

0.96 0.97 0.98 0.96 0.96

0.99 0.97 0.95 0.99 0.96

0.98 0.81 0.57 0.89 0.40

0.77 0.09 0.39 0.41 0.11

No single local center: T3, T4, SLCP and two moons. On T3 and T4 every single-center estimator degrades with n; the free spline estimator degrades least (C2ST 0.70 and 0.61 at n = 100 against 0.79–0.95 for the others, Table 10) but is not accurate either. On the two symmetric sbibm benchmarks (Table 13) every single-center estimator fails at every size, including inside the training range; the free spline estimator retains far more of the structure (MMD an order of magnitude lower) yet its C2ST at n = 100 is still 0.91 on SLCP and 0.97 on two moons, so it does not recover the posterior either. 23

F

M ECHANISM ABLATIONS AND MATCHED CONTROLS

Each subsection answers one question about the main-text comparison: whether the free baselines are under-trained, whether the ordering depends on the training maximum, whether the imposed rate matters once the center is held fixed, whether the exponent can be estimated from the training range, whether the residual network needs the size input, and how a compositional estimator compares. F.1

B EST- EFFORT FREE BASELINES

The free baselines of the main text share the training budget of every other estimator. To check whether they lose because they are under-trained rather than because they cannot extrapolate, Free-G and Free-NSF were retrained on all six regular benchmarks under a best-effort protocol: the same architecture, prior-standardized parameters, the learning rate selected from {10−3 , 3 × 10−4 , 10−4 } by in-range validation NLL, three times the step budget and a longer early-stopping patience. No out-of-range data enter training, selection or calibration. Table 14 places the result beside the standard estimators at the training boundary and at 5×. Within the training range, the extra effort narrows the gap on four benchmarks: on both Gaussian linear tasks, the Bernoulli GLM and SIR the best-effort baselines fit markedly better at the boundary, and on Gaussian linear with a uniform prior the best-effort Free-NSF matches the CE estimators there. On the Gaussian mixture and Lotka–Volterra, by contrast, the large in-range errors remain. Beyond the range the picture is the same on every benchmark: the best-effort baselines still widen and drift at 5×, their C2ST never falls below 0.83 while CE-NSF stays at or below 0.76, and on Gaussian linear the tuned Free-G loses most of its joint coverage. Extra training effort therefore narrows the in-range gap but does not change the out-of-range ordering. F.2

A SECOND TRAINING MAXIMUM ON SIR

To check that the ordering of the main text is a statement about the relative size n/N rather than about the particular boundary N = 20, every estimator (and the two best-effort free baselines) was retrained on SIR with training sizes drawn uniformly from {20, . . . , 80} and paired evaluation at n ∈ {40, 80, 200, 400}, with fresh reference posteriors at those sizes and three seeds. Table 15 reports both training maxima on the common n/N axis. This study reuses the CE arms of the development runs, whose center is the likelihood-based Fisher-scoring estimator, at both training maxima; their N = 20 rows therefore differ slightly from the main-text SIR rows, which use the likelihood-free moment center of Appendix C. The study was not rerun with the moment center because the ordering at N = 20 is the same under either center. CE retains the best C2ST at 5N at both training maxima, but the absolute errors differ. At N = 80, CE-G has a width ratio of 1.16 and joint coverage of 0.86 at 5N , compared with 1.83 and 0.96 at N = 20; CE-NSF undercovers at 5N in both settings. The best-effort free baselines improve in-range fitting but still deteriorate out of range. Increasing the training maximum therefore preserves the observed C2ST ordering without making either accuracy or coverage invariant to N . F.3

D OES THE IMPOSED RATE MATTER AT A FIXED CENTER ?

This study compares, on the Gaussian mixture, two fixed-center controls built on the same ten-step EM center, the full CE-G model, and a standardized, tuned Free-G baseline. The two controls keep the posterior mean exactly at the EM estimate (residual mean constrained to zero, no correction network) and differ only in whether the covariance carries the imposed n−1/2 factor; CE-G is the full method; Free-G (standardized, tuned) is the free baseline with prior-standardized parameters and a learning rate selected by in-range validation. The study uses the main-text protocol (training sizes up to 20, evaluation at 20 and 100; Table 16). With the center held fixed, imposing the contraction factor improves C2ST and KL at 5× over letting the covariance scale freely, while joint coverage moves from 0.974 to 0.909 against a nominal 0.95; the rate parameterization therefore matters on its own for agreement with the reference, not for coverage. In this comparison the fixed-center, imposed-scale control has lower C2ST and KL at 5× than the full CE-G model, so the study isolates a benefit of scaling at a fixed center but does not show

24

Table 14: Best-effort free baselines beside the standard estimators at the training boundary (n = 20) and at 5× (n = 100); seed means over three seeds. Bold marks the best entry in each column and block (closest to the target where the metric has one). n = 20 (training boundary) estimator

n = 100 (5×)

SD ratio

center err.

cov.

C2ST

SD ratio

center err.

cov.

C2ST

Gaussian linear Free-G Free-G (best effort) Free-NSF Free-NSF (best effort) CE-G CE-NSF

1.09 1.01 1.06 1.03 1.00 1.01

0.11 0.07 0.11 0.08 0.05 0.04

0.97 0.96 0.98 0.96 0.96 0.96

0.61 0.54 0.58 0.54 0.51 0.51

1.41 0.98 1.53 1.42 1.06 0.98

0.59 0.81 0.53 0.44 0.12 0.10

0.86 0.39 1.00 1.00 0.97 0.93

0.89 0.93 0.87 0.83 0.57 0.56

Gaussian linear (uniform) Free-G Free-G (best effort) Free-NSF Free-NSF (best effort) CE-G CE-NSF

1.97 1.05 1.12 1.04 1.01 1.02

1.03 0.13 0.15 0.09 0.10 0.08

0.98 0.96 0.98 0.96 0.94 0.94

0.88 0.64 0.62 0.56 0.61 0.59

3.62 1.13 1.62 1.37 1.04 0.98

2.42 0.78 0.74 0.48 0.10 0.08

0.55 0.61 0.99 0.98 0.97 0.93

0.98 0.92 0.92 0.83 0.57 0.56

Gaussian mixture Free-G Free-G (best effort) Free-NSF Free-NSF (best effort) CE-G CE-NSF

3.42 3.14 3.76 3.75 1.14 1.20

2.48 2.36 2.48 2.34 0.07 0.06

0.93 0.94 0.94 0.95 0.94 0.94

0.90 0.89 0.89 0.89 0.55 0.52

5.88 5.10 5.41 6.82 0.87 1.04

3.67 3.20 4.60 3.68 0.06 0.05

0.98 0.98 0.86 0.98 0.88 0.93

0.95 0.95 0.95 0.96 0.54 0.51

Bernoulli GLM Free-G Free-G (best effort) Free-NSF Free-NSF (best effort) CE-G CE-NSF

8.88 2.84 2.97 1.71 1.04 1.04

7.00 1.95 1.94 0.88 0.13 0.13

0.94 0.95 0.95 0.96 0.94 0.95

1.00 0.99 0.99 0.92 0.67 0.64

19.49 5.84 6.02 3.17 1.10 1.08

16 4.82 4.70 3.14 0.12 0.08

0.88 0.77 0.91 0.87 0.97 0.95

1.00 1.00 1.00 1.00 0.68 0.65

SIR Free-G Free-G (best effort) Free-NSF Free-NSF (best effort) CE-G CE-NSF

3.07 1.41 1.80 1.32 1.08 1.30

2.07 0.47 0.70 0.37 0.17 0.16

0.96 0.98 0.97 0.98 0.93 0.93

0.92 0.71 0.79 0.68 0.58 0.56

7.43 2.91 3.78 2.71 2.07 2.31

18 5.27 8.07 4.80 0.41 0.33

0.36 0.47 0.38 0.44 0.97 0.76

0.99 0.94 0.96 0.94 0.77 0.74

Lotka–Volterra Free-G Free-G (best effort) Free-NSF Free-NSF (best effort) CE-G CE-NSF

34.30 29.73 3.84 2.14 1.11 1.43

22 20 1.83 1.10 0.19 0.19

0.96 0.97 0.97 0.97 0.93 0.94

0.99 0.99 0.98 0.95 0.71 0.66

75.06 65.84 7.45 3.70 0.90 1.03

47 45 9.39 3.05 0.13 0.19

0.26 0.85 0.56 0.64 0.79 0.66

1.00 1.00 1.00 0.99 0.74 0.76

25

Table 15: SIR at two training maxima against the relative size n/N . C2ST

95% joint coverage

SD ratio

N

estimator

0.5N

1N

2.5N

5N

0.5N

1N

2.5N

5N

0.5N

1N

2.5N

5N

20 20 20 20 20

Free-G Scale-G CE-G Free-NSF CE-NSF

0.88 0.88 0.57 0.74 0.55

0.92 0.90 0.55 0.79 0.53

0.97 0.96 0.64 0.92 0.61

0.99 0.99 0.76 0.96 0.75

2.27 2.30 1.02 1.46 1.47

3.07 2.60 1.00 1.80 1.26

4.37 3.09 1.14 2.44 1.77

7.43 4.46 1.83 3.78 2.40

0.94 0.95 0.94 0.95 0.94

0.96 0.92 0.92 0.97 0.93

0.61 0.46 0.95 0.63 0.93

0.36 0.15 0.96 0.38 0.70

80 80 80 80 80 80 80

Free-G Free-G (best effort) Scale-G CE-G Free-NSF Free-NSF (best effort) CE-NSF

0.96 0.78 0.95 0.58 0.89 0.76 0.54

0.97 0.82 0.96 0.58 0.91 0.80 0.55

0.99 0.92 0.98 0.64 0.97 0.91 0.64

0.99 0.96 0.99 0.72 0.99 0.96 0.72

5.30 1.40 4.36 1.04 2.69 1.43 1.60

6.79 1.84 4.51 1.03 3.50 1.79 1.47

11.63 3.05 5.90 1.02 5.95 2.69 1.94

17.81 5.07 7.79 1.16 10.05 4.20 2.59

0.95 0.94 0.95 0.94 0.96 0.94 0.94

0.97 0.97 0.92 0.94 0.97 0.97 0.94

0.73 0.73 0.61 0.92 0.71 0.72 0.92

0.56 0.58 0.31 0.86 0.47 0.52 0.74

an additional gain from the learned location adjustment: once the center is accurate and the rate is right, most of the gain is already obtained, which is the scope stated in Section 6. Table 16: Fixed-center controls on the Gaussian mixture at the training boundary (n = 20) and 5× (n = 100); means across three seeds. method

C2ST at n = 20

KL at n = 20

C2ST at n = 100

KL at n = 100

cov. at n = 100

0.898 0.552 0.574 0.558

2.317 0.121 0.174 0.137

0.942 0.537 0.623 0.513

2.980 0.096 0.328 0.033

0.966 0.870 0.974 0.909

Free-G (standardized, tuned) CE-G (log n residual) Fixed EM center, free scale Fixed EM center, imposed scale

The same two controls were run on SIR and Lotka–Volterra with the simulator-based moment center, under the main-text protocol and on the paired cases of Appendix D (Table 17, C2ST only). Imposing the contraction factor at a fixed center is again what matters: with the covariance free, C2ST rises with sample size on both tasks, whereas with the factor imposed it stays close to CE-G at every size and is slightly below it at 5× (SIR 0.750 against 0.772, Lotka–Volterra 0.704 against 0.737). The learned location adjustment therefore adds nothing measurable on these two tasks either. The spline head does: CE-NSF has the lowest C2ST at every in-range size and on SIR at every size, but at 5× on Lotka–Volterra it falls behind the fixed-center Gaussian control, and on both tasks its joint coverage at n = 100 is well below nominal (Appendix D), so the in-range benefit of residual flexibility does not carry over intact to the largest extrapolation. Table 17: Fixed-center controls on SIR and Lotka–Volterra: C2ST at every evaluation size (training sizes 5 to 20); means across three seeds, lowest per column in bold. SIR estimator Fixed center, free scale Fixed center, imposed scale CE-G CE-NSF

F.4

Lotka–Volterra

n = 10

n = 20

n = 50

n = 100

n = 10

n = 20

n = 50

n = 100

0.653 0.634 0.640 0.602

0.642 0.572 0.579 0.559

0.789 0.647 0.663 0.614

0.897 0.750 0.772 0.740

0.768 0.729 0.738 0.707

0.806 0.704 0.713 0.665

0.886 0.656 0.671 0.655

0.940 0.704 0.737 0.757

C AN THE EXPONENT BE ESTIMATED FROM THE TRAINING RANGE ?

CE-NPE imposes the contraction rate. When no asymptotic argument supplies it, the only data available at training time are simulations at the training sizes, and the natural estimate is the rate at which the statistical center itself contracts there: for each of J = 64 parameters drawn from the prior and each training size n ∈ {5, 8, 12, 16, 20} we simulate R = 32 replicate datasets, compute the spread of cstat across replicates, and fit log SD = a − α̂ log n by least squares. This needs no likelihood and no reference posterior. Table 18 reports α̂ on every benchmark next to the rate CE-NPE imposes (1/2 throughout) and the true asymptotic rate. 26

Three readings follow. On the regular simulators other than the mixture, the estimate is close to the regular rate. On T1 it gives 0.80, which is the finite-n contraction of the maximum over n ∈ [5, 20] rather than the asymptotic 1; imposing 0.80 is exactly the estimated-rate row of Table 11, which repairs most of the T1 failure, so the estimate is useful where the asymptotic rate is unknown and wrong. On the Gaussian mixture, the estimate is 1.38, the narrow mixture component is identified between n = 5 and n = 8, and the center spread contracts much faster than n−1/2 in this transient regime before settling to the regular rate. An estimator that trusted the training range would impose an exponent that over-contracts out of range. A training-range estimate is therefore a check on the model-class rate, not a substitute for it: when the two disagree, the disagreement itself indicates that the training range is not yet in the regime whose rate one wants to extrapolate. Table 18: Contraction exponent estimated from the training range. benchmark

imposed rate

true asymptotic rate

estimated α

1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2

1/2 1/2 1/2 1/2 1/2 1/2 1 1/2 ridge 1/2

0.50 0.48 1.38 0.46 0.56 0.48 0.80 0.43 0.44 0.47

Gaussian linear Gaussian linear (uniform) Gaussian mixture Bernoulli GLM SIR Lotka–Volterra T1 (support) T2 (boundary) T3 (singular) T4 (bimodal)

F.5

D OES THE RESIDUAL NETWORK STILL NEED n?

Gamma model and protocol. Training sizes are {20, 50, 100, 200, 500, 1000} and evaluation sizes are {1000, 2000, 5000, 10000}. The encoder is a Deep Sets network over the log-transformed observations. The structured center is the log of the maximum-likelihood (equivalently methodof-moments) rate estimate plus a learned correction applied in rescaled units, cn = log θ̂MLE + n−1/2 ∆ϕ (sn , log n); the residual density on zn is a one-dimensional rational-quadratic spline flow with linear tails (Durkan et al., 2019). Each diagnostic uses five seeds. √ We hold the structured center and n transformation fixed and change only whether the residual network receives log n, 1/n, or no size input. As Table 19 shows, all three variants have a fitted exponent of 0.500 and nearly identical SD ratios, coverage, center error, and Kolmogorov–Smirnov (KS) statistics for simulation-based calibration (SBC) (Talts et al., 2018) at 10×. The log n context improves KL from 0.0023 without context to 0.0018; 1/n gives 0.0012. These differences concern finite-sample density fitting, while the contraction exponents remain nearly identical. Sample size is not unnecessary: n still enters the analytic transformation and the structured center. In this Gamma ablation, explicit neural context becomes secondary once the dominant contraction law is encoded. Table 19: Matched residual-context comparison at n = 10,000 (means over five seeds; coverage is the nominal 90% region).

residual context log n 1/n none F.6

KL

SD ratio

90% cov.

center error

SBC KS

exponent

0.0018 0.0012 0.0023

0.991 0.998 0.992

0.886 0.889 0.888

0.061 0.059 0.065

0.039 0.038 0.041

0.500 0.500 0.500

A SCORE - BASED COMPOSITIONAL COMPARATOR ON THE SIX REGULAR BENCHMARKS

Compositional estimators handle a large n by learning the posterior score for a block of observations and summing block scores at test time, so every n is supported without any single estimator having

27

Table 20: Score-based compositional estimator on the six regular benchmarks. PF-NPSE (three seeds) with the published annealed Langevin sampler and with the GAUSS sampler of Linhart et al. (2026). Coverage requires a density and is not available for PF-NPSE. C2ST estimator Gaussian linear Free-NSF PF-NPSE (Langevin) PF-NPSE (GAUSS) CE-NSF

SD ratio

n = 10 n = 20 n = 50 n = 100 n = 10 n = 20 n = 50 n = 100 0.55 0.56 0.51 0.51

0.58 0.57 0.52 0.51

0.72 0.57 0.55 0.54

0.87 0.62 0.59 0.56

1.01 1.11 0.99 1.02

1.06 1.12 0.99 1.01

1.26 1.03 1.00 0.99

1.53 0.98 0.98 0.98

Gaussian linear (uniform) Free-NSF 0.58 PF-NPSE (Langevin) 0.60 PF-NPSE (GAUSS) 0.57 CE-NSF 0.59

0.62 0.61 0.59 0.59

0.78 0.64 0.68 0.56

0.92 0.72 0.74 0.56

1.03 1.10 0.94 1.04

1.12 1.14 0.96 1.02

1.35 1.06 1.00 1.00

1.62 1.04 1.04 0.98

Gaussian mixture Free-NSF PF-NPSE (Langevin) PF-NPSE (GAUSS) CE-NSF

0.87 0.94 0.93 0.56

0.89 0.95 0.94 0.52

0.93 0.96 0.95 0.51

0.95 0.96 0.96 0.51

3.14 5.65 4.24 1.33

3.44 6.06 4.27 1.15

4.08 5.62 5.19 1.05

4.76 5.83 6.10 1.01

Bernoulli GLM Free-NSF PF-NPSE (Langevin) PF-NPSE (GAUSS) CE-NSF

0.98 0.75 0.71 0.67

0.99 0.83 0.79 0.64

1.00 0.89 0.88 0.63

1.00 0.94 0.93 0.65

2.24 1.18 1.02 1.05

2.97 1.33 1.15 1.04

4.42 1.14 1.09 1.05

6.02 1.10 1.07 1.08

SIR Free-NSF PF-NPSE (Langevin) PF-NPSE (GAUSS) CE-NSF

0.74 0.78 0.65 0.60

0.79 0.83 0.68 0.56

0.92 0.84 0.75 0.61

0.96 0.91 0.81 0.74

1.46 1.70 0.98 1.58

1.80 2.09 1.09 1.30

2.44 2.19 1.32 1.50

3.78 3.76 1.53 2.31

Lotka–Volterra Free-NSF PF-NPSE (Langevin) PF-NPSE (GAUSS) CE-NSF

0.97 0.97 0.94 0.71

0.98 0.98 0.95 0.66

0.99 0.98 0.97 0.66

1.00 0.99 0.98 0.76

2.93 2.12 1.54 1.57

3.84 2.51 1.73 1.43

5.43 2.17 1.85 1.19

7.45 2.51 2.12 1.03

to extrapolate. The score network of PF-NPSE (Geffner et al., 2023) conditions on a block of up to 20 observations; at test time the n observations are partitioned into blocks of 20 plus a remainder, each block is scored by the same network given its block size, and the block scores are summed with the diffused prior score subtracted once per additional block. For conditionally independent observations, as on all six benchmarks, the factorization of the posterior into block posteriors is exact at noise level zero, but the sum of diffused block scores is an approximation at positive noise levels (Geffner et al., 2023; Linhart et al., 2026); at n ≤ 20 the dataset forms a single block and no composition takes place, so only the two out-of-range sizes exercise it. We run the authors’ public implementation under the main-text protocol: the same number of training datasets as one main-text estimator, each simulated with 20 replicates and randomly masked to a block size in {5, . . . , 20} (so about 1.6 times the replicates a main-text estimator sees), the analytic diffused prior of the authors’ code, which ignores the mild prior truncation of SIR and Lotka–Volterra, in-range validation only, three seeds, and the first 200 of the paired evaluation cases at every size. Two samplers draw from the composed score: the published annealed Langevin sampler and the GAUSS sampler of Linhart et al. (2026), which corrects the composition with a Gaussian approximation. Neither provides a density, so coverage is not reported. Table 20 places both beside Free-NSF and CE-NSF.

28

Composition extrapolates well where the block posterior is easy to learn. On Gaussian linear, PFNPSE keeps the correct width at every size and its C2ST at n = 100 is close to that of CE-NSF and far better than that of Free-NSF. On the uniform-prior variant and on SIR it lies between the free baseline and the center-estimated estimator. On the Gaussian mixture, Bernoulli GLM, and Lotka–Volterra it fails at every size, inside the training range as well as outside, as the free baselines do. The mechanism differs from that of free conditioning: on Bernoulli GLM the composed width stays within a third of the √ reference at every size while the center error, in posterior standard deviations, grows roughly as n, so the composed posterior has approximately the right width around a center whose error does not contract. The GAUSS sampler improves width and C2ST over the published sampler on most benchmarks and sizes, most clearly on SIR, without changing the ordering. Among the estimators in Table 20, CE-NSF has the best C2ST at n = 100 on all six benchmarks, and PF-NPSE is the closest alternative only on Gaussian linear. Configuration and per-case results accompany the code.

G

R EAL - DATA APPLICATION : MUTATION RATE AND GROWTH FROM 1000 G ENOMES

Population genetics supplies a real-data instance of the paper’s setting in which the center is a classical moment construction rather than one of ours, and in which neither a ground-truth parameter nor a tractable reference posterior exists: in-model evaluation is truth-based, and real-data evaluation uses only reference-free criteria. The five main-text estimators are compared under the main-text protocol; nothing else changes. G.1

TASK AND CENTER

Under the neutral coalescent (Kingman, 1982) with exponential population growth, genomic window i yields the summary triple xi = (Si , πi , ηi ) for k haplotypes: segregating sites, mean pairwise differences, and folded singletons. All three means are linear in the per-kb scaled mutation rate θ, while growth raises the rare-variant fraction η/S and lowers π/S, so a growth rate is identified jointly with θ from the moments. The growth parameter is the raw rate divided by four, which keeps the rescaled residual inside the spline head’s core interval (Appendix A). Prior Simulator

Dimensions Center

Reference posterior References G.2

θ ∼ U [0.2, 6] per kb; growth parameter ∼ U [0, 2] (raw rate U [0, 8]) coalescent tree of k = 20 under exponential growth (time-changed merging); mutations Poisson on
branches, class counts ξb , then S = P P k ξ , π = b(k − b)ξ / b b b b 2 , η = ξ1 + ξk−1 θ ∈ R2 , xi ∈ R3 simulator method of moments in the Watterson tradition (Watterson, 1975): the growth rate from the two observed ratios (π̄/S̄, η̄/S̄) against moment curves precomputed once by simulation, then θ from the scale S̄; a grid lookup with quadratic refinement at query time, likelihoodfree and contracting at n−1/2 none: the likelihood is intractable Kingman (1982); Watterson (1975)

R EAL DATA

Summaries are computed from 1000 Genomes phase-3 chromosome 22 (The 1000 Genomes Project Consortium, 2015): k = 20 haplotypes (ten YRI individuals), biallelic PASS SNVs, non-overlapping 10 kb windows with strict-mask accessibility at least 0.9, giving 404 windows (150,223 segregating sites). An evaluation dataset of size n subsamples n windows at least 60 kb apart; 400 such overlapping datasets are drawn per size with a fixed seed shared by every estimator. The real ratios (π̄/S̄ = 0.258 against the constant-size 0.282; η̄/S̄ = 0.355 against 0.297) show the rare-variant excess of a growing population. Window independence is approximate and the neutral model is misspecified for real genomes; the criteria below do not condition on the model being correct.

29

G.3

R EFERENCE - FREE CRITERIA

With no ground truth, each (estimator, n) is scored over the 400 subsampled datasets by: • Consistency: the fraction of datasets whose two marginal central 95% intervals both contain the full-data moment fit (θ = 1.70, raw growth 0.75, from all 404 windows), and the rms standardized distance to it (radius). Because the fit is estimated from windows that overlap the evaluation subsets and uses the same moment functional as the CE center, this is a consistency diagnostic that favours CE-NPE, not a test of nominal coverage. • Predictive: the energy distance (Székely & Rizzo, 2013) between posterior-predictive window summaries (one simulated window per posterior draw, draws clipped to the prior box) and held-out real windows, on standardized summaries, over the first 64 datasets per size; this criterion is independent of the anchor. • Posterior readout: the typical posterior of (θ, γ) at each n (Figure 4). G.4

R ESULTS

In-model SBC and coverage give complementary diagnostics. Using simulation-based calibration (Talts et al., 2018), we assess posterior ranks alongside joint coverage (Table 22). At 5× the largest training size the free heads are over-wide or rank-miscalibrated, Scale-G’s learned center has substantial location error, and CE-G and CE-NSF have the lowest SBC KS statistics although their joint coverage is below nominal; Free-NSF has a radius and coverage closer to their targets but a larger SBC KS statistic. No single diagnostic therefore establishes posterior calibration (Section 4.2). On the real genome CE-NPE maintains consistency with the full-data moment fit. The consistency of CE-G and CE-NSF with the full-data fit is flat from n = 10 to n = 100 while every other estimator’s decays, their standardized distance to the fit stays constant while the others’ grows, and their predictive energy distance matches the best free baseline within seed variability (Table 23). At n = 100 their posteriors sit on the full-data fit with a positive growth rate, as expected for an expanding population, whereas Scale-G halves the growth rate and shifts θ by more than a posterior SD (Figure 4). Table 22: In-model protocol, every size (seed means): rms standardized truth distance (radius; target 1), joint 95% HPD coverage, and maximum marginal SBC KS; bold: within one seed SD of the best per column. n = 10 Free-G Free-NSF Scale-G CE-G CE-NSF

n = 20

n = 50

n = 100

rad.

cov.

KS

rad.

cov.

KS

rad.

cov.

KS

rad.

cov.

KS

1.04 0.98 1.00 0.97 1.00

0.94 0.95 0.95 0.95 0.95

0.052 0.033 0.042 0.046 0.046

0.87 0.92 1.06 0.98 0.95

0.98 0.98 0.96 0.96 0.95

0.063 0.062 0.077 0.053 0.048

0.82 0.92 1.34 1.14 1.04

1.00 0.97 0.87 0.94 0.91

0.087 0.110 0.167 0.091 0.039

0.82 1.01 1.64 1.24 1.10

0.98 0.92 0.75 0.87 0.86

0.131 0.224 0.181 0.075 0.068

Table 23: 1000 Genomes windows, every size (seed means): consistency with and standardized distance to the full-data moment fit, and posterior-predictive energy distance to held-out windows; bold: within one seed SD of the best per column. n = 10 Free-G Free-NSF Scale-G CE-G CE-NSF

n = 20

n = 50

n = 100

cons.

rad.

energy

cons.

rad.

energy

cons.

rad.

energy

cons.

rad.

energy

0.79 0.80 0.83 0.83 0.83

1.32 1.28 1.21 1.23 1.18

0.162 0.158 0.158 0.153 0.160

0.81 0.79 0.77 0.81 0.83

1.32 1.52 1.44 1.32 1.28

0.132 0.123 0.129 0.122 0.121

0.77 0.72 0.58 0.85 0.86

1.32 1.63 1.87 1.22 1.22

0.094 0.090 0.095 0.087 0.088

0.66 0.58 0.28 0.83 0.85

1.59 1.92 2.79 1.23 1.28

0.085 0.077 0.085 0.079 0.082

30

Related documents

• Sample-Size Extrapolation in Simulation-Based Inference via Statistical Scaling#1341703
• Model-assisted estimation with a training subsample: a two-phase sampling approach with design-based variance estimation#926925
• Neural Estimation of Summary Statistics for Simulation-Based Inference with Non-Asymptotic Guarantees#1297383
• Two-way Guidance: Simultaneous Learning of Summary Statistics for Generative Bayes#1296857
• Spectra: Exact Component Transport for Test-Time Prior Adaptation in Simulation-Based Inference#1217756