Arbitrary Marginal Neural Ratio Estimationfor Simulation-based Inference
François Rozet
University of Liège
francois.rozet@uliege.be

Gilles Louppe
University of Liègeg.louppe@uliege.be

Abstract
In many areas of science, complex phenomena are modeled by stochastic parametricsimulators, often featuring high-dimensional parameter spaces and intractable
likelihoods. In this context, performing Bayesian inference can be challenging. In
this work, we present a novel method that enables amortized inference over arbitrary
subsets of the parameters, without resorting to numerical integration, which makes
interpretation of the posterior more convenient. Our method is efficient and canbe implemented with arbitrary neural network architectures. We demonstrate theapplicability of the method on parameter inference of binary black hole systems
from gravitational waves observations.

1

Introduction

Formally, a simulator is a stochastic forward model that takes a vector of parameters θ ∈ Θ as input,
samples internally a series z ∈ Z ∼ p(z|θ) of latent variables and finally produces an observationx ∈ X ∼ p(x|θ, z) as output, thereby defining an´ implicit likelihood p(x|θ). This likelihood typicallyis intractable as it corresponds to p(x|θ) = Z p(x, z|θ) dz, the integral of the joint likelihoodp(x, z|θ) over all possible trajectories through the latent space Z. In Bayesian inference, we areinterested in the posteriorp(θ|x) =p(x|θ)p(θ)
p(x|θ)p(θ)=´p(x)
p(x|θ0 )p(θ0 ) dθ0Θ(1)

for some observation(s) x and a prior p(θ), which not only involves the intractable likelihood p(x|θ)
but also an intractable integral over the parameter space Θ. The omnipresence of this problem gave
rise to a rapidly expanding field of research [1] commonly referred to as simulation-based inference.
Recent approaches [2–5] are to learn a surrogate model p̂(θ|x) of the posterior and, then, proceed asif the latter was tractable.
However, domain scientists are not always interested in the full set of simulator parameters at once.
In particular, when interpreting posterior predictions, they generally study several small parameter
subsets, like singletons or pairs,´ while ignoring the others. This task corresponds to estimating themarginal posterior p(θa |x) = Θb p(θ|x) dθb over parameter subspaces Θa ≤ Θ of interest, while
the complement subspaces Θb : Θa × Θb = Θ are unobserved. To this end, most applications [4, 5]
resort to numerical integration of a surrogate p̂(θ|x) of the full posterior, which is computationallyexpensive if Θb is large.
A solution to get rid of numerical integration is to learn directly a surrogate p̂(θa |x) by considering
θb as part of the latent variables. If we are interested in a single or a few predetermined subspaces,
this approach is reasonable and leads to accurate estimation of marginal posteriors [6, 7]. However, if
we need to choose arbitrarily the subspace Θa at inference time, this solution is not viable anymore
as there exists an exponential number (2dim(Θ) − 1) of marginal posteriors.
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

Contribution We build upon neural ratio estimation (NRE) [3, 8] to enable integration-lessmarginal posterior estimation over arbitrary parameter subspaces. The key idea is to introduce,as input of the ratio estimator, a binary mask a ∈ {0, 1}dim(Θ) indicating the current subspace Θa ofinterest. Intuitively, this allows the network to distinguish the subspaces and, thereby, to learn a different ratio for each of them. Our method, dubbed arbitrary marginal neural ratio estimation (AMNRE),
can be implemented with arbitrary neural network architectures, including multi-layer perceptrons(MLPs) [9] and residual networks [10]. AMNRE is an amortized method, meaning that inference issimulation-free and can be repeated several times with distinct observations, without retraining. The
counterpart is that AMNRE could require a lot of training simulations to produce accurate predictions.
The implementation is available at https://github.com/francois-rozet/amnre.
Related work Imputation methods [11–13] were the first to introduce a binary mask to condition
networks with respect to which features are missing. This trick allowed to train a single generative
network for all combinations of missing features. Our method differs in that it does not generate likely
replacements for the missing features but evaluates the likeliness, conditionally to an observation, of
those that are provided.

