Efficient identification of informative features in simulation-based inference
Authors: Jonas Beck
Michael Deistler
Yves Bernaerts
Jakob H. Macke
Philipp Berens
Venue: NeurIPS 2022 Accept (2022)
Keywords: Simulation-based Inference
Feature Selection
Hodgkin-Huxley Model
Mechanistic Models
Neuroscience

Abstract: Simulation-based Bayesian inference (SBI) can be used to estimate the parameters of complex mechanistic models given observed model outputs without requiring access to explicit likelihood evaluations. A prime example for the application of SBI in neuroscience involves estimating the parameters governing the response dynamics of Hodgkin-Huxley (HH) models from electrophysiological measurements, by inferring a posterior over the parameters that is consistent with a set of observations. To this end, many SBI methods employ a set of summary statistics or scientifically interpretable features to estimate a surrogate likelihood or posterior. However, currently, there is no way to identify how much each summary statistic or feature contributes to reducing posterior uncertainty. To address this challenge, one could simply compare the posteriors with and without a given feature included in the inference process. However, for large or nested feature sets, this would necessitate repeatedly estimating the posterior, which is computationally expensive or even prohibitive. Here, we provide a more efficient approach based on the SBI method neural likelihood estimation (NLE): We show that one can marginalize the trained surrogate likelihood post-hoc before inferring the posterior to assess the contribution of a feature. We demonstrate the usefulness of our method by identifying the most important features for inferring parameters of an example HH neuron model. Beyond neuroscience, our method is generally applicable to SBI workflows that rely on data features for inference used in other scientific fields.

Efficient identification of informative features in simulation-based inference
Jonas Beck
Michael Deistler
University of Tübingen
University of Tübingen
j jonas beckbuni www.bingen.ndee
mi chael. deint www.rrruuntttttnnnnnn..mm wwwebingeennn.
Yves Bernaerts
Jakob H. Macke
University of Tübingen
University of Tübingen
yves. www.aertsuunii -wweiinttubbinnnnccmm
Max Planck Institute for Intelligent Systems de
jakob
Philipp Berens
University of Tübingen
philipp berens@uni- bereens@un-tuuebingen twwbingen.ddeee de
Abstract
Simulation basod Bayesian inference (SBI) can be used to estimate the parameters complex mechanistic models given observed model outputs without requiring access explicit likelihood evaluations. A prime example for the application of SBI in neuroscience involves estimating the parameters governing the response dynamics of Bodgkin--Huxhae ments, by inferring a posterior over the parameters that is consistent with a set observations. To this end, many SBI methods employ a of summary statisties scientifically interpretable features to estimate surrogate likelihood or posterior. However, currently, there is no way to identify how much each summary statistic or feature contributes to reducing posterior uncertainty. To address this challenge, one could simply compare the posteriors with and without a given feature included in the inference process. However, for large or nested feature sets, this would necessitate repeatedly estimating the posterior, which is computationally expensive or even prohibitive. Here, we provide a more efficient approach based on the SBI method neural likelihood estimation (NLE): We show that one can marginalize the trained surrogate likelihood post-hoc before inferring the posterior assess feature. We demonstrate the usefulness of our method by identifying the most important features for inferring parameters of example HH neuron model. Beyond neuroscience, our method is generally applicable to SBI workflows that rely data features for inference used in other scientific fields.
(HH) models from electrophysiological measure-
the contribution
Introduction
Mechanistic models elegant way to encode scientific knowledge about the world in the form of numerical simulations. They include models such as Kepler's laws of planetary motion (T. the SEIR model (2) for describing the spread of infectious diseases or the Boodgkin-Huxiee (HH) model for the dynamics of action potentials in neurons (3). Efficiently constraining the parameters of such models by measurements a key problem in many disciplines (477). Since these models give rise to intractable likelihood functions, however, classical likelihood Bayesian methods such as variational inference or Markov Chain Monte Carlo (9) cannot be used directly. There are algorithms,
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

MDN
(marginal) likelhood
posterior
prior
simulated data
FSLM
B
peramee
mechanistic model
observation
Poramerr
NLE
2
Figure 1: Feature Selection Through Likelihood Marginalization (FSLM) is method to identify informative features in simulation- inference (SBI). It builds on Neural Likelihood Estimation (NLE) with mixture density networks (MDN). This requires prior mechanistic model to simulate data and an observed data point. 1. Parameters are sampled from the prior and used to simulate synthetic dataset. 2. A MDN learns the probabilistic relationship between data (or data features) and underlying parameters, in the form of a tractable surrogate likelihood. 3. After the MDN has learned approximate the likelihood, it can be conditioned on the observation to yield a likelihood estimate that is parameterised as a mixture of Gaussians. 4. This surrogate likelihood is then combined with the prior distribution to obtain the posterior distribution, i.e. the space of parameters consistent with both prior knowledge and data. Consistent parameters are assigned high, inconsistent parameters a low probability. The naive approach for identifying the importance of features would be to repeat 2. for different sets of data features (grey). With FSLM, the likelihood estimate can be marginalized post-hoc and a single estimate therefore sufficient to obtain multiple posterior distributions (black).
the parameter space,
such as ABC-MCMC [10] or pseudo-marginal methods [1] that deal with this problem, however, they have potentially slow convergence rates and are computationally expensive. To overcome this problem, several techniques collectively known as simulation based inference (SBI) have recently been developed [12]. Leveraging the ability to simulate the model, these techniques obtain estimates of the likelihood posterior from simulated data. The most recent of these algorithms such as neural likelihood estimation (NLE) [33 or neural posterior estimation (NPE) 16) employ state-of-the-as neural density estimators to learn tractable surrogates of these functions to be evaluated instead of the real quantities. To this end, summary statisties or features capturing essential aspects the model's dynamics are defined by the scientist to reduce the high dimensional model output to manageable scale, turn decreasing the problem complexity and computation costs. This often done by hand to emphasize specific aspects of the data or to aid scientific interpretation [12] but can also be automated
23
In neuroscience, many different approaches have been developed in order to find suitable parameters of models of neural activity [4, 24 29]. SBI approaches have been used infer parameters in biophysical neuron models from measurements of neural activity in a Bayesian way [15] 30
3T). For inference of HH models from electrophysiological data, features such as action potential threshold or width and resting membrane potential spike count can be used, reflecting measures that electrophysiologists use to characterize recorded neurons.
For scientific interpretation of SBI results, it is often of interest which features, or combinations of features, have the biggest impact on the posterior and which parameters they affect specifically. To this end, one could compare the posterior uncertainty estimated with and without including a specific feature in the SBI method the increase in uncertainty resulting from not relying on that feature can be used as a measure of its importance and on average is equivalent to its mutual information. Of course, this approach can also be applied to whole subsets of features. To evaluate an entire set of features exhaustively would require re- re-estimating the posterior many times (Fig. [ grey), which would scale prohibitively with the number of features. To address this issue, we here introduce a method called Feature Selection Through Likelihood Marginalization (FSLM) to compute posteriors for arbitrary subsets of features, without the need for repeated training (Fig. T black). To achieve this, we here use NLE 32] and exploit the marignalization properties of mixture density networks

