stat.MLOct 12, 2026

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

Authors: Shiyi Sun, Lei Fang, Geoff K Nicholls

Organizations: Department of Statistics, University of Oxford, Oxford, United Kingdom. · Mathematical Institute, University of Oxford, Oxford, United Kingdom.

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

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 simulator-based 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.

Figures & tables

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)z_{n}=n^{\alpha}(\theta-c_{n}) 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.
Figure 2: The same location error becomes more serious as uncertainty shrinks.
Figure 3: Results on the six regular benchmarks. Seed means with ±1\pm 1 seed-SD bands. The shaded band is the training range; dotted lines mark the targets ( 11 , 0.950.95 , and 0.50.5 ). SD ratio and center error are on log axes.
Figure 4: Posterior contraction on real 1000 Genomes windows. (a) The typical posterior of each estimator as nn grows from 1010 to 100100 , with the moment fit to all windows as the cross. (b) Against nn : consistency , the fraction of datasets whose 95%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 .
Appendix figures & tables23 assets

Supplementary material from the paper’s appendix.

Appendix

Priorθ∼N(0, 0.1 I10)\theta\sim\mathcal{N}(0,\,0.1\,I_{10})
Simulatorxi∣θ∼N(θ, 0.1 I10)x_{i}\mid\theta\sim\mathcal{N}(\theta,\,0.1\,I_{10})
Dimensionsθ∈R10\theta\in\mathbb{R}^{10} , xi∈R10x_{i}\in\mathbb{R}^{10}
Centersample mean xˉ\bar{x}
Reference posteriorconjugate Gaussian posterior
ReferencesLueckmann et al. (2021)
Appendix
Table 5
Priorθ∼U[−1,1]10\theta\sim U[-1,1]^{10}
Simulatorxi∣θ∼N(θ, 0.1 I10)x_{i}\mid\theta\sim\mathcal{N}(\theta,\,0.1\,I_{10})
Dimensionsθ∈R10\theta\in\mathbb{R}^{10} , xi∈R10x_{i}\in\mathbb{R}^{10}
Centersample mean xˉ\bar{x}
Reference posteriorproduct of truncated normals
ReferencesLueckmann et al. (2021)
Appendix
Table 6
Priorθ∼U[−10,10]2\theta\sim U[-10,10]^{2}
Simulatorxi∣θ∼12N(θ,I)+12N(θ,0.01 I)x_{i}\mid\theta\sim\tfrac{1}{2}\mathcal{N}(\theta,I)+\tfrac{1}{2}\mathcal{N}(\theta,0.01\,I)
Dimensionsθ∈R2\theta\in\mathbb{R}^{2} , xi∈R2x_{i}\in\mathbb{R}^{2}
Centerten EM steps from the coordinate-wise median
Reference posteriorgrid integration, 201×201201\times 201 points around the MAP
ReferencesSisson et al. (2007) ; Beaumont et al. (2009)
Appendix
Table 7
Priorθ=(β,f)\theta=(\beta,f) : β∼N(0,2)\beta\sim\mathcal{N}(0,2) , f∼N(0,(F⊤F)−1)f\sim\mathcal{N}\bigl(0,(F^{\top}F)^{-1}\bigr) , with Fi,i−2=1F_{i,i-2}=1 , Fi,i−1=−2F_{i,i-1}=-2 , Fi,i=1+(i−1)/9F_{i,i}=1+\sqrt{(i-1)/9} , Fi,j=0F_{i,j}=0 otherwise, 1≤i,j≤91\leq i,j\leq 9
Simulatorxi=(xi,1,…,xi,10)x_{i}=(x_{i,1},\dots,x_{i,10}) , xi,1=∑t=1Tztx_{i,1}=\sum_{t=1}^{T}z_{t} , xi,2:10=1xi,1Vzx_{i,2:10}=\frac{1}{x_{i,1}}Vz , zt∼Bern(η(vt⊤f+β))z_{t}\sim\mathrm{Bern}\bigl(\eta(v_{t}^{\top}f+\beta)\bigr) , η(⋅)=exp⁡(⋅)/(1+exp⁡(⋅))\eta(\cdot)=\exp(\cdot)/(1+\exp(\cdot)) , frozen input between time bins t−8t-8 and tt : vt∼N(0,I)v_{t}\sim\mathcal{N}(0,I) , V=[v1,…,vT]V=[v_{1},\dots,v_{T}]
Dimensionsθ∈R10\theta\in\mathbb{R}^{10} , xi∈R10x_{i}\in\mathbb{R}^{10}
Fixed parametersduration T=100T=100 time bins
Centersimulator method of moments: damped Gauss–Newton steps matching X⊤σ(Xθ)X^{\top}\sigma(X\theta) to the replicate mean
Reference posterioradaptive importance sampling from the Newton MAP
ReferencesLueckmann et al. (2017) ; Gonçalves et al. (2020)
Appendix
Table 8
Priorθ∼U[−3,3]5\theta\sim U[-3,3]^{5}
Simulatorxi=(xi,1,…,xi,4)x_{i}=(x_{i,1},\dots,x_{i,4}) , xi,j∼N(mθ,Sθ)x_{i,j}\sim\mathcal{N}(m_{\theta},S_{\theta}) , where mθ=(θ1,θ2)⊤m_{\theta}=(\theta_{1},\theta_{2})^{\top} , Sθ=(s12ρs1s2ρs1s2s22)S_{\theta}=\bigl(\begin{smallmatrix}s_{1}^{2}&\rho s_{1}s_{2}\\ \rho s_{1}s_{2}&s_{2}^{2}\end{smallmatrix}\bigr) , s1=θ32s_{1}=\theta_{3}^{2} , s2=θ42s_{2}=\theta_{4}^{2} , ρ=tanh⁡θ5\rho=\tanh\theta_{5}
Dimensionsθ∈R5\theta\in\mathbb{R}^{5} , xi∈R8x_{i}\in\mathbb{R}^{8}
Centermoment estimator on the branch θ3,θ4>0\theta_{3},\theta_{4}>0
Reference posteriorimportance sampling on a fundamental domain
ReferencesPapamakarios et al. (2019) ; Greenberg et al. (2019) ; Hermans et al. (2020) ; Durkan et al. (2020)
Appendix
Table 9
Priorθ∼U[−1,1]2\theta\sim U[-1,1]^{2}
Simulatorxi∣θ=[rcos⁡α+0.25rsin⁡α]+[−∣θ1+θ2∣/2(−θ1+θ2)/2]x_{i}\mid\theta=\begin{bmatrix}r\cos\alpha+0.25\\ r\sin\alpha\end{bmatrix}+\begin{bmatrix}-|\theta_{1}+\theta_{2}|/\sqrt{2}\\ (-\theta_{1}+\theta_{2})/\sqrt{2}\end{bmatrix} , where α∼U(−π/2,π/2)\alpha\sim U(-\pi/2,\pi/2) , r∼N(0.1, 0.012)r\sim\mathcal{N}(0.1,\,0.01^{2})
Dimensionsθ∈R2\theta\in\mathbb{R}^{2} , xi∈R2x_{i}\in\mathbb{R}^{2}
Centermoment estimator on the branch θ1+θ2≥0\theta_{1}+\theta_{2}\geq 0
Reference posteriorgrid integration on a fundamental domain, expanded by the symmetry group
ReferencesGreenberg et al. (2019)
Appendix
Table 10
Priorβ∼LogNormal(log⁡0.4, 0.5)\beta\sim\mathrm{LogNormal}(\log 0.4,\,0.5) , γ∼LogNormal(log⁡0.125, 0.2)\gamma\sim\mathrm{LogNormal}(\log 0.125,\,0.2) , truncated to R0=β/γ∈[1.5,8]R_{0}=\beta/\gamma\in[1.5,8] (about 12%12\% of prior mass removed; below 1.51.5 no epidemic occurs within the horizon, above 88 the peak precedes the first observation and β\beta is weakly identified). Inference is in (log⁡β,log⁡γ)(\log\beta,\log\gamma) .
Simulatorxi=(xi,1,…,xi,10)x_{i}=(x_{i,1},\dots,x_{i,10}) , xi,t∼Binomial(1000, I(t)N)x_{i,t}\sim\mathrm{Binomial}\bigl(1000,\,\frac{I(t)}{N}\bigr) stored as fractions, where II is simulated from dSdt=−βSIN\frac{dS}{dt}=-\beta\frac{SI}{N} dIdt=βSIN−γI\frac{dI}{dt}=\beta\frac{SI}{N}-\gamma I dRdt=γI\frac{dR}{dt}=\gamma I
Dimensionsθ∈R2\theta\in\mathbb{R}^{2} , xi∈R10x_{i}\in\mathbb{R}^{10}
Fixed parameterspopulation N=106N=10^{6} , duration T=160T=160 , initial conditions (S(0),I(0),R(0))=(N−1,1,0)(S(0),I(0),R(0))=(N-1,1,0)
Centersimulator method of moments: coarse grid, then three damped Gauss–Newton steps matching the simulator’s mean output to the replicate mean
Reference posteriorgrid integration in whitened Laplace coordinates; Gaussian to high accuracy, fitted contraction exponent 0.510.51
ReferencesKermack and McKendrick (1927)
Appendix
Table 11
Priorα∼LogNormal(−0.125, 0.5)\alpha\sim\mathrm{LogNormal}(-0.125,\,0.5) , β∼LogNormal(−3, 0.5)\beta\sim\mathrm{LogNormal}(-3,\,0.5) , γ∼LogNormal(−0.125, 0.5)\gamma\sim\mathrm{LogNormal}(-0.125,\,0.5) , δ∼LogNormal(−3, 0.5)\delta\sim\mathrm{LogNormal}(-3,\,0.5) , truncated to parameters for which neither species falls below e−5e^{-5} individuals at any observation time (about 8%8\% of prior mass removed; an effectively extinct species carries no information about its rates). Inference is in log-rate coordinates.
Simulatorxi,1,t∼LogNormal(log⁡X(t), 0.1)x_{i,1,t}\sim\mathrm{LogNormal}(\log X(t),\,0.1) , xi,2,t∼LogNormal(log⁡Y(t), 0.1)x_{i,2,t}\sim\mathrm{LogNormal}(\log Y(t),\,0.1) stored as log-counts, where X,YX,Y are simulated from dXdt=αX−βXY\frac{dX}{dt}=\alpha X-\beta XY dYdt=−γY+δXY\frac{dY}{dt}=-\gamma Y+\delta XY
Dimensionsθ∈R4\theta\in\mathbb{R}^{4} , xi∈R20x_{i}\in\mathbb{R}^{20}
Fixed parametersduration T=20T=20 time units, initial conditions (X(0),Y(0))=(30,1)(X(0),Y(0))=(30,1)
Centersimulator method of moments: 13413^{4} grid, then four damped Gauss–Newton steps with multi-start
Reference posterioradaptive importance sampling with a Student- tt proposal at the Gauss–Newton MAP; Gaussian to high accuracy, fitted contraction exponent 0.500.50
ReferencesLotka (1920)
Appendix
Table 12
 C2STSD ratio95%95\% joint coverage