2

Arbitrary marginal neural ratio estimation

NRE The principle of NRE [3] is to train a classifier network dφ : Θ × X 7→ [0, 1] to discriminate
between pairs (θ, x) equally sampled from the joint distribution p(θ, x) and the product of themarginals p(θ)p(x). Formally, the optimization problem isφ∗ = arg min
E 0 L(dφ (θ, x)) + L(1 − dφ (θ0, x)) ,(2)
φp(θ,x)p(θ )where L(p) = − log p is the negative log-likelihood. For this task, the decision function modelingthe Bayes optimal classifier [3] isp(θ, x)d(θ, x) =,(3)p(θ, x) + p(θ)p(x)thereby defining the likelihood-to-evidence (LTE) ratiop(θ, x)p(x|θ)
p(θ|x)d(θ, x)===.
(4)r(θ, x) =1 − d(θ, x)p(θ)p(x)p(x)p(θ)

Consequently, NRE gives access to an estimator log rφ (θ, x) = logit(dφ (θ, x)) of the LTE log-ratio
and a surrogate p̂(θ|x) = rφ (θ, x)p(θ) for the posterior density.
AMNRE With the additional binary mask a ∈ {0, 1}dim(Θ) , the classifier takes the formdφ (θa , x, a) and the optimization problem becomesφ∗ = arg min
E 0 E L(dφ (θa , x, a)) + L(1 − dφ (θa0 , x, a)) ,(5)
φp(θ,x)p(θ ) p(a)where θa = (θi : ai = 1) and p(a) is a mask distribution. In this context, the Bayes optimal classifier(see Appendix A) isp(θa , x)d(θa , x, a) =,(6)p(θa , x) + p(θa )p(x)meaning that AMNRE gives access to an estimator log rφ (θa , x, a) = logit(dφ (θa , x, a)) of all
marginal LTE log-ratios and a surrogate p̂(θa |x) = rφ (θa , x, a)p(θa ) for all marginal posteriors.
AMNRE does not have any particular architectural requirements, with the notable exception of thevariable input size of θa . To make the method more convenient, θa is replaced by the element-wiseproduct θ · a (θa · 1 and θb · 0), carrying the same information at fixed size. The mask a is still required
as input since a zero in θ · a does not unambiguously indicate a zero in a. To prevent numericalstability issues when dφ (θa , x, a) → 1, the approximate log-ratio log rφ (θa , x, a) is extracted fromthe neural network and the class prediction is recovered by application of the sigmoid function.
The mask distribution is an important part of AMNRE’s training. If some masks a have a small
probability p(a) to be selected, it is likely that the estimator will not model their respective marginal
posteriors as well as other, more frequent masks. In our experiments, we adopt a uniform mask
distribution p(a) = (2dim(Θ) − 1)−1 , leaving the study of this aspect to future work.
2

φlog rφ (θa , x, a)

a
x

σ

θ·a

dφ (θa , x, a)Figure 1: Illustration of AMNRE’s classifier architecture.

3

Experiments and results

3.1

Simple likelihood and complex posteriorθ∗
99.7 %
95.5 %
68.3 %
Ground-truth (MCMC)
AMNRE (1d, 2d)

θ2

2
0
−2
2

θ3

Papamakarios et al. [14] introduce a toy
simulator with a 5-dimensional parameter space Θ ⊆ R5 , for which the likelihood is tractable. Despite its simple likelihood, the simulator has a complex posterior (SLCP) with four symmetric modes.
Hence, SLCP is a non-trivial posterior estimation benchmark that allows to retrieve
the ground-truth posterior through Markov
chain Monte Carlo (MCMC) sampling [15,16] of the likelihood.

0
−2
2

3.2

2

0

2

2

−

0

2

2

−

0

2

2

−

0

2

0

2

−

−

2

θ5

θ4