---

which can be used with NLE as a density estimator: instead of re-estimating the surrogate likelihood from seratch, we marginalize it analytically with respect to a given (set of) features before applying Bayes rule to obtain the posterior estimate. This way, we can efficiently compare the posterior uncertainty with and without including a certain feature.
(MDNs) [3].
For a simple linear Gaussian model and non-linear HH models, we show that the obtained posterior estimates are accurate and robust as when repeatedly retraining density estimators from scratch, despite being much faster to compute. We then apply our algorithm to study which features are the most useful for constraining posteriors over parameters of the HH model. Our tool allows us also to study which model parameter is affected by which summary feature. Finally, we suggest a greedy feature selection strategy select useful features from a predefined set.
2 Methods
2.1 Neural Likelihood Estimation (NLE)
Neural Likelihood Estimation (32 is SBI method which approximates the likelihood p(æ: () of the data æ given the model parameters 0 by training a conditional neural density estimator on data generated from simulations of mechanistic model. Using Bayes rule, the approximate likelihood can then be used to obtain an estimate of the posterior (Fig. D.. Unlike NPE [15]16), NLE requires an additional sampling or inference step [4] to obtain the posterior distribution. However, it allows access to the intermediate likelihood approximation. a property we will later exploit to develop our efficient feature selection algorithm (see Sec. 2.3
First, a set of N parameters are sampled from the prior distribution 6, ~ p(0). n (L....,N), NY. With these parameters, the simulator is run to implicitly sample from the model's likelihood function according to Fr ~ p()) Here, xx = (xI..... NN,)) are feature vectors that are usually taken to be a function of the simulator output & with the individual features r; = fi(s), where i E {1. (1.....NJ) }. rather than the output directly. For mechanistic models whose output time series, also reduces the dimensionality of the data to a set of lower dimensional data features. The resulting training data {I,,. On h:N p(x; 0) can then be used to train conditional density estimator (x 0). parameterized by Φ. to approximate the likelihood function. Here, we use a Mixture Density Network for this task (33), the output density of which is parameterised by a mixture of Gaussians, which can be marginalized analytically. Thus, p(x: 0) = Σx TNN ((u,, Σx) with the parameters (14; Σx: л) being non-linear functions of the inputs 0 and the network parameters The parameters & are optimized by maximizing the log-likelihood [(x, 0) * En log ( O) over the training data with respect to d. As the number of training samples goes to infinity, this is equivalent to maximizing the negative Kullback-Leibler (KL) divergence between the true and approximate posterior for every æ E supp(æ)) (3]:
Ep(0..) [log (x 0)] - p()) [Dxz (p())))))) qe (x)))) + const.
After obtaining such a tractable likelihood surrogate, Bayes rule can be used to obtain an estimate of the posterior conditioned observation v, (see Eq. 2 for instance via Markov Chain Monte Carlo (MCMC):
p(0 x.) x (at ()p(O)
2.2 A naive algorithm for quantifying feature importance
Given this posterior estimate, we would now like to answer the following question: for which feature I; from vector of features x does the uncertainty of the posterior estimate p(0 x) increase the most, when it is ignored? For a single feature, a naive and costly algorithm does the following: iterate over 2: to obtain 2 (#1,-. "0;1-...;-11, 11 ENy) NNn) train a total of N+1 density estimators obtain the likelihoods (0) and p(x (0), sample their associated posteriors and compare ((0) to the reference posterior p(0 x) with every feature present. The same procedure can also be applied to quantify the contribution of any arbitrary subset of features. As this procedure requires estimating the likelihood based on the reduced feature set from scratch for each feature (set) (Fig. D grey arrows), is computationally costly. To quantify the contribution of different features more efficiently, we would thus need a way to avoid re- re-estimating the posterior with each feature left out.
3

---

a
p& Ho Kg)
Di xg
Modell peremteere
IONAA (BRp)
FSLM
10
Ground buth NLE-MON
6,
FSIM
Removed feature
Figure 2: FSLM can accurately compute posteriors with one feature marginalized out for linear model described in Eq. a. ID and 2D marginals of the full posterior estimated by NLE and the ground truth, which be analytically computed. b. Posterior distribution when Ty removed from the feature set. FSLM is as accurate as re-estimating NLE on the reduced feature set. c. Increase in the uncertainty of the marginal posterior distributions measured using the IQR ratio (as defined in Sec. (2.4) for each feature, for the analytical posteriors and those computed using FSLM. The joint influence of 2i and 22 02, as per our model, is clearly visible.
2.3 Efficient quantification of feature importance through post-hoc likelihood marginalization
To more efficiently estimate posteriors for any subset of features, we propose to marginalize the likelihood estimate obtained via NLE post-hoc, as opposed training a separate density estimator
for every new set of features to obtain the marginal estimates (Sec. 2.2 Fig. T We refer to this algorithm as efficient Feature Selection through Likelihood Marginalization (FSLM).
Suppose we want to evaluate the joint increase in posterior uncertainty for a subset of features. We partition our feature vector æ - 222.. such that æ2 contains the features be removed from the inference. To avoid the need for training the surrogate likelihood for such a feature subset, we can rewrite the posterior with respect to æi as
p(0 x1) (x) 0)p(0) - =
f g0 g0( 22 0 p(0)
Since we parameterize the approximate likelihood qa as an MDN, we can perform this marginalization analytically: For a K-component MDN, the parameters of each mixture component can be similarly partitioned [35] Chapter 2.3.2]:
Σ
Σ2,
H4 = | Me2l,
Σx
k E
(4)
Σx
Σx
Then the marginalization with respect to æ2 results in the following posterior distribution
p(0 x1) (x 0))(())
N Mx,1 Ex,11) p(0).
(5)
where (0), (()) and Σ = Σ(e). FSLM can thus sample from the posterior given arbitrary feature subsets without estimating the surrogate likelihood do from scratch.
We implement FSLM using python 3.8, building on top of the public sbi library (36). All the code
available at wigght....oooo...... wtthuu.commbberesslaa//a/c running Intel(R) Xeon(R) Gold 6226R CPUs @ 2.90GHz.
repo All computations were done internal cluster
2.4 Measures of posterior uncertainty
Assuming the NLE procedure converges and the density estimator is sufficiently flexible, differences in posterior uncertainty can be ascribed to differences in information content of features or noise the data. Given that that noise in the data remains constant, relative differences in posterior uncertainty can be used to assess the contributions of individual summary features. Hence, we can identify informative summary statistics by quantifying the change in uncertainty of p(0) x.)))))) with respect to p(0) æ). For this two metrics:

