# On the accuracy and precision of correlation functions and field-level inference in cosmology

Florent Leclercq, Alan Heavens

<figure><img src="image_1.png" /></figure>

## To cite this version:

Florent Leclercq, Alan Heavens. On the accuracy and precision of correlation functions and field-level inference in cosmology. Monthly Notices of the Royal Astronomical Society: Letters, 2021, 506 (1), pp.L85-L90. 10.1093/mnrasl/slab081. hal-05347495

HAL Id: hal-05347495

https://hal.science/hal-05347495v1

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# On the accuracy and precision of correlation functions and field-level inference in cosmology

Florent Leclercq* and Alan Heavens

Imperial Centre for Inference and Cosmology (ICIC) & Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK

Accepted 2021 July 12. Received 2021 July 12; in original form 2021 March 12

## ABSTRACT

We present a comparative study of the accuracy and precision of correlation function methods and full-field inference in cosmological data analysis. To do so, we examine a Bayesian hierarchical model that predicts lognormal (LN) fields and their two-point correlation function. Although a simplified analytic model, the LN model produces fields that share many of the essential characteristics of the present-day non-Gaussian cosmological density fields. We use three different statistical techniques: (i) a standard likelihood-based analysis of the two-point correlation function; (ii) a likelihood-free (simulation-based) analysis of the two-point correlation function; (iii) a field-level analysis, made possible by the more sophisticated data assimilation technique. We find that (a) standard assumptions made to write down a likelihood for correlation functions can cause significant biases, a problem that is alleviated with simulation-based inference; and (b) analysing the entire field offers considerable advantages over correlation functions, through higher accuracy, higher precision, or both. The gains depend on the degree of non-Gaussianity, but in all cases, including for weak non-Gaussianity, the advantage of analysing the full field is substantial.

**Key words:** methods: statistical – cosmological parameters – large-scale structure of Universe.

## 1 INTRODUCTION

As cosmological surveys increase in size, the statistical errors decrease and the possibility for new discoveries increases. This increase in precision needs to be accompanied by equally high accuracy in order to avoid the wrong inferences being drawn. With the imminent arrival of the Euclid satellite (Laureijs et al. 2011) and the Legacy Survey of Space and Time (LSST Science Collaboration 2012), the time is right to assess whether the standard methods and more advanced tools for statistical analysis of late-time, non-Gaussian, cosmological fields are likely to be adequate. In this letter, we compare methods by analysing an analytic model, namely the lognormal (LN) model, that has many of the salient features of the late-time fields such as the galaxy density field or the weak lensing cosmic shear field.<fnref n="1" />

We focus on three analysis methods: (i) the most common standard technique, i.e. a likelihood-based analysis (LBA) of the two-point correlation function (2PCF), assuming a Gaussian distribution with fixed covariance matrix; (ii) the relatively new technique (to cosmology) of simulation-based inference (SBI), also known as likelihood-free inference, based on the 2PCF; (iii) and the computationally expensive but powerful data assimilation (DA) technique, which allows field-level inference without any compression of the data. We investigate fields of varying levels of non-Gaussianity, and find

that field-level inference offers advantages in accuracy and precision over the 2PCF that can be extremely large. For significantly non-Gaussian fields, the precision of LBA and SBI are similar, but SBI is typically more accurate (less ‘biased’). DA, which uses all the field values and captures all of the information, outperforms both and provides highly accurate and precise results. We also study weakly non-Gaussian fields, where one might expect the precision of 2PCF analyses to equal that of field-level inference, as the common lore is that the statistical properties in the limit of vanishing non-Gaussianity are completely captured by the 2PCF. This turns out not to be the case, and, although we find comparable statistical errors in the shape parameter of field correlations, DA gives much more precise measurements of the degree of non-Gaussianity.

## 2 MODEL