We apply AMNRE on SLCP and compare0
the learned surrogates with the ground−2
truth posterior. Training details are pro2
vided in Appendix C. In Figure 2, we observe that AMNRE 1d and 2d surrogates0
are in close agreement with the ground−2
truth. The structure of the distribution, represented by the credible regions, is modθ1
θ2
θ3
θ4
θ5
eled correctly, even in low density regions.
We also note that the four symmetric modes Figure 2: Ground-truth posterior against AMNRE 1d(see θ3 and θ4 ) are properly recovered, and 2d marginal surrogates, for an observation x∗ of thewhich is sometimes challenging for tradi- SLCP testing set. Density is averaged over three training
tional sampling methods. Concerning the instances. Contours represent the 68.3 %, 95.5 % andparameter θ5 , we observe that the network 99.7 % highest posterior density regions. Stars represent
is slightly underconfident around the mode, the true parameters θ∗ of the observation.
which could indicate that, among the fiveparameters of SLCP, θ5 is the hardest to infer. Finally, in Figure 4, we see that AMNRE is also able
to recover the full 5d posterior and the predictions are very consistent with the 1d and 2d surrogates.
Gravitational waves

In recent years, the observations of gravitational waves (GWs) from compact binary coalescences
systems have had a massive impact on our understanding of the Universe, partly thanks to inferenceof the systems’ parameters. To obtain posterior samples, the LIGO/Virgo collaboration currently
applies MCMC [15, 16] or nested sampling [17, 18] algorithms to involved physical models of thelikelihood of emitted waves [19, 20]. With these approaches, posterior calculation typically takes
days for binary black hole (BBH) mergers and has to be repeated from scratch for each observation.
As a proof of concept, we employ AMNRE to infer the full 15-dimensional set of precessing quasicircular BBH parameters, given GW observations from the LIGO/Virgo detectors. The simulator
details are provided in Appendix B. After training (see Appendix C for details), we evaluate thelearned surrogate model on data surrounding GW150914, the first recorded GW event [21]. As
reference, we use the posterior samples produced by Bilby [20] with the dynesty [18, 22] nested3sampler (MIT License), which leverages the true likelihood. It takes 3 days for Bilby to complete theposterior inference of GW150914, while our network builds histograms (100 bins per dimension) of
all 1d and 2d marginal posteriors in about 1 second, on a single 1080Ti GPU.
95.5 %
68.3 %
Bilby (15d)
AMNRE (1d, 2d)

m2

75
50
1.5

25
0.08

tc

δ

0.0

0.00

α

a2

a1

300

0

5

0.

1.

5
2.

dL

600

5
− .0
1.
5

−1.5

−0.08
900

δ

0.6
0.3
0.6
0.3

θ1

3.0
1.5

θ2

3.0
1.5

6

5

0

5

0

5

0

1.

3.

1.

3.

1.

3.

6

3

3

0.

a1

dL

0.

tc

0.

m2

0.

50
− 7
0. 5
08
0.
00
0.
08
30
0
60
0
90
0

m1

75
25

50

1.5

25

θJN

3.0

a2

θ1

θ2

θJN

Figure 3: AMNRE 1d and 2d surrogate marginal posteriors against marginalized Bilby posterior
samples over a subset of the parameters, for the GW150914 observation. Density is averaged overthree training instances.
As can be seen in Figure 3, the surrogate marginal posteriors of AMNRE share the same structure as
the marginalized posterior inferred by Bilby. For some parameter subsets, especially those containing
the masses m1 and m2 and the inclination angle θJN , the predictions present significant inaccuracies.
For the masses, the surrogates are underconfident but predict the correct modes. For other parameters,
including the coalescence time tc , luminosity distance dL and sky location (α, δ), the surrogates are
in close agreement with Bilby.

4

Conclusions

