Nuisance hardened data compression for fast likelihood-free inference

TL;DR

This paper introduces a method to perform likelihood-free inference that efficiently marginalizes over nuisance parameters by re-casting them as latent variables and using optimal data compression, reducing the number of simulations needed.

Contribution

The authors develop a novel approach to marginalize nuisance parameters directly within likelihood-free inference, using data compression that preserves Fisher information and is insensitive to nuisances.

Findings

Efficient inference of cosmological parameters with fewer simulations.

Nuisance parameters can be marginalized without increasing simulation cost.

Method validated on supernovae and weak lensing cosmology problems.

Abstract

In this paper we show how nuisance parameter marginalized posteriors can be inferred directly from simulations in a likelihood-free setting, without having to jointly infer the higher-dimensional interesting and nuisance parameter posterior first and marginalize a posteriori. The result is that for an inference task with a given number of interesting parameters, the number of simulations required to perform likelihood-free inference can be kept (roughly) the same irrespective of the number of additional nuisances to be marginalized over. To achieve this we introduce two extensions to the standard likelihood-free inference set-up. Firstly we show how nuisance parameters can be re-cast as latent variables and hence automatically marginalized over in the likelihood-free framework. Secondly, we derive an asymptotically optimal compression from $N$ data down to $n$ summaries -- one per interesting parameter -- such that the Fisher information is (asymptotically) preserved, but the summaries are insensitive (to leading order) to the nuisance parameters. This means that the nuisance marginalized inference task involves learning $n$ interesting parameters from $n$ "nuisance hardened" data summaries, regardless of the presence or number of additional nuisance parameters to be marginalized over. We validate our approach on two examples from cosmology: supernovae and weak lensing data analyses with nuisance parameterized systematics. For the supernova problem, high-fidelity posterior inference of $\Omega_m$ and $w_0$ (marginalized over systematics) can be obtained from just a few hundred data simulations. For the weak lensing problem, six cosmological parameters can be inferred from $\mathcal{O}(10^3)$ simulations, irrespective of whether ten additional nuisance parameters are included in the problem or not.

Figures

Equations

$$\mathbf{\varphi }\leftarrow q(\mathbf{\varphi }),$$

$$\mathbf{d}\leftarrow \mathrm{simulator}(\mathbf{d}|\mathbf{\varphi }),$$

$$\mathbf{t}=\mathrm{compression}(\mathbf{d}),$$

$$p(\mathbf{\varphi }|{\mathbf{t}}_o)\propto p({\mathbf{t}}_o|\mathbf{\varphi })p(\mathbf{\varphi }).$$

$$\mathbf{t}=\nabla {\mathcal{L}}_{\ast }$$

$$\mathbf{\varphi }\leftarrow q(\mathbf{\varphi }),$$

$$\mathbf{d}\leftarrow \mathrm{simulator}(\mathbf{d}|\mathbf{\varphi }),$$

$$\mathbf{t}=\mathrm{compression}(\mathbf{d}),$$

$$\mathrm{input}\mathrm{parameters}\mathbf{\varphi },$$

$${\mathbf{z}}_1\leftarrow p({\mathbf{z}}_1|\mathbf{\varphi }),$$

$${\mathbf{z}}_2\leftarrow p({\mathbf{z}}_2|{\mathbf{z}}_1,\mathbf{\varphi }),$$

$$⋮$$

$${\mathbf{z}}_n\leftarrow p({\mathbf{z}}_n|{\mathbf{z}}_1,\dots ,{\mathbf{z}}_{n-1},\mathbf{\varphi }),$$

$$\mathbf{d}\leftarrow p(\mathbf{d}|{\mathbf{z}}_1,\dots ,{\mathbf{z}}_n,\mathbf{\varphi }),$$

$$p(\mathbf{t}|\mathbf{\varphi })\equiv \int p({\mathbf{z}}_1|\mathbf{\varphi })$$

$$\times p(\mathbf{t}|{\mathbf{z}}_1,\dots ,{\mathbf{z}}_n,\mathbf{\varphi })d{\mathbf{z}}_{1}d{\mathbf{z}}_2\dots d{\mathbf{z}}_n.$$

$$\mathbf{\theta }\leftarrow q(\mathbf{\theta }),$$

$$\mathbf{\eta }\leftarrow p(\mathbf{\eta }|\mathbf{\theta }),$$

$$\mathbf{d}\leftarrow \mathrm{simulator}(\mathbf{d}|\mathbf{\varphi },\mathbf{\eta }),$$

$$\mathbf{t}=\mathrm{compression}(\mathbf{d}).$$

$$\mathcal{L}(\mathbf{\varphi })$$

$$\approx {\mathcal{L}}_{\ast }+{\mathbf{t}}^{\mathrm{T}}\delta \mathbf{\varphi }-\frac{1}{2}\delta {\mathbf{\varphi }}^{\mathrm{T}}\mathbf{F}\delta \mathbf{\varphi },$$

$$\bar{\mathcal{L}}(\mathbf{\theta })$$

$$={\mathcal{L}}_{\ast }+{\mathbf{t}}_{\mathbf{\theta }}^{\mathrm{T}}\delta \mathbf{\theta }-\frac{1}{2}\delta {\mathbf{\theta }}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{\theta }\mathbf{\theta }}\delta \mathbf{\theta }$$

$$+\frac{1}{2}({\mathbf{F}}_{\mathbf{\eta }\mathbf{\theta }}\delta \mathbf{\theta }-{\mathbf{t}}_{\mathbf{\eta }}{)}^{\mathrm{T}}{\mathbf{F}}_{\mathbf{\eta }\mathbf{\eta }}^{-1}({\mathbf{F}}_{\mathbf{\eta }\mathbf{\theta }}\delta \mathbf{\theta }-{\mathbf{t}}_{\mathbf{\eta }})$$

$$\displaystyle\mathbf{F}\equiv\left(\begin{array}[]{@{}cc}\mathbf{F}_{\boldsymbol{\theta}\boldsymbol{\theta}}&\mathbf{F}_{\boldsymbol{\theta}\boldsymbol{\eta}}\\ \mathbf{F}_{\boldsymbol{\eta}\boldsymbol{\theta}}&\mathbf{F}_{\boldsymbol{\eta}\boldsymbol{\eta}}\end{array}\right).$$

$${\bar{\mathbf{t}}}_{\mathbf{\theta }}={\mathbf{t}}_{\mathbf{\theta }}-{\mathbf{F}}_{\mathbf{\theta }\mathbf{\eta }}{\mathbf{F}}_{\mathbf{\eta }\mathbf{\eta }}^{-1}{\mathbf{t}}_{\mathbf{\eta }}.$$