We consider discrete two-dimensional LN random fields (Coles & Jones 1991) of size $N_{pix}$. Such fields can be seen as variables of a Bayesian hierarchical model (BHM) represented in Fig. 1. Specifically, the starting point is a Gaussian random field (GRF) g (see Peacock & Heavens 1985; Bardeen et al. 1986) of size $N_{pix}$, with zero mean and a covariance matrix $\Xi_G$. We further assume that the 2PCF $\xi_G(r)$, giving the coefficients of $\Xi_G$, is parametrized by a single variable $\beta$, such that

$$ \mathcal{A}(g|\beta) = G[0, \Xi_G] \quad \text{with} \quad \xi_G(r) = \exp\left(-\frac{1}{4}\frac{r^2}{\beta^2}\right), \quad (1) $$

<footnote>* E-mail: florent.leclercq@polytechnique.org</footnote>

<footnote n="1">¹This letter is accompanied by supplementary appendices, which provide supporting information on several aspects of our model and study.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Graphical representation of the Bayesian hierarchical model used in this work. The rounded yellow boxes represent probability distributions and the purple squares represent deterministic functions. Three different statistical methods, as indicated by dashed rectangles, are considered.</figcaption></figure>

where *r* is the separation between two grid points, in units of pixels. The GRF *g* is made LN via the following transform:

$$ \mathrm{LN}(\mathbf{g}, \alpha) = \frac{1}{\alpha} \left[ \exp \left( \alpha \mathbf{g} - \frac{1}{2} \alpha^2 \right) - 1 \right], \qquad (2) $$

which involves one free parameter *α*. The final field *f* is assumed to be a noisy realisation, the expectation value being given by equation (2). For simplicity, we assume a Gaussian noise model, uniform in real space, with standard deviation *σ*. Therefore,

$$ \mathcal{P}(\mathbf{f}|\mathbf{g}, \alpha) = \mathcal{G}[\mathrm{LN}(\mathbf{g}, \alpha), \mathrm{diag}(\sigma^2)]. \qquad (3) $$

We note that in equation (2), the term $-\frac{1}{2}\alpha^2$ in the exponential ensures that the mean of LN($\mathbf{g}, \alpha$) is zero,<fnref n="2" /> and the prefactor $\frac{1}{\alpha}$ ensures that the signal to noise ratio in *f* is independent of *α*, in the limit $\alpha \to 0$.

Remarkably, in the absence of noise, LN fields have analytic correlation functions at all orders (Coles & Jones 1991). In particular, their 2PCF is given by

$$ \xi_{\mathrm{LN}}(r) \equiv \frac{1}{\alpha^2} \left\{ \exp \left[ \alpha^2 \exp \left( -\frac{1}{4} \frac{r^2}{\beta^2} \right) \right] - 1 \right\}, \qquad (4) $$

whose gradients with respect to *α* and *β* are readily obtained (see appendix A of the supplementary material).

The equations above can be interpreted as a simple model for observations of the matter distribution in the Universe. As can be observed in equation (1), *β* captures the strength of correlations in *g*, and can therefore be seen as a proxy for the initial matter power spectrum. From equation (2), it appears that *α* characterizes the strength of non-Gaussianities in the final field, and therefore it can been seen as a proxy for the non-linear evolution of the matter field. Note that the signal is arbitrarily close to the original Gaussian

random field *g* at early times, i.e. in the limit $\alpha \to 0$. All other physical and observational processes are modelled stochastically by the additive Gaussian noise.

Throughout this paper, we adopt $N_{\mathrm{pix}} = 20^2$, $\sigma = 0.01$, and uniform (hyper-)priors $\mathcal{P}(\alpha)$ and $\mathcal{P}(\beta)$ on the intervals [0,2] and [0.2,0.8], respectively. We aim at inferring $\alpha$ and $\beta$ from observations of one or several realizations of *f*, using different techniques.

### 3 LIKELIHOOD-BASED ANALYSIS OF THE 2PCF