---

Model
Method
Training time
Sampling time
Total time
KL
NLE + MDN
3.57 0.79
0.33 + 0.04
3.89 0.80
0.06 0.10
LGM [min]
NLE MAF
11,40 + 2.22
0.57 + 0.28
11.98 + 2.30
0.05 + 0.12
FSLM
0.64 0.21
0.66 + 0.10
1.30 0.26
0.07
0.16
NLE + MDN
84,20 + 9.28
10.34 + 0.54
94.53 + 9.24
29.04 + 8.13
HH [h]
NLE MAF
130.54 = 29.18
19.58 + 4.64
150.12 = 30.86
27.96 10.15
FSLM
7.21 + 1.49
13.27 + 1.04
20.48 = 2.08
28.27 + 7.40
Table 1: Performance of FSLM (mean + standard deviation across 10 runs), compared repeated runs of NLE employing MDN or MAF. For the LGM example, all times are minutes [min] and rejection sampling was used to sample the posterior. For the HH model, all times are hours [h] and MCMC was used. The last column (KL) measures the KL divergence (as defined in Sec. 2.4). For the LGM KL is measured between the posterior and the analytic ground truth. Since no ground truth data is available for the HH model, we use a NLE+MIN estimate of the posterior trained on substantially larger set of summary statistics (see Appendix, Tab. 3 The mean and standard deviation were obtained across all 10 runs and for all subsets that exclude one feature. The timings are single threaded.
Ratio of IQR of the marginals: We use the ratio of inter quantile ranges (IQR) of the DDmmaginnals of p(0 æ1) and p(0 z). This lets us deduce if specific features reduce the uncertainty of some measure of spread as it is much less susceptible to
parameters 8 more than others. We use IQRs outliers compared to the standard deviation. The resulting visualization allows to precisely pin down which parameter is affected by which summary feature.
Mon-parametric estimate of KL. divergence: As a second metric, we use a non-parametric estimate of the KL divergence between p(0 21) and p x). Since posterior estimate, on average, cannot become more constrained if features are removed, increases in the KL estimate indicate an increase in uncertainty overall. To estimate the KL-divergence between p(O xxi) and
p(0 æ), we employ a purely sample based estimator (37 38 This is necessary because NLE only allows to evaluate unnormalized posteriors which circumvents having evaluate them explicitly. The estimator makes use of 1-nearest neighbor search and is implemented using k-d trees [39] to reduce the 2N comparison operations to a manageable number. With samples X
{X+] 1. and
Y = (Y)),M == from p(0) | æ and p(0 æ)) respectively, the KL-divergence (L(P( x p) x))
is estimated as
min, X Y| min X X jf
M
DKL
log
+ log
(6)
N --1'
where d the dimensionality of the samples, and || denotes the fn-norm. We note that the Maxi-
mum Mean Discrepancy (MMD) 40 41 could be used an alternative measure, but preliminary
experiments indicated that the KI. estimate was more consistent.
3 Results
3.1
Example 1: Linear Gaussian model
We illustrate our algorithm in simple toy problem, where by construction, we know which features should have which effect. We used a Linear Gaussian Model (LGM) with a three dimensional [ T1, F2, 23]. 0 is
parameter space e (o. oz 02] and four dimensional feature space
drawn from a uniform prior 0; ~ 2d-5, 5). i E {0, 2), with the features æ being noisy linear transformations of the parameter vector O. The mean (0)) and covariance are defined by
1 0
0
1
0
m()) p4 + LO, == L
0
1 1
0_
With scale and shift Mg- We construct L such that we can easily verify that the posterior reacts to and interacts with marginalized features as expected. our example, we let feature 2i solely depend on 8g 2i on both 0, and 0g x2 on 0g. and 23 without influence any model parameter.
5

---