$$p(\mathbf{\eta }|\mathbf{d})$$

$$\approx \int e^{-\frac{1}{2}(\mathbf{\eta }-{\mathbf{\mu }}_{\mathbf{\eta }}(\mathbf{\theta }){)}^{†}{\left({\mathbf{F}}_{\mathbf{\eta }\mathbf{\eta }}^{-1}-{\mathbf{F}}_{\mathbf{\eta }\mathbf{\theta }}^{-1}{\left[{\mathbf{F}}_{\mathbf{\theta }\mathbf{\theta }}^{-1}\right]}^{-1}{\mathbf{F}}_{\mathbf{\theta }\mathbf{\eta }}^{-1}\right)}^{-1}(\mathbf{\eta }-{\mathbf{\mu }}_{\mathbf{\eta }}(\mathbf{\theta }))}$$

$$\times p(\mathbf{\theta }|{\bar{\mathbf{t}}}_{\mathbf{\theta }})d\mathbf{\theta },$$

$${\mathbf{\mu }}_{\mathbf{\eta }}(\mathbf{\theta })=\hat{\mathbf{\eta }}+{\mathbf{F}}_{\mathbf{\eta }\mathbf{\theta }}^{-1}{\left[{\mathbf{F}}_{\mathbf{\theta }\mathbf{\theta }}^{-1}\right]}^{-1}(\mathbf{\theta }-\hat{\mathbf{\theta }}),$$

$$m_{\mathrm{B}}=5{\log }_{10}\left[\frac{D_{\mathrm{L}}^{\ast }(z;\mathbf{\theta })}{10\mathrm{pc}}\right]$$

$$+M_{\mathrm{B}}+\delta M\Theta (M_{\mathrm{stellar}}-10^{10}{M}_{\odot }),$$

$$D_{\mathrm{L}}^{\ast }(z;\mathbf{\theta })=\frac{(1+z)c}{100}{\int }_0^z\frac{d{z}^{\prime }}{\sqrt{{\Omega }_{\mathrm{m}}(1+z^{\prime }{)}^3+(1-{\Omega }_{\mathrm{m}})(1+z^{\prime }{)}^{3(w_0+1)}}},$$

$$\ln p(\mathbf{d}|\mathbf{\varphi })=-\frac{1}{2}(\mathbf{d}-\mathbf{\mu }(\mathbf{\varphi }){)}^T{\mathbf{C}}^{-1}(\mathbf{d}-\mathbf{\mu }(\mathbf{\varphi }))-\frac{1}{2}\ln |\mathbf{C}|,$$

$$\mathbf{t}\equiv {\nabla }_{\mathbf{\varphi }}{\mathcal{L}}_{\ast }={\nabla }_{\mathbf{\varphi }}^T{\mathbf{\mu }}_{\ast }{\mathbf{C}}^{-1}(\mathbf{d}-{\mathbf{\mu }}_{\ast }),$$

$${\bar{\mathbf{t}}}_{\mathbf{\theta }}={\mathbf{t}}_{\mathbf{\theta }}-{\mathbf{F}}_{\mathbf{\theta }\mathbf{\eta }}{\mathbf{F}}_{\mathbf{\eta }\mathbf{\eta }}^{-1}{\mathbf{t}}_{\mathbf{\eta }},$$

$$\displaystyle\boldsymbol{\mu}_{\mathrm{P,\boldsymbol{\theta}}}=(0.3,\;-0.75),\;\mathbf{C}_{\mathrm{P,\boldsymbol{\theta}}}=\left({\begin{array}[]{cc}0.4^{2}&-0.24\\ -0.24&0.75^{2}\\ \end{array}}\right),$$

$$C_{\ell ,\alpha \beta }^{\gamma \gamma }$$

$$w_{\alpha }(\chi )=\frac{3{\Omega }_{\mathrm{m}}{H}_0^2}{2}\chi {\int }_{\chi }^{{\chi }_{\mathrm{H}}}d{\chi }^{\prime }{n}_{\alpha }({\chi }^{\prime })\frac{{\chi }^{\prime }-\chi }{{\chi }^{\prime }},$$

Full text

††pubyear: 2018††pagerange: Nuisance hardened data compression for fast likelihood-free inference–A

1 Introduction

Likelihood-free inference is emerging as a new paradigm for performing principled inference from cosmological data using forward simulations only, incorporating exactly all known effects that can be successfully simulated without relying on approximate likelihoods (Schafer & Freeman, 2012; Cameron & Pettitt, 2012; Weyant et al., 2013; Robin et al., 2014; Lin & Kilbinger, 2015; Ishida et al., 2015; Akeret et al., 2015; Jennings et al., 2016; Hahn et al., 2017; Kacprzak et al., 2017; Carassou et al., 2017; Davies et al., 2017; Alsing et al., 2018b, 2019).

In Alsing et al. (2018b, 2019) we showed that likelihood-free inference is achievable from just $\mathcal{O}(10^{3})$ forward simulations for typical problems in cosmology, with $n\sim 6$ or so cosmological parameters of interest. However, while many problems have a small number of interesting parameters, they often come with a large number of additional nuisance parameters that ultimately need marginalizing over. In the standard Bayesian approach, interesting and nuisance parameters are first inferred together and the resulting joint posterior is then marginalized over the nuisances a posteriori. For likelihood-free methods, solving the initially higher-dimensional inference task (over interesting and nuisance parameters together) would typically demand a much larger number of forward simulations, even though the nuisance parameters are later marginalized out anyway. When forward simulations are expensive, this may pose a critical impasse for likelihood-free analyses.

In this paper we show that in the likelihood-free setting, nuisance marginalized posteriors can be inferred directly from forward simulations, bypassing the higher-dimensional “interesting plus nuisance parameters" inference task entirely. The result is that for problems with a given number of interesting parameters, the number of simulations required for performing likelihood-free inference is (roughly) the same irrespective of the number of additional nuisance parameters that need marginalizing over. For problems in cosmology and astrophysics where often a modest number of interesting parameters are accompanied by a large number of nuisances, this is a major step toward enabling likelihood-free analyses where simulations are expensive.

To enable direct inference of nuisance marginalized posteriors we introduce two extensions to the standard likelihood-free inference set-up. Firstly, we re-cast nuisance parameters as local latent variables in the forward simulations so they are automatically marginalized over with no additional approximations or loss of information. This way, unknown nuisances are effectively treated as an additional source of noise in the problem, allowing nuisance marginalized posteriors to be inferred directly.