The first technique we consider is usual in cosmology, and consists in a LBA of the estimated 2PCF $\xi^{\mathrm{obs}}$. Given a realization $f^{\mathrm{obs}}$ of the model, an estimator of the 2PCF is constructed by taking the outer product of $f^{\mathrm{obs}}$ and itself, and then averaging these two-point correlations in different *r*-bins. This is a deterministic process, denoted *S*, that compresses any field *f* into a new data vector $\xi$ of size $N_{\xi}$, where $N_{\xi}$ is the number of *r*-bins. Formally, this can be seen as adding a layer at the bottom of the BHM, shown in green in Fig. 1.

Without any prior knowledge of the sampling distribution $\mathcal{P}(\xi|\alpha, \beta)$, we make the common assumption that it is a multivariate Gaussian distribution,

$$ \mathcal{P}(\xi|\alpha, \beta, \mathbf{\Sigma}) = \mathcal{G}[\xi_{\mathrm{LN}}(\alpha, \beta), \mathbf{\Sigma}], \qquad (5) $$

where the expectation value is given by the analytic form of equation (4) and the covariance matrix is denoted $\mathbf{\Sigma}$. With the LN model, it would be possible to compute $\mathbf{\Sigma}$ analytically using the expression of the four-point function. However, in this paper we adopt the more general scenario where $\mathbf{\Sigma}$ is unknown and can only be evaluated through the use of simulations. We further make the common assumption that $\mathbf{\Sigma}$ does not depend on $\alpha$ and $\beta$. Given $N_{\mathrm{sims}}$ realizations of $\xi$, denoted $\xi_i$ ($1 \le i \le N_{\mathrm{sims}}$), $\mathbf{\Sigma}$ is evaluated using the usual unbiased estimator

$$ \hat{\mathbf{\Sigma}} \equiv \frac{1}{N_{\mathrm{sims}} - 1} \sum_{i=1}^{N_{\mathrm{sims}}} (\xi_i - \bar{\xi})(\xi_i - \bar{\xi})^{\mathrm{T}}, \\ \text{where } \bar{\xi} \equiv \frac{1}{N_{\mathrm{sims}}} \sum_{i=1}^{N_{\mathrm{sims}}} \xi_i. \qquad (6) $$

When estimating parameters from Gaussian-distributed data in situations where the covariance matrix is unknown and estimated from simulations, the Gaussian likelihood should be replaced by a multivariate *t*-distribution (Sellentin & Heavens 2016). Thereby, the final likelihood that we use satisfies

$$ \mathcal{P}(\xi|\alpha, \beta, \hat{\mathbf{\Sigma}}) \propto \left[ 1 + \frac{(\xi - \xi_{\mathrm{LN}})^{\mathrm{T}} \hat{\mathbf{\Sigma}}^{-1} (\xi - \xi_{\mathrm{LN}})}{N_{\mathrm{sims}} - 1} \right]^{-\frac{N_{\mathrm{sims}}}{2}}. \qquad (7) $$

With the estimated 2PCF $\xi^{\mathrm{obs}}$ and uniform priors, the LBA posterior is thus $\mathcal{P}_{\mathrm{LBA}}(\alpha, \beta|\xi^{\mathrm{obs}}) \propto \mathcal{P}(\xi^{\mathrm{obs}}|\alpha, \beta, \hat{\mathbf{\Sigma}})$.

### 4 SIMULATION-BASED INFERENCE USING THE 2PCF

The second technique we consider is SBI, that is to say, a likelihood-free analysis of the estimated 2PCF $\xi^{\mathrm{obs}}$. Since we aim at comparing different inference techniques, the model is the same as in Sections 2 and 3 for the generation of LN fields and the estimation of $\xi$. Only the statistical assumptions differ.

SBI is known to be difficult when the size of the data vector is large. To solve this issue, we add a compression step, denoted by *C*.

<footnote n="2">²Normally this term is multiplied by the variance of the Gaussian field, which is unity by construction.</footnote>
Resulting compressed variables correspond to an additional layer at the bottom of the BHM, represented in blue in Fig. 1. In this work, we use a compression using the score function, a generalization of the MOPED algorithm (Heavens, Jimenez & Lahav 2000; Alsing & Wandelt 2018). Specifically, to compress the data, we assume that the sampling distribution of $\xi$ is Gaussian-distributed, i.e.