This work introduces AMNRE, a novel simulation-based inference method that enables integrationless marginal posterior estimation over arbitrary parameter subspaces. Through our experiment with
the SLCP toy simulator, we demonstrate that the proposed algorithm is indeed able to recover the
ground-truth posterior and marginalize it arbitrarily. This experiment also highlights the capacity of
AMNRE to model multi-modal distributions, even using a very basic MLP architecture.
The second experiment consists in applying AMNRE to the problem of BBH parameter inference
from GW observations. This proof of concept demonstrates that AMNRE is able to analyze GW
events several order of magnitude faster (seconds instead of days) than traditional sampling methods.
However, if most of the surrogate marginal posteriors seem accurate, some present significant
inaccuracies. Possible causes are a lack of estimator expressiveness or insufficient simulation budget;
aspects we do not properly study in this work. Still, we believe these results to be a promisingdemonstration of the applicability of AMNRE for convenient interpretation of the posterior inchallenging scientific settings.
4

Acknowledgments
The authors would like to thank Antoine Wehenkel, Arnaud Delaunoy and Joeri Hermans for theinsightful discussions and comments. Gilles Louppe is recipient of the ULiège - NRB Chair on Big
Data and is thankful for the support of the NRB.

References
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]

[15]
[16]
[17]
[18]
[19]

[20]
[21]
[22]

Kyle Cranmer et al. “The frontier of simulation-based inference”. In: Proceedings of the
National Academy of Sciences 117.48 (2020), pp. 30055–30062.
David Greenberg et al. “Automatic posterior transformation for likelihood-free inference”. In:
International Conference on Machine Learning. PMLR. 2019, pp. 2404–2414.
Joeri Hermans et al. “Likelihood-free mcmc with amortized approximate ratio estimators”. In:
International Conference on Machine Learning. PMLR. 2020, pp. 4239–4248.
Pedro J Gonçalves et al. “Training deep neural density estimators to identify mechanistic
models of neural dynamics”. In: Elife 9 (2020), e56261.
Stephen R Green et al. “Complete parameter inference for GW150914 using deep learning”.
In: Machine Learning: Science and Technology 2.3 (2021), 03LT01.
Arnaud Delaunoy et al. “Lightning-Fast Gravitational Wave Parameter Inference through
Neural Amortization”. In: arXiv preprint arXiv:2010.12931 (2020).
Benjamin Kurt Miller et al. “Truncated Marginal Neural Ratio Estimation”. In: arXiv preprint
arXiv:2107.01214 (2021).
Kyle Cranmer et al. “Approximating likelihood ratios with calibrated discriminative classifiers”.
In: arXiv preprint arXiv:1506.02169 (2015).
Kurt Hornik et al. “Multilayer feedforward networks are universal approximators”. In: Neuralnetworks 2.5 (1989), pp. 359–366.
Kaiming He et al. “Deep residual learning for image recognition”. In: Proceedings of the IEEE
conference on computer vision and pattern recognition. 2016, pp. 770–778.
Jinsung Yoon et al. “Gain: Missing data imputation using generative adversarial nets”. In:
International Conference on Machine Learning. PMLR. 2018, pp. 5689–5698.
Mohamed Ishmael Belghazi et al. “Learning about an exponential amount of conditional
distributions”. In: arXiv preprint arXiv:1902.08401 (2019).
Yang Li et al. “ACFlow: Flow models for arbitrary conditional likelihoods”. In: International
Conference on Machine Learning. PMLR. 2020, pp. 5831–5841.
George Papamakarios et al. “Sequential neural likelihood: Fast likelihood-free inference withautoregressive flows”. In: The 22nd International Conference on Artificial Intelligence andStatistics. PMLR. 2019, pp. 837–848.
W Keith Hastings. “Monte Carlo sampling methods using Markov chains and their applications”. In: (1970).
Ming-Hui Chen et al. Monte Carlo methods in Bayesian computation. Springer Science &Business Media, 2012.
John Skilling. “Nested sampling for general Bayesian computation”. In: Bayesian analysis 1.4(2006), pp. 833–859.
Edward Higson et al. “Dynamic nested sampling: an improved algorithm for parameter
estimation and evidence calculation”. In: Statistics and Computing 29.5 (2019), pp. 891–913.
John Veitch et al. “Parameter estimation for compact binaries with ground-based gravitationalwave observations using the LALInference software library”. In: Physical Review D 91.4
(2015), p. 042003.
Gregory Ashton et al. “BILBY: a user-friendly Bayesian inference library for gravitationalwave astronomy”. In: The Astrophysical Journal Supplement Series 241.2 (2019), p. 27.
Benjamin P Abbott et al. “Observation of gravitational waves from a binary black hole merger”.
In: Physical review letters 116.6 (2016), p. 061102.
Joshua S Speagle. “dynesty: a dynamic nested sampling package for estimating Bayesianposteriors and evidences”. In: Monthly Notices of the Royal Astronomical Society 493.3 (2020),
pp. 3132–3158.
5[23] Sebastian Khan et al. “Frequency-domain gravitational waves from nonprecessing black-holebinaries. II. A phenomenological model for the advanced detector era”. In: Physical Review D
93.4 (2016), p. 044007.
[24] Alejandro Bohé et al. “PhenomPv2–technical notes for the LAL implementation”. In: LIGO
Technical Document, LIGO-T1500602-v4 (2016).
[25] Djork-Arné Clevert et al. “Fast and accurate deep network learning by exponential linear units(elus)”. In: arXiv preprint arXiv:1511.07289 (2015).
[26] Sergey Ioffe et al. “Batch normalization: Accelerating deep network training by reducing
internal covariate shift”. In: International conference on machine learning. PMLR. 2015,
pp. 448–456.
[27] Diederik P Kingma et al. “Adam: A method for stochastic optimization”. In: arXiv preprint
arXiv:1412.6980 (2014).
[28] Ilya Loshchilov et al. “Decoupled weight decay regularization”. In: arXiv preprint
arXiv:1711.05101 (2017).
[29] Ilya Loshchilov et al. “Sgdr: Stochastic gradient descent with warm restarts”. In: arXiv preprint
arXiv:1608.03983 (2016).

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes] In Sections 3 and 4.
(c) Did you discuss any potential negative societal impacts of your work? [No](d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [No] Some assumptions are implicit.
(b) Did you include complete proofs of all theoretical results? [No] An informal proof ofthe correctness of AMNRE is provided in Appendix A.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] The
implementation is available as a public GitHub repository.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [Yes] In Appendix C.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [N/A](d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes] In the caption of Figure 5 (a single
1080TI GPU).
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] In Appendix B.
(b) Did you mention the license of the assets? [Yes](c) Did you include any new assets either in the supplemental material or as a URL? [Yes](d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A]
6

