HNPE: Leveraging Global Parameters for Neural Posterior Estimation
Authors: Pedro L. C. Rodrigues
Thomas Moreau
Gilles Louppe
Alexandre Gramfort
Venue: NeurIPS 2021 Poster (2021)
Keywords: Simulation-based Inference
Bayesian Inference
Hierarchical Models
Normalizing Flows
Computational Neuroscience
Approximate Inference

Abstract: Inferring the parameters of a stochastic model based on experimental observations is central to the scientific method. A particularly challenging setting is when the model is strongly indeterminate, i.e. when distinct sets of parameters yield identical observations. This arises in many practical situations, such as when inferring the distance and power of a radio source (is the source close and weak or far and strong?) or when estimating the amplifier gain and underlying brain activity of an electrophysiological experiment. In this work, we present hierarchical neural posterior estimation (HNPE), a novel method for cracking such indeterminacy by exploiting additional information conveyed by an auxiliary set of observations sharing global parameters. Our method extends recent developments in simulation-based inference (SBI) based on normalizing flows to Bayesian hierarchical models. We validate quantitatively our proposal on a motivating example amenable to analytical solutions and then apply it to invert a well known non-linear model from computational neuroscience, using both simulated and real EEG data.

HNPE: Leveraging Global Parameters for Neural Posterior Estimation
Pedro L. C. Rodrigues Inria, CEA, Université Paris-Saclay, France
Thomas Moreau Inria, CEA, Université Paris-Saclay, France
pedro centoorodriguee wwwriiueeeennia.....
thomas Mhomas.moreau@inge
Gilles Louppe
Alexandre Gramfort
University of Liège, Belgium g s.loupp@uulieee.ce
Inria, CEA, Université Paris-Saclay, France
alexandre, gramfort@inria.ft
Abstract
Inferring the parameters of stochastic model based on experimental observations is central to the scientific method. A particularly challenging setting is when the model is strongly indeterminate, i.c. when distinet of parameters yield identical observations. This arises in many practical situations, such as when inferring the distance and power of radio source (is the source close and weak far and strong?) or when estimating the amplifier gain and underlying brain activity an clectrophysiologcal experiment. In this work, we present hierarchical neural posterior estimation (HNPE), novel method for cracking such indeterminacy by exploiting additional information conveyed by an auxiliary set of observations sharing global parameters. Our method extends recent developments in simulation- based inference (SBI) based normalizing flows Bayesian hierarchical models. We validate quantitatively our proposal on a motivating example amenable to analytical solutions and then apply it to invert well known non-linear model from computational neuroscience.
1 Introduction
Simulation- based inference (SBI) has the potential to revolutionize experimental science as it opens the door to the inversion of arbitrary complex non-linear computer models, such as those found in physics, biology, or neuroscience (Cranmer et al., 2020). The only requirement is to have access mathematical model implemented simulators in neuroscience (e.g. Leon et al. 2013), it could estimate properties the brain closer to the cellular level, thus closing the gap between the neuroimaging and computational neuroscience communities. Grounded in Bayesian statistics, recent SBI techniques profit from recent advances in deep generative modeling to approximate the posterior distributions over the full simulator parameters. Their intrinsic quantification of uncertainties reveals whether certain parameters are worth (or not) scientific interpretation given some experimental observation.
simulator. When applied biophysical models and
SBI concerned with the estimation a conditional distribution over parameters interest O. Given some observation æp. the goal is to compute the posterior p(0 (0)))) It generally happens that some of these parameters are strongly coupled, leading very structured posteriors with low dimensional sets of equally likely parameters values. For example, this happens when the data generative process depends only the products of some parameters: multiplying one of such parameters by constant and another by its inverse will not affect the output. Performing Bayesian inference on such models naturally leads to "ridge" or "banana shape" in the posterior landscape, as seen e.g. in Figure 4 of Gonçalves al. (2020). More formally the present challenge posed as soon as the model likelihood
35th Conference Neural Information Processing Systems (NeurIPS 2021).

---