estimatorn=10n=10n=20n=20n=50n=50n=100n=100n=10n=10n=20n=20n=50n=50n=100n=100n=10n=10n=20n=20n=50n=50n=100n=100
Gaussian linear
Free-G0.570.610.750.891.021.091.241.410.960.970.970.86
Scale-G0.550.570.650.711.011.001.001.000.960.950.910.87
CE-G0.510.510.550.571.001.001.031.060.960.960.960.97
Free-NSF0.550.580.720.871.011.061.261.530.960.980.991.00
Appendix
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}n\in\{5,\dots,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).
 T1: parameter-dependent supportT2: boundary parameter
methodKLSD ratiocenter err.C2STKLSD ratiocenter err.C2ST
Free-G12.271.094.720.920.121.000.590.61
Scale-G1.312.692.020.820.240.830.390.61
CE-G1.112.710.540.780.100.810.180.58
Free-NSF0.652.020.870.700.171.020.520.61
CE-NSF0.492.821.810.680.180.840.240.59
Appendix
Table 10: Nonregular prototypes at n=100n=100 ( 5×5\times the largest training size; ten seeds). Each block violates a different assumption of the structured representation.
 CE-GCE-NSF
imposed α\alphaSD ratiocenter err.cov. 90%C2STKLSD ratiocenter err.cov. 90%C2STKL
1/21/2 (default)2.710.540.970.781.102.821.810.970.680.49
0.80.8 (estimated)1.400.250.910.720.671.420.590.920.590.14
11 (oracle ablation)0.920.150.860.690.480.920.230.850.560.03
Appendix
Table 11: Result on T1 at n=100n=100 .
 T2, θ0=0\theta_{0}=0 stratumT3T4