$$ \mathcal{L} \equiv \ln \mathcal{P}(\xi|\alpha, \beta, \mathbf{\Sigma}) = -\frac{1}{2}(\xi - \xi_{\text{LN}})^{\top}\mathbf{\Sigma}^{-1}(\xi - \xi_{\text{LN}}) - \frac{1}{2}\ln|\mathbf{\Sigma}|. \quad (8) $$

As we further assume that $\mathbf{\Sigma}$ is independent of parameters ($\nabla \mathbf{\Sigma} = 0$), the Fisher information matrix takes a simple form:

$$ \mathbf{F} = -E[\nabla\nabla\mathcal{L}] = \nabla\xi_{\text{LN}}^{\top}\mathbf{\Sigma}^{-1}\nabla\xi_{\text{LN}}. \quad (9) $$

Using a fiducial point $\theta_* = (\alpha_*, \beta_*)$ in parameter space, a quasi maximum-likelihood estimator for the parameters is $\tilde{\theta} = \theta_* + \mathbf{F}_*^{-1}\nabla\mathcal{L}_*$, where the inverse of the Fisher matrix and the gradient of the log-likelihood are evaluated at the fiducial point. In our case, the expression is therefore (Alsing & Wandelt 2018)

$$ \tilde{\theta} = \theta_* + \mathbf{F}_*^{-1} [\nabla\xi_{\text{LN}*}^{\top} \hat{\mathbf{\Sigma}}^{-1} (\xi - \xi_{\text{LN}*})], \quad (10) $$

with

$$ \mathbf{F}_* = \nabla\xi_{\text{LN}*}^{\top} \hat{\mathbf{\Sigma}}^{-1} \nabla\xi_{\text{LN}*}. \quad (11) $$

The quasi maximum-likelihood estimator $\tilde{\theta}$ compresses a data vector of arbitrary size $p$ summaries, where $p$ is the number of target parameters of the problem (two in our case). It is an optimal compression, in the sense that it saturates the lower bound of the Fisher information inequality. Compressing simulated data $\xi$ via equations (10) and (11) requires only the estimated covariance matrix $\hat{\mathbf{\Sigma}}$ already used for the LBA [equation (6)], and the gradients of the expectation value $\nabla\xi_{\text{LN}}^{\top}$, for which we use the analytic expressions [equations (A.1) and (A.2)]. It yields two numbers that we note $\tilde{\alpha}$ and $\tilde{\beta}$. Applying the same compressor $\mathcal{C}$ to $\xi^{\text{obs}}$ yields $\tilde{\alpha}^{\text{obs}}$ and $\tilde{\beta}^{\text{obs}}$. We note that, in this section, the Gaussian assumption of equation (8) is used only for compression and not for subsequent inference. In SBI, compression can be performed under an approximate likelihood without introducing biases in the result.

By definition, SBI shall only rely on forward evaluations of the data model. In the last few years, sophisticated algorithms have been developed to drastically reduce the number of simulations needed and/or to scale to high dimension (Leclercq 2018; Alsing, Wandelt & Feeney 2018; Alsing et al. 2019; Leclercq et al. 2019). Given that we are only interested in a two-parameter problem and that our simulations are computationally very cheap, in this paper we rely on the simplest solution for SBI, namely likelihood-free rejection sampling (sometimes also known as Approximate Bayesian Computation). To do so, we replace the parametric assumption for the likelihood $\mathcal{P}(\xi|\alpha, \beta)$ [equation (7)] by a measurement of the discrepancy between simulated $\{\tilde{\alpha}, \tilde{\beta}\}$ and observed $\{\tilde{\alpha}^{\text{obs}}, \tilde{\beta}^{\text{obs}}\}$ compressed data. We choose the Euclidean distance,

$$ \Delta(\tilde{\alpha}, \tilde{\beta}) = \sqrt{(\tilde{\alpha} - \tilde{\alpha}^{\text{obs}})^2 + (\tilde{\beta} - \tilde{\beta}^{\text{obs}})^2}. \quad (12) $$