a
NLE MAF
NLE MDN
FSLM
20
96e
Model pormmtor
12960) 13R
APA
ANO
AIP
APC
AFT
APC
JPN
AP7
APM
Ant
Removed feeture
b
Uee-t0
Uanatt
100 mA
10 ms
120
APA of oostentt SXTD pAmTpees (M1))
115
110
NOT of postentes
105
0.0001
0.8
0.0001
0.8
Bunl
-0- 9ur -980.)
(Ry--19K))
Poull(()
-00
APT
APA
Figure 3: Removal individual features affects posterior uncertainty and samples in a HH model: a. The ratio of IQRs for the marginal posterior distributions (as defined in Sec. 2.4). when individual features are left out. Values are large (dark) when the feature has a large effect on the marginal posterior for the specified model parameter, and small otherwise (light). b. Simulations from the HH model with the original true parameters (0., blue), parameters drawn from the full posterior (0 full. orange) and two posterior distributions when leaving the features APT and APA (0 (APT green and red, 0, APA. respectively). e. Posterior distribution with the full feature set (orange) and after removal of the highly constraining feature APT (green) shows the effect some of the 1D and 2D marginals (see Fig/6 for complete posteriors) compared to the full posterior (orange). d. As c. for the less informative feature APW. The effect of the removal of informative uninformative features from the posterior estimate on the feature values of the simulated parameter samples.

---

The resulting posterior is well concentrated on 8, and the NLE+MDN estimate matches the ana- lytically computed posterior well (Fig. 21. To show that FSLM can correctly identify individual feature importance, while being more efficient than re-running NLE for each subset, we compared samples from the FSI. -Serived posterior estimates for every three feature subset to those obtained using the analytical marginal likelihood p(z 0), E {0, 1, 2, 3} and brute-force NLE estimates using separate ( ()) for each subset (Fig. 2b, Appendix Fig 5 a-c). The brute-force baseline was computed with both a MDN and a Masked
TaaskeeAAttoggeessss Autoregressive Flow (MAF) [42]. The MDNs had three hidden layers and 10 mixture components. The MAF consisted of five 3 layer deep MADES 43). Training was done on 10,000 samples, then 500 samples were drawn from each posterior, for a total of 10 different initializations. Training was terminated when the log- fog-likelihoods had stopped improving on validation set that was split from the training set at 10% over 20 epochs.
The posterior estimates obtained from FSLM matched those from the analytic ground truth very accurately (Fig. 2 Appendix Fig. 5 and Tab. DT while being several times faster to compute than both NLE alternatives (see Tab. T. As expected from the construction of our model, we observed that 5 2)0o) was not only less constrained compared p(0 2) (Fig. 22). but that p(0n 22.0) = p(0)) 5). We also found that FSLM correctly identified the dependency structure between the other parameters and features (Appendix Fig. ()-)), with FSLM matching the analytic ground truth while being much faster than both NLE approaches (Tab.
3.2 Example 2: Hodgkin-Huxley neuron model
Next, we turned to a more realistic complex mechanistic model used for simulations of neurons. For
experiments used single compartment Hodgkin-Huxkee model (3), Eq. 8 derived from Pospischil et al. [44]. In addition to the standard sodium, potassium and leak currents to generate spiking this model also considers slow non-inacie non-inactivating K current to model spike-frequency adaptation and a high-threshold Ca current to generate bursting.
dV.
I,- = C. Asmaa dt
+ m N(V, - ENa) + OV - Ek)
(8)
+ [Ve - Fieak) DMP(V, - Ex) + 9u F.ca)
Eq. 8 models the evolution of the membrane voltage V-2 V(t) given stimulus I(t) I. Di. i E {Na, K,1, M, L} are the maximum conductances of the sodium, potassium, leak, adaptive potassium and calcium ion channels, respectively, and E, are the associated reversal potentials. I denotes the current per unit area, C' the membrane capacitance, Asoma the compartment area and n, m, h, g, r and p represent the fraction of independent gates in the open state, based Hodgkin and Huxley (3] For the remaining equations and parameter settings, see Appendix A.3
We then compared the performance of FSLM and the two versions of brute-force NLE to measure how commonly used summary features constrain the posterior of an exemplary voltage trace. Specifically, used set of 10 summary statistics commonly of interest to neurophysiologists, which capture the response dynamics of the HH model and can be used to differentiate cortical cell types. They briefly explained in Table 3 To extract these features from the simulated traces use the Allen SDK
To train the density estimator we simulated I million observations with parameters drawn from Uniform prior over biologically plausible parameter ranges A.3). However, because substantial fraction of the generated observations did not produce spiking behavior, some features such spike rate were not defined. Hence, these observations could not be used to train the density estimator. We employed two different techniques introduced by Lueckmann et al. [15) to reduce the amount of invalid samples returned by the prior and to compensate for any remaining invalid observations in the training data (for details see Appendix A.2). This yielded accurate posterior estimates, despite incomplete data. The change in uncertainty was then assessed on 3000 samples from the MDN estimates of p(0) æ)) and p(0 x..)). Architectures and training were kept the same as in Sec 3.1
We found that all three algorithms identified the same two features as yielding the highest increase in posterior uncertainty (Fig. Ba): the action potential threshold (APT) and the action potential count (APC), while a MAF additionally identified the mean resting potential (μ(Vrest)) as a constraining feature. For the HDN-bascd density estimates, APT strongly affected the threshold voltage Vr, the
7

---