estimatorSD ratiocenter err.cov. 90%KLridge dist.mode dist.
Exact posterior1.00 →\to 1.000.00 →\to 0.000.96 →\to 0.940.00 →\to 0.00––
Free-G0.85 →\to 0.740.19 →\to 0.290.70 →\to 0.930.11 →\to 0.113.4 →\to 7.99.5 →\to 21.2
Scale-G0.83 →\to 0.860.19 →\to 0.610.53 →\to 0.030.13 →\to 0.433.4 →\to 7.99.5 →\to 21.2
CE-G0.83 →\to 0.790.19 →\to 0.220.60 →\to 0.580.11 →\to 0.163.4 →\to 6.39.2 →\to 12.7
Free-NSF1.01 →\to 1.140.03 →\to 0.580.77 →\to 0.260.01 →\to 0.233.4 →\to 7.89.5 →\to 21.0
Appendix
Table 12: Fixed- θ0\theta_{0} strata, n=20→100n=20\to 100 (ten seeds). T2 at the boundary ( θ0=0\theta_{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 nn for any single-center summary, the exact posterior mean included. Accuracy on T3 and T4 is read from Table 10 .
 C2STMMD95%95\% joint coverage
estimatorn=10n=10n=20n=20n=50n=50n=100n=100n=10n=10n=20n=20n=50n=50n=100n=100n=10n=10n=20n=20n=50n=50n=100n=100
SLCP
Free-G0.980.980.990.990.0640.0730.0840.0930.960.990.940.87
Scale-G0.980.981.001.000.0580.1420.3630.5320.970.920.200.00
CE-G0.980.980.980.990.0400.0720.1440.2080.950.940.870.69
Free-NSF0.750.790.830.910.0040.0040.0070.0130.920.980.960.82
Appendix
Table 13: The two symmetric multimodal sbibm benchmarks.
 n=20n=20 (training boundary)n=100n=100 ( 5×5\times )
estimatorSD ratiocenter err.cov.C2STSD ratiocenter err.cov.C2ST
Gaussian linear
Free-G1.090.110.970.611.410.590.860.89
Free-G (best effort)1.010.070.960.540.980.810.390.93
Free-NSF1.060.110.980.581.530.531.000.87
Free-NSF (best effort)1.030.080.960.541.420.441.000.83
Appendix
Table 14: Best-effort free baselines beside the standard estimators at the training boundary ( n=20n=20 ) and at 5×5\times ( n=100n=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).
  C2STSD ratio95%95\% joint coverage
NNestimator0.5N0.5N1N1N2.5N2.5N5N5N0.5N0.5N1N1N2.5N2.5N5N5N0.5N0.5N1N1N2.5N2.5N5N5N
20Free-G0.880.920.970.992.273.074.377.430.940.960.610.36
20Scale-G0.880.900.960.992.302.603.094.460.950.920.460.15
20CE-G0.570.550.640.761.021.001.141.830.940.920.950.96
20Free-NSF0.740.790.920.961.461.802.443.780.950.970.630.38
20CE-NSF0.550.530.610.751.471.261.772.400.940.930.930.70
Appendix
Table 15: SIR at two training maxima against the relative size n/Nn/N .
methodC2ST at n=20n=20KL at n=20n=20C2ST at n=100n=100KL at n=100n=100cov. at n=100n=100
Free-G (standardized, tuned)0.8982.3170.9422.9800.966
CE-G ( log⁡n\log n residual)0.5520.1210.5370.0960.870
Fixed EM center, free scale0.5740.1740.6230.3280.974
Fixed EM center, imposed scale0.5580.1370.5130.0330.909
Appendix
Table 16: Fixed-center controls on the Gaussian mixture at the training boundary ( n=20n=20 ) and 5×5\times ( n=100n=100 ); means across three seeds.
 SIRLotka–Volterra
estimatorn=10n=10n=20n=20n=50n=50n=100n=100n=10n=10n=20n=20n=50n=50n=100n=100
Fixed center, free scale0.6530.6420.7890.8970.7680.8060.8860.940
Fixed center, imposed scale0.6340.5720.6470.7500.7290.7040.6560.704
CE-G0.6400.5790.6630.7720.7380.7130.6710.737
CE-NSF0.6020.5590.6140.7400.7070.6650.6550.757
Appendix
Table 17: Fixed-center controls on SIR and Lotka–Volterra: C2ST at every evaluation size (training sizes 55 to 2020 ); means across three seeds, lowest per column in bold.
benchmarkimposed ratetrue asymptotic rateestimated α\alpha
Gaussian linear1/21/21/21/20.500.50
Gaussian linear (uniform)1/21/21/21/20.480.48
Gaussian mixture1/21/21/21/21.381.38
Bernoulli GLM1/21/21/21/20.460.46
SIR1/21/21/21/20.560.56
Lotka–Volterra1/21/21/21/20.480.48
Appendix
Table 18: Contraction exponent estimated from the training range.
residual contextKLSD ratio90% cov.center errorSBC KSexponent
log⁡n\log n0.00180.00180.9910.9910.8860.8860.0610.0610.0390.0390.5000.500
1/n1/n0.00120.00120.9980.9980.8890.8890.0590.0590.0380.0380.5000.500
none0.00230.00230.9920.9920.8880.8880.0650.0650.0410.0410.5000.500
Appendix
Table 19: Matched residual-context comparison at n=10,000n=10{,}000 (means over five seeds; coverage is the nominal 90%90\% region).
 C2STSD ratio
estimatorn=10n=10n=20n=20n=50n=50n=100n=100n=10n=10n=20n=20n=50n=50n=100n=100
Gaussian linear
Free-NSF0.550.580.720.871.011.061.261.53
PF-NPSE (Langevin)0.560.570.570.621.111.121.030.98
PF-NPSE (GAUSS)0.510.520.550.590.990.991.000.98
CE-NSF0.510.510.540.561.021.010.990.98
Appendix
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.
Priorθ∼U[0.2, 6]\theta\sim U[0.2,\,6] per kb; growth parameter ∼U[0, 2]\sim U[0,\,2] (raw rate U[0,8]U[0,8] )
Simulatorcoalescent tree of k=20k=20 under exponential growth (time-changed merging); mutations Poisson on branches, class counts ξb\xi_{b} , then S=∑bξbS=\sum_{b}\xi_{b} , π=∑bb(k−b)ξb/(k2)\pi=\sum_{b}b(k-b)\xi_{b}/\binom{k}{2} , η=ξ1+ξk−1\eta=\xi_{1}+\xi_{k-1}
Dimensionsθ∈R2\theta\in\mathbb{R}^{2} , xi∈R3x_{i}\in\mathbb{R}^{3}
Centersimulator method of moments in the Watterson tradition ( Watterson, 1975 ) : the growth rate from the two observed ratios (πˉ/Sˉ,ηˉ/Sˉ)(\bar{\pi}/\bar{S},\bar{\eta}/\bar{S}) against moment curves precomputed once by simulation, then θ\theta from the scale Sˉ\bar{S} ; a grid lookup with quadratic refinement at query time, likelihood-free and contracting at n−1/2n^{-1/2}
Reference posteriornone: the likelihood is intractable
ReferencesKingman (1982) ; Watterson (1975)
Appendix
Table 25
 n=10n=10n=20n=20n=50n=50n=100n=100
 rad.cov.KSrad.cov.KSrad.cov.KSrad.cov.KS
Free-G1.040.940.0520.870.980.0630.821.000.0870.820.980.131
Free-NSF0.980.950.0330.920.980.0620.920.970.1101.010.920.224
Scale-G1.000.950.0421.060.960.0771.340.870.1671.640.750.181
CE-G0.970.950.0460.980.960.0531.140.940.0911.240.870.075
CE-NSF1.000.950.0460.950.950.0481.040.910.0391.100.860.068
Appendix
Table 22: In-model protocol, every size (seed means): rms standardized truth distance (radius; target 11 ), joint 95%95\% HPD coverage, and maximum marginal SBC KS; bold: within one seed SD of the best per column.
 n=10n=10n=20n=20n=50n=50n=100n=100
 cons.rad.energycons.rad.energycons.rad.energycons.rad.energy
Free-G0.791.320.1620.811.320.1320.771.320.0940.661.590.085
Free-NSF0.801.280.1580.791.520.1230.721.630.0900.581.920.077
Scale-G0.831.210.1580.771.440.1290.581.870.0950.282.790.085
CE-G0.831.230.1530.811.320.1220.851.220.0870.831.230.079
CE-NSF0.831.180.1600.831.280.1210.861.220.0880.851.280.082
Appendix
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.

Explore similar work

CardsList

  1. Conservative neural posterior estimation via distributionally robust training

    May 27, 2026William Laplante, Yuga Hikida, Charita Dellaporta +2Simulation-Based InferenceNeural Posterior Estimation

    Department of Statistical Science, University College London, UK · The Alan Turing Institute, UK · Department of Computer Science, Aalto University, Finland

  2. Mixed neural posterior estimation for simulators with discrete and continuous parameters

    May 13, 2026Jan Boelts, Cornelius Schröder, Jonas Beck +3Simulation-Based InferenceParameter Estimation

    appliedAI Institute for Europe · Machine Learning in Science, University of Tübingen · Tübingen AI Center +3

  3. High-Dimensional Simulation-Based Inference in Latent Spaces

    Sep 29, 2026Lars Kühmichel, Stefan T. Radev, Bhanu Prasanna Koppolu +3Simulation-Based InferenceNeural Posterior Estimation

    Department of Statistics, TU Dortmund University, Germany · Department of Cognitive Science, Rensselaer Polytechnic Institute