Given a small threshold $\epsilon$, the algorithm for likelihood-free rejection sampling is then to iterate many times the following procedure: draw $(\alpha, \beta)$ from a proposal distribution (in our case, the uniform priors), simulate $\tilde{\alpha}$ and $\tilde{\beta}$ using the full BHM represented in Fig. 1, compute the discrepancy $\Delta(\tilde{\alpha}, \tilde{\beta})$, then accept $(\alpha, \beta)$ as a sample of the approximate posterior if $\Delta(\tilde{\alpha}, \tilde{\beta}) \le \epsilon$, and reject it otherwise. Therefore, given $N_{\text{sbi}}$ such tries with discrepancies $\Delta_j(\tilde{\alpha}, \tilde{\beta}) (1 \le j$

$\le N_{\text{sbi}}$), the resulting SBI posterior satisfies

$$ \mathcal{P}_{\text{SBI}}(\alpha, \beta | \tilde{\alpha}^{\text{obs}}, \tilde{\beta}^{\text{obs}}, \epsilon) \propto \sum_{j=1}^{N_{\text{sbi}}} \mathbf{I}_{[0, \epsilon]} [\Delta_j(\tilde{\alpha}, \tilde{\beta})], \quad (13) $$

where $\mathbf{I}_{[0, \epsilon]}$ is the indicator function of the interval $[0, \epsilon]$.

# 5 FIELD-LEVEL INFERENCE VIA DATA ASSIMILATION

The third and last technique we consider is a field-level analysis, where we infer $\alpha$ and $\beta$ from the full observed field $f^{\text{obs}}$, rather than its estimated 2PCF. Field-level inference uses all of the data to hand. This is done by sampling all the variables appearing above $f$ in the BHM (i.e. $\alpha, \beta$ and $g$), conditional on $f^{\text{obs}}$. This means that the $N_{\text{pix}}$ variables of $g$, which are latent variables of the problem in Sections 3 and 4, now become target parameters of the problem. The full joint posterior satisfies

$$ \begin{aligned} \mathcal{P}(\alpha, \beta, g | f^{\text{obs}}) &\propto \mathcal{P}(f^{\text{obs}} | g, \alpha, \beta) \mathcal{P}(g, \alpha, \beta) \\ &= \mathcal{P}(f^{\text{obs}} | g, \alpha) \mathcal{P}(g | \beta) \mathcal{P}(\beta) \mathcal{P}(\alpha). \end{aligned} \quad (14) $$

The first term is given by equation (3), the second by equation (1), and the third and fourth are the uniform priors. Only the probabilistic equations of the LN model (those given in Section 2) appear here, without any additional compression or statistical assumption. Therefore, from a mathematical point of view, field-level inference provides the exact solution to the inference problem.

Given its high dimensionality, sampling from the posterior written down in equation (14) requires advanced statistical techniques such as Hamiltonian Monte Carlo (Duane et al. 1987) or more recent variants such as a No-U Turn Sampler (NUTS, Hoffman & Gelman 2014). We refer to such techniques as DA techniques, as they permit assimilation of the full observed field values $f^{\text{obs}}$ into the field-level LN model [equations (1)–(3)]. Once samples of $\mathcal{P}(\alpha, \beta, g | f^{\text{obs}})$ are obtained, the DA posterior on $\alpha$ and $\beta$ is simply given by marginalizing over the field values $g$,

$$ \mathcal{P}_{\text{DA}}(\alpha, \beta | f^{\text{obs}}) = \int \mathcal{P}(\alpha, \beta, g | f^{\text{obs}}) dg, \quad (15) $$

which is trivially obtained from the sampled values of $\alpha$ and $\beta$.

# 6 NUMERICAL RESULTS

In this section, we present some numerical results of our comparative study.

We first analyse four different LN field realizations with ground truth values $\alpha = 1$ and $\beta = 0.5$ (see appendix B). Corresponding posteriors on $\alpha$ and $\beta$, obtained with the three different methods, are shown in Fig. 2. For the LBA of $\xi$, the posteriors $\mathcal{P}_{\text{LBA}}(\alpha, \beta | \xi^{\text{obs}})$ [equation (7)] are evaluated on a regular grid of 80 × 80 points and shown by the green contours. We used $N_{\text{sims}} = 500$ simulations run at the ground truth values of $\alpha$ and $\beta$, well above the number of non-empty r-bins ($N_{\xi} = 68$), so as to get an accurate estimate of $\mathbf{\Sigma}$ (see appendix C). For SBI, we use the ground truth as fiducial point $(\alpha_* = 1, \beta_* = 0.5)$ to compress $\xi$. We generated a large set of $N_{\text{sbi}} \gtrsim 600 000$ simulations. Such a pool size is sufficient to ensure that at least 3000 samples are accepted when doing likelihood-free rejection sampling with $\epsilon = 0.05$, for each of the four fields. The value of $\epsilon = 0.05$ was chosen as a compromise between reducing the marginal variance of $\alpha$ and $\beta$ (by decreasing $\epsilon$), and keeping a sensible acceptance rate (which goes to zero as $\epsilon \to 0$, see appendix
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of different analysis techniques for four different LN fields. The posteriors are obtained using a LBA of the 2PCF [equation (7), green], SBI using the 2PCF [equation (13), blue], and DA at the field level [equation (15), orange]. 1-, 2- and 3σ contours are shown. Ground truth values are $α = 1$ and $β = 0.5$.</figcaption></figure>

D). The SBI posteriors $\mathcal{P}_{\text{SBI}}(\alpha, \beta | \tilde{\alpha}^{\text{obs}}, \tilde{\beta}^{\text{obs}}, \varepsilon)$ [equation (13)] are shown in blue. For DA, we choose a NUTS initialized with automatic differentiation variational inference (ADVI, Kucukelbir et al. 2017). During the initialization phase, ADVI automatically determines a variational family for the posterior and optimizes the variational objective. Quantities required for NUTS, such as gradients of the data model with respect to the parameters $(\alpha, \beta, g)$ and the mass matrix are also computed automatically. For each of the four fields $f^{\text{obs}}$, we run 16 independent chains. We determine the length of the residual burn-in phase (usually a few thousand samples), after ADVI initialization, by examining trace plots of the parameters. We drop this burn-in phase, and obtain samples of the high-density regions of the posteriors. We enforce sufficient convergence by making sure that, for each of the 402 parameters $(\alpha, \beta$ and the $20^2$ field values of $g$), the Gelman–Rubin statistic $\hat{R}$ is smaller than 1.05. After marginalizing over $g$, the final posteriors $\mathcal{P}_{\text{DA}}(\alpha, \beta | f^{\text{obs}})$ [equation (15)] are shown in orange in Fig. 2. We checked that the marginal standard deviations of $\alpha$ and $\beta$, obtained via LBA, SBI, and DA, are consistent with Cramér–Rao bounds, which can be computed analytically for a 2PCF analysis (see appendix E) and for a noise-free field-level analysis (see appendix F).

Since the four fields analysed in Fig. 2 are independent, the final posterior for each method is obtained by multiplying the four likelihoods, and is shown in Fig. 3. From Figs 2 and 3, it appears that the incorrect assumptions made in Section 3 for the LBA of the 2PCF can cause strong biases, particularly when the field contains outlying values (e.g. field 4, see Fig. B.1 of the supplementary material). Combining different fields does not entirely solve the issue. On the contrary, SBI using the same 2PCF provides accurate results, at a comparable level of precision. As the exact mathematical solution to the problem, field-level inference via DA achieves unbiased and far more precise results. The value of $α = 1$ is rather large, so one might expect the 2PCF to lose information, but the relative performance of the field-level inference is strikingly good.