(b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]

A

Correctness of AMNRE

Reformulating (5), we have˚L=p(θ, x)p(θ0 ) E L(dφ (θa , x, a) + L(1 − dφ (θa0 , x, a)) dθ dx dθ0p(a)¨ Θ×X ×Θ=
E p(θ, x)L(dφ (θa , x, a)) + p(θ)p(x)L(1 − dφ (θa , x, a)) dθ dx
Θ×X p(a)¨= Ep(θa , x)L(dφ (θa , x, a)) + p(θa )p(x)L(1 − dφ (θa , x, a)) dθa dx,
p(a)
{z
}
Θa ×X |`(dφ (θa ,x,a))which is minimized only if each term `(dφ (θa , x, a)) is itself minimized. Assuming p(θa )p(x) > 0,
if p(θa , x) = 0, ` is uniquely minimized by the value 0. Otherwise, if p(θa , x) > 0, the minimum isonly reached by a value q such that
0=

d`(q)
dq

dL(q)
dL(1 − q)+ p(θa )p(x)
dq
dq
−1
1
= p(θa , x)
+ p(θa )p(x)q1−qp(θa , x)q== d(θa , x, a).
p(θa , x) + p(θa )p(x)= p(θa , x)⇔

Importantly, if p(θa , x) = 0, d(θa , x, a) = 0 and still minimizes `. Therefore, as long as p(θa )p(x) >0, d(θa , x, a) is the optimal classifier.

B

Simulators

B.1

Simple likelihood and complex posterior

In this toy simulator, θ ∈ R5 parametrizes a 2d multivariate Gaussian from which four points are
independently sampled to construct an observation x. The generative process [14] isθi ∼ U(−3, 3)

s1 = θ32 ,i = 1, . . . , 5

fors2 = θ42 ,µ = (θ1 , θ2 ),ρ = tanh(θ5 ) 2ρs1 s2s1
Σ=ρs1 s2s22

x = (z1 , . . . , z4 ) where zj ∼ N (µ, Σ),Q
for which the likelihood p(x|θ) = j p(zj |θ) is tractable.
B.2

Gravitational waves

As we do not have enough knowledge in the domain, we borrow the BBH simulator implemented by
Green et al. [5] (MIT License). We succinctly describe the generative process in this section. For
more information, please refer to the original paper or the implementation.
7

Prior We perform inference over the full 15-dimensional set of precessing quasi-circular BBH parameters: component masses (m1 , m2 ), reference phase φc , coalescence time tc , luminosity distancedL , spin magnitudes (a1 , a2 ), spin angles (θ1 , θ2 , φ12 , φJL ), inclination angle θJN , polarization
angle ψ, and sky location (α, δ). To analyze GW150914, we take a prior uniform over
10 M ≤ mi ≤ 80 M
−0.1 s ≤ tc ≤ 0.1 s
100 Mpc ≤ dL ≤ 1000 Mpc
0 ≤ ai ≤ 0.88
and standard over the remaining quantities. We take tc = 0 to be the trigger time of GW150914 andconstraint m1 ≥ m2 .
Waveform generation The simulator generates waveforms using the IMRphenomPv2 frequencydomain processing model [23, 24] and assumes stationary Gaussian noise with respect to the noise
power spectral density (PSD) estimated from 1024 s of detector data prior to GW150914. The
frequency ranges from 20 to 1024 Hz and each waveform has a duration of 8 s. The waveforms arewhitened with respect to the estimated PSD.
Waveform processing The observations are quite large (16 384 features) and, thereby, impracticalto store on disk and feed to a neural network. To alleviate this problem, the waveforms are compressed
to a reduced-order basis corresponding to the first 128 components of a singular value decomposition(SVD). Using more SVD components did not help producing better predictions, likely due to thehigher ratio of noise in less significant components.
Since an observation corresponds to two waveforms from two geographically distant detectors (H1and L1) and frequency-domain signals are represented by complex-valued vectors, each processed
observation is a vector of 128 × 2 × 2 = 512 real-valued numbers.
Noise For the training set, the noise of the detectors is not added to the stored waveforms. Instead,noise realizations are sampled with respect to the PSD in real time during training, which effectively
increases the size of the training set.

C

Experimental details

Datasets For each simulator, we use three fixed datasets of pairs (θ, x) ∼ p(θ, x) to train, validate
and test AMNRE, respectively. The sizes of the datasets are provided in Table 1.
Table 1: Dataset sizes for each simulator.
Simulator

Training set

Validation set

Testing set

SLCP
GW

1 048 576
4 194 304

131 072
131 072

131 072
131 072

Architectures For SLCP, we use an MLP with 7 hidden layers of 256 neurons and ELU [25]activation functions. For GW, the classifier is a residual residual network [10] consisting of 17
residual blocks of 2 linear layers with 512 neurons and ELU [25] activation functions. In the blocks,
we insert batch normalization layers [26] before the activation functions.
Training All networks are optimized with the AdamW [27, 28] stochastic optimization algorithm.
At each epoch, the batches are built by sampling without replacement from the training set. The
independent parameters θ0 are obtained by shifting circularly (i ← i + 1 and n ← 1) the batchof parameters θ. Each element in the batch has a different mask, sampled from the uniform mask
distribution.
For SLCP, we apply a “reduce on plateau” scheduling to the learning rate, that is, we divide thelearning rate by a factor 2 each time the loss on the validation set has not decreased for 7 consecutive
8

epochs. The training stops when the learning rate reaches 10−6 or lower. For GW, we apply a learning
rate cosine annealing [29] over 512 epochs. Other hyperparameters are provided in Table 2.
Table 2: Training hyperparameters.

D

Hyperparameter

SCLP

GW

Optimizer
Weight decay
Batch size
Batches per epoch
EpochsScheduling
Initial learning rate
Final learning rate

AdamW
10−4
1024
256
reduce on plateau
10−3
10−6

AdamW
10−4
1024
1024
512
cosine annealing
2 × 10−4
10−6

Additional figuresθ∗
99.7 %
95.5 %
68.3 %
5d (MCMC)
1d and 2d

θ2

2
0
−2

θ3

2
0
−2

θ4

2
0
−2

θ5

2
0

θ1

θ2

θ3

2

0

2

2

−

0

2

2

−

0

2

2

−

0

2

2

−

0

−

2

−2

θ4

θ5

Figure 4: AMNRE full 5d surrogate posterior against 1d and 2d marginal surrogates, for an observation of the SLCP testing set. The predictions of AMNRE for the marginal posteriors are consistent
with the predictions for the full posterior, marginalized onto the 1d and 2d subspaces.

9

1.2
1.0

validation
training

1.0

Loss

Loss

1.2

validation
training

0.8

0.8
0.6

0.6

0.4

0.4
0

25

50

75

100

125

0

100

Epoch

200

300

400

500

Epoch

Figure 5: Mean training and validation losses of AMNRE surrogate models for SLCP (left) and GW(right) simulators. Each color corresponds to a different training instance. All instances converge
without signs of overfitting. Training takes around 5 minutes for SLCP and 8 hours for GW, on asingle 1080Ti GPU.

MNRE
AMNRE

θ1
θ2
θ3
θ4
θ5
θ1 , θ2
θ1 , θ3
θ2 , θ3
θ1 , θ4
θ2 , θ4
θ3 , θ4
θ1 , θ5
θ2 , θ5
θ3 , θ5
θ4 , θ5
−0.5

0.0

0.5

1.0

KL to ground-truth

Figure 6: KL divergence to the marginalized ground-truth posterior of 1d and 2d surrogate marginal
posterior histograms. The bars represent the mean and standard deviation over 64 observations fromthe SLCP testing set. AMNRE does not diverges more from the ground-truth than MNRE [3, 7],
despite using a single network for all subspaces.

10

1.0

0.6

0.4

0.6

0.4

0.2

0.2

0.0

0.0
0.0

0.2

m1
m2
tc
dL
a1
a2
θ1
θ2
θJN
α
δ

0.8

CDF(p)

0.8

CDF(p)

1.0

θ1
θ2
θ3
θ4
θ5

0.4

0.6

0.8

1.0

0.0

p

0.2

0.4

0.6

0.8

1.0

p

Figure 7: Cumulative distribution function (CDF) of the percentiles p of 8192 parameters θ∗ in the
one-dimensional surrogate marginal posteriors p̂(θi |x∗ ) for pairs (θ∗ , x∗ ) of the SLCP (left)and GW(right) testing sets. If the surrogate posterior is consistent with the prior, i.e. if Ep(x) p̂(θ|x) ≈ p(θ),
the percentiles should be distributed uniformly between 0 and 1. Since the CDFs lie close to the
diagonal, we conclude that the surrogates are consistent with the prior.
Importantly, one cannot conclude that the network models properly the posterior from this result, as
any distribution consistent with the prior, including the prior itself, would present diagonal CDFs.
Green et al. [5] inadvertently draw this erroneous conclusion.

11