Secondly, we derive an asymptotically optimal compression of $N$ data down to $n$ summaries – one per interesting parameter – such that the Fisher information is preserved, but the summaries are insensitive (to leading order) to the nuisance parameters. This allows us to keep the number of informative “nuisance hardened" summary statistics at one per interesting parameter, irrespective of the number of nuisance parameters in the problem.

With nuisance hardened compressed summaries and nuisances cast as latent parameters, the nuisance marginalized likelihood-free inference task involves learning $n$ interesting parameters from data compressed down to $n$ informative summaries, regardless of the presence or number of additional nuisance parameters. The complexity and hence number of simulations required will then be (roughly) independent of the number of nuisances, expanding the range of problems that can be feasibly tackled with likelihood-free methods.

The plan of the paper is as follows. In §2 and 3 we review the key ideas of likelihood-free inference and data-compression. In §4 we show how to re-cast nuisances as latent variables so they are automatically marginalized over in the likelihood-free setting. In §5 we derive an asymptotically optimal scheme for compressing data down to summary statistics that are hardened to nuisance parameters. In §6–7 we validate the likelihood-free inference of nuisance marginalized posteriors on two case studies from cosmology: cosmological parameter inference from supernovae and cosmic shear data, with nuisance parameterized systematics. We conclude in §8.

2 likelihood-free inference

The ideas presented in this paper are very general and can be applied to any likelihood-free inference approach. However, for simplicity we will frame the discussion in terms of density-estimation likelihood-free inference (DELFI; Papamakarios et al., 2018; Lueckmann et al., 2018; Alsing et al., 2018b, 2019), which we briefly review below.

Suppose we have $p$ parameters $\boldsymbol{\phi}$ that we want to learn from data $\mathbf{d}$, that we have compressed down to some small number of informative summaries $\mathbf{t}\equiv\mathbf{t}(\mathbf{d})$. In §3 we will review two key methods for compressing large datasets down to $p$ summaries (that we will later extend to nuisance hardened data compression schemes), so throughout the paper we will also assume that $\mathbf{t}\in\mathbb{R}^{p}$.

Density-estimation likelihood-free inference proceeds in three steps:

Generate parameter-data summary pairs $\{\boldsymbol{\phi},\mathbf{t}\}$ by running forward simulations,

[TABLE]

where $q(\boldsymbol{\phi})$ is some parameter proposal distribution for choosing where to run simulations111The parameter proposal may be adaptively modified to optimize acquisition of relevant simulated parameter-data pairs (Papamakarios et al., 2018), or Bayesian optimization-style acquisition rules may be adopted for even more optimal simulation acquisition (Lueckmann et al., 2018). and left-arrows indicate drawing random realizations. 2. 2.

Fit to these simulations $\{\boldsymbol{\phi},\mathbf{t}\}$ a flexible parameterized model for the conditional density, $p(\mathbf{t}|\boldsymbol{\phi})$ – the sampling distribution of the summary statistics as a function of the parameters222There are three choices for fitting a density to $\{\boldsymbol{\phi},\mathbf{t}\}$: a conditional density estimator for the posterior $p(\boldsymbol{\phi}|\mathbf{t})$ (Papamakarios & Murray, 2016; Lueckmann et al., 2017), a model for the joint density $p(\boldsymbol{\phi},\mathbf{t})$ (Alsing et al., 2018b), or a conditional density estimator for the likelihood $p(\boldsymbol{\phi}|\mathbf{t})$ (Papamakarios et al., 2018; Lueckmann et al., 2018). In this paper we will focus on the latter; see Alsing et al. (2019) and Papamakarios et al. (2018) for the relative merits of these three approaches.. Mixture density networks (MDN) and masked autoregressive flows (MAF) have proven useful conditional density estimators in this context (Papamakarios & Murray, 2016; Papamakarios et al., 2018; Alsing et al., 2019). 3. 3.

Use the learned sampling density $p(\mathbf{t}|\boldsymbol{\phi})$ evaluated at the observed data $\mathbf{t}_{o}$ in Bayes’ theorem to obtain the (target) posterior:

[TABLE]

The full inference task is hence reduced to learning a $p$-dimensional density $p(\mathbf{t}|\boldsymbol{\phi})$, as a function of $p$ parameters $\boldsymbol{\phi}$, from a set of simulated parameter-data pairs $\{\boldsymbol{\phi},\mathbf{t}\}$.

In order to enable fast likelihood-free inference, requiring the fewest simulations possible, we want to keep the dimensionality on the left and right and sides of the condition in $p(\mathbf{t}|\boldsymbol{\phi})$ as low as possible for the problem we are trying to solve.

3 Data compression and summary statistic choice

When performing likelihood-free inference, it is usually necessary to compress large datasets down to a manageable number of informative summary statistics. This data compression should ideally be done so that is preserves as much of the information content of the full data (with respect to the parameters of interest) as possible. In this section we review three key data-compression schemes that summarize $N$ down to $p$ numbers – one per parameter – with the aim of preserving the (Fisher) information content of the data.

3.1 Approximate score compression

When the likelihood is known, the score function – the gradient of the log-likelihood with respect to the parameters – yields asymptotically optimal compression of $N$ data down to $p$ summaries (Alsing & Wandelt, 2018). Hence,

[TABLE]

represent a natural choice of summaries for parameters $\boldsymbol{\phi}$, where $\mathcal{L}$ denotes the log-likelihood, $*$ denotes evaluation about some fiducial parameters $\boldsymbol{\phi}_{*}$ and $\nabla\equiv\frac{\partial}{\partial\boldsymbol{\phi}}$. When the gradient is taken close to the maximum-likelihood parameters, the summaries $\mathbf{t}$ saturate the information inequality, preserving the Fisher information content of the data. If a good expansion point is not available a priori, $\boldsymbol{\phi}_{*}$ can be iterated toward to the maximum likelihood for improved optimality if needed (Alsing & Wandelt, 2018).

Whilst the score presents a natural summary statistic choice for likelihood-free inference, it is only available if the likelihood-function is known. For likelihood-free applications where the likelihood is not known a priori, the score may still provide a guiding principle for choosing informative summaries. For many problems, an approximate likelihood (eg., Gaussian) may be used for the purposes of data compression, the only price being some loss of optimality. If no obvious likelihood approximation presents itself, once can learn the conditional density $p(\mathbf{d}|\boldsymbol{\theta})$ from simulations locally in a small neighborhood around $\boldsymbol{\phi}_{*}$, which can then be used to define an approximate score. Alternatively, the score-function may be regressed directly from simulations Brehmer et al. (2018a, b, c).

Note that for Gaussian data where the sensitivity to the model parameters is only in the mean, score-compression is equivalent to moped (Heavens et al., 2000), and where the sensitivity to the parameters is only in the covariance it is equivalent to the optimal quadratic estimator (Tegmark et al., 1997).

3.2 Neural network parameter estimators

An emerging trend in cosmology is to find efficient parameter estimators from complex datasets by training deep neural networks to regress cosmological parameters from simulated data (Ravanbakhsh et al., 2016; Schmelzle et al., 2017; Gupta et al., 2018; Ribli et al., 2018; Fluri et al., 2018a). The resulting trained neural networks can be viewed as radical data compression schemes, compressing large data sets down to a set of parameter estimators, whose sampling distributions are unknown. These neural network parameter estimators can be straightforwardly used in a subsequent likelihood-free analysis to obtain Bayesian posteriors from neural network summarized data sets.

3.3 Information maximizing neural networks

Combining the ideas of Fisher information preserving score-compression and neural network parameter estimators, Information Maximizing Neural Networks (imnn; Charnock et al., 2018) parameterize the data compression function $\mathbf{t}(\mathbf{d}):\mathbb{R}^{N}\rightarrow\mathbb{R}^{p}$ as a neural network, and train the network on a set of forward simulations of the data such that it maximizes the retained Fisher information content of the compressed summaries. The result is a likelihood-free Fisher optimal data compression scheme that is trained only on forward simulations, without any further (optimality-reducing) assumptions. In contrast to neural network parameter estimators, imnns only need to run simulations around some fiducial model. This simplifies the learning task and has empirically been shown in test cases to give optimal compression with simple network architectures (Charnock et al., 2018).

In maximizing the Fisher information content of the summaries, imnns also yield an estimate of the Fisher information matrix as a byproduct; this Fisher matrix estimate is essential for projecting nuisance parameters out of nuisance hardened summaries as we will see in §5.

4 Likelihood-free inference of nuisance marginalized posteriors

Suppose now that our $p$ parameters $\boldsymbol{\phi}=(\boldsymbol{\theta},\boldsymbol{\eta})$ constitute $n$ interesting parameters, $\boldsymbol{\theta}$, and $m$ nuisance parameters, $\boldsymbol{\eta}$.

When performing likelihood-free inference, parameter-data pairs $\{\boldsymbol{\phi},\mathbf{t}\}$ are generated by running forward simulations:

[TABLE]

where $q(\boldsymbol{\phi})$ is some proposal distribution. The “simulator" is just a set of probabilistic statements that generate a realization of data $\mathbf{d}$ from some input (hyper-) parameters $\boldsymbol{\phi}$, via some intermediate local latent variable layers $\{\mathbf{z}\}$:

[TABLE]

where each latent layer depends on some or all of the parameters that appear upstream in the simulation. Typical examples of latent variables appearing in cosmological forward simulations might be the amplitudes and phases of the initial potential perturbations in the Universe, the (true) redshifts and other physical properties of galaxies in a survey, realizations of stochastic foregrounds or detector artifacts, etc.

In the likelihood-free set-up we aim to fit a flexible model for $p(\mathbf{t}|\boldsymbol{\phi})$ to a set of simulated pairs $\{\boldsymbol{\phi},\mathbf{t}\}$. Crucially, the target $p(\mathbf{t}|\boldsymbol{\phi})$ is implicitly (by construction) marginalized over all of the local latent variables, since

[TABLE]

This means that if we just re-cast the nuisance parameters as an additional intermediate latent layer in the forward simulations, they will implicitly be marginalized over in the subsequent likelihood-free inference. This is easy: factorizing the prior over interesting and nuisance parameters as $p(\boldsymbol{\phi})=p(\boldsymbol{\theta})p(\boldsymbol{\eta}|\boldsymbol{\theta})$, then data-interesting parameter pairs $\{\boldsymbol{\theta},\mathbf{t}\}$ can be generated from simulations as,

[TABLE]

These data-interesting parameter pairs $\{\boldsymbol{\theta},\mathbf{t}\}$ can then be fit with a conditional density estimator for $p(\mathbf{t}|\boldsymbol{\theta})$, implicitly marginalizing over the now local latent nuisance parameters $\boldsymbol{\eta}$. Note that no further approximations have been made in changing the target density from $p(\mathbf{t}|\boldsymbol{\theta},\boldsymbol{\eta})$ to the nuisance-marginalized $p(\mathbf{t}|\boldsymbol{\theta})$. Uncertain nuisance parameters and their degeneracies with the interesting parameters are fully accounted for in the forward simulations Eq. (4), so the estimated $p(\mathbf{t}|\boldsymbol{\theta})$ will converge to the marginalized likelihood $\int p(\mathbf{t}|\boldsymbol{\theta},\boldsymbol{\eta})p(\boldsymbol{\eta}|\boldsymbol{\theta})d\boldsymbol{\eta}$.

The nuisance-marginalized inference task has hence been reduced from inferring a $p$-dimensional density333Assuming there are $p$ compressed summaries $\mathbf{t}\in\mathbb{R}^{p}$, ie., one per parameter. This will be the case if either the score of an approximate likelihood or an information maximizing neural network is used for the summary statistics (see §3; Alsing et al., 2018b; Charnock et al., 2018). as a function of $p$ parameters $p(\mathbf{t}|\boldsymbol{\theta},\boldsymbol{\eta})$, to inferring a $p$-dimensional density as a function of only the $n=p-m$ interesting parameters $p(\mathbf{t}|\boldsymbol{\theta})$. When the number of nuisance parameters is even modestly large, this already gives a major reduction in the computational complexity of the likelihood-free inference task.

5 Nuisance parameter hardened data compression

Having successfully reduced the number of parameters that must be fit for in the inference step (§4), the goal of this section is to reduce the number of data summaries that need to be considered, from $\mathbf{t}\in\mathbb{R}^{p}$ to some smaller set $\bar{\mathbf{t}}_{\theta}\in\mathbb{R}^{n}$ – one per interesting parameter.

Taking either the (approximate) score or a trained imnn as the compression scheme, each compressed summary $t_{\phi_{i}}$ is constructed to be as informative as possible about the corresponding parameter $\phi_{i}$. However, degeneracies between parameters mean that each summary will also in general be sensitive to all other parameters. Our goal is to find a reduced set of “nuisance hardened" summary statistics $\bar{\mathbf{t}}_{\theta}\in\mathbb{R}^{n}$ that are asymptotically optimal summaries for the $n$ interesting parameters (preserving the Fisher information) whilst having their sensitivity to the nuisance parameters projected out.