function is non-injective w.r.t 0. and is not strictly due to the presence of noise on the output observations. In statisties and econometrics literature, such models are called partially identified models (Gustafson, 2014).
To alleviate the ill- ill-posedness of the estimation problem, one may consider a hierarchical Bayesian model (Gelman and Hill, 2007) where certain parameters are shared among different observations. In other words, the model's ' parameters 0, for an observation o; are partitioned into 0, = {ax, (9), where a;; is set of sample specific (or local) parameters, and B corresponds to shared (or global) parameters. For this broad class of hicrarchical models, the posterior distribution for a set X {a1,..... 2N of N observations can be written as (Tran et al., 2017):
N
p(a B12) o ((3)
II 3)p(a,l1)) F=1
(1)
Hierarchical models share statistical strength across observations, hence resulting in sharper posteriors and more reliable estimates of the (global and local) parameters and their uncertainty. Examples of applications of hierarchical models are topic models (Blei et al., 2003), matrix factorization algorithms hSalakhutdinov et al., 2013), including Bayesian non-parametrics strategies (Teh and Jordan, 2010), and population genetics (Bazin et al., 2010).
In this work, we further assume that the likelihood function p(z; (3) is implicit and intractable, which implies that traditional MCMC methods not be used to estimate the posterior distribution. This setup leads to so-called likelihood-free inference (I.FI) problems and many algorithms (Papa- makarios and Murray, 2016; Greenberg et al., 2019; Hermans et al., 2020; Durkan et al., 2020b) have recently been developed to earry out inference under this scenario. These methods all operate by learning parts of the Bayes' rule, such as the likelihood function, the seeeihood-to-cviddee or the posterior itself. Approaches for LFI in hierarchical models exist, but are limited. Bazin et al. (2010) extend approximate Bayesian computation (ABC) into a two-step procedure in which local and global variables are estimated. Tran et al. (2017) adapt variational inference to hicrarchical implicit models using a GAN-like training approach, while Brehmer et al. (2019) and Hermans al. (2020) extend amortized likelihood ratios to deal with global parameters, but cannot do inference on local parameters. Motivated by the posterior estimates of individual samples, we consider sequential neural posterior estimation approach derived from SNPE-C (Greenberg et al., 2019).
ratio,
The paper is organized as follows. First, we formalize our estimation problem by introducing the notion of global and local parameters, and instantiate it on a motivating example amenable to analytic posterior estimates allowing for quantitative evaluation. Then, we propose neural posterior estimation technique based on a pair of normalizing flows and a deepset architecture (Zaheer et al., 2017) for conditioning on the set X' of observations sharing the global parameters; we call our method 'hierarchical neural posterior estimation', or simply HNPE. Results on an application with time series produced by a non-linear model from computational neuroscience (Ableidinger et al., 2017) demonstrate the gain in statistical power of our approach thanks to the use of auxiliary observations. We also use this model to analyse real brain signals, giving a full demonstration of the power of LFI to relate parameters from theoretical models to real experimental recordings.
2 Hierarchical models with global parameters
Motivating example. Consider a stochastic model with two parameters, and B, that generates as output x a, + E, where
(O, o²). We assume that both parameters have uniform prior distribution a, 8 ~ [1, 1] and that a is known and small. Our goal is to obtain the posterior distribution of (a, :3) for a given observation za = ag/o + This simple example describes common situations where indeterminacy emerges. For instance, Zn could be the radiation power measured by sensor, the intensity of the emitting source, and / the inverse squared distance of the sensor to the source. In this case, given measurement may have been due to either close weak sources (a J and /3 1) or far strong ones (a t and 3 L). Using Bayes' rule and considering a small we can write (see Appendix A for more details)
e (no (000-0))/220 /20² pa)
p(a, B|zo) ~0
(2)
22
log(1/2o)
where 1 10,y) (a:) is indicator function that equals for E [a,, bj and zero elsewhere. Note that the first term in the product converges to 5(1ro a:9) as + 0 and that the joint posterior distribution has

---

an infinite number of pairs (a, 3) with the same probability, revealing the parameter indeterminacy
of this example. Indeed, for [0, 1] and E [s, 11, all pairs of parameters (5. ()) yield the same observations and the likelihood function p( () constant. Thus, the posterior distribution has level sets with ridge "banana shape" along these solutions. The top row of Figure I on Page 4 portrays the joint and the marginal posterior distributions when (ao, So) = (0.5, 0.5) and a = 0.
Exploiting the additional information in X. Our motivating example illustrates a situation where two parameters related in such way that one may not be known without the other. In practice, however, it possible that one these parameters is shared with other observations. For instance, this is the case when a single source of radiation is measured with multiple sensors located at different unknown distances. The power of the source fixed across multiple measurements and its posterior can be better inferred by aggregating the information from all sensors. Our goal in this section is to formalize such setting so to leverage this additional information and obtain a posterior distribution that 'breaks' parameter indeterminacy. Note that the root cause of the statistical challenge here is not the presence of noise, but rather the intrinsic structure of the observation model.
To tackle the inverse problem of determining the posterior distribution of parameters (oo, 33 given observation zn of stochastic model, we consider the following scenario. We assume that the model's 's structure is such that On is a parameter specific to each observation (local), while , is shared among different observations (global). Yet. both unknown. We consider having access to set X' { NN) ) additional observations generated with the same B as 20-
Taking the model's hierarchical structure into account we use Bayes' rule to write
p( 3 X²)
p(axa X))((
x p( 3. (())
< p(ac0 G. 20 20) p(r; (3)
< p(ann, 3.0
II p((]x))
which shows how the initial posterior distribution p(ac: B| is modified by additional observations from sharing the same 3 as To. In Section 3, we present a strategy for approximating such posterior distribution when the likelihood function of the stochastic model of interest is intractable and, therefore, the posterior distributions p(am) z0) and (B]ro, X)) have be approximated with conditional density estimators.
Motivating example with multiple observations. We now detail the effect of X' on the posterior distribution of our motivating example. The N observations in {ao] U X such that x; a,,0 + for é O,... N with as ~ (4(0, 1] drawn from the same prior. The posterior distribution may be written as (see Appendix A)
I (ma((B) N Nlog(1/2a) BN
p(oa, Blaoo X) Rs pl
(4)
(1/p -1)
where j = max( [2a) U X). This expression shows how the initial full posterior distribution (2) changes with the extra information conveyed by X. It can be also shown that as N + 0 (no additional observations) the posterior distribution converges back to p(0o, 3 (. Figure portrays the joint and marginal posterior distributions with = 10 and = 100.
3 HNPE: neural posterior estimation on Bayesian hierarchical models
When the likelihood function of the stochastic model is intractable, MCMC methods commonly used for posterior estimation are not applicable, since they depend on the evaluation of likelihood ratios, which are not available analytically nor numerically. To bypass such difficulty, we employ tools from eikeihood-free inference (LFT) to directly estimate approximation to the posterior distribution using conditional neural density estimator trained over simulations the model. In what follows, we present novel neural network architecture for approximating the posterior distribution of hierarchical model with global parameters based normalizing flows. We also describe the training procedure for learning the parameters of the network using multi-round procedure known as sequential neural posterior estimation or SNPE-C (Greenberg et al., 2019).

---

plcr. (110)
p(o.m)) mi)
718 mg)
Round
Extra observations V
Nt. X)
p/e Xt)
8 p(88ks.. X2)
10
W(e, (0 ((ii)
16-
in
0.
pla, A X)
pleri .7)
pi8 x
100
Number of simulation round
Figure I: Results on the motivating example from Section 2 with o = 0. (Left) Plots of the analytic (black) and upproximated (red) posterior distributions; the ground truth values on and 3o which generate zi are indicated with dashed lines. Approximations are oblained using the strategy described in Section 3 with R = 5 rounds of 10 simulations from the model. We observe that adding N 10 and N 100 observations X' significantly reduces uncertainty over the estimates oo and Po. (Right) Evolution of the Sinkhorn divergence W, between the analytic posterior distribution p(0|zo. X)) and our approximation d6 trained with increasing number of simulations per round. The larger the simulation budget, closer the learned posterior the analytic one for any number of extra observations. Sequential refinement of the posterior (dash) improves the approximation.
3.1 Approximating the posterior distribution with two normalizing flows
We approximate p( B|ro, X) based on its factorization (3) follows:
p(3 X) gø. (B)xx; fos (X))
(5)
p(ao0 xo) Ja. (cro) z0)
where Ge, and le, are normalizing flows, i.e,, invertible neural networks capable of transforming data points sampled from simple distribution, c.g. Gaussian, to approximate any probability density
function (Papamakarios et al., 2019). The function fø, is a deepset neural network (Zaheer et al.,
2017) structured as fø, (X) 9 EM h ((g)) where h a neural network parametrized
representation for the data points in X and is network parametrized by (i(2) that processes the average value of the embeddings. Note that this aggregation step is crucial for imposing the invariance to permutation of the neural network. It would also be possible to choose other permutation invariant operations, such the maximum value of the set the sum of its elements, but we have observed more stable performance on our experiments when aggregating the observations by their average. It is possible to show that fex is an universal approximator invariant the ordering of its inputs (Zaheer al., 2017). Such property is important for our setting because the ordering of the extra observations in X' should not influence the approximation of the posterior distribution. We refer to approximation either by its factors 9d, and go., by d with Φ (Φ1 Φ2, Φ3)-
by ( that generates
Estimating We estimate ái by minimizing the average Kullback-Leible divergence between the posterior distribution p(ac, Bl X) and 9s (as, Bl X') for different Zo and X::
min. +
Ex(or)))) KL(p(an, 3 No. x²) 1a (oo. B1xo,
where (L(pllga) = if, and only if, p(o 3 Ir, X) ( 2)) We may rewrite the optimization problem terms of each of its parameters get
min.
(a(()) + L (2)1)))))))
(6)
01.00

---

with
(a( = Eplxg a.. a..9) (og((pz((a)) mo) L O3) = - -Bl, a..3) [log(1>, (3) f (2)))
Training from simulated data. In practice, we minimize the objective function in (6) using a Monte-Carke Carlo approximation with data points generated using the factorization p( X, 0x., 3)
p(3) p(z, p() where p(a;, 3) = pla, (3)(3)))((3) is prior distribution describing our initial knowledge about the parameters, and p( Or4, (3) related to the stochastic output of the simulator for a given pair of parameters (as, (). More concretely, the training dataset is generated as follows: First, sample set of parameters from the prior distribution such that (3) ~ p(a;, (3) with = I...., and i = O.....N. Then, for each (i,j)-pain generate observation from the stochastic simulator xf p(x)oe,, (1) so that each observation ag is accompanied by its corresponding N extra observations X "of....,."" xxy) The losses La and La are then approximated by
C # Eド= log(ees (11Bi,23)) and 23 = - Eド= log(es, (8 fes (X*)))
3.2 Refining the approximation with multiple rounds
The optimization strategy above minimizes the KL divergence between the true posterior distribution p and the approximation 9d. on average, for all possible values of zo and X. This is sometimes called amortization, since the posterior distribution is expected to be well approximated for every possible observation. However, when the observed data is scarce and/or difficult to obtain or simulations of the model are costly, it might be useful to focus the capacity gs to better estimate the posterior distribution for a specific choice of 7g and X.
We target the approximation lo to Th and X using adaptation to SNPE-C (Greenberg et al., 2019). This algorithm uses a multiround strategy in which the data points used for minimizing the loss function L and obtaining parameters ((r) round are obtained from simulations with cn. 3 ~ 96t--1 (a.. Bloo. &). At round r = O, parameters Org and / are generated from their prior distributions, which boils down to the procedure described in Section 3.1. Note that an important point is that for the different rounds, the extra observations X' should be simulated with the parameters af drawn from the original prior distribution p(a, since the posterior distribution returned by the multi-round procedure is only largeted for observation Zg. We refer the reader to Greenberg et al. (2019) for further details on the usual SNPE-C procedure, notably a proof of convergence (which extends our case) of the targeted version of 9 the correct posterior density p( p(arp, B| X)) as the number of simulations per round tends to infinity. Algorithm describes the procedure for obtaining a( 3 @0 X') after R rounds simulations.
Algorithm 1: Sequential posterior estimation for hicrarchical models with global parameters
Input :observation To. X, prior (0)), simulator for round r = 1 to R do for sample
2
to n do S(a 3) for (agg, 3i) 1r-1)); Draw x
Draw a set of extra observations Xi = (S(a),))) for af (1133)
s
Train Qgte) to minimize L + Cj:
Set next proposal p! = gatr) R);
7 return posterior 9gon (20, X))
4 Experiments
All experiments described next are implemented with Python (Python Software Fundation, 2017) and the sbi package Trejero-Cantent et al., 2020) combined with PyTorch (Paszke et al., 2019). Pyro (Bingham al., 2018) and nflows (Durkan et al., 2020a) for posterior estimation' In all
Code is available in the supplementary materials.

---

W)()))
Wi.3
10-1
10
Figure 2: Results on the mo- tivating example described in Section 2 (a = 0.05). We see that the marginal pos- teriors get sharper around the ground truth parameter when more observations are available.
10-2
10-2
10-1
HNPE
20-4
h.ACC
LFVI
10-
10
10a
10?
102 10a
10
Number of observations N
Number of observations N
experiments, we use the Adam optimizer (Kingma and Ba, 2014) with default parameters, learning rate of 5.10 and batch size of 100. The code required for reproducing most of the results presented in the paper is available at https:/ //ggtttt...... cwwmccmm.pcc..... httppplloorrggge gueii/NueaHHPP /HNPE
4.1 Results on the motivating example
To evaluate the impact of leveraging multiple observations when estimating the parameters of a hierarchical model. we the model presented in Section 2, where the observation 2o is obtained the product of two parameters og and in with independent uniform prior distributions in [0, 1] (we consider the case where σ = 0). The set of extra observations X' { is obtained by fixing the same global parameter 8g for all æ; and sampling local parameters og from the prior distribution.
Our approximation to the posterior distribution consists of two conditional neural spline flows of linear order (Durkan et al., 2019), (d, and Jæh. both conditioned by dense neural networks with one layer and 20 hidden units. We use neural spline flows because of the highly non-Gaussial aspect of the analytic marginal posterior distributions, which can well captured by this class of normalizing flows. In general. however, the true posterior distribution is not available, of normalizing flows might be justifiable, especially if one's main goal is simply identify set of parameters generating a given observation. We set the function fø, to be simply an averaging operation over the elements of X' the observations in this case are scalar, the only parameters to be learned Algorithm are Φ1 and Φb2
using other classes
We first illustrate in Figure the analytic posterior distribution p(a, (,, x() and approximation (a, B Fo X²) with an increasing number of extra observations (N = 0, 10, 100). For N = 0. i.c. only zo is available, we observe ridge shape in the joint posterior distribution, which is typical of situations with indeterminacies where all solutions (5. .) have the same probability. The addition of few extra observations resolves this indeterminacy and concentrates the analytic postcrior distribution on reduced support (zo, min(1, #e) )] x [M; 11, where μ = max([zo) U X). Moreover, on this support, the solutions are no longer equally probable due the 3 - factor that increases the probability of solutions close to p.. Also note that our estimated posterior is close to the analytic one in all cases.
To have a quantitative evaluation of the quality of approximations 4s. in Figure/ we plot the Sinkhorn divergence (Feydy et al., 2019) W. for = 0.05 between the analytical posterior p(a, 9/o0, X and our approximation for different numbers of simulations per round (cf. Algorithm I). The curves display the median value for nine different choices of (ag, 86) and the transparent area represent the first and the third quartiles. As expected, we note that as the number of simulations per round increases, the approximation gets closer to the analytic solution. The figure also confirms the intuition that, in general, the sequential refinement of multiple rounds Icads to better approximations of the true posterior distribution for fixed observation.
Our next analysis assesses how the posterior approximation concentrates around given point in the (oo. 3) space as the number extra observations N increases. In Figure 2, we display the Wasserstein distances between the marginals of the learned posterior distribution and a Dirac at the ground truth values generating the observation Zo; we consider the results nine different choices of (ao, 3o) but display only the median results. We see that the distance to the Dirac for the global parameter 3 decreases as more observations are added to X, but for the local parameter ag it stabilizes

---

Figure 3: Posterior esti- mates for the parameters of the neural mass model ob- tained on 8 of data sampled at 128 Hz and simulated us- ing C 135, / 220, a= 2000, and g = (rep- resented with black dots in the figure). We observe that the concentration of the pos- terior distribution around Dirac converges to a lower bound when N grows, re- flecting intrinsic uncer- tainty in the parameter es- timation that cannot be duced with more observa- tions.
Normaliee: Wils. Re)
1:25
1.00
400
200
Namber of extra observations N
4300
2000
20
100)
200
2301
1 20000 4000
400
-20
20
on lower bound. This happens because / is observed several times and, therefore, expected to be well estimated, whereas the local parameter is obtained "through the lens" of the estimated / with information from a single observation Hn corrupted by additive noise (og z za/ (). We compare our method (HNPE) with three other approaches: a naive posterior estimation using single normalizing flow with the same capacity as the approximation with tor and lp2- i.e. two layers with 20 hidden a units each, in which we stack the observations from To and the average from those in X' as context variables, the hierarchical ABC (h-ABC) proposed in Bazin al. (2010) and the likelihood-frim variationaa inference (LFVI) presented Tran et al. (2017). The flow-based approaches are trained with R = rounds of 10 simulations and A-ABC has the same simulation budget with acceptance rate of %. We see that the naive approach has very similar performance HNPE, mainly due the low dimensionality of the example being considered (in Section 4.2 we show an example where the naive architecture has similar performance HNPE well but taking much longer train). LFVI captures well the global parameter as increases, but il performs poorly for local parameter. Indeed, we have not found any evidence in the literature showing that LFVI could well estimate local parameters. For instance, all examples in Tran el al. (2017) involve only global variables. The posterior estimated with h-ABC does not concentrate for any of the parameters, which indicates that it would probably need a larger simulation budget to atlain results comparable to the other methods.
4.2 Inverting non-linear model from neuroscience
We consider a class of non-linear models from computational neuroscience known as neural mass modely (Jansen and Rit, 1995) (NMM). These models of cortical columns consist of set of physi- ologically motivated stochastic differential equations able to replicate oscillatory electrical signals observed with electroencephalograpp Such models are used in large-scale simulators (Sanz Leon al., 2013) to generate realistic neural signals oscillating at different frequencies and serve as building blocks for several simulation studies in cognitive and clinical neuroscience (Aerts el al., 2018) In what follows, version of such models presented Ableidinger et al. (2017) and use the C++ implementation in the supporting code of Buck war el al. (2019). In simple terms, the NMM that we consider may be seen a generative model taking input set of four parameters and generating output time series æ. The parameters the neural mass model are:
(EEG) or using intracranial electrodes (Deco el al., 2008).
focus in the stochastic
C, which represents the degree of connectivity between excitatory and inhibitory neurons in the cortical column modelled by the NMM. This connectivity is at the root of the temporal behavior of x and only certain ranges of values generate oscillations.
μ and model the statistical properties of the incoming oscillations from other neighbouring cortical columns. They drive the oscillations of the NMM and their amplitudes have direct effect on the amplitude of x.

---

9 represents a gain factor relating the amplitude of the physiological signal s generated by the system of differential equations for given set (C, (), and the electrophysiology measurements X, expressed in Volts.
The reader is referred to Appendix B for the full description of the stochastic differential equations defining the neural mass model.
Note that the NMM described above suffers from indeterminacy: the same observed signal zo could be generated with larger (smaller) values of and smaller (larger) values of Jr and a. Fortunately, other auxiliary signals
it is common to record several chunks of signals within an experiment, T1, F1,..... N obtained with the same instrument setup (and, therefore, the same gain g) can be exploited. Using the formalism presented in Section 3, we have that a = (C, ()) and [
In what follows, we describe the results obtained when approximating the posterior distribution p(C. μ. o, g X') with Algorithm 1 using R = 2 rounds and n = 50000 simulations per round. Each simulation corresponds to 8 seconds of signal sampled at 128 Hz, so each simulation out- puts vector of 1024 samples. The prior distributions of the parameters are independent uniform distributions defined as:
C' ~ (~, 250) /4 ~ (50, 500) a ~ ((,0 5000) ~ U1( 4(-30, 30, +-30)
where the intervals were chosen based review of the literature neural mass models (Jansen and Rit, 1995; David and Friston, 2003; Deco al., 2008). Note that the gain parameter g is given in decibels (dB), which is standard scale when describing amplifiers in experimental setups. We have, therefore, that at) = 109/10s(t) s(t).
It is standard practice in likelihood-Fie inference to extract summary features from both simulated and observed data in order reduce its dimensionality while describing sufficiently well the statistical behavior the observations. In the present experiment, the summary features consist of the logarithm of the power spectral density (PSD) of each observed time series (Percival and Walden, 1993). The PSD is evaluated in 33 frequency hins between zero and 64 Hz (half of the sampling rate). This leads to a setting with 4 parameters to estimate given observations defined in a 33- 33-dimensional space.
The normalizing flows qΦb and qo, used in our approximations are masked autoregressive flows (MAF) (Papamakarios al., 2017) consisting of three stacked masked autoencoders (MADE) (Ger- main et al., 2015), each with two hidden layers of 50) units, and a standard normal base distribution as input to the normalizing flow. This choice of architecture provides sufficiently flexible functions capable of approximating complex posterior distributions. We refer the reader to Papamakarios al. (2019) for more information on the different types of normalizing flows. We fix function fo, to be a simple averaging operation over the elements of X, so only parameters Φ1 and Φ2 are learned from data.
Results on simulated data. We first consider case in which the observed time series To simulated by the neural mass model with a particular choice of inpul parameters. In the lower left part of Figure 3. we display the smoothed histograms of the posterior approximation 94 obtained when conditioning on just To (N- = 0) or zo and X with N = 9. We see that when N O, parameters j and a have large variances and that some of the pairwise joint posterior distributions have a ridge shape that reveals the previously described indeterminacy relation linking g with μ and a. When N= 9, the variances of the parameters decrease and we obtain a posterior distribution that is more concentrated around the true parameters generating Zp. This concentration is explained by the sharper estimation of the g parameter, which is obtained using Tig and ten auxiliary observations.
In the upper right part of Figure 3, we evaluate how HNPE concentrates around the true parameters when N increases and plot the results using two other architectures: a "naive" architecture taking context variables a stacking of 2n and the elements in X, and an "aggregation" " architecture which stacks xn and the average of X' as context variables; both architectures use a normalizing flow with 10 layers of two hidden layers and 50 units each. We evaluate the concentration of the posterior distributions via its Wasserstein distance a Dirac located at the ground truth parameter. For each parameter in the model, i.e. j E (C. p. a. g). we have three curves, di duuive daggegg and diNeE which describe how the posterior marginal of j converges to a Dirac when N increases for each architecture. Since the parameters have very different scaling, standard value across different ground truth parameters and values of N. We then take the mean along j so to obtain three final curves duaie, daggreg. and dHNPE. We consider ten choices of ground
normalize the curves by dividing them by their
8

---

N 0
N
eyes state 4000
eys state
dooodd
2000
2000
open
-20
14
Figure 4: Posterior estimates for the parameters of the neural mass model computed human EEG signals. Data were collected in two different experimental conditions: eyes closed (in blue) or eyes open (in orange). All signals are 8 long and recorded at 128 Hz. We see that when N = 9 the posterior distributions concentrates, and that the global gain parameter gets similar in both eyes conditions. We observe that the posterior on the 3 parameters of the neural mass model clearly separate between the 2 conditions when N
truth parameters and show the curves with the normalized median distances. For N = O, the posterior the fact that the curves do not start
distribution is supposed to be indeterminate ("banana-shape") at the same point has no proper interpretation. For N > the curve HNPE uniformly below the other methods, and they converge to a plateau. This demonstrates the existence of a lower bound for the concentration of the posterior approximation, which can be interpreted as irreducible variance on the estimation of the parameters. The rather good performance for the "aggregation" architecture as compared to HNPE is likely due to the fact that, for the example considered here, taking the average of the elements in X leads sufficient statistic of the observations.
We have also considered setting which the summary statistics of the observed time series are learned from the data instead of being fixed to the log power spectral densities, i.c. when fox is learned. We have used the YuleNet architecture proposed by Rodrigues and Gramfort (2020) the example with neural mass models and report the results in Appendix B. In all our experiments, we did not see significant changes in the performance of our model so we did not include it in our evaluation increased the complexity of the model and its computational burden.
Results EEG data. One of the most commonly observed oscillations in EEG are known a waves (Lopes Silva, 1991). These waves are characterized by their frequency around 10 Hz and are typically strengthened when closing our eyes. To relate this phenomenon to the underlying biophysical parameters of the NMM model, we estimated the posterior distribution over the model parameters on EEG signals recorded during short periods of eyes open or eyes closed. Data consists of recordings taken from public datasct (Cattan al., 2018) in which subjects were asked to keep their eyes open or closed during periods of 8s (sampling frequency of 128 Hz). Results for one subject of the dataset are presented in Figure with zp being either recording with eyes closed (in blue) eyes open (in orange). We consider situations in which extra-observations posterior approximation (N = 0) or when = 9 additional observations from both eyes-closed and available. When N 9,
used for the
eyes-open conditions concentrates both eyes conditions. More interestingly, we observe that the posterior on the 3 parameters of the neural mass model clearly separate between the conditions when N = 9. Looking at parameter C, see that it concentrates around 130 for the eyes closed data while it peaks around for eyes open. This finding is perfectly inline with previous analysis of the model (Jansen and Rit, 1995). Signals used in this experiment are presented in Appendix C.
that the gain parameter, which is global,

---

Discussion
In this work, propose HNPE, likelihood-free inference approach able to leverage set of additional observations to boost the estimation of the posterior. This improvement made possible by hicrarchical model where all available obscrvations share certain global parameters. A dedicated neural network architecture based on normalizing flows is proposed, opposed to the usual approach of LFI practitioners that often choose a "one size fits all" neural density estimator. We also provide training procedure based on simulations from the model and based on the sequential approach from Greenberg et al. (2019). Although the number of additional observations (N) was fixed in our analysis and experiments, this parameter could be randomized and amortized during learning and enable the posterior approximation be fed with sets of auxiliary observations of varying sizes, making it more flexible for upplications. To do so, it would necessary to simulate datasets with varying sizes of N so to ensure that the several simulations are IID between them; note, however, that this would have a significant computational cost. Our approach could be extended multi-level models using a similar factorized architecture; we did not consider such generic hierarchical models to keep the presentation clear and because motivating examples did not require such complexity. Note, also, that HNPE could implemented with other lypes of conditional density estimators apart from normalizing flows, as long as the hierarchical structure of the global parameters is embedded into the structure of the approximator.
It well known that methods for likelihood-free inference are often difficult to validate: our method is no exception. We have considered toy models for which the analytic form of the target posterior are available to have a precise way of assessing the quality of approximation and avoiding such difficulties. Nevertheless, further research on validation schemes for L.FI methods remain great interest, specially for more general settings for which the analytic posterior is unknown. Note, also, that I.FI methods can require large number of simulations in order approximate the posterior distribution and may, therefore, lead to a non-negligible carbon footprint. This can be mitigated with the development of new methods for optimizing the number of simulations required for a given error tolerance, e.g. choosing the sampled parameters 0, for which the simulations æ; are the most useful for training the posterior approximation.
We demonstrated that HNPE could be reliably applied to neuroscience considering stochastic model with non-linear differential equations. Very encouraging results on human EEG data open the door to more biologically informed descriptions and quantitative analysis of such non-invasive recordings.
Acknowledgments and Disclosure of Funding
This work was granted access to the HPC resources of IDRIS under allocations DA1-AD01111111220 made by GENCI. GL is recipient of the ULiège- NRB Chair on Big Data and is thankful for the support of the NRB, AG thanks the support of the ERC-StG SLAB (ID:676943) and the ANR BrAIN (ANR-20-CHIA0001) grants.
References
Markus Ableidinger, Evelyn Buckwar, and Harald Hinterleitner. A stochastic version of the Jansen and Rit neural mass model: Analysis and numerics. The Journal of Mathematical Neuroscience, 7 (1), August 2017. doi:
Hannelore Aerts, Michael Schirner, Ben Jeurissen, Dirk Van Roost, Eric Achten, Petra Ritter, and Daniele Marinazzo. Modeling brain dynamics in brain tumor patients using the virtual brain. eNeuro, 5(3), June 2018. ISSN 2373-2822. Society for Neuroscience.
Erie Bazin, Kevin J Dawson, and Mark A Beaumont. Likelihood-freee inference of population structure and local adaptation in bayesian hierarchical model. Genetics, 185(2):587-602, June 2010. doi: 10. h1t34///eeeiics119999 109.11239.11111
Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman. Pyro: Deep Universal Probabilistic Programming. Journal of Machine Learning Research, 2018.
David M. Blei, Andrew Y. Ng, and Michael I. Jordan. Latent dirichlet allocation. J. Mach. Learn. Res., 3(mul)):933-1022 March 2003. ISSN 1532-4435.
10

---

Johann Brehmer, Siddharth Mishra-Sharma, Joeri Hermans, Gilles Louppe, and Kyle Cranmer. Mining for dark matter substructure: Inferring subhalo population properties from strong lenses with machine learning. The Astrophysical Journal, 886(1):49, 2019.
Evelyn Buckwar, Massimiliano Tamborrino, and Irene Tubikanec. Spectral density- density-based and measure-preserving ABC for partially observed diffusion processes. an illustration on hamil- tonian SDEs. Statistics and Computing, 30(3):627-644 November 2019. doi: 10.1007/
411222-01--09909-6.
Grégoire Cattan, Pedro L. C. Rodrigues, and Marco Congedo. EEG alpha waves dataset. December 2018. doi: 0t.5221/zenndo.220 2348892..
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-base inference. Proceedings of the National Academy of Sciences, (17748):3005533006 2020. ISSN 0027-8424.
Olivier David and Karl J. Friston. neural mass model for MEG/EEG: Neurolmage, 20(3):
1743-1755, November 2003, doi:
Gustavo Deco, Viktor Jirsa, Peter A. Robinson, Michael Breakspear, and Karl Friston. The dynamic brain: From spiking neurons to neural masses and cortical fields. PLOS Computational Biology, 4(8):1-35, 08 2008. doi: 10.137 ici.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché Buc, E. Fox, and R. Garnett, editors, Advances Neural Information Processing Systems, volume 32, pages 7511-7522, 2019.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. nflows: normalizing flows in PyTorch. November 2020a. doi: 10.528 h/weeodo. 42962887.
Conor Durkan, lain Murray, and George Papamakarios. On contrastive learning for likelihood-free inference. In Hal Daumé IIII and Aarti Singh, editors. Proceedings of the 37th International Conference on Machine Learning, volume 119 Proceedings of Machine Learning Research, pages 2771-2781. PMLR, 13-18 Jul 2020b.
Jean Feydy, Thibault Séjoumé, François-Xave Xavier Vialard, Shun-ichi Amari, Alain Trouve, and Gabriel Peyré. Interpolating between optimal transport and mmd using sinkhorn divergences. In Kamalika Chaudhuri and Masashi Sugiyama, editors, Proceedings of Machine Learning Research, volume 89, pages 2681- PMLR, 16-18 Apr 2019.
Andrew Gelman and Jennifer Hill. Data analysis using regression and multilevelhierarccical models, volume Analytical methods for social research. Cambridge University Press, New York, 2007.
Mathieu Germain, Karol Gregor, Jain Murray, and Hugo Larochelle. Made: Masked autoencoder for distribution estimation. In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 881-889, Lille, France, 07-09 Jul 2015. PMLR.
Pedro Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacherr Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F Podlaski, Sara A Haddad, Tim P Vogels, David Greenberg, and Jakob H Macke. Training deep neural density estimators to identify mechanistic models of neural dynamics. elife, 9ce56261, sep 2020. ISSN 2050-084X.
David Greenberg, Marcel Nonnenmacher, and Jakob Macke. Automatic posterior transformation for likelihood- free inference. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning. volume 97 of Proceedings of Machine Learning Revearch, pages 2404-2414 14, PMLR, 09-15 Jun 2019.
Paul Gustafson. Bayesian inference in partially identified models: Is the shape of the posterior distribution useful? Electronic Journal of Statistics, 8(1), January 2014. doi: 10. 10.1214/14eejs811.
Joeri Hermans, Volodimir Begy, and Gilles Louppe. Likelihood-fre MCMC with amortized ap- proximate ratio estimators. In Hal Daumé IIII and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume I 119, pages 4239 4248. PMLR, 13-18 Jul 2020.
Ben H. Jansen and Vincent G. Rit. Electreenceehallogamm and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics, 73(4):357-36 September 1995. doi: 10. 10.1007/bf00199441..
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 141.6994, 2014.
11

---

Paula Sanz Leon, Stuart A. Knock, M. Marmaduke Woodman, Lia Domide, Jochen Mersmann, Anthony R. McIntosh, and Viktor Jirsa. The virtual brain: simulator of primate hrain network dynamics. Frontiers in Neuroinformattis, 7, 2013. doi: 10.3389/fninf:20) 20..00000.......000000000...
Fernando Lopes da Silva. Neural mechanisms underlying brain waves: from neural membranes networks. Electreencephalograppy and Clinical Neurophysiology, 79(2):81 - 93, 1991. ISSN 0013-4694. doi: http:////oooorr/1100 0001000013349999999999999999
George Papamakarios and lain Murray. Fast e-free inference of simulation models with bayesian conditional density estimation. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances Neural Information Processing Systems, volume 29, pages 1028-1036, 2016.
George Papamakarios, Theo Pavlakou, and lain Murray, Masked autoregressive flow for density estimation. In I. Guyon, U. V. Luxburg, S. Bengio, Wallach, R. Fergus. S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 2338-234 2347. Curran Associates, Inc., 2017.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. arXiv preprint arXiv:19/2.02762, 2019.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PoTochh: An Imperative Style, Righ-Perffrmance Deep Learning Library. In Advances in Neural Information Processing Systems (NeurIPS), page 12, Vancouver, BC, Canada, 2019.
Donald B. Percival and Andrew T. Walden. Spectral Analysis for Physical Applications. Cambridge University Press, 1993. doi: 10.007//C0978051 1622762.
Python Software Fundation. Python Language Reference, version 3.6, 2017.
Pedro L. Rodrigues and Alexandre Gramfort. Learning summary features of time series for likelihood free inference. arXiv preprint arXiv:2012.0280) 2020.
Salakhutdinov, J. B. Tenenbaum, and A. Torralba. Learning with hierarchical-deey models.
IEEE Transactions on Pattern Analysis and Machine Intelligence, 1958-11951 2013. doi: 11110//TAAM..012..2200
Paula Sanz Leon, Stuart Knock, M. Woodman, Lia Domide, Jochen Mersmann, Anthony McIntosh, and Viktor Jirsa. The virtual brain: simulator of primate brain network dynamics. Frontiers in Neuroinformatiis, 7:10, 2013. ISSN 1662-5196. doi: 10.3389/fninf..213..000
Yee Whye Teh and Michael Jordan. Hierarchical bayesian nonparametric models with applications. Bayesian nonparametrics, 1: 1:158-207, 2010.
Alvaro Tejero- Caneteo,, Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pedro Gonçalves. David S. Greenberg, and Jakob H. Macke. sbi: A toolkit for simulation based inference. Journal of Open Source Software, 5(52):2505, 2020. doi: 10.21 (05/soss.02505.
Dustin Tran, Rajesh Ranganath, and David Blei. Hierarchical implicit models and likelihood-free variational inference. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vish- wanathan, and Garnett, editors, Advances Neural Information Processing Systems, volume 30, pages 5523-5533. Curran Associates, Inc., 2017.
Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ Salakhutdinov, and Alexander Smola. Deep sets. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30, pages 3391 1-301. Curran Associates, Inc., 2017.
12


Reviews:
[{"review_id": "bsNUpuE9VUa", "reviewer": "Reviewer_tXAy", "summary": "The paper extends SNPE-C to the situation where multiple sets of parameters yield identical observations, but where the parameters can be split into local and global ones, which are specific to, or shared across all observations, respectively. It is shown how this additional information can be used to build a flow-based posterior estimator yielding sharper distributions.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper considers simulation-based inference in the presence of strongly coupled parameters leading to a non-injective likelihood function. This parameter estimation problem for such a system is ill-posed, but under additional assumptions this can be resolved with a hierarchical Bayesian model in which some parameters are shared across all observations.\n\nThe authors propose a practical realization of this solution in the form of a \"h-flow\" model, consisting of two normalizing flows modeling the posterior for the local and global parameters, and utilizing a deepset architecture for the global parameters to enforce permutation invariance. The network is optimized with KL divergence, and the authors show how to train it within the multi-round SNPE-C scheme.\n\nh-flow is then experimentally tested on a simple model dependent on the product of two parameters, and on a 4d neural-mass model (NMM) of a cortical column with 3 local parameters and 1 global one. Performance is comparable to a naive normalizing flow in the simple case, but a significant improvement is seen for the NMM. To further show practical applicability, the authors then use h-flow to estimate NMM parameters on real EEG signals recorded under eyes open/closed conditions, and show that this results in separated distributions for the two conditions, in line with expectations.\n\nThe paper is clearly written, provides simple and illustrative motivating examples, and keeps the most relevant information in the main text. The proposed method is practically useful, and easy to implement, which should lead to adoption among practitioners already using SNPE-C or similar SBI methods. The experimental results could be made a bit more convincing with additional models or more baseline methods in the case of the NMM.\n\n----\nScore revised +1 based on the authors' responses and the promise to report improved baselines in the NMM test case.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No concerns seen here. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "bm12xJWpyP", "reviewer": "Reviewer_fnEF", "summary": "This paper proposes an approach to perform statistical inference for hierarchical models with intractable likelihoods. It extends an existing approach for simulation-based inference (SNPE-C) in order to use models that have global and local random variables. \n\nThe authors propose a scheme in which two conditional flows a trained: The first flow models a distribution over global parameters given an observation and $N$ permutation invariant auxiliary observations. Permutation invariance could e.g. be accounted for with a deepset architecture. The experiments do not make use of this flexibility, relying on an average operation instead. The second flow then models a distribution over local parameters given the observation and global parameters. They refer to the resulting inference scheme jointly as h-Flow, for hierarchical flow.\n\nA motivating toy example is provided, along with a relevant, more realistic simulator from Neuroscience (a Neural Mass Model; NMM). h-Flow is compared against a naive approach on both of these examples. The comparison turns out about equal on the toy example and favorably for h-Flow on the NMM. On the toy example, h-Flow is additionally compared against a hierarchical ABC variant and LFVI, comparing favorably to both of them.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "# Main Review\n\n## Significance and Originality\n\nOverall, the paper is well motivated and introduces a strategy which can potentially be useful for many applications of SBI. The proposed solution might come at no big surprise, but in fact this can be viewed as a strength: As long as the proposed strategy of training two flows (see summary above) performs robustly and demonstrably better than what one would naively do with SNPE-C, and given that it outperforms existing methods, I would view this paper as being a valuable addition to the conference that should be presented. \n\n\n## Clarity\n\nThe paper is generally well-written and easy to follow. I especially enjoyed that the authors introduced a motivating example that then finds its structure reflected in the applied setting.\n\n\n## Quality\n\nMy reasons for not giving a higher score, and currently tending to rather reject the paper, are many open questions and some specific requests I have regarding the comparisons. I detail these concerns below and hope that the authors will address them, in which case I will revise my evaluation accordingly.\n\n### Comparisons to Naive Approaches\n\nMy primary concern is that the paper does not yet clearly demonstrate that the proposed method is better than naive approaches using SNPE-C with a single flow (instead of using two flows, as h-Flow does). In the toy experiment, no difference between h-Flow and the naive approach is found (Figure 2). On the NMM, a difference is clearly visible in Figure 3, but as of now I am unsure how to interpret this result exactly and how it was obtained. Some specific questions and requests that can hopefully help clarify this issue follow.\n\n**Q1**: In Figure 3, the difference between the naive approach and h-Flow _is biggest_ for $N=0$. I was surprised by this, since my understanding is that we expect improvements in resolving indeterminacies by additional observations, i.e. when $N>0$. In the toy example, which served as a motivation for the NMM, as $N \\rightarrow 0$ the posterior distribution converges to $p(\\alpha, \\beta | x_0)$. In other words, for the toy example, we do not expect a difference (and empirically don't find) a difference in performance when $N=0$. Could the authors comment on the difference in performance for $N=0$ on the NMM? \n\n**Q2**: Regarding the results presented in Figure 3 the authors write: \"We see that the factorization scheme yields uniformly better results than the naive architecture. This is mainly because the dimensionality of the context variables in the naive case gets too big when more observations become available, making the training of the normalizing flow in the non-factorized setting harder.“. Why exactly do the context variables in the naive case get too big with more observations? I would have thought that for the naive scheme we simply aggregate the additional observations before passing them as context to the single flow. If we do not aggregate them, we assume that the ordering matters, which is would put the naive implementation at a big disadvantage and is not how one would naively apply SNPE-C with a single flow in this context. It would be important to comment and clarify this.\n\n**R1**: As a request for new experiments: The results in Figure 3 are reported in terms of Wasserstein distance between a dirac delta on ground truth parameters and the inferred posterior (median without error bars for 5 repeats). It would be important to extend this analysis for a much larger number of ground truth parameters and report results such that the variability can be evaluated (e.g. mean and standard errors). Ideally, a secondary metric to the Wasserstein distance was reported in the supplement, such as the more commonly used negative log probability of true parameters _averaged_ over as many ground truth parameters as feasible (this is crucial, otherwise this metric would be hard to interpret).\n\n### Results for h-ABC and LFVI on NMM\n\n**R2**: As a second request, it would be good to provide results for h-ABC and LFVI on the NMM model for completeness.\n\n\n## Additional questions\n\n**Q3**: In Figure 1 (right) we observe that the $N=10$ case is worse relative to $N=0$. To me this seemed counterintuitive at first and I think it would be good if the authors commented on it in the manuscript. My interpretation is that the density estimation problem gets more difficult when we have more observations (posterior gets more concentrated). So resolving indeterminacies leads to less accurate inference (as compared to the analytic solution). Would the authors agree with this interpretation?\n\n**Q4**: \"We did not carry out more experiments on this data-driven setting because of difficulties due to numerical instabilities in the training procedure when N increases and for certain choices of ground truth parameters.\". It would be important to investigate this issue further, if possible for the current publication. Why would this only occur for certain ground truth parameters and not for others? \n\n**Q5**: For the toy example the authors set $\\sigma=0$, so that $\\epsilon = 0$. Usually SBI is only applied in the context of stochastic simulators. Would results change in a meaningful way if the experiments were repeated for a small value of $\\sigma$?\n\n**Q6**: The NSF in Figure 1 does not seem able to handle the sharp transitions that occur when $N=10$ or $N=100$. Do the authors have any suggestions to remedy this? \n\n**Q7**: On the NMM, MAFs instead of the more expressive NSFs are used (whereas the toy example does use NSFs). What are the reasons behind this choice and how do results change for h-Flow and the naive strategy when using them instead?\n\n**Q8**: \"Although the number of additional observations (N) was fixed in our analysis and experiments, this parameter could be randomized and amortized during learning and enable the posterior approximation to be fed with sets of auxiliary observations of varying sizes making it more flexible for applications.“ $\\rightarrow$ My understanding is that the aggregation is simply done through the average in the experiments. What exactly would not work/require changes to have flexible N?\n\n**Q9**: Have the authors done any internal experiments in which $\\phi_3$ was learned as well? If so, were there any problems that they encountered that restricted the experiments to learning $\\phi_1$ and $\\phi_2$ only?\n\n\n## Name of the algorithm\n\nI see two potential disadvantages associated with with naming the algorithm h-Flow. They are:\n1. One of the strengths of SNPE is that it can be used with any conditional density estimator. In the future, perhaps entirely different classes of CDEs will be invented. I might be misunderstanding something (please correct me if I'm wrong) but to me it seems there is nothing flow-specific about the proposed algorithm. In particular I assume that we can use it with all kinds of CDEs as long as they work with embeddings/can be designed to be permutation invariant and we can evaluate log probabilities and draw samples. For example, we could just as well use this with MDNs, which might be more suitable depending on the application/knowledge we have about a problem.\n2. h-Flow in the paper title sounded to me as if a single architecture was trained, i.e. that this paper was about an entirely novel flow architecture. Instead, two flows are trained, and if the approach gets extended to multi-level architectures, training would involve more and more flows, not a single one. \n\n**Q10**: What are the authors thoughts on these remarks regarding their algorithm's name?\n\n## Minor\n\n- Fig. 3: There are display issues with this figure when opening the PDF using Preview on macOS. Most marginals are missing and the x-axis of the inset figure is cut off (3 instead of 30)\n- L11: validate quantitatively our proposal $\\rightarrow$ validate our proposal quantitatively\n- L103: have to be approximated $\\rightarrow$ we decide to approximate them\n- L178: $5 \\times 10^{−4}$ rather than $5.10^{−4}$?\n- L227: captures well the global parameter $\\rightarrow$ captures the global parameter well \n- L229: could well estimate local parameters $\\rightarrow$ could estimate local parameters well\n\n## Update\n\nI thank the authors for their reply and have updated my score.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors state in their checklist that potential negative societal impacts of their work are not applicable. I agree that an extended discussion of negative impacts might not be strictly necessary, given the nature of this work is about statistical inference applied to a certain class of models (with hierarchical structure and intractable likelihoods). Thinking more broadly, the societal impact and potentially adverse effects will depend on the setting in which the algorithm gets deployed. However, I would suggest to at least briefly touch upon the general difficulties associated with diagnostics/criticism in the SBI context.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "7"}, "extra_text": {}}, {"review_id": "3AyBEaIU6a", "reviewer": "Reviewer_u4up", "summary": "The paper addresses the inference of illposed models, where the likelihood is intractable, and proposes to utilize hierarchical structures in combination with neural density estimators (normalizing flows) to ameliorate identifiability of the parameters. It motivates the approach with a simple analytical tractable example, and validates the approach on modelled and real EEG time-series.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The problem the authors consider is interesting and important for scenarios, where models are parametrized ambiguously. However, I consider the present work as too immature for publication, because to me, the proposed approach (making use of additional observations and separating the parameters in global and local ones) is in my opinion not sufficiently justified explored. Some questions I would like to see answered are\n+ When is the separation off variables into local and global possible or desirable. E.g. when you remove the assumption that alpha and beta are upper bounded by 1 in the motivating example, the separation wouldn’t help any more in my opinion. I do not say, that the motivating example isn’t interesting, but I would like to have seen more discussion, on when it makes sense to have this separation. Hence, I would appreciate some general discussion, when such the proposed separation makes sense.\n+ The problems, that are considered, are very low dimensional, which is making me wonder whether the neural density estimators are really necessary here. Just demonstrating that the separation ameliorates the inference by using Gaussians (or mixtures of them, or other simple estimators) for the density approximation in Eq. (5), might be more insightful, and easier to train.\n+ Overall, I was missing details about how many parameters are used in the different models, how long it takes them to train, etc. So I couldn’t make myself a picture, whether comparisons are fair or not.\n+ For the fact, that the paper doesn’t provide any theory (except for a specific toy example), how this separation is helping model inference, I think the empirical results are too little, to consider the proposed method as helpful in general.\n\n### Update\n\nReading the comments of the other reviewers, and the response of the authors, I am reconsidered my score and raised it to 6 (though decreasing my confidence score). This is because,\n\n+ Reading the reviewers' comments, it seems I underestimated the contribution of the paper, leveraging information of the hierarchical structure of models in the context of likelihood-free inference.\n+ I mentioned originally, that I found the manuscript written too complicated and mentioned as examples invariance of the averaging operation to permutations, and that $\\phi_3$ is introduced, but seems to be never learnt. I understand for the first point, that though trivial, it is a crucial point for the construction of the neural posterior, so it makes sense to give the fact more space than needed on first sight. The fact, that $\\phi_3$ is not needed in the shown experiments (as I got from the response to reviewer fnEF), still poses the question, whether it is necessary to introduce it, or just to mention the possibility of having, but simplify notation.\n+ Experiments seem to be computational heavy to provide new baselines for the NMM example, so it is maybe too much to ask.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "+ As mentioned above, I’d like to see some additional and rigour discussion on when such separation of variables helps the inference.\n+ Please report number of parameters of models, what device was used to train them etc. to allow the reader to evaluate the showed results.\n+ In general, more experiments, that demonstrate the improvement achieved by the proposed method \n+ I sometimes feel, the paper is written too complicated. For example, in l134 it is written that $f_{\\phi^{(3)}}$ is “invariant to the ordering of the input”, which is obvious, since the function in question is a sum over inputs (which is commutative). I thought quite a bit about this comment, and I am not sure whether I am not seeing something, or why it is important to mention here. Another example is, that $f_{\\phi^{(3)}}$ seems never to be used in the experiments, which makes me wonder, why is it introduced in the first place.\n+ Some typos:\n++ L129f: indices of phi are mixed\n++ There are several approximation signs in the appendix that should be equal, and some equal signs that should be approximations.\n++ L107: It should be \\beta not \\beta_0\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5h"}, "extra_text": {}}, {"review_id": "2ZXNB2KHj5D", "reviewer": "Reviewer_2NmZ", "summary": "The authors propose a new likelihood-free inference method for Bayesian hierarchical models, called h-Flow, that exploits additional information provided by a set of auxiliary observations with shared, global parameters. h-Flow works by approximating two relevant posterior distributions, i.e. that of the global parameter beta and the local parameters alpha given beta, with separate normalizing flows, which are then multiplied to yield the joint posterior distribution. The normalizing flows are trained by minimising the KL-divergence between the true and approximated posterior distribution and the training data is refined over several rounds by means of SNPE-C (also known as APT).", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The authors provide an approach to likelihood-free inference, which is a known and important task, for Bayesian hierarchical models. While the method is new, to me it appears to be a slight variation of SNPE-C that includes a second normalizing flow for global parameters that is conditioned on a set of auxiliary variables. I am thus mildly worried about novelty and significance. I would encourage the authors to emphasise their contributions and the difference to related work further.\n\nThe derivations are generally clear and technically sound. While the examples in their experiments are very interesting, I don't think they highlight the method very well. It is still unclear to me whether or not any improvements are due to the proposed method, or a) an increased flexibility or b) simply having more auxiliary observations that are correlated with new observations. The authors attempt to clarify this by using a naive approach that ignores the hierarchical structure, but this only consists of one normalizing flow with a lower capacity than the two flows of h-Flow. Similarly, as the authors said, N>0 generally yields better posteriors but this is naturally not very surprising. In my opinion, a fair baseline for e.g. the EEG experiments would have been a comparison of N=0 vs N=9 with both naive architectures and h-Flow, but a setting where both approaches have the same capacity.\n\nThe submission is generally clearly-written, at least until Section 3. I believe the experimental section could have been more concise with some explanations about the models and approaches being moved to the appendix, leaving space for, perhaps, one more experiment or more elaborate baseline comparisons.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors did not comment on limitations of their method and the potential societal impact it may have. \n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Dear authors, 

congratulations on your paper being accepted at Neurips. As you know, the reviewers gave extensive comments and feedback on the submission, and we urge you to ensure that this feedback is incorporated into the final version of the article. In particular, reviewers felt that the name of the method, 'h-flow', is misleading, as the paper does not present a new flow, but rather a new SBI method (which will likely, but not necessarily, be applied to flows). We would ask you to take this advice seriously. 

With best regards and congratulations, your AC
Decision: Accept (Poster)