In the limit $α \to 0$, the field becomes Gaussian, so one would expect the 2PCF (closer to a sufficient statistical summary of the field), to fare better. We investigated this hypothesis for a weakly non-Gaussian field with $α = 0.2$ (field 5 in Fig. B.1), repeating the same treatment as described above. The results are shown in Fig. 4.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** The posterior from the set of four independent LN fields. Using the same information (the 2PCF), SBI (blue) gives more accurate results than a LBA (green), which presents a bias of ~2σ. Using all the data to hand, field-level inference achieves unbiased and far more precise results.</figcaption></figure>

In this regime, the 2PCF does capture well the spatial correlations of the field, as demonstrated by a comparable marginal variance on $β$ for the three methods. However, field-level inference remains the only technique that gives precise results for the non-Gaussianity parameter $α$. This result can be interpreted from the model: as the LN transform [equation (2)] is a pixel-by-pixel operation, $α$ is essentially one-point information, which is not captured well by the 2PCF, which is constant up to small corrections of order $α^2$ when $α \to 0$, as is readily seen from equation (4). Consistently, from equations (A.1)
<figure><img src="image_7.png" /><figcaption>**Figure 4.** As in Fig. 2, except for a single near-Gaussian field with $\alpha = 0.2$. Note that even for this almost Gaussian field, the 2PCF does not capture all the information; both a LBA and SBI fare well for the parameter $\beta$, but field-level inference is the only technique that gives precise and accurate results for $\alpha$.</figcaption></figure>

and (A.2), it appears that $\partial\xi_{LN}/\partial\alpha$ goes more rapidly to zero than $\partial\xi_{LN}/\partial\beta$ as $r$ increases.

# 7 DISCUSSION AND SUMMARY

In this letter, we have compared the posteriors obtained from applying three data analysis techniques to LN fields. The LN field is occasionally used as an approximation to the late-time cosmological density field since its one-point distribution agrees approximately with that of the evolved matter density field, and by choice of the Gaussian two-point function, it can be made to agree at two-point level as well. The non-Gaussianity parameter $\alpha$ is a proxy for non-linear growth. The three techniques we have investigated are (i) LBA: the standard technique of using correlation functions, with an assumption that these summary statistics are Gaussian-distributed; (ii) SBI, or likelihood-free inference, where (optimally) compressed arbitrary summary statistics are used without statistical assumptions about their distribution; and (iii) field-level inference, where the BHM is inferred in a mathematically exact way via DA, and all the data are used. We find that for fields with significant non-Gaussianity, field-level inference gives far more precise and accurate posteriors than the standard LBA method, which can sometimes give highly inaccurate answers dependent on the data realization. It is also far more precise than SBI using the 2PCF. Even for nearly Gaussian fields, field-level inference outperforms both in precision.

We conclude that field-level inference should always be preferred to correlation functions when it is possible, particularly for non-Gaussian fields for which it gives the largest gain in precision. With specialized software, field-level inference with tens of millions of parameters is possible (e.g. Lavaux, Jasche & Leclercq 2019 for galaxy clustering, Porqueres et al. 2021 for weak lensing). When field-level inference is not possible for computational considerations, SBI gives more accurate answers than the standard LBA, which must arise from a breakdown of the Gaussian likelihood assumption (see appendix G, where we investigate the influence of our setup parameters on the accuracy of LBA).

The standard approach in weak lensing uses summary statistics assuming a Gaussian likelihood function, but this assumption is not strictly correct (Sellentin & Heavens 2018). For correlation functions, Hartlap et al. (2009) and Sellentin, Heymans & Harnois-Déraps (2018) claim significant errors in posterior values. Lin et al. (2020) used principal component analysis to decorrelate the correlation function, showing small biases, but making a further assumption of independence. For pseudo-power spectra, a Gaussian likelihood may be accurate enough (Taylor et al. 2019; Upham, Brown & Whittaker 2021), but the method needs the covariance matrix, which may be hard to calculate because of non-Gaussianities and super-sample covariance. Generally, likelihood-based analyses suffer from the difficulty of supplying an accurate covariance matrix, for which the parameter dependence, non-Gaussianity of the underlying field, and inclusion of super-sample covariance are challenging.