For a given likelihood $\mathcal{L}$, the Fisher information maximizing summaries for $n$ interesting parameters that are insensitive to some $m$ nuisances is just the score of the likelihood marginalized over the nuisance parameters, evaluated at some fiducial parameter values. We use this as our guiding principle for deriving nuisance hardened summary statistics. In the short derivation below we will work in the regime where the likelihood-function is (assumed) known: for likelihood-free applications, the derived result can be applied to the approximate methods for estimating the score described in §3, or imnn compression.

Nuisance hardened summaries can be obtained given some likelihood-function $\mathcal{L}$ as follows. Firstly, we Taylor expand the log-likelihood to second order in the parameters about some expansion point $\boldsymbol{\phi}_{*}$,

[TABLE]

where $\mathbf{t}\equiv\nabla\mathcal{L}_{*}$ is the score and $\mathbf{J}\equiv-\nabla\nabla^{\mathrm{T}}\mathcal{L}_{*}$ the observation matrix, which we have replaced by its expectation value – the Fisher matrix $\mathbf{F}=-\langle\nabla\nabla^{\mathrm{T}}\mathcal{L}_{*}\rangle$ – in the second line. Gradients are defined as $\nabla\equiv(\nabla_{\boldsymbol{\theta}},\nabla_{\boldsymbol{\eta}})$.

Next, we marginalize the approximate likelihood over the nuisance parameters. In the quadratic expansion of the log-likelihood this is just a Gaussian integral, which after some algebra (dropping $\boldsymbol{\theta}$-independent terms) gives:

[TABLE]

where $\mathbf{t}_{\boldsymbol{\theta}}\equiv\nabla_{\boldsymbol{\theta}}\mathcal{L}_{*}$, $\mathbf{t}_{\boldsymbol{\eta}}\equiv\nabla_{\boldsymbol{\eta}}\mathcal{L}_{*}$, and the various blocks of the Fisher matrix are defined from:

[TABLE]

Finally, we obtain the nuisance-hardened summaries as the score of the approximate nuisance-marginalized likelihood $\bar{\mathbf{t}}_{\boldsymbol{\theta}}\equiv\nabla_{\boldsymbol{\theta}}\bar{\mathcal{L}}_{*}$, giving:

[TABLE]

This is the main result of this section: Eq. (13) tells us how to form nuisance-hardened summaries $\bar{\mathbf{t}}_{\boldsymbol{\theta}}$ from score statistics $\mathbf{t}$, where the only extra ingredient required for projecting out the nuisance parameter sensitivities is the Fisher matrix.

The projection in Eq. (13) makes good intuitive sense; expected covariances $\mathbf{F}_{\boldsymbol{\theta}\boldsymbol{\eta}}$ between interesting and nuisance parameters are projected out, weighted by the expected marginal variances of the nuisance parameters $\mathbf{F}^{-1}_{\boldsymbol{\eta}\boldsymbol{\eta}}$.

How successfully the projection in Eq. (13) removes the nuisance parameter dependence of the summaries depends on how well the approximations in Eq. (5) are satisfied. In the asymptotic (or linear-model) limit where the likelihood is Gaussian in the parameters, the projection and hence nuisance-hardening will be exact. When the likelihood is non-Gaussian in the parameters or there are non-linear degeneracies between interesting and nuisance parameters, the projection will be approximate. In this case, as always with likelihood-free inference, the only price to pay is some loss of optimality; the nuisance-hardened summaries will be sub-optimal, but provided the same compression scheme is applied to the data and simulations no biases are introduced.

For likelihood-free applications, Eq. (13) can be applied to approximate score summaries (see §3). Alternatively, the projection in Eq. (13) can also be straightforwardly applied to information maximizing neural networks. In maximizing the Fisher information, imnns implicitly learn the score of the (true) likelihood some (non-linear) transformed data vector. The trained imnn gives the compressed summaries $\mathbf{t}$ and the Fisher information matrix $\mathbf{F}$ as output, to which Eq. (13) can be applied directly to project out the nuisance parameters.

Note that the result in Eq. (13) was also obtained by a different route in Zablocki & Dodelson (2016). In that work they sought the linear combinations of (Gaussian) data that maximize the Fisher information for each parameter in turn, whilst constrained to be insensitive (to leading order in the expectation) to all other parameters. Our results show that the Zablocki & Dodelson (2016) result can be interpreted as the score-compression of any likelihood, marginalized over nuisance parameters in the Laplace approximation.

With the nuisance-hardened summary statistics in hand and nuisances re-cast as local latent variables (§4), the nuisance marginalized inference task is now reduced to learning an $n$-dimensional density $p(\bar{\mathbf{t}}_{\boldsymbol{\theta}}|\boldsymbol{\theta})$ as a function of the $n$ interesting parameters only, irrespective of the presence or number of additional nuisance parameters to be marginalized over. For even a modest number of nuisance parameters this represents a significant reduction in computational complexity compared to inferring the higher-dimensional $p(\mathbf{t}|\boldsymbol{\theta},\boldsymbol{\eta})$ and marginalizing over $\boldsymbol{\eta}$ a posteriori.

5.1 Recovering an approximate nuisance parameter posterior

It is possible to construct an approximate marginal nuisance parameter posterior a posteriori, from the nuisance hardened analysis. The nuisance posterior can be approximated as follows:

[TABLE]

where in the second line we replace the marginal interesting-parameter posterior by the nuisance hardened version, and take a Laplace approximation for the nuisance conditional $p(\boldsymbol{\theta},\boldsymbol{\eta}|\mathbf{d})$ about $(\boldsymbol{\theta}_{*},\boldsymbol{\eta}_{*})$ (cf., Eq. 5). The conditional mean in the Laplace approximation term is given by

[TABLE]

where $(\hat{\boldsymbol{\theta}},\hat{\boldsymbol{\eta}})=(\boldsymbol{\theta}_{*},\boldsymbol{\eta}_{*})+\mathbf{F}^{-1}\mathbf{t}$ is one Newton iteration towards the maximum of the approximate likelihood used to obtain the Fisher matrix and score. The integral in Eq. (5.1) can then be approximated as a sum over samples from the interesting-parameter posterior, giving a pseudo Blackwell-Rao estimator for the nuisance parameter posterior.

While likelihood-free inference using nuisance-hardened summaries is robust in the sense that approximations made in the compression/nuisance projection steps cannot bias the parameter inferences, the same cannot be said of the above approximate nuisance parameter posterior in Eq. (5.1). The location and scale of the approximate nuisance posterior is determined (partly) by the approximate likelihood used for performing the compression, ie., the first term in the integral Eq. (5.1). If the approximate likelihood used for the compression is poor, the location and scale of the nuisance posterior will be biased, and in the limit where the nuisance and interesting parameters are independent, Eq. (5.1) just gives a Gaussian approximation of the approximate likelihood. Nevertheless, Eq. (5.1) gives a rough inference of the nuisances when the approximate likelihood is reasonable, and will improve as correlation between the interesting and nuisance parameters becomes more important.

6 validation case I: supernovae data analysis

In this section we validate likelihood-free inference of nuisance marginalized posteriors on a simple case: inferring cosmological parameters from supernovae data in the presence of systematics.

Supernova data provide a great opportunity for likelihood-free methods, since the data are impacted by a large number of systematic biases and selection effects that need to be carefully accounted for to obtain robust cosmological parameter inferences.

Here, for the purpose of validation, we perform a simple analysis of the JLA data (Betoule et al., 2014) under assumptions that allow us to compare the likelihood-free results against an exact (known) likelihood. The set-up is identical to Alsing et al. (2018b), which we review briefly below.

6.1 JLA data and model

The JLA sample is comprised of $740$ type Ia supernovae for which we have estimates for the apparent magnitudes $m_{\mathrm{B}}$, redshifts $z$, color at maximum-brightness $C$ and stretch $X_{1}$ parameters characterizing the lightcurves. We take the data vector to be the vector of estimated apparent magnitudes $\mathbf{d}=(\hat{m}_{\mathrm{B}}^{1},\hat{m}_{\mathrm{B}}^{2},\dots,\hat{m}_{\mathrm{B}}^{M})$, with uncertainties in the redshift, color and stretch approximately accounted for in the covariance matrix of the observed apparent magnitudes (Betoule et al., 2014).

We assume that the apparent magnitudes of type Ia supernovae depend on the luminosity distance to the source at the given redshift $D^{*}_{\mathrm{L}}(z)$, a reference absolute magnitude for type Ia supernovae, and calibration corrections for the stretch $X_{1}$ and color at maximum-brightness $C$ (Tripp, 1998),

[TABLE]

where $\boldsymbol{\theta}$ are the cosmological parameters (see below), $\alpha$ and $\beta$ are calibration parameters for the stretch and color, and $\tilde{M}_{\mathrm{B}}$ and $\delta M$ characterize the host stellar-mass ($M_{\mathrm{stellar}}$) dependent reference absolute magnitude. $\Theta$ is the Heaviside function.

The cosmological model enters via the luminosity distance-redshift relation: we will assume a flat $w$CDM universe with cold dark matter (with total matter density parameter $\Omega_{\mathrm{m}}$) and dark energy characterized by equation-of-state $p/\rho=w_{0}$. The luminosity distance-redshift relation is given by

[TABLE]

where $c$ is the speed of light in vacuum. The cosmological parameters of interest are $\boldsymbol{\theta}=(\Omega_{\mathrm{m}},w_{0})$, and we treat the remaining parameters $\boldsymbol{\eta}=(\alpha,\beta,M_{\mathrm{B}},\delta M)$ as nuisances.

6.2 Likelihood and data compression

For this validation case we assume the data are Gaussian,

[TABLE]

with mean given by Eq. (6.1), and we assume a fixed covariance matrix from Betoule et al. (2014) (see Alsing et al., 2018b for details of the covariance matrix).

We take compressed data summaries to be the score of the Gaussian likelihood, ie.,

[TABLE]

where we take fiducial parameters for the expansion point $\boldsymbol{\theta}_{*}=(0.202,-0.748,-19.04,0.126,2.644,-0.0525)$444Found in a few iterations of $\hat{\boldsymbol{\theta}}_{k+1}=\hat{\boldsymbol{\theta}}_{k}+\mathbf{F}^{-1}_{k}\mathbf{t}_{k}$, and ‘$*$’ indicated evaluation at the fiducial parameters.

Projection of the nuisance parameters is performed following Eq. (13), giving nuisance hardened summary statistics for the interesting (cosmological) parameters:

[TABLE]

where the Fisher information matrix is given by $\mathbf{F}=\nabla_{\boldsymbol{\phi}}\boldsymbol{\mu}^{T}\mathbf{C}^{-1}\nabla_{\boldsymbol{\phi}}^{T}\boldsymbol{\mu}$.

For a full sophisticated implementation of likelihood-free inference to supernova data analysis, in practise, may use an approximate likelihood or information maximizing neural network for performing data compression.

6.3 Simulations

For this validation case, simulations are just draws from the (exact) assumed sampling distribution of the data, ie., drawing Gaussian data from Eq. (18) (given parameters).

6.4 Priors

We assume broad Gaussian priors on the cosmological parameters $\boldsymbol{\theta}=(\Omega_{\mathrm{m}},w_{0})$ with mean and covariance (following Alsing et al., 2018b):

[TABLE]

with additional hard prior boundaries $\Omega_{\mathrm{m}}\in[0,0.6]$ and $w_{0}\in[-1.5,0]$.

We take independent Gaussian priors on the nuisance parameters $\boldsymbol{\eta}=(\alpha,\beta,M_{\mathrm{B}},\delta M)$, with means $\boldsymbol{\mu}_{\mathrm{P,\boldsymbol{\eta}}}=(-19.05,\;0.125,\;2.6,\;-0.05)$ and standard deviations $\boldsymbol{\sigma}_{\mathrm{P,\boldsymbol{\eta}}}=(0.1,\;0.25,\;0.025,\;0.05)$.

6.5 DELFI set-up

The density-estimation likelihood-free inference (DELFI) algorithm learns the sampling distribution of the data summaries as a function of the parameters: we parameterize the data sampling distribution by a mixture density network with three Gaussian components, and simulations are ran in sequential batches with an adaptive proposal density (following Papamakarios et al., 2018; Alsing et al., 2019). Further technical details of the network architectures and DELFI algorithm set-up used here can be found in Appendix A. We used the public DELFI implementation pydelfi555https://github.com/justinalsing/pydelfi (Alsing et al., 2019).

6.6 Results

We ran the JLA inference problem in two scenarios. First, we inferred the posterior for the six interesting plus nuisance parameters jointly, ie., learning the $6$-dimensional density $p(\mathbf{t}_{\boldsymbol{\theta}},\mathbf{t}_{\boldsymbol{\eta}}|\boldsymbol{\theta},\boldsymbol{\eta})$ as a function of the six parameters. Secondly, we inferred the posterior for the two interesting parameters directly marginalized over the nuisances, ie., learning the 2-dimensional density $p(\bar{\mathbf{t}}_{\boldsymbol{\theta}}|\boldsymbol{\theta})$ as a function of the two cosmological parameters.

In Figure 1 we see that the nuisance marginalized posterior obtained in the likelihood-free setting is in excellent agreement with exact marginal posterior (obtained by sampling the full 6-parameter likelihood with MCMC). This means that for this particular example, the nuisance hardened data compression described in §5 is effectively lossless. This is expected if the Laplace approximation is a good approximation for integrating out the nuisance parameters, the likelihood used for the score compression is close to the true likelihood (in this validation case it is exact), and the expansion is performed close to the peak. The further one deviates from those conditions, the more lossy we can expect the compression and projection to be.

Figure 2 shows the convergence of the negative log-likelihood of the simulated parameter-data summary samples under the learned neural density estimator. For this problem, inferring the nuisance marginalized posterior directly gives roughly an order of magnitude reduction in the number of simulations required: the two-parameter problem converges after a few hundred simulations, while the full six-parameter problem takes a few thousand. We note also that for the two-parameter problem, the number of simulations required here is an order of magnitude fewer than was reported for the same problem using Bayesian Optimization Likelihood-free Inference (BOLFI) in Leclercq (2018) (in that work, the nuisance parameters were kept fixed rather than marginalized over as they are here). Note also that this small number of simulations is also in spite of assuming a flexible model for the sampling distribution of the data (summaries); a three component Gaussian mixture, where in this case the data are actually Gaussian.

Figure 3 shows the recovered Blackwell-Rao estimator for the approximate nuisance parameter posterior. In this specific case, since the likelihood used to perform the compression was exact, the assumptions underpinning the Blackwell-Rao estimated nuisance posterior are satisfied and the recovered nuisance posterior approximation is hence very reasonable.

7 validation case II: cosmic shear with photo-z systematics

In this section we demonstrate the utility of inferring nuisance marginalized posteriors in the likelihood-free setting on a second validation example: inferring cosmological parameters from cosmic shear observations, with nuisance parameters describing systematic uncertainties in the photometric redshift distributions.

Cosmic shear data are a good candidate for likelihood-free analysis, since they are impacted by a number of effects that can be included (to varying degrees) in forward simulations but are hard to build into a high-fidelity likelihood function. These include non-linear physics and baryonic feedback on small scales (Rudd et al., 2008; Harnois-Déraps et al., 2015), intrinsic alignments (Joachimi et al., 2015), shape and photo-$z$ measurement systematics (Massey et al., 2012; Mandelbaum, 2018; Salvato et al., 2018), blending (Mandelbaum, 2018), reduced shear corrections (Krause & Hirata, 2010), redshift-dependent galaxy-property distributions (Kannawadi et al., 2018), non-trivial sampling distributions for common summary statistics (Sellentin et al., 2018), etc., all of which can result in biased inferences if not carefully accounted for. Likelihood-free methods may also allow us to extract extra information from non-standard lensing observables (eg., magnification Hildebrandt et al., 2009; van Waerbeke, 2010; Hildebrandt et al., 2013; Duncan et al., 2013; Heavens et al., 2013; Alsing et al., 2015), and summary statistics such as convergence peaks (Kratochvil et al., 2010; Fluri et al., 2018b), bispectra (Cooray & Hu, 2001) etc.

Here, for the purpose of validation, we construct a simplified toy model for cosmic shear data where we can validate the likelihood-free results against an exact (known) likelihood.

7.1 Tomographic cosmic shear data and model

The tomographic angular shear power spectra contain much of the cosmological information contained in the weak lensing distortion fields and are straightforwardly predicted from theory. For a given set of cosmological parameters $\boldsymbol{\theta}$, the predicted angular power spectra between tomographic bins $\alpha$ and $\beta$ are given by666In the Limber approximation, (Limber, 1954). (Kaiser, 1992, 1998; Hu, 1999, 2002; Takada & Jain, 2004; Kitching et al., 2017),

[TABLE]

where cosmological parameters $\boldsymbol{\theta}$ enter via the comoving distance-redshift relation $\chi(z)$, the 3D power spectrum of matter fluctuations $P_{\delta}(k;\chi)$, and the lensing weight functions which are given by

[TABLE]

where $n_{\alpha}(\chi)d\chi=p_{\alpha}(z)dz$ is the redshift distribution for galaxies in redshift bin $\alpha$.

The tomographic redshift distributions $p_{\alpha}(z)$ must be estimated from photometric data and are subject to a range of systematic biases that typically must be parameterized and marginalized over (see eg., Salvato et al., 2018 for a review). For this demonstration we parameterize systematic uncertainties in the tomographic redshift distributions as an unknown shift of the estimated distributions, ie., setting

[TABLE]

where $\hat{p}_{\alpha}(z)$ is the estimated redshift distribution for tomographic bin $\alpha$ and the systematic shift parameters $b_{\alpha}$ – one per tomographic bin – constitute the nuisance parameters $\boldsymbol{\eta}$ that we ultimately want to marginalize over.

For this demonstration we will assume the data vector to be a set of estimated (noise biased) band powers,

[TABLE]

We consider a flat-$w$CDM cosmology parameterized by $\boldsymbol{\theta}=(\sigma_{8},\Omega_{\mathrm{m}},\Omega_{\mathrm{b}}h^{2},h,n_{s},w_{0})$.

7.2 Likelihood and data compression

For this validation case study, we assume that the data (band powers) are independently Wishart distributed,

[TABLE]

where the elements of the expected band powers are given by

[TABLE]

$N_{\ell,ij}$ is the shape noise power spectrum (from intrinsic random galaxy ellipticities) and $\nu_{\mathcal{B}_{k}}$ is the total number of modes contributing to band $k$. We will assume isotropic shape noise $N_{\ell,ij}=\sigma_{e}^{2}/\bar{n}_{i}\delta_{ij}$ with intrinsic ellipticity variance $\sigma_{e}^{2}$ and mean galaxy number density $\bar{n}_{i}$ in tomographic bin $i$. We will assume $\nu_{\mathcal{B}_{k}}=f_{\mathrm{sky}}\sum_{\ell\in\mathcal{B}_{k}}(2\ell+1)$ modes per band to approximately mimick partial sky-fraction coverage $f_{\mathrm{sky}}$.

We take our compressed summaries as the score of the Wishart likelihood in Eq. (28),

[TABLE]

where ‘$*$’ indicates evaluation about fiducial parameters $\boldsymbol{\theta}_{*}$ which we take to be $\boldsymbol{\theta}_{*}=(0.8,0.3,0.05,0.7,0.96,-1)$.

Projection of the nuisance parameters is performed following Eq. (13) as usual, giving nuisance hardened summary statistics:

[TABLE]

where the Fisher matrix in this case is given by,

[TABLE]

For a full sophisticated implementation of likelihood-free inference to cosmic shear on real data, an approximate likelihood or information maximizing neural network may be used for performing data compression.