10
15-
29RE
33% 1016
Fosture
APA
50000000 05
APA
APT
A22
APC
FAAMAA
510 m
Seaturea
Run
Figure 4: Greedy selection of the most informative features from candidate set for the HH model: Example trajectory a single run using FSLM for greedy search (black). Parameters are added to the feature set from left right (right most point includes all features). Medians and confidence intervals for 50 runs, when parameters are selected at random (blue) vs. when selected by the greedy-search procedure (orange). b. Order in which the features chosen for the first 10 feature searches.
leak conductance Gleak as well as the sodium conductance INa. APC, influenced the leak reversal potential Eteak. the leak conductance gleak well as C and @Na (Fig. 31.. When flow was used, A did not seem to be constrained by either APT or APC and (Vrrst)) additionally constrained the leak reversal potential (Eteak), as well as the adaptive Potassium current 9Ms. In contrast, the action potential amplitude (APA), as well as other features, did not have a sizable impact on the posterior in this example (Fig. 3a. While the results for both MDN- algorithms
equivalent and have strong overlap with the flow-based predictions, FSLM was significantly faster than both brute-force alternatives (= 5-fold and 7-fold speedups respectively, Table I and overall almost accurate, when measured against a common baseline [
As consequence, simulations based on samples from p(0 x)) and p(0 (AAAAA) were comparable to the observed voltage trace Vo((), while the sample drawn from p(0 a APT showed different spiking behavior, with much longer inter spike intervals in the simulation (Fig. 3D.. This resulted from higher values for the sodium conductance INa and threshold voltage Vr, compared to 6,, for the particular simulation shown (Fig. Ba,c). The measured increase uncertainty is also directly visible in the two-dimensional marginal distributions of the posterior, where the removal of APT lead to a much less constrained posterior estimate (Fig Br), while the removal of APA did not affect the posterior uncertainty much 31). Also, the measured summary statistics from samples from p(0 x (PP) and p(0 x APA) showed the expected behavior: the APT feature was much more variable in samples from ş(0 APT))) than p(O a:), while the estimated value for APA did not vary as drastically even for samples from $(0) (AAAAA) (Fig. 3e).
3.3 Efficient greedy feature selection using FSLM
In the previous section, we showed that FSLM can be used to study how the removal of individual features increases the uncertainty of the posterior. Next show how FSLM can be used to study how the addition of features decreases uncertainty and provide a method to efficiently search for informative features by greedy selection of new features.
First, we used NLE to obtain a posterior estimate given a candidate set of features. We then used FSLM to obtain all the ID marginal likelihoods, sampled their associated posteriors p(O 2x)) and computed their KL divergences to the full posterior. We subsequently selected the feature 2; that minimizes the KL, before repeated the procedure for all 2D marginals 5(0 (z,, (z), which include the selected feature z;. We can iterate this procedure until a desired number of features is identified. This is equivalent combinations of features the scheme is easily extensible to larger beam sizes.
seam-search (45) with beam size equal to one. To investigate
To demonstrate this, we turned again to the HH example already studied in Sec. 3.2 million samples for training and sampled every posterior 1000 times. To check for consistency, we repeated the procedure 50 times. Compared to selecting features at random, selecting features based
We used 1

---

on the results of our greedy algorithm consistently yielded more constrained posterior estimates (see Fig. Fa). Furthermore, we found the three most informative features to be (V APT, AHP. respectively, which were consistently chosen in this order (see Fig. 4D).
4 Discussion
We set out develop an efficient method characterize the relative contribution of different features to constraining the posterior in simulation-based inference. To this end, we used NLE with a mixture density network to estimate the surrogate likelihood and marginalize it analytically. We showed that our method, FSLM, correctly identifies influential features both in simple linear Gaussian model, as well as in a complex mechanistic model, namely the HH model in neuroscience. In both cases FSLM was as accurate as repeated application of NLE. but many times faster, independent of whether the likelihood was approximated by MDNs or MAFs. While the different architectures lead to subtle differences in their predictions of important features of HH models, we found them still to be mostly in agreement. FSLM can thus provide guidance on how to specify, prepare and select appropriate summary statistics, while also allowing to study which parameters are constrained by which features to gain insight into the scientific structure of the model. In contrast to methods developed for automatically deriving optimal (combinations of) summary features 17 21 23), FSLM focuses on making the evaluation and selection of individual scientifically relevant Teature sets interpretable. Our visualization method allowed us to identify which features influenced the posterior uncertainty over which model parameters, and for which model parameters the posterior was not influenced by any feature. A similar approach could also be implemented for Approximate Bayesian
Computation (ABC) 46 47
Compared to approaches that aim to directly quantify the mutual information (MI) between sets of parameters and features (23), our method aims to quantify the relationship between model parameters and features for individual observations. Computing estimates of the MI for simulator models is difficult without access to closed form likelihood and purely sample based estimators [48] 49 do
not allow for inspection of individual or small sets of observations. However, this is of interest when different cells, cell types or cell families are constrained by different feature sets. Nonetheless, FSLM still retains the option to obtain single ranking of features actoss a whole dataset by averaging over all observations in the set, which would then be equivalent to the MI between posterior and prior (the trained neural network is amortized and, thus, requires no retraining for different observations).
For the case of Hodgkin Huxkey model, frequently used in neuroscience as a mechanistic model of neural activity, we showed that action potential threshold and action potential count were the most informative features, jointly constraining the action potential threshold, the resting potential, membrane capacitance as well as model conductances (except for that of potassium). Removing these features from the inference procedure led to much more poorly constrained posteriors and simulations that did not match the original trace.
We note some limitations of FSLM and the metrics introduced: First, the inherent dependence on MDNs for the conditional neural density estimation limits overall accuracy for highly complex posteriors or likelihoods, compared to e.g. normalizing flows, as it only works well for likelihoods that are well approximated by mixtures of Gaussians. Additionally, the use of full covariance matrices can be problematic in nigher-dimensiona. feature spaces due the curse of dimensionality and one would have to evaluate whether the use of e.g. diagonal covariances is more suitable. Second, since our metrics are relative to the overall level of uncertainty in a given posterior marginal, for very sharp distributions small increase in variance can already lead to large increase in our metrics, which would not be the case for a broad distribution. While this may be the desired behavior in some cases, it might be misleading in others. Third, the effect of feature may be judged differently depending on whether the effect of removing it from a large set of features is being assessed, or whether is being added early in a greedy feature selection run. Finally, the approximation quality of NLE, thus also FSLM, depends crucially on the choice of prior (50). However, posterior predictive checks [51
[2] and coverage tests [53) applicable to NLE can also be applied to FLSM to verify the correctness of the inference results.
Our method opens the door for careful studies on the relationship between data features and model distributions. neuroscience, for example, we could explore whether commonly used features such as those used here constrain the posterior equally across different cellular families types for

---

example, is action potential threshold equally important for modeling fast-spiking interneurons and for cellular families with more complex firing behavior, or are additional features needed in this case? Given the efficient exploration of the relationship between features and posterior uncertainty, we believe that our method will also be applicable to wide range of inference problems beyond neuroscience.
Acknowledgments
We thank the German Research Foundation for support through the Excellence Cluster "Machine Learning New Perspectives for Science" (EXC 2064, 390727645) and Heisenberg Professorship to PB (BE5601/8-1) 1). JHM and PB are also members of the Tübingen AI Center funded by the We thank the International Max
German Ministry for Science and Education (FKZ Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Michael Deistler. We also thank our colleagues, in particular Jonathan Oesterle and Lisa Koch, for their valuable comments and fruitful discussions.
10

---

References
[1] J. L. Russell. "Kepler's Laws of Planetary Motion: 1609-1666". The British Journal for the History of Science (1964), pp. 1-24.
[2] Herbert W. Hethcote. "The Mathematics of Infectious Diseases".. .IAMM Review (2000), pp. 599
653.
[3] A. L. Hodgkin and A. F. Huxley. "A quantitative description of membrane current and its
application to conduction and excitation nerve". The Journal of Physiology (1952), pp. 500- 544.
[4] Astrid A. Prinz, Cyrus P. Billimoria, and Eve Marder. "Alternative to Hand-Tuning Conductance- Based Models: Construction and Analysis of Databases of Model Neurons". Journal of Neurophysiology (2003), pp. 3998-4011
[5] Mark S. Goldman et al. "Global Structure, Robustness, and Modulation Neuronal Models". The Journal of Neuroscience (2001), pp. 5229-5238.
[6] Johann Brehmer and Kyle Cranmer. "Flows for simultaneous manifold learning and density estimation".. Advances in neural information processing systems 33: Annual conference on neural information processing systems, NeurIPS 2020. 2020.
[7] Sergey Kravtsov. "Dynamies and Predictability of Hemispheric-Scalee Multidecadal Climate Variability in an Observationally Constrained Mechanistic Model". Journal of Climate (2020),
pp. 4599 4620.
[8] Martin Wainwright and Michael Jordan. "Graphical Models, Exponential Families, and Variational Inference". Foundations and Trends@ in Machine Learning (2007), pp. 1-305.
[9] Radford M Neal. "Probabilistic Inference Using Markov Chain Monte Carlo Methods" (1993),
p. 144.
[10] Paul Marjoram et al. "Markov chain Monte Carlo without likelihoods". Proceedings of the National Academy of Sciences (2003), pp. 15324-15328.
[11] Christophe Andrieu and Gareth O. Roberts. "The Monte Carlo computations". The Annals of Statistics (2009), pp. 697-725.
approach for efficient
[12] Kyle Cranmer, Johann Brehmer, and Gilles Louppe. "The frontier of simulation-based in- ference". Proceedings of the National Academy of Sciences of the United States of America (2020), pp. 30055-30063
[13] George Papamakarios. "Neural Density Estimation and Likelihood-free Inference". arXiv:1910.1323: (2019).
[14] George Papamakarios and lain Murray. "Fast e-free inference of simulation models with bayesian conditional density estimation". Advances in neural information processing systems 29: Annual conference neural information processing systems. 2016, pp. 1028-1036.
[15] Jan-Matthis Lueckmann et "Flexible statistical inference for mechanistic models of neural dynamics" Advances in neural information processing systems 30: Annual conference on neural information processing systems. 2017, pp. 1289-1299.
[16] David S. Greenberg, Marcel Nonnenmacher, and Jakob H. Macke. "Automatic posterior transformation for likelihood-free inference". Proceedings of the 36th international conference on machine learning, ICML (CML.2019 2019, pp. 2404-2414.
[17] Tian Qi Chen et al. "Neural ordinary differential equations". Advances in neural information processing systems 31: Annual conference on neural information processing systems, NeurlPS 2018. 2018, pp. 6572-6583.
[18] Michael GB Blum et al. "A comparative review of dimension reduction methods in approximate Bayesian computation". Statistical Science (2013), pp. 189-208.
[19] Bai Jiang et al. "Learring Summary Statistic for Approximate Bayesian Computation via Deep Neural Network". Statistica Sinica (2018).
[20] Justin Alsing, Benjamin Wandelt, and Stephen Feeney. "Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology". MMonthly Notices of the Royal Astronomicul Society (2018), pp. 2874-2885.
[21] Rafael Izbicki, Ann B. Lee, and Taylor Pospisil. "ABC-CDE: Towards Approximate Bayesian Computation with Complex High-Dimensional Data and Limited Simulations". Journal of Computational and Graphical Statistics (2019), pp. 481-492.

---

[22]
Justin Alsing and Benjamin Wandelt. "Nuisance hardened data compression for fast likelihood- free inference". Monthly Notices of the Royal Astronomical Society (2019), pp. 5093-5103. Yanzhi Chen et al. "Neural Approximate Sufficient Statistics for Implicit Models". Interna- tional Conference Learning Representations. 2021.
[23]
[24] Astrid A. Prinz, Dirk Bucher, and Eve Marder. "Similar network activity from disparate circuit parameters". Nature Neuroscience (2004), pp. 1345-1352.
[25] Hillel Ori, Eve Marder, and Shimon Marom. "Cellular function given parametric variation in the Hodgkin and Huxley model of excitability". Proceedings of the National Academy of Sciences (2018), E821 1-8821...
[26] Shaul Druckmann et al. "A novel multiple objective optimization framework for constraining conductance-basee neuron models by experimental data". Frontiers in Neuroscience (2007),
p.
[27] Roy Ben-Shalom et "Optimizing ion channel models using a parallel genetic algorithm on graphical processors". Journal of Neuroscience Methods (2012), pp. 183-194.
[28] Pablo Achard and Erik De Schutter. "Complex parameter landscape for complex neuron model". PLoS computational biology (2006), e94.
[29] A. L. Taylor, J.-M. Goaillard, and E. Marder. "How Multiple Conductances Determine Elec- trophysiological Properties in a Multicompartment Model". Journal of Neuroscience (2009), pp. 5573-5586.
[30] Pedro J Gonçalves et al. "Training deep neural density estimators to identify mechanistic models of neural dynamics". eLife (2020), e56261.
[31] Jonathan Oesterle et al. "Bayesian inference for biophysical neuron models enables stimulus optimization for retinal neuroprosthetics" eLife (2020), e54997.
[32] George Papamakarios, David C. Sterratt, and lain Murray. "Sequential neural likelihood: Fast likelihood- free inference with autoregressive flows". The 22nd international conference on artificial intelligence and statistics, AISTATS 2019. 2019, pp. 837-848.
[33] Christopher M. Bishop. "Mixture density networks". Technical Report (1994).
[34] Manuel Glöckler, Michael Deistler, and Jakob H. Macke. "Variational methods for simulation-
based inference". International Conference on Learning Representations. 2021.
[35] Christopher Bishop. Pattern Recognition and Machine Learning. New York, 2006. ISBN:
978-0-387-31073-2.
[36] Alvaro Tejero-Cantero et al. "sbi: A toolkit for simulation-based Source Software (2020), p. 2505.
Journall of Open
[37] Bai Jiang. "Approximate bayesian computation with kullback-leibee divergence as data discrepancy". International conference on artificial intelligence and statistics, AISTATS 2018.
2018, pp. 1711-1721
[38] Fernando Perez-Cruz. "Kullback-Leibe divergence estimation of continuous distributions". 2008 IEEE International Symposium on Information Theory. Toronto, ON, Canada, 2008, pp. 1666-1670. ISBN: 978-1-4244-2256-2.
[39] Jon Louis Bentley. "Multidimensional binary search trees used for associative searching". Communications of the ACM (1975), pp. 509-517.
[40] Arthur Gretton et al. "A kernel two-sample test". The Journal of Machine Learning Research
(2012), pp. 723-773.
[41] James Lloyd and Zoubin Ghahramani. "Statistical Model Criticism using Kernel Two Sample Tests". Advances in Neural Information Processing Systems. 2015.
[42] George Papamakarios, lain Murray, and Theo Pavlakou. "Masked autoregressive flow for density estimation". Advances in neural information processing systems 30: Annual conference on neural information processing systems 2017. 2017, pp. 2338-2347.
[43] Mathieu Germain et al. "MADE: Masked autoencoder for distribution estimation". "stimation..Procedding of the 32nd international conference machine learning, ICML 2015. 2015, pp. 881.
[44] Martin Pospischil et al. "Minimal Hodgkin type models for different classes of cortical and thalamic neurons"'. Biological Cybernetics (2008), pp. 427-441.
[45] David Furcy and Sven Koenig. "Limited discrepancy beam search".. Proceedings of the 19th international joint conference on Artificial intelligence. 2005, pp. 125-131.
[46]
Handbook of Approximate Bayesian Computation. New York, 2018. ISBN: 978- 978-1-315-11719.5
12

---

[47] M. A. Beaumont et al. "Adaptive approximate Bayesian computation". Biometrika (2009).
pp. 983-990.
[48] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger. "Estimating mutual information".
Physical Review E (2004), p. 066138.
[49] Nicholas Carrara and Jesse Ernst. "On the Estimation of Mutual Information". Proceedings
(2020), p. 31.
[50] Joeri Hermans et al. "Averting A Crisis In Simulation-Based Inference". arXiv:2110.06581 (2021).
[51] Jonah Gabry et al. "Visualization in Bayesian workflow". Journal of the Royal Statistical Society: Series A (Statistics in Society) (2019), pp. 389-402.
[52] Sean Talts al. "Validating Bayesian Inference Algorithms with Simulation Based Calibra- tion". arXiv:1804.067: (2020).
[53] David Zhao al. "Diagnostics for conditional density models and Bayesian inference algo- rithms". R Proceedings of the Thirty-Seventh Conference on Uncertainty Artificial Intelligence 2021, pp. 1830-1844 1840.
[54] John Charles Butcher and Nicolette Goodwin. Numerical methods for ordinary differential equations. 2016.
13

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Sec.
(c) Did you discuss any potential negative societal impacts of your work? [N/A] We do not see negative social impact of our work.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A]
(b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main ex- perimental results (either in the supplemental material or URL)? [Yes] See Sec.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Sec. 3.1 and 3.2
(c) Did you report error bars (e-g., with respect the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Sec. 2.3
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you eite the creators? [Yes]
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material or as a URL? [Yes]
(d) Did you discuss whether and how consent was obtained from people whose data you're 're using/curating? [Yes] No data from third parties used.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [Yes]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "ZPTAVPrDxAA", "reviewer": "Reviewer_9rkA", "summary": "The paper proposed a method to determine the \"informativeness\" of data features by comparing the width of posteriors with some data features marginalized out. The authors compare a naive solution where many models are trained to their proposed technique which takes advantage of the properties of Mixture Density Networks (DMN) as the architecture for Neural Likelihood Estimation (also known as synthetic likelihood) to trivially marginalize over data features without retraining.\n\nTheir method Feature Selection Through Likelihood Marginalization (FSLM):\n- allows for closed-form marginalization of data features, thereby enabling access to posteriors which are conditioned on arbitrarily marginalized data random variables without retraining.\n- comes with a greedy feature selection strategy to select useful features from a predefined set.\n- is compared to ground truth posteriors from simple linear Gaussian model and against a non-linear neuronal model. ", "questions": "- Why is this method superior to a method which estimates the mutual information between parameters and data with some data features marginalized out?\n", "limitations": "- I think this was also made clear.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Originality\n\n### Strengths\n- The closed-form marginalization is simple and easily seen as a result of using a MDN. I've not seen it used before and I find it interesting for their purposes.\n\n### Weaknesses\n- The paper does not reference the work of Benjamin Wandelt and coauthors which aims to transform simulation data into a compressed format which maximizes the information shared by the data and its compressed representation for Neural Likelihood Estimation. This is not strictly the same subject as Wandelt's work does not quantify the value of specific dimensions of input data; however, it is a highly relevant dimension reduction technique. Please consider citing all of the papers on the [PyDelfi](https://pydelfi.readthedocs.io/en/latest/intro.html) website (there are three).\n- There is other tangentially related work regarding marginalizing over parameters. In particular, [Truncated Marginal Neural Ratio Estimation](https://arxiv.org/abs/2107.01214), [Arbitrary Marginal Neural Ratio Estimation for Simulation-based Inference](https://arxiv.org/abs/2110.00449), and (from pydelfi) [Nuisance hardened data compression for fast likelihood-free inference](https://arxiv.org/abs/1903.01473). Perhaps it would be interesting to include for the purposes of dimensionality reduction.\n\n## Quality\n### Strengths\n- The Gaussian linear model example is chosen well to show the relevance of independence between parameters and data and the corresponding effect on the estimated posterior.\n- The posterior visualization, tables which indicate computational cost and estimated KL divergence, and search result visualizations are all quite clear.\n- The limitations of using MDNs are made clear and the benefit for feature search is as well.\n\n### Weaknesses\n- Although the method seems nice, I still don't understand why the practitioner would not instead want to compute the mutual information. We do not know what sort of observations we will have in the future and the mutual information is likely more robust to future data (rather than overfitting feature selection to x_o.)\n\n## Clarity\n- Clarity was satisfactory in my opinion.\n\n\n## Significance\n### Strengths\n- They're certainly a nice feature of MDNs and NLE.\n\n### Weaknesses\n- MDNs seem to be of limited value in the times of normalizing flows, etc.\n- I also don't see why setting the feature importance based off a single x_o is really so interesting. Selected features should have some sort of predictive power for future data, right? If so, shouldn't we then integrate over all possible data, i.e. estimate the mutual information between different feature sets and the parameters?\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "QLjQkDGsl16", "reviewer": "Reviewer_Pf14", "summary": "The authors propose a new method for determining informative features in simulation based inference. Recently proposed simulation based inference methods proceed by (1) generating parameter samples from the prior, (2) simulating to generate conditional samples of the data given the parameter, (3) fitting an approximation to the likelihood of data features under the simulator model consisting of a mixture model parameterized by neural networks and (4) sampling from the posterior over parameters given the data, using the conditional mixture model as an approximation to the true likelihood. The authors observe that to obtain a posterior conditional on a subset of features, one can use the property of multivariate normals that their marginals are analytically tractable, and replace the full approximate likelihood with the marginal approximate likelihood.  ", "questions": "There are some serious typos and lack of clarity in Eqn. 6 - there's min_j, but j does not appear; lowercase m and n are undefined; and it's unclear what min^N_j means (min raised to the Nth power??).\n\nMore generally I'm confused by the motivation for the KL estimator. KL estimation from paired samples is notoriously difficult and unreliable in high dimensions. I think it would make more sense to use maximum mean discrepancy or related techniques, which are well-established for moderately high dimensional data, and  which offer the further advantage of providing a witness function that can help determine where the two distributions differ. (see https://jmlr.csail.mit.edu/papers/v13/gretton12a.html and https://papers.nips.cc/paper/2015/hash/0fcbc61acd0479dc77e3cccc0f5ffca7-Abstract.html)\n\n\n\n", "limitations": "The authors discuss some limitations. One additional limitation that I think is important is that the success of NLE in general depends crucially on the choice of prior; if the true parameters are far out in the tails of the prior, the likelihood approximation is likely to be very poor, particularly as it restricts the tails of the likelihood to be Gaussian (and will thus fail if the true likelihood has, for instance, very fat tails). I worry that the proposed method does not take into account uncertainty or systematic errors in the approximation quality", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The method itself is very straightforward, enough so that it seems like a trick any researcher in the field would be able to figure out, as it relies on well-known properties of multivariate normal distributions. However, the method is presented very clearly. Identifying informative features in simulation based inference is an important and poorly studied problem (to my knowledge); this paper offers smart and clear problem statement and setup. Moreover, writing a paper about the method helps motivate the particular advantages of neural likelihood estimation in contrast to competing techniques (such as NPE). \n\nThe evaluations are clear and compelling, especially the biological application. Hodgkin Huxley models are notoriously hard to fit, and the connection between the features that electrophysiologists focus on and actual underlying ion channel mechanisms are notoriously difficult to understand. This paper thus offers a compelling advance for the computational neuroscience/electrophysiology field in particular.\n\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "IvAHl64OJE", "reviewer": "Reviewer_GQHb", "summary": "This paper develops an efficient feature selection method in simulation-based inference, where each evaluation of the likelihood is costly. Conventionally, to test the sensitivities of the posterior distribution to specific features, one would need to re-evaluate the posterior with different combinations of features many times. The key idea of the method is to estimate the likelihood as an intermediate quantity instead of directly estimating the posterior, and analytically marginalize the joint likelihood for each feature instead of re-evaluating from scratch each time. The marginalization is made possible by the use of mixed density networks with gaussian properties where different parameters are easily partitioned. Using a Hodgkin-Huxley model as an example, the paper demonstrates how this method can be used for selecting features at a much smaller computational cost.", "questions": "i) It would be nice to have a better \"chance level\" curve for the random selection in Fig 4a, i.e., average over multiple random sequences of features.\n\nii) The speedup by the proposed method is evident, but 20 hours is still a long time if one wants to repeat. Can the authors comment on how/where the algorithm can use parallelization to reduce computation time?\n\niii) I had two questions, and they both boil down to understanding the computation time cost of considering more features:\n\n- There is still a human curated step required in the procedure, namely the pre-selected feature list. Can the authors comment on how this new level of efficiency could change the constraints related to the choice of the pre-selected feature list? That is, can one now run through a very long list of candidate features, including all the silly things that one would not expect to be important a priori?\n- Regarding the HH model experiment: I agree that the greedy features selection is a completely reasonable thing to do. But looking at Fig 4a made me wonder if This method also provides easy ways to explore combinations of features. (Specifically, from Fig 4a, I was wondering if the combination of APT, AHP, and APW were important, based on the drop in the blue curve at the addition of APT, and the drop in the black curve at the addition of APW.) How much additional cost would it be to investigate combinations of features?", "limitations": "The authors have adequately addressed the limitations of their methods.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This work seems to be a small but useful technical improvement to simulation-based inference (SBI). I am not an expert of the SBI method space in particular, so I cannot speak much on the originality or the novelty of the idea. But I have seen usages of SBI and surrogate likelihood/posterior approaches in neuroscience as well as in other domains (computational chemistry), and I have not seen an analytical marginalization approach like this before. Whenever the mixed gaussian model (with the convenient partitioning properties that allow marginalization) is a good fit to the system, this method will provide an immediate practical advantage. The paper is mostly clear and well written.\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "-MmtLhBBrLw", "reviewer": "Reviewer_STJL", "summary": "The paper introduces a method to efficiently estimate the importance of several features / outcome variables of a simulation model in a Bayesian setting. The method is based on mixture density networks and the fact that marginalization of Gaussian distributions is analytically tractable and relatively simple.\nThe authors use their method to demonstrate in 2 practical examples that the method reliable identifies feature importance and demonstrate the method's superior speed. ", "questions": "- What kind of density estimator was used for NLE? Just an MDN? Why not a flow? Have alternatives different from the default be tried? This should be mentioned. (Could any of those choices reasonable lead to NLE working better than your method, e.g. a flow instead of an MDN.)\n- How is time for training compared? Is there a predefined number of iterations or are you using a stopping criterion? Or is this just to say that you train one MDN and marginalize instead of retraining an MDN on the subset? (I.e. if you used a flow for NLE it's not clear how to compare time.)\n- Are you modelling the MDN with a full covariance matrix? If so, what would you suggest for higher dimensional data spaces?\n", "limitations": "The authors discuss some limitations in the discussion section. While there are some elements on the limitations of MDNs, I'd like to know more about the covariance matrix that has to be estimated and hence the scalability. (A reasonable argument to be made would be that scalability isn't that important since handcrafted summary statistics are usually of limited dimensionality, but I think this should be discussed.)", "rating": 5, "confidence": 5, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "__Strengths__\n\n- Reasoning about the importance of _interpretable_ summary statistics is commonly (and unjustly) ignored in the era of neural networks. The authors therefore tackle an issue that is highly relevant for reliable and trustworthy inference. \n- Marginalisation using the MDN is indeed an efficient approach and using the so computed marginal distributions to compute the feature important makes sense to me.\n- The method is easy to understand and implement\n- The writing is coherent and what is said is easy to follow (but the text does lack crucial information, see below)\n\n__Weaknesses__\n\n- The paper seems rushed: there is very little information on many relevant details. For example, what kind of density estimator was used for NLE? A flow? If so, what kind of flow? How many layers? I have many questions, see below. If any of those are answered in the text and I missed them, could the authors please let me know?\n- The paper contains only 2 examples, one of which is just a simple Gaussian. I understand the value of the Gaussian example, but the experiments section then only has one \"real\" model, which puts this section on the lighter side. (Perhaps the relatively long description of the HH model could have gone to the appendix in favour of another example.)\n- The literature review is a bit biased towards more recent developments and perhaps a bit unfair towards earlier contributions. For example:\n\n> likelihood-based Bayesian methods such as variational inference [8] or Markov Chain Monte Carlo [9] cannot be used.\n\nThat's a bit imprecise. Of course, MCMC can be used, in the form of ABC-MCMC [1] or pseudo-marginal methods [2]. The argument against those methods is more the computational cost and potentially slow convergence, not however, that it's impossible.\n\n- Equation 11 in the appendix is confusing. The RHS has a variable $valid$ which the LHS does not have. The reference for proofs is to imprecise, I could not find the relevant part. The reference should point to the theorem that establishes the result. \n\n[1] Marjoram, P., Molitor, J., Plagnol, V., & Tavaré, S. (2003). Markov chain Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences, 100(26), 15324-15328.\n\n[2] Andrieu, C., & Roberts, G. O. (2009). The pseudo-marginal approach for efficient Monte Carlo computations. The Annals of Statistics, 37(2), 697-725.\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept