HNPE: Leveraging Global Parameters for Neural Posterior Estimation
Year: 2022

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 simulationbased 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.

heading: 1 Introduction
text: 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 to a mathematical model implemented as a simulator. When applied to biophysical models and simulators in neuroscience (e.g. Leon et al. 2013), it could estimate properties of 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.
SBI is concerned with the estimation of a conditional distribution over parameters of interest ✓. Given some observation x0, the goal is to compute the posterior p(✓|x0). It generally happens that some of these parameters are strongly coupled, leading to very structured posteriors with low dimensional sets of equally likely parameters values. For example, this happens when the data generative process depends only on the products of some parameters: multiplying one of such parameters by a constant and another by its inverse will not affect the output. Performing Bayesian inference on such models naturally leads to a “ridge” or “banana shape” in the posterior landscape, as seen e.g. in Figure 4 of Gonçalves et al. (2020). More formally the present challenge is posed as soon as the model likelihood
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
function is non-injective w.r.t. ✓, and is not strictly due to the presence of noise on the output observations. In statistics and econometrics literature, such models are called partially identified models (Gustafson, 2014).
To alleviate the 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 ✓i for an observation xi are partitioned into ✓i = {↵i, }, where ↵i is a set of sample specific (or local) parameters, and corresponds to shared (or global) parameters. For this broad class of hierarchical models, the posterior distribution for a set X = {x1, . . . , xN} of N observations can be written as (Tran et al., 2017):
p(↵1, . . . ,↵N , |X ) / p( ) NY
i=1
p(xi|↵i, )p(↵i| ). (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 (Salakhutdinov 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(xi|↵i, ) is implicit and intractable, which implies that traditional MCMC methods can not be used to estimate the posterior distribution. This setup leads to so-called likelihood-free inference (LFI) problems and many algorithms (Papamakarios and Murray, 2016; Greenberg et al., 2019; Hermans et al., 2020; Durkan et al., 2020b) have recently been developed to carry out inference under this scenario. These methods all operate by learning parts of the Bayes’ rule, such as the likelihood function, the likelihood-to-evidence ratio, 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 hierarchical implicit models using a GAN-like training approach, while Brehmer et al. (2019) and Hermans et 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 a sequential neural posterior estimation approach derived from SNPE-C (Greenberg et al., 2019).
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 a 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.

heading: 2 Hierarchical models with global parameters
text: Motivating example. Consider a stochastic model with two parameters, ↵ and , that generates as output x = ↵ + ", where " ⇠ N (0, 2). We assume that both parameters have uniform prior distribution ↵, ⇠ U [0, 1] and that is known and small. Our goal is to obtain the posterior distribution of (↵, ) for a given observation x0 = ↵0 0+". This simple example describes common situations where indeterminacy emerges. For instance, x0 could be the radiation power measured by a sensor, ↵ the intensity of the emitting source, and the inverse squared distance of the sensor to the source. In this case, a given measurement may have been due to either close weak sources (↵ # and ") or far strong ones (↵ " and #). Using Bayes’ rule and considering small we can write (see Appendix A for more details)
p(↵, |x0) ⇡ e (x0 ↵ )
2/2 2
p 2⇡ 2
1[0,1](↵)1[0,1]( )
log(1/x0) , (2)
where 1[a,b](x) is an indicator function that equals one for x 2 [a, b] and zero elsewhere. Note that the first term in the product converges to (x0 ↵ ) as ! 0 and that the joint posterior distribution has
an infinite number of pairs (↵, ) with the same probability, revealing the parameter indeterminacy of this example. Indeed, for x 2 [0, 1] and 2 [x, 1], all pairs of parameters ( x , ) yield the same observations and the likelihood function p(·| x , ) is constant. Thus, the posterior distribution has level sets with a ridge or “banana shape” along these solutions. The top row of Figure 1 on Page 4 portrays the joint and the marginal posterior distributions when (↵0, 0) = (0.5, 0.5) and = 0.
Exploiting the additional information in X . Our motivating example illustrates a situation where two parameters are related in such a way that one may not be known without the other. In practice, however, it is possible that one of 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 is 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 as 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 (↵0, ) given an observation x0 of a stochastic model, we consider the following scenario. We assume that the model’s structure is such that ↵0 is a parameter specific to each observation (local), while is shared among different observations (global). Yet, both are unknown. We consider having access to a set X = {x1, . . . , xN} of additional observations generated with the same as x0. Taking the model’s hierarchical structure into account we use Bayes’ rule to write
p(↵0, |x0, X ) = p(↵0| , x0, X )p( |x0, X ) / p(↵0| , x0)p(x0, X| )p( )
/ p(↵0| , x0)p( ) QN i=0 p(xi| ) / p(↵0, |x0)p( ) N QN i=1 p( |xi)
(3)
which shows how the initial posterior distribution p(↵0, |x0) is modified by additional observations from X sharing the same as x0. 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(↵0| , x0) and p( |x0, X ) have to 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 + 1 observations in {x0} [ X are such that xi = ↵i 0 + " for i = 0, . . . , N with ↵i ⇠ U [0, 1] drawn from the same prior. The posterior distribution may be written as (see Appendix A)
p(↵0, |x0, X ) ⇡ p(↵0, |x0) 1[µ,1]( ) N N log(1/x0)
(1/µN 1) , (4)
where µ = max({x0} [ 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(↵0, |x0). Figure 1 portrays the joint and marginal posterior distributions with N = 10 and N = 100.

heading: 3 HNPE : neural posterior estimation on Bayesian hierarchical models
text: 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 likelihood-free inference (LFI) to directly estimate an approximation to the posterior distribution using a conditional neural density estimator trained over simulations of the model. In what follows, we present a novel neural network architecture for approximating the posterior distribution of a hierarchical model with global parameters based on normalizing flows. We also describe the training procedure for learning the parameters of the network using a multi-round procedure known as sequential neural posterior estimation or SNPE-C (Greenberg et al., 2019).

heading: 3.1 Approximating the posterior distribution with two normalizing flows
text: We approximate p(↵0, |x0, X ) based on its factorization (3) as follows: p( |x0, X ) ⇡ q 1( |x0, f 3(X ))
p(↵0| , x0) ⇡ q 2(↵0| , x0) (5)
where q 1 and q 2 are normalizing flows, i.e., invertible neural networks capable of transforming data points sampled from a simple distribution, e.g. Gaussian, to approximate any probability density function (Papamakarios et al., 2019). The function f 3 is a deepset neural network (Zaheer et al., 2017) structured as f 3(X ) = g (1)3 ⇣ 1 N PN i=1 h (2)3 (xi) ⌘ , where h is a neural network parametrized by (1)3 that generates a new representation for the data points in X and g is a network parametrized by (2)3 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 as the maximum value of the set or 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 f 3 is an universal approximator invariant to the ordering of its inputs (Zaheer et 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 our approximation either by its factors q 1 and q 2 or by q with = { 1, 2, 3}.
Estimating . We estimate by minimizing the average Kullback-Leibler divergence between the posterior distribution p(↵0, |x0, X ) and q (↵0, |x0, X ) for different x0 and X :
min. Ep(x0,X ) h KL(p(↵0, |x0, X )kq (↵0, |x0, X )) i ,
where KL(pkq ) = 0 if, and only if, p(↵0, |x0, X ) = q (↵0, |x0, X ). We may rewrite the optimization problem in terms of each of its parameters to get
min. 1, 2, 3 L↵( 2) + L ( 1, 3) (6)
with L↵( 2) = Ep(x0,X ,↵0, ) [log(q 2(↵0| , x0)] ,
L ( 1, 3) = Ep(x0,X ,↵0, ) [log(q 1( |x0, f 3(X )))] .
Training from simulated data. In practice, we minimize the objective function in (6) using a Monte-Carlo approximation with data points generated using the factorization p(x0, X ,↵0, ) = p( ) QN i=0 p(xi|↵i, )p(↵i| ) where p(↵i, ) = p(↵i| )p( ) is a prior distribution describing our initial knowledge about the parameters, and p(xi|↵i, ) is related to the stochastic output of the simulator for a given pair of parameters (↵i, ). More concretely, the training dataset is generated as follows: First, sample a set of parameters from the prior distribution such that (↵ji ,
j) ⇠ p(↵i, ) with j = 1, . . . , n and i = 0, . . . , N . Then, for each (i, j)-pair, generate an observation from the stochastic simulator xji ⇠ p(x|↵ j i ,
j) so that each observation xj0 is accompanied by its corresponding N extra observations X j = {xj1, . . . , x j N}. The losses L↵ and L are then approximated by
Ln↵ = 1n Pn j=1 log(q 2(↵ j 0| j , x j 0)) and Ln = 1n Pn j=1 log(q 1( j |xj0, f 3(X j))) .

heading: 3.2 Refining the approximation with multiple rounds
text: The optimization strategy above minimizes the KL divergence between the true posterior distribution p and the approximation q , on average, for all possible values of x0 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 of q to better estimate the posterior distribution for a specific choice of x̃0 and X̃ .
We target the approximation q to x̃0 and X̃ using an 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) at round r are obtained from simulations with ↵0, ⇠ q (r 1)(↵0, |x̃0, X̃ ). At round r = 0, parameters ↵0 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 ↵ji drawn from the original prior distribution p(↵i| ), since the posterior distribution returned by the multi-round procedure is only targeted for observation x̃0. 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 to our case) of the targeted version of q to the correct posterior density p(↵0, |x̃0, X̃ ) as the number of simulations per round tends to infinity. Algorithm 1 describes the procedure for obtaining q(↵0, |x̃0, X̃ ) after R rounds of n simulations.
Algorithm 1: Sequential posterior estimation for hierarchical models with global parameters
Input :observation x̃0, X̃ , prior p(0), simulator S 1 for round r = 1 to R do 2 for sample j = 1 to n do 3 Draw xj0 = S(↵ j 0, ) for (↵ j 0,
j) ⇠ p(r 1); 4 Draw a set of extra observations X j =
S(↵ji , j) for ↵ j i ⇠ p(0)(·| j) N i=1
; 5 Train q (r) to minimize Ln↵ + Ln ; 6 Set next proposal p(r) = q (r)(·|x̃0, X̃ ); 7 return posterior q (R)(·|x̃0, X̃ )

heading: 4 Experiments
text: All experiments described next are implemented with Python (Python Software Fundation, 2017) and the sbi package (Tejero-Cantero et al., 2020) combined with PyTorch (Paszke et al., 2019), Pyro (Bingham et al., 2018) and nflows (Durkan et al., 2020a) for posterior estimation1. In all
1Code is available in the supplementary materials.
experiments, we use the Adam optimizer (Kingma and Ba, 2014) with default parameters, a learning rate of 5.10 4 and a batch size of 100. The code required for reproducing most of the results presented in the paper is available at https://github.com/plcrodrigues/HNPE

heading: 4.1 Results on the motivating example
text: To evaluate the impact of leveraging multiple observations when estimating the parameters of a hierarchical model, we use the model presented in Section 2, where the observation x0 is obtained as the product of two parameters ↵0 and 0 with independent uniform prior distributions in [0, 1] (we consider the case where = 0). The set of extra observations X = {xi}Ni=1 is obtained by fixing the same global parameter 0 for all xi and sampling local parameters ↵i from the prior distribution.
Our approximation to the posterior distribution consists of two conditional neural spline flows of linear order (Durkan et al., 2019), q 1 and q 2 , both conditioned by dense neural networks with one layer and 20 hidden units. We use neural spline flows because of the highly non-Gaussian aspect of the analytic marginal posterior distributions, which can be well captured by this class of normalizing flows. In general, however, the true posterior distribution is not available, so using other classes of normalizing flows might be justifiable, especially if one’s main goal is simply to identify a set of parameters generating a given observation. We set the function f 3 to be simply an averaging operation over the elements of X as the observations in this case are scalar, so the only parameters to be learned in Algorithm 1 are 1 and 2.
We first illustrate in Figure 1 the analytic posterior distribution p(↵, |x0, X ) and the approximation q (↵, |x0, X ) with an increasing number of extra observations (N = 0, 10, 100). For N = 0, i.e. only x0 is available, we observe a ridge shape in the joint posterior distribution, which is typical of situations with indeterminacies where all solutions ( x , ) have the same probability. The addition of a few extra observations resolves this indeterminacy and concentrates the analytic posterior distribution on a reduced support [x0,min(1, x0µ )] ⇥ [µ; 1], where µ = max({x0} [ X ). Moreover, on this support, the solutions are no longer equally probable due to the N factor that increases the probability of solutions close to µ. Also note that our estimated posterior is close to the analytic one in all cases.
To have a quantitative evaluation of the quality of our approximations q , in Figure 1 we plot the Sinkhorn divergence (Feydy et al., 2019) W✏ for ✏ = 0.05 between the analytical posterior p(↵, |x0, X ) and our approximation for different numbers of simulations per round (cf. Algorithm 1). The curves display the median value for nine different choices of (↵0, 0) 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 leads to better approximations of the true posterior distribution for a fixed observation.
Our next analysis assesses how the posterior approximation concentrates around a given point in the (↵0, ) space as the number of 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 x0; we consider the results on nine different choices of (↵0, 0) but display only the median results. We see that the distance to the Dirac for the global parameter decreases as more observations are added to X , but for the local parameter ↵0 it stabilizes
on a 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 x0 corrupted by additive noise (↵0 ⇡ x0/ ). We compare our method (HNPE) with three other approaches: a naive posterior estimation using a single normalizing flow with the same capacity as the approximation with q 1 and q 2 , i.e. two layers with 20 hidden units each, in which we stack the observations from x0 and the average from those in X as context variables, the hierarchical ABC (h-ABC) proposed in Bazin et al. (2010) and the likelihood-free variational inference (LFVI) presented in Tran et al. (2017). The flow-based approaches are trained with R = 5 rounds of 104 simulations and h-ABC has the same simulation budget with acceptance rate of 1%. We see that the naive approach has very similar performance to HNPE, mainly due to the low dimensionality of the example being considered (in Section 4.2 we show an example where the naive architecture has similar performance to HNPE as well but taking much longer to train). LFVI captures well the global parameter as N increases, but it performs poorly for the 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 et 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 attain results comparable to the other methods.

heading: 4.2 Inverting a non-linear model from neuroscience
text: We consider a class of non-linear models from computational neuroscience known as neural mass models (Jansen and Rit, 1995) (NMM). These models of cortical columns consist of a set of physiologically motivated stochastic differential equations able to replicate oscillatory electrical signals observed with electroencephalography (EEG) or using intracranial electrodes (Deco et al., 2008). Such models are used in large-scale simulators (Sanz Leon et 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 et al., 2018). In what follows, we focus in the stochastic version of such models presented in Ableidinger et al. (2017) and use the C++ implementation in the supporting code of Buckwar et al. (2019). In simple terms, the NMM that we consider may be seen as a generative model taking as input a set of four parameters and generating as output a time series x. The parameters of the neural mass model are:
• 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 a direct effect on the amplitude of x.
• g represents a gain factor relating the amplitude of the physiological signal s generated by the system of differential equations for a 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 x0 could be generated with larger (smaller) values of g and smaller (larger) values of µ and . Fortunately, it is common to record several chunks of signals within an experiment, so other auxiliary signals x1, . . . , xN 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 ↵ = (C, µ, ) and = g.
In what follows, we describe the results obtained when approximating the posterior distribution p(C, µ, , g|x0, X ) with Algorithm 1 using R = 2 rounds and n = 50000 simulations per round. Each simulation corresponds to 8 seconds of a signal sampled at 128 Hz, so each simulation outputs a vector of 1024 samples. The prior distributions of the parameters are independent uniform distributions defined as:
C ⇠ U(10, 250) µ ⇠ U(50, 500) ⇠ U(0, 5000) g ⇠ U( 30,+30)
where the intervals were chosen based on a review of the literature on neural mass models (Jansen and Rit, 1995; David and Friston, 2003; Deco et al., 2008). Note that the gain parameter g is given in decibels (dB), which is a standard scale when describing amplifiers in experimental setups. We have, therefore, that x(t) = 10g/10s(t).
It is standard practice in likelihood-free inference to extract summary features from both simulated and observed data in order to reduce its dimensionality while describing sufficiently well the statistical behavior of 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 bins 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-dimensional space.
The normalizing flows q 1 and q 2 used in our approximations are masked autoregressive flows (MAF) (Papamakarios et al., 2017) consisting of three stacked masked autoencoders (MADE) (Germain 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 et al. (2019) for more information on the different types of normalizing flows. We fix function f 3 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 a case in which the observed time series x0 is simulated by the neural mass model with a particular choice of input parameters. In the lower left part of Figure 3, we display the smoothed histograms of the posterior approximation q obtained when conditioning on just x0 (N = 0) or x0 and X with N = 9. We see that when N = 0, parameters µ and 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 . When N = 9, the variances of the parameters decrease and we obtain a posterior distribution that is more concentrated around the true parameters generating x0. This concentration is explained by the sharper estimation of the g parameter, which is obtained using x0 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 as context variables a stacking of x0 and the elements in X , and an “aggregation" architecture which stacks x0 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 to a Dirac located at the ground truth parameter. For each parameter in the model, i.e. j 2 {C, µ, , g}, we have three curves, djnaive, d j aggreg, and djHNPE, 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, we normalize the curves by dividing them by their standard value across different ground truth parameters and values of N . We then take the mean along j so to obtain three final curves dnaive, daggreg, and dHNPE. We consider ten choices of ground
truth parameters and show the curves with the normalized median distances. For N = 0, the posterior distribution is supposed to be indeterminate (“banana-shape"), so the fact that the curves do not start at the same point has no proper interpretation. For N > 0, the curve HNPE is 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 an 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 to a sufficient statistic of the observations. We have also considered a setting in 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.e. when f 3 is learned. We have used the YuleNet architecture proposed by Rodrigues and Gramfort (2020) on 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 as it increased the complexity of the model and its computational burden.
Results on EEG data. One of the most commonly observed oscillations in EEG are known as ↵ waves (Lopes da 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 4 model parameters on EEG signals recorded during short periods of eyes open or eyes closed. Data consists of recordings taken from a public dataset (Cattan et al., 2018) in which subjects were asked to keep their eyes open or closed during periods of 8 s (sampling frequency of 128 Hz). Results for one subject of the dataset are presented in Figure 4 with x0 being either a recording with eyes closed (in blue) or eyes open (in orange). We consider situations in which no extra-observations are used for the posterior approximation (N = 0) or when N = 9 additional observations from both eyes-closed and eyes-open conditions are available. When N = 9, we see that the gain parameter, which is global, concentrates for both eyes conditions. More interestingly, we observe that the posterior on the 3 parameters of the neural mass model clearly separate between the 2 conditions when N = 9. Looking at parameter C, we see that it concentrates around 130 for the eyes closed data while it peaks around 70 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.
Discussion
In this work, we propose HNPE, a likelihood-free inference approach able to leverage a set of additional observations to boost the estimation of the posterior. This improvement is made possible by a hierarchical model where all available observations share certain global parameters. A dedicated neural network architecture based on normalizing flows is proposed, as opposed to the usual approach of LFI practitioners that often choose a “one size fits all” neural density estimator. We also provide a 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 to be fed with sets of auxiliary observations of varying sizes, making it more flexible for applications. To do so, it would be 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 to multi-level models using a similar factorized architecture; we did not consider such generic hierarchical models to keep the presentation clear and because our motivating examples did not require such complexity. Note, also, that HNPE could implemented with other types 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 is 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 so to have a precise way of assessing the quality of our approximation and avoiding such difficulties. Nevertheless, further research on validation schemes for LFI methods remain of great interest, specially for more general settings for which the analytic posterior is unknown. Note, also, that LFI methods can require a large number of simulations in order to 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 ✓i for which the simulations xi are the most useful for training the posterior approximation.
We demonstrated that HNPE could be reliably applied to neuroscience considering a 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 2021-AD011011172R1 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-CHIA0016) grants.

Review summaries:
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).

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.   The 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.  A 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.