7.3 Simulations

In this simple validation case we simulate data (given parameters) by drawing band powers as Wishart random variates, following Eq. (28).

We assume a survey set-up similar to the upcoming ESA Euclid survey: $15,000$ square degrees with a mean galaxy number density of $\bar{n}=30\,\mathrm{arcmin}^{-2}$, an overall galaxy redshift distribution $n(z)\propto z^{2}\mathrm{exp}\left[-(1.41z/z_{m})^{1.5}\right]$ with a median redshift $z_{m}=0.9$, Gaussian photo-$z$ uncertainties with standard deviation $\sigma_{z}=0.05*(1+z)$, and ten tomographic bins with equal mean galaxy number density per bin. Modes are binned into ten logarithmically spaced bands between $\ell=10$ and $\ell=3000$.

The data for this demonstration are simulated following Eq. (28) assuming values for the cosmological parameters (based on Aghanim et al., 2018) $\sigma_{8}=0.811$, $\Omega_{\mathrm{m}}=0.315$, $\Omega_{\mathrm{b}}h^{2}=0.0224$, $h=0.674$, $n_{s}=0.965$, $w_{0}=-1.03$, and systematics parameters $b_{i}=0\;\forall i$.

7.4 Priors

For the interesting parameters $\boldsymbol{\theta}=(\sigma_{8},\Omega_{\mathrm{m}},\Omega_{\mathrm{b}}h^{2},h,n_{s},w_{0})$ we assume Gaussian priors with means $\boldsymbol{\mu}_{\boldsymbol{\theta}}=(0.8,0.3,0.0224,0.674,0.96,-1.0)$, standard deviations $\boldsymbol{\sigma}_{\boldsymbol{\theta}}=(0.1,0.1,0.00015,0.005,0.3,0.3)$, and hard parameter limits $\sigma_{8}\in\left[0.4,1.2\right]$, $\Omega_{\mathrm{m}}\in\left[0,1\right]$, $\Omega_{\mathrm{b}}h^{2}\in\left[0,0.1\right]$, $h\in\left[0.4,1\right]$, $n_{\mathrm{s}}\in\left[0.7,1.3\right]$. The tight priors on $\Omega_{\mathrm{b}}h^{2}$ and $h$ (which are poorly constrained by weak lensing data) are taken from (Aghanim et al., 2018).

For the nuisance parameters $\{b_{i}\}$ we assume independent zero-mean Gaussian priors with standard deviation $0.05$ and hard limits $b_{i}\in\left[-0.1,0.1\right]$.

7.5 DELFI set-up

As in §6, the DELFI algorithm learns the sampling distribution of the data summaries as a function of the parameters. As before we parameterize the sampling distribution by a mixture density network with three Gaussian components, and simulations are ran in sequential batches with an adaptive proposal density (using the public code pydelfi Alsing et al., 2019; see also Appendix A).

7.6 Results

We ran the cosmic shear inference task described above, inferring the six interesting cosmological parameters, implicitly marginalizing over the ten additional nuisances describing the photo-$z$ systematics.

Figure 4 shows the convergence of the negative log-likelihood of the simulated parameter-data summary samples under the learned neural density estimator, as a function of the number of forward simulations. We see that convergence is achieved after $\mathcal{O}(10^{3})$ forward simulations, irrespective of the presence of ten additional nuisance parameters that are implicitly marginalized over here.

Figure 5 shows the recovered nuisance marginalized posterior for the six cosmological parameters using likelihood-free inference (red) versus MCMC sampling of the exact nuisance plus interesting parameter problem and marginalizing a posteriori. We see that again in this case, the nuisance hardened data compression has successfully retained essentially all of the information content of the data with respect to the interesting parameters, yielding a high-fidelity nuisance marginalized (likelihood-free) posterior inference.

8 Conclusions

We have shown that nuisance marginalized posteriors can be inferred directly in the likelihood-free paradigm, massively reducing the number of simulations required to perform likelihood-free inference in the presence of nuisance parameters. This opens up likelihood-free methods to problems with expensive forward simulations, a small number of parameters of interest, but a large number of additional nuisance parameters; a common scenario in cosmological data analysis.

To achieve fast direct inference of nuisance marginalized posteriors, we do the following: first, we showed that nuisance parameters can easily be re-cast as local latent variables and hence implicitly marginalized over in the likelihood-free inference framework (§4). Secondly, we derived a scheme for compressing $N$ data down to $n$ informative summaries – one per parameter of interest – such that the compressed summaries (asymptotically) retain the Fisher information content of the data for the interesting parameters whilst being “hardened" to the nuisance parameters. This means that just $n$ data summaries can be used for likelihood-free inference, regardless of the number of nuisance parameters.

The result is that the complexity of the inference task – and hence the number of simulations required for performing inference – is dictated only by the number of interesting parameters considered in the problem.

Acknowledgements

This work is supported by the Simons Foundation. Justin Alsing was partially supported by the research project grant “Fundamental Physics from Cosmological Surveys" funded by the Swedish Research Council (VR) under Dnr 2017-04212. Benjamin Wandelt acknowledges support by the Labex Institut Lagrange de Paris (ILP) (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme Investissements d’avenir under the reference ANR-11-IDEX-0004-02.

Appendix A DELFI algorithm set-up

We parameterize the sampling distribution of the data summaries $p(\mathbf{t}|\boldsymbol{\theta})$ as a mixture density network. Simulations are run in batches and the mixture network is re-trained after each new batch of simulations. In the first cycle, parameters for performing forward simulations are drawn from a broad Gaussian centered on the prior mean and with covariance equal to nine times the (estimated) Fisher matrix. For subsequent batches of simulations, parameters are drawn from the geometric mean of the current posterior estimate from the mixture network and the prior777This proposal density is inspired by the optimal parameter proposal for Approximate Bayesian Computation (Alsing et al., 2018a)..

The neural networks are trained using the stochastic gradient optimizer adam (Kingma & Ba, 2014), with a batch-size of one tenth of the training set during each re-training cycle and a learning rate of $0.001$. Over-fitting is prevented during training using early-stopping; during each re-training cycle, 10% of the training set is set aside for validation and training is terminated when the validation-loss does not improve after $20$ epochs.

We set the size of the mixture network architecture and number of new simulations per training cycle depending on the number of parameters to be inferred. We use the following rule of thumb for scaling the architecture and simulation batches: all of the networks have two hidden layers with $5\times$ as many units as parameters to be inferred, and simulations are run in batches of $50\times$ the number of parameters in the problem. Of course, this is a very rough heuristic and careful application of neural density estimators should involve a cross-validation search over network architectures.