We note that there are prospects for improving the precision of SBI results by including more summary statistics, such as the field skewness, its three-point function, or summaries defined via machine learning. This is also the case for likelihood-based inference, but the problems of supplying an accurate covariance matrix, and the assumption of a Gaussian likelihood, still remain.

# ACKNOWLEDGEMENTS

We thank Andrew Jaffe, Guilhem Lavaux, and Benjamin Wandelt for useful discussions. This work made use of the ELFI (Lintusaari et al. 2018), pyDELFI (Alsing et al. 2019), and pyMC3 (Salvatier, Wiecki & Fonnesbeck 2016) codes. This work was done within the Aquila Consortium.

# DATA AVAILABILITY

The code and data underlying this letter are publicly available on Github, at https://github.com/florent-leclercq/correlations_versus-field.

# REFERENCES

* Alsing J., Wandelt B., 2018, MNRAS, 476, L60

* Alsing J., Wandelt B., Feeney S., 2018, MNRAS, 477, 2874

* Alsing J., Charnock T., Feeney S., Wandelt B., 2019, MNRAS, 488, 4440

* Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15

* Coles P., Jones B., 1991, MNRAS, 248, 1

* Duane S., Kennedy A. D., Pendleton B. J., Roweth D., 1987, Phys. Lett. B, 195, 216

* Hartlap J., Schrabback T., Simon P., Schneider P., 2009, A&A, 504, 689

* Heavens A. F., Jimenez R., Lahav O., 2000, MNRAS, 317, 965

* Hoffman M. D., Gelman A., 2014, J. Mach. Learn. Res., 15, 1593

* Kucukelbir A., Tran D., Ranganath R., Gelman A., Blei D. M., 2017, J. Mach. Learn. Res., 18, 1

* Laureijs R. et al., 2011, Euclid Definition Study Report. preprint (arXiv:1110.3193)

* Lavaux G., Jasche J., Leclercq F., 2019, preprint (arXiv:1909.06396)

* Leclercq F., 2018, Phys. Rev. D, 98, 063511

* Leclercq F., Enzi W., Jasche J., Heavens A., 2019, MNRAS, 490, 4237

* Lin C.-H., Harnois-Déraps J., Eifler T., Pospisil T., Mandelbaum R., Lee A. B., Singh S., LSST Dark Energy Science Collaboration, 2020, MNRAS, 499, 2977

* Lintusaari J. et al., 2018, J. Mach. Learn. Res., 19, 1

* LSST Science Collaboration, 2012, Large Synoptic Survey Telescope: Dark Energy Science Collaboration. preprint (arXiv:1211.0310)

* Peacock J. A., Heavens A. F., 1985, MNRAS, 217, 805

* Porqueres N., Heavens A., Mortlock D., Lavaux G., 2021, MNRAS, 502, 3035
Salvatier J., Wiecki T. V., Fonnesbeck C., 2016, PeerJ Comput. Sci., 2, e55

Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132

Sellentin E., Heavens A. F., 2018, MNRAS, 473, 2355

Sellentin E., Heymans C., Harnois-Déraps J., 2018, MNRAS, 477, 4879

Taylor P. L., Kitching T. D., Alsing J., Wandelt B. D., Feeney S. M., McEwen
J. D., 2019, Phys. Rev. D, 100, 023519

Upham R. E., Brown M. L., Whittaker L., 2021, MNRAS, 503, 1999

SUPPORTING INFORMATION

Supplementary data are available at *MNRASL* online.

**Figure B.1.** The different LN fields analysed in this paper, with ground truth values of $\alpha$ and $\beta$ indicated above the panels.

**Figure C.1.** Simulations of the 2PCF of LN fields run at ground truth values $\alpha = 1, \beta = 0.5$.

**Figure D.1.** SBI analysis of the 2PCF of LN fields.

**Figure D.2.** Effect of varying $\varepsilon$ on the SBI posterior.

**Figure G.1.** Box plots showing the distribution of $\ln(\mathcal{L}_{\text{gt}}/\mathcal{L}_{\text{max}})$, a proxy for the bias of LBA, as a function of setup parameters.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.