Non-Gaussianity in the Weak Lensing Correlation Function Likelihood -- Implications for Cosmological Parameter Biases

TL;DR

This study investigates the non-Gaussian features of weak lensing shear correlation functions, assesses their impact on cosmological parameter estimation, and proposes PCA-based data compression methods to mitigate biases and retain information.

Contribution

It demonstrates that non-Gaussianity has minimal bias on parameters and introduces a PCA-based likelihood modeling approach for efficient data analysis.

Findings

Non-Gaussianity shows significant skewness and kurtosis in shear correlation functions.

Parameter biases due to non-Gaussian likelihoods are negligible for 100 deg² mock data.

PCA reduces data dimensionality by a factor of ~5 while preserving most cosmological information.

Abstract

We study the significance of non-Gaussianity in the likelihood of weak lensing shear two-point correlation functions, detecting significantly non-zero skewness and kurtosis in one-dimensional marginal distributions of shear two-point correlation functions in simulated weak lensing data. We examine the implications in the context of future surveys, in particular LSST, with derivations of how the non-Gaussianity scales with survey area. We show that there is no significant bias in one-dimensional posteriors of $\Omega_{\rm m}$ and $\sigma_{\rm 8}$ due to the non-Gaussian likelihood distributions of shear correlations functions using the mock data ($100$ deg$^{2}$). We also present a systematic approach to constructing approximate multivariate likelihoods with one-dimensional parametric functions by assuming independence or more flexible non-parametric multivariate methods after decorrelating the data points using principal component analysis (PCA). While the use of PCA does not modify the non-Gaussianity of the multivariate likelihood, we find empirically that the one-dimensional marginal sampling distributions of the PCA components exhibit less skewness and kurtosis than the original shear correlation functions.Modeling the likelihood with marginal parametric functions based on the assumption of independence between PCA components thus gives a lower limit for the biases. We further demonstrate that the difference in cosmological parameter constraints between the multivariate Gaussian likelihood model and more complex non-Gaussian likelihood models would be even smaller for an LSST-like survey. In addition, the PCA approach automatically serves as a data compression method, enabling the retention of the majority of the cosmological information while reducing the dimensionality of the data vector by a factor of $\sim$5.

Figures

Equations

$$n(z)\propto z^{\alpha }\exp \left[-{\left(\frac{z}{z_0}\right)}^{\beta }\right],$$

$${ϵ}^{\mathrm{obs}}=\frac{\gamma +{ϵ}_{\mathrm{int}}}{1+\gamma {ϵ}_{\mathrm{int}}^{\ast }}.$$

$$\mathcal{A}=\left(\begin{array}[]{cc}1-\kappa-\gamma_{1}&-\gamma_{2}\\ -\gamma_{2}&1-\kappa+\gamma_{1}\end{array}\right),$$

$${\gamma }_t=-\mathbf{Re}(\gamma e^{-2\mathrm{i}\varphi }),{\gamma }_{\times }=-\mathbf{Im}(\gamma e^{-2\mathrm{i}\varphi }),$$

$${\xi }_{\pm }^{ij}(\theta )=⟨{\gamma }_t^i{\gamma }_t^j⟩(\theta )\pm ⟨{\gamma }_{\times }^i{\gamma }_{\times }^j⟩(\theta ).$$

$${\hat{\mathbf{C}}}^{ij}=\frac{1}{\nu }\sum _k^{N_s}({\xi }_k^{ij}-{\bar{\xi }}^{ij})({\xi }_k^{ij}-{\bar{\xi }}^{ij}{)}^T,$$

$${\xi }_{\text{theory}}^{ij}(\theta )/⟨{\xi }_{\text{mock}}^{ij}(\theta )⟩.$$

$$C_{\kappa \kappa }^{ij}(l)=\frac{9H_0^4{\Omega }_{\mathrm{m}}^2}{4c^4}{\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi \frac{g^i(\chi )g^j(\chi )}{a^2(\chi )}{P}_{\delta }\left(\frac{l}{f_K(\chi )},\chi \right),$$

$$g^i(\chi )={\int }_{\chi }^{{\chi }_{\mathrm{h}}}\mathrm{d}{\chi }^{\prime }{n}^i({\chi }^{\prime })\frac{f_K({\chi }^{\prime }-\chi )}{f_K({\chi }^{\prime })}.$$

$${\xi }_{\pm }^{ij}(\theta )=\int \frac{dll}{2\pi }{J}_{0/4}(l\theta )C_{\kappa \kappa }^{ij}(l),$$

$$L(\xi \mid \pi )\propto \exp [-\frac{1}{2}(\xi -{\xi }_{\pi }{)}^T{\mathbf{C}}^{-1}(\xi -{\xi }_{\pi })],$$

$$\hat{\Psi }=\frac{\nu -N_d-1}{\nu }(\hat{\mathbf{C}}{)}^{-1}$$

$$L(\xi \mid \pi )\propto \exp [-\frac{1}{2}(\xi -{\xi }_{\pi }{)}^T\hat{\Psi }(\xi -{\xi }_{\pi })].$$

$$\hat{\Lambda }=Q^T\hat{\mathbf{C}}Q$$

$$\begin{array}{rl}\mathrm{Edgeworth}(x)& =\frac{1}{\sqrt{2\pi }}\exp \left[-\frac{(\Delta x{)}^2}{2}\right]\\ & \cdot [1+\frac{{\kappa }_3}{3!}H{e}_3\left(\Delta x\right)+\frac{{\kappa }_4}{4!}H{e}_4\left(\Delta x\right)\\ & +\frac{10{\kappa }_3^2}{6!}H{e}_6\left(\Delta x\right)+\cdots ],\end{array}$$

$$H{e}_n(x)=(-1{)}^n{e}^{x^2/2}\frac{d^n}{d{x}^n}{e}^{-x^2/2},$$

$$\tilde{f}(k)={\int }_{-\infty }^{\infty }{e}^{\mathrm{i}kx}f(x)\mathrm{d}x.$$

$$\ln \tilde{f}(k)=\sum _{n=1}^{\infty }{\kappa }_n\frac{(\mathrm{i}k{)}^n}{n!}.$$

$$\beta (x)=\sum {\alpha }_i{\psi }_i(x).$$

$$\text{Skew}[X]=\frac{\mathbb{E}\left[(X-\mu {)}^3\right]}{{\sigma }^3}$$

$$\text{Kurt}[X]=\frac{\mathbb{E}\left[(X-\mu {)}^4\right]}{{\sigma }^4}-3$$

$$D^2(p,q)=2\mathbb{E}\Vert X-Y\Vert -\mathbb{E}\Vert X-X^{\prime }\Vert -\mathbb{E}\Vert Y-Y^{\prime }\Vert$$

$$\begin{array}{rl}T& =\frac{1}{m^2}\sum _{i=1}^m\sum _{j=1}^n{K}_h(X_i,X_j)-\frac{2}{mn}\sum _{i=1}^m\sum _{j=1}^n{K}_h(X_i,Y_j)\\ & +\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n{K}_h(Y_i,Y_j).\end{array}$$

$${\hat{\xi }}_{XY}(\theta )=$$

$${\mathcal{A}}_W(\theta )=\int d^2{\theta }^{\prime }W({\theta }^{\prime }+\theta )W({\theta }^{\prime })=\int \frac{d^2\mathbf{\ell }}{(2\pi {)}^2}{e}^{-\mathrm{i}\mathbf{\ell }\cdot \theta }\tilde{W}(\mathbf{\ell })\tilde{W}(-\mathbf{\ell })=\int \frac{d\mathbf{\ell }\ell }{2\pi }{J}_0(\ell \theta )\tilde{W}(\ell )\tilde{W}(-\ell )$$

$${\text{S}}_3({\hat{\xi }}_{{\kappa }_1{\kappa }_2}({\theta }_i){\hat{\xi }}_{{\kappa }_3{\kappa }_4}({\theta }_j){\hat{\xi }}_{{\kappa }_5{\kappa }_6}({\theta }_l))$$

$$=⟨{\hat{\xi }}_{{\kappa }_1{\kappa }_2}({\theta }_i){\hat{\xi }}_{{\kappa }_3{\kappa }_4}({\theta }_j){\hat{\xi }}_{{\kappa }_5{\kappa }_6}({\theta }_l)⟩+2⟨{\hat{\xi }}_{{\kappa }_1{\kappa }_2}({\theta }_i)⟩⟨{\hat{\xi }}_{{\kappa }_3{\kappa }_4}({\theta }_j)⟩⟨{\hat{\xi }}_{{\kappa }_5{\kappa }_6}({\theta }_l)⟩$$

$$-\left\{⟨{\hat{\xi }}_{{\kappa }_1{\kappa }_2}({\theta }_i)⟩\text{Cov}({\hat{\xi }}_{{\kappa }_3{\kappa }_4}({\theta }_j),{\hat{\xi }}_{{\kappa }_5{\kappa }_6}({\theta }_l))+\text{perms}\right\}$$

$${⟨{\kappa }_1{\kappa }_2{\kappa }_3{\kappa }_4{\kappa }_5{\kappa }_6⟩}_{ij}=$$

$$W(\theta )W({\theta }^{\prime })W(\theta +{\theta }_i)W({\theta }^{\prime }+{\theta }_j)W({\theta }^{\prime \prime })W({\theta }^{\prime \prime }+{\theta }_l)$$

$${⟨{\kappa }_1{\kappa }_2{\kappa }_3{\kappa }_4{\kappa }_5{\kappa }_6⟩}_{ij}=$$

$$\times e^{i({\mathbf{\ell }}_1-{\mathbf{q}}_1)\cdot \theta }{e}^{i({\mathbf{\ell }}_2-{\mathbf{q}}_2)\cdot (\theta +{\mathbf{r}}_{\mathbf{i}})}{e}^{i({\mathbf{\ell }}_3-{\mathbf{q}}_3)\cdot {\theta }^{\prime }}{e}^{i({\mathbf{\ell }}_4-{\mathbf{q}}_4)\cdot ({\theta }^{\prime }+{\theta }_j)}{e}^{i({\mathbf{\ell }}_5-{\mathbf{q}}_5)\cdot {\theta }^{\prime \prime }}{e}^{i({\mathbf{\ell }}_6-{\mathbf{q}}_6)\cdot ({\theta }^{\prime \prime }+{\theta }_l)}$$

Full text

1 Introduction

Weak gravitational lensing is the deflection of light by the gravitational field of large-scale structure, which leads to minute distortions of the observed galaxy images compared to their original shapes in the galaxy source plane. Measuring the correlation functions of the galaxy shapes is therefore a way to measure the growth of structure and the geometry of the Universe (e.g., Bartelmann & Schneider, 2001; Hoekstra & Jain, 2008; Kilbinger, 2015; Mandelbaum, 2018) and hence a promising avenue to constrain cosmology (Huff et al., 2014; Jee et al., 2016; Hildebrandt et al., 2018; Troxel et al., 2018a; Hikage et al., 2019).

With the next generation of weak lensing surveys, such as the Large Synoptic Survey Telescope (LSST111http://www.lsst.org/lsst; Ivezić et al. 2019), the Wide-Field Infrared Survey Telescope (WFIRST222https://wfirst.gsfc.nasa.gov) and Euclid333http://sci.esa.int/euclid/, we expect data sets that are both wider and deeper compared to current surveys (e.g. KiDS444http://kids.strw.leidenuniv.nl, DES555https://www.darkenergysurvey.org, HSC666http://hsc.mtk.nao.ac.jp/ssp/) in the near future. The high-quality data from upcoming lensing surveys is expected to reduce the statistical uncertainty in weak lensing measurements compared to current surveys by an order of magnitude. In order to fully exploit the cosmological constraining power of weak lensing surveys, much effort has been committed to understand observational systematic effects, such as uncertainties in the measurement of galaxy shapes (e.g., Massey et al., 2013; Mandelbaum et al., 2015) and photometric redshifts. On the astrophysical side, major progress has been made in modeling systematics that affect the interpretation of the weak lensing signal, such as the intrinsic alignment of galaxies (e.g., Mandelbaum et al., 2011; Joachimi et al., 2015; Troxel & Ishak, 2015; Krause et al., 2016), the nonlinear evolution of the dark matter density field (e.g., Takahashi et al., 2012; Heitmann et al., 2014), and baryonic effects that modify the latter (e.g., Eifler et al., 2015; Mead et al., 2015; Chisari et al., 2018).

However, inaccuracies in the last step of the analysis, the inference of cosmological parameters from measurements of observables such as $\xi_{\pm}$ and some model for the likelihood function, are less well-explored in the weak lensing community. Uncertainties related to parameter space sampling, discrepancy metrics, and the likelihood of the summary statistics itself are important aspects of the cosmological interpretation that can lead to potential biases in the analysis.

While likelihood-free approaches such as Approximate Bayesian Computation (e.g., Akeret et al., 2015; Peel et al., 2017) are beginning to emerge as a tool for cosmological inference, most analyses still assume a likelihood function to transition from observations to cosmological parameters. Among the possible choices, the multivariate Gaussian likelihood function is the simplest and most commonly used.

Even though the Cosmic Microwave Background (CMB) temperature field is close to a Gaussian and the non-Gaussian features in CMB power spectra are small, the CMB data analyses (Hinshaw et al., 2013; Planck Collaboration et al., 2014, 2016) have made some progress beyond this simple Gaussian assumption.

They use second-order statistics to capture the cosmological information in the underlying temperature and polarization field, which are (very close to) Gaussian. For idealized CMB observations (full-sky, isotropic beam, spatially uniform noise) the empirical power spectra of the underlying Gaussian temperature and polarization fields has a Wishart distribution given the model power spectra. At sufficiently high $\ell$, the likelihood function of power spectra approaches a multivariate Gaussian following the Central Limit Theorem. When going beyond the idealized case, the inclusion of potentially non-Gaussian foreground distributions such as galactic dust emission, Cosmic Infrared Background, and radio point sources breaks the initial assumption that the measured field is Gaussian and consequently breaks the conclusion that the likelihood of the power spectra beyond a certain $\ell$ is well-approximated by a multivariate Gaussian. As further detailed in Planck Collaboration et al. (2016), the Planck analysis masks foreground contaminants and assumes that the non-Gaussian features are subdominant outside of the masked regions.

The situation is different for weak lensing. Due to non-linear structure evolution at late times, the shear field itself is non-Gaussian, which invalidates the premise of the CMB argument. Nevertheless, weak lensing analyses have assumed a multivariate Gaussian likelihood function to be the underlying distribution of shear two-point statistics (Fu et al., 2008; Semboloni et al., 2011; Huff et al., 2014; Hildebrandt et al., 2016; Troxel et al., 2018a). In these analyses, the non-Gaussianity of the shear field enters only via non-linear matter density power spectra (Takahashi et al., 2012; Heitmann et al., 2014) that are used to model the observed two-point statistics and via so-called non-Gaussian covariances, which indicate that a non-vanishing four-point function is included in the covariance computation.

Non-Gaussian shear covariances due to non-linear clustering have been studied in Takada & Jain (2009); Sato et al. (2009); Sato et al. (2011b); Harnois-Déraps & van Waerbeke (2015). These studies use the multivariate Gaussian shear likelihood with the contribution of covariance from non-Gaussian fields included. The impact on cosmological constraints depends on the scales considered and on depth and area of the survey. Computing non-Gaussian covariances is essential in cosmic shear analyses, and most recent cosmic shear measurements have opted for different strategies, such as analytical computation of Gaussian and non-Gaussian covariances (e.g., Jee et al., 2016; Hildebrandt et al., 2017; Krause et al., 2017), or for covariance estimated through numerical simulations (e.g. Heymans et al., 2013). For future surveys, several studies indicate the high computational costs of a brute-force numerical simulations approach (Dodelson & Schneider, 2013; Taylor & Joachimi, 2014), which led to the development of new covariance estimators (Joachimi, 2017; Friedrich & Eifler, 2018). For cosmic shear, Barreira et al. (2018) have shown that the Gaussian covariance plus Super Sample Covariance terms are sufficient, both of which can be easily implemented analytically.

When going beyond the Gaussian likelihood function for the convergence power spectrum, previous studies have considered a lognormal distribution and the copula method for describing the non-Gaussian distribution (Taruya et al., 2002; Hilbert et al., 2011; Sato et al., 2011a, 2010). In configuration space, Hartlap et al. (2009) revisited the assumption of Gaussianity for the two-point correlation function and tested the non-Gaussianity with independent component analysis (ICA). The authors measure the distribution of the correlation function in 9600 realizations of ray-tracing simulations resembling the Chandra Deep Field South lensing analysis and perform three full likelihood analyses: using ICA, a standard multivariate Gaussian with a ray-tracing covariance, and with a Gaussian covariance. This paper triggered several attempts to build analytical expressions for the likelihood function of the shear two-point correlation function that improve over the standard Multivariate Gaussian approximation. For example, Schneider & Hartlap (2009) have shown that two-point correlation functions of Gaussian fields cannot take arbitrary values since this would violate the constraint of non-negativity of power spectrum. As a consequence the sampling distribution of the correlation functions cannot be an exact multivariate Gaussian. Keitel & Schneider (2011) employ Fourier mode expansion and characteristic functions of a Gaussian random field to derive an analytical expression for the likelihood function of its uni- and bi-variate correlation functions. In Wilking & Schneider (2013); Wilking et al. (2015), the authors transform the correlation functions such that a quasi-Gaussian approximation of the likelihood function is justified and tested its performance with simulations. A recent paper by Sellentin & Heavens (2018) explored the high-order correlations between the data points of the CFHTLenS cosmic shear correlation functions in search for non-Gaussianity. In Sellentin et al. (2018), the authors measure the skewed distributions of weak lensing shear correlation functions in simulations and follow the CMB literature in developing an analytical expression for the correlation function likelihood.

It is well-established in the literature that the sampling distribution of shear two-point correlation functions is not strictly Gaussian. Despite this, the Gaussian likelihood model is still the standard in weak lensing likelihood analysis for current surveys. If the likelihood is not Gaussian, then analysing the data with a Gaussian assumption could bias the cosmological parameter constraints. It is therefore important to quantify potential cosmological parameter biases that may arise due to un-modeled aspects of the likelihood function, so as to determine whether an alternative way of modeling the likelihood function is needed. This effect was discussed in the literature (e.g. Hartlap et al. 2009; Sellentin et al. 2018; Taylor et al. 2019), but agreement on whether the Gaussian approximation cause a significant bias for cosmological weak lensing analyses has not yet been reached. It is difficult to fully demonstrate the impact of the non-Gaussian likelihood on cosmological parameter constraints. Part of the difficulty comes from the dimensionality of the problem, with reconstructing in detail the full high-dimensional non-Gaussian likelihood currently being an unsolved problem. In this work, we do not address the full problem, but make some advance upon previous work by modelling the non-Gaussian likelihood with a large ensemble of weak lensing simulations, and further estimate the biases on cosmological parameters.

This paper is structured as follows: In Sect. 2, we describe the details of our simulated weak lensing data. In Sect. 3, after showing how the data vectors are modeled theoretically, we describe the likelihood analysis and likelihood models. Section 4 expands on the systematic approach of assessing the performance of likelihood models and data compression. In Sect. 5, we show the results for non-Gaussianity of weak lensing observables and the impact on cosmological parameter estimates. Section 6 contains our discussion and conclusions.

2 Simulations

The simulated (mock) weak lensing data that are used in this paper are based on the Scinet LIght Cone Simulations777https://slics.roe.ac.uk/ (Harnois-Déraps & van Waerbeke, 2015; Harnois-Déraps et al., 2018, SLICS hereafter), which are specifically tailored for statistical studies of weak lensing analyses. They consist of a series of lines-of-sight (LOS) of 100 deg2 each, constructed by ray-tracing in their own independent realization. In the simulations, no masks are applied, and hence the patches are 10 by 10 deg2. The underlying $N$-body simulations evolved $1536^{3}$ dark matter particles in a box length of 505 $h^{-1}$Mpc, and produced 18 mass planes between redshift 0.0 and 3.0, which are then converted into shear maps using the Born approximation. We used 932 such independent realizations in which the initial random seeds changed prior to the $N$-body run, with the assumed cosmology fixed to that of WMAP9+SN+BAO flat $\Lambda$CDM cosmology (Hinshaw et al., 2013): $\Omega_{\rm m}=0.2905$, $\Omega_{\Lambda}=0.7095$, $\Omega_{\rm b}=0.0473$, $h=0.6898$, $\sigma_{8}=0.826$ and $n_{s}=0.969$.

The mock galaxy catalogues are then created in a way that is meant to reproduce the redshift distributions of weak lensing source galaxies in LSST (Chang et al., 2013):

[TABLE]

with $\{\alpha,\beta,z_{0}\}=\{1.21,1.0,0.5\}$, assuming a source number density of 26 gal/arcmin2. The mocks are split in 10 tomographic redshift bins $n_{i}(z)$, each containing the same number of galaxies. These distributions are further smoothed by a Gaussian kernel of width $\sigma=(1+z)\sigma_{z}$ and $\sigma_{z}=0.02$. For the detailed redshift distributions of the 10 LSST-like source bins, see Fig. A1 in Harnois-Déraps et al. (2018).

Besides the cosmological shear $\gamma$ (see Sect. 3.1), the observed ellipticity $\epsilon^{\rm obs}$ includes the intrinsic shapes of galaxies $\epsilon_{\rm int}$ through the shear addition formula:

[TABLE]

In the above expression, shear and ellipticities are written as complex variables, and the $\epsilon^{\rm obs}_{1/2}$ components are recovered from the real and imaginary parts respectively. The two components of the intrinsic galaxy shapes are each drawn from a Gaussian distribution with zero mean, a standard deviation of 0.29 inspired by the KiDS-450 (Hildebrandt et al., 2016) data, and the constraint that $|\epsilon_{\rm int}|^{2}\leq 1$.

In our analysis, we measure the non-Gaussian shapes of the likelihood function of shear correlation functions with the SLICS simulation, and then use them to estimate the parameter biases for large weak lensing surveys such as LSST. Due to the differences in survey areas, different sources of uncertainty dominate: the 100 deg2 SLICS simulations are shape-noise dominated while LSST data will be cosmic-variance dominated. Therefore, in this work we sometimes switch off the intrinsic shape noise in order to separately understand the contributions of cosmic variance and shape noise to the shape of the likelihood function. For all results that are presented, we refer to the results as ‘without shape noise’ or ‘with shape noise’.

3 Method

3.1 Cosmic shear correlation function data vector

The weak lensing effect is mathematically approximated as a linear transformation that maps the unlensed location to the lensed location. The transformation matrix $\mathcal{A}$, which connects the shape of a source with the observed images, can be written as

[TABLE]

where $\kappa$ is the convergence and $\gamma_{1},\gamma_{2}$ are the two components of the spin-two shear $\gamma=\gamma_{1}+\mathrm{i}\gamma_{2}$ in Cartesian coordinates. Conventionally the coordinate system is rotated so that the separation vector is parallel to the $x$-axis. The shear components are decomposed into the tangential direction ($\gamma_{t}$) and the cross direction ($\gamma_{\times}$) in the rotated coordinates:

[TABLE]

where $\phi$ is the polar angle of the separation vector of the two galaxies. From the shear components we can write the two shear correlation functions as a function of angular separation

[TABLE]

Here $i$ and $j$ are indices of tomographic redshift bins, and the angle brackets refer to ensemble average. The correlation functions are computed from the SLICS simulations with the package TreeCorr888https://github.com/rmjarvis/TreeCorr (Jarvis et al., 2004).

The angular bins for measuring shear correlation functions initially divide the scales from 0.32 to 400 arcmin into logarithmically spaced bins of width $\Delta({\rm ln}\theta)=0.23$; however, we further apply angular selections to minimize the impact of known limitations in the simulations (Harnois-Déraps & van Waerbeke, 2015). We require $\theta>0.8$ arcmin for $\xi_{+}$, to avoid resolution effects in the simulations at small scales. Since $\xi_{-}$ is more sensitive to small-scale uncertainties, we apply a more aggressive constraint: $\theta>6.5$ arcmin. We also require $\theta<160$ arcmin for $\xi_{+}$ to avoid scales with significant power loss due to the box size. On the other hand, $\xi_{-}$ is less affected by the box size effect within 400 arcmin; therefore we do not introduce extra constraints on the large scale for $\xi_{-}$.

These scales are similar to those used in the recent KiDS-450 cosmic shear analysis by Hildebrandt et al. (2017), and presume that the analysis has separately accounted for uncertainties at small scales due to baryon feedback, modeling of the non-linear power spectrum and of the intrinsic alignments of galaxies999Note that these effects could contribute to the level of non-Gaussianity in the data, but that is beyond the scope of this paper.. Modeling these uncertainties at the precision required by LSST will be challenging. In the recent DES analysis (Troxel et al., 2018b) conservative scale cuts were implemented in order to avoid biases due to imperfect modeling of said astrophysical effects. For this paper, we assume an optimistic scenario in which the LSST weak lensing measurement pipeline includes these small angular scales.

The complete data vector is the concatenation of all $\xi_{+}$ and $\xi_{-}$ values across all tomographic bins and $\theta$ bins. Each simulated realization has 55 correlation functions, each with 42 angular bins (24 for $\xi_{+}$ and 18 for $\xi_{-}$).

The estimator of the data covariance from the many simulation ensemble is defined as

[TABLE]

where $i$ and $j$ indicate the tomographic redshift bins, $\xi_{k}$ is the data vector of the $k$-th realization, $\bar{\xi}$ is the mean data vector across all simulated realizations, and $\nu=N_{s}-1$ is the number of degrees of freedom given that the mean is estimated from the data. If the number of data-points, $N_{d}$, exceeds the number of realizations, $N_{s}$, we can neither ensure that the data covariance matrix is positive definite nor control the error in the data covariance matrix and its inverse. Therefore, we rebinned the data vector to reduce the number of points from 2310 to 770 by combining the angular bins in groups of three. Our final $\theta$ binning was chosen such that $N_{s}=932$, $N_{d}=770$ and hence $N_{d}<N_{s}$.

An illustration of the data vector for a particular set of tomographic bins, the diagonal covariance matrix elements, and a comparison with analytic theory predictions is given in Fig. 1. The theoretical predictions are based on the HALOFIT method (Takahashi et al., 2012), which models the nonlinear power spectrum with a fitting function. The HALOFIT model agrees well in general with the measurements from the simulations on the scales of interest, but there are still percent-level errors compared to the simulations (Harnois-Déraps & van Waerbeke, 2015). The finite box effect of the simulations also introduces power drops with respect to the theoretical predictions (Harnois-Déraps & van Waerbeke, 2015). The mismatch between the simulations and the theory will be compensated by a correction factor

[TABLE]

The correction factor above is not intended to fix these or other limitations of the simulations (which have some simplifications compared to reality, e.g., lack of baryonic physics). Rather, the correction is applied in order to assure that any bias we observe in our likelihood analysis comes from the covariance instead of the very small but nonzero intrinsic mismatch between the simulations and the theory. Note that the correction factor is applied before the computation of the covariance, skewness and kurtosis, and is therefore included in these estimators.

3.2 Modeling of observables

For the simulated likelihood analyses in Sect. 5 we employ the CosmoLike analysis and forecasting software package101010www.cosmolike.info. CosmoLike has been used in several forecasts exploring joint analyses of multiple cosmological probes (Eifler et al., 2014; Krause & Eifler, 2017; Schaan et al., 2017) and systematics mitigation strategies, such as the impact of baryons and intrinsic alignment (Eifler et al., 2015; Krause et al., 2016). On the observational side, the code was used in a weak lensing analysis of Sloan Digital Sky Data (Huff et al., 2014), the analysis of DES science verification data (Becker et al., 2016) and the recent DES Year 1 analysis (Krause et al., 2017; Troxel et al., 2018a; Abbott et al., 2018).

We compute the linear power spectrum of the best fitting flat $\Lambda$CDM cosmology for WMAP9 + BAO + SN (Hinshaw et al., 2013) using the Eisenstein & Hu (1999) transfer function and model the non-linear evolution of the density field as described in Takahashi et al. (2012). From the density power spectrum $P_{\delta}(k,z)$, we compute the shear power spectrum using the Limber and the flat sky approximations as

[TABLE]

with $l$ being the 2D wave vector perpendicular to the line of sight, $\chi$ denoting the comoving coordinate, $\chi_{\mathrm{h}}$ is the comoving distance to the horizon, $a(\chi)$ is the scale factor, and $f_{K}(\chi)$ the comoving angular diameter distance (throughout set to $\chi$ since we assume a flat Universe). The lens efficiency $g^{i}$ is defined as an integral over the redshift distribution of source galaxies $n(\chi(z))$ in the $i^{\mathrm{th}}$ tomographic interval

[TABLE]

We compute the cosmic shear two-point functions $\xi_{\pm}$ using the flat-sky approximation

[TABLE]

with $J_{n}(x)$ the $n$-th order Bessel function of the first kind.

3.3 Likelihood functions and data covariance matrix

From Bayes’ theorem, we can compute the posterior distribution of the cosmological parameters. In the standard likelihood analysis, the likelihood function is parametrized as a multivariate Gaussian function

[TABLE]

where $\mathbf{\xi}$, $\mathbf{\pi}$ and $\mathbf{C}$ denote the data vector, cosmological parameters and the covariance matrix respectively. The covariance matrix is fixed throughout the analysis.

Analytically computed covariance matrices are noise-free and can be factorized or inverted without complications. However, our covariance matrix estimated from numerical simulations is inherently noisy, and the noise level is affected by the number of realizations ($N_{s}=932$ in our case) and the size of the data vector ($N_{d}=770$).

To compute the likelihood function in the form of Eq. (11), we need the inverse covariance matrix $\Psi$, also called precision matrix. An unbiased estimator of the precision matrix is given by (Anderson, 2003; Hartlap et al., 2007; Taylor et al., 2013)

[TABLE]

in the case that the noise is Gaussian-distributed and the data points are statistically independent with $\nu=N_{s}-1$.

If the observed data is drawn from a multivariate Gaussian, Sellentin & Heavens (2016) show that after marginalizing over the noisy covariance matrix, the likelihood measured from the simulation realizations follows a multivariate t-distribution rather than a multivariate Gaussian. An earlier correction proposed by Hartlap et al. (2007) uses the unbiased inverse covariance matrix and keeps the multivariate Gaussian distribution to construct the likelihood function when the underlying distribution is Gaussian and the covariance is noisy. Thus, the Hartlap’s correction method leads to

[TABLE]

The Sellentin-Heavens likelihood shows improvement over Hartlap’s correction in terms of parameter inference with the marginalization over noise in covariance matrices. But it is still unclear how to extend the t-distributed likelihood function to cases where the underlying distribution of the data is non-Gaussian. In this paper, we adopt Eq. (12) instead, since it is easier to apply to near-Gaussian distributions. Since we quantify the difference between the Gaussian likelihood and the non-Gaussian likelihood through bias in cosmological parameter space, the bias depends more on the asymmetry of the likelihood distributions and thus is less sensitive to the difference between these two methods.

3.4 PCA transformation

For the shear correlation functions $\xi_{\pm}$, the multivariate Gaussian likelihood function with the form described in Eq. (11) is the most commonly used likelihood model in the literature. Since it is not trivial to build robust multivariate non-Gaussian likelihood functions, we perform the principal component analysis (PCA) transformations first on the data vector to remove the correlation between the data points. PCA is an orthogonal transformation that transforms data points into coordinates without the linear correlations. It can be described as:

[TABLE]

with columns of the transformation matrix $Q$ containing the eigenvectors of the covariance matrix $\hat{\mathbf{C}}$ estimated from the simulations. After the PCA transformation, the matrix $\hat{\Lambda}$ is diagonal and the diagonal elements are the eigenvalues of the covariance matrix $\hat{\mathbf{C}}$. The new coordinates are usually referred to as principal components, and the components are sorted according to decreasing variance. In addition, PCA automatically concentrates the information into a smaller number of principal components. For instance, Harnois-Déraps & Pen (2013) studied the noisy covariance matrix of the matter power spectrum and found that the information content of the leading principal components remains stable as they raised the noise level by reducing the number of measurements from 200 to only 4.

Note that in Eq. (14) we use the covariance instead of its inverse in the PCA. It may sound counterintuitive that the data points with the highest variance, i.e. the highest uncertainty, contain more information. However, the principal components with highest variance are also those with the largest signal-to-noise ratio in our data (this was also found in Harnois-Déraps & Pen (2013)), and thus contain more information. To quantify the information content on the training data contained in the first N components, we perform a Markov Chain Monte Carlo (MCMC) analysis on the mock data with different numbers of principal components, and define the information content as the square root of the determinant of the inverse parameter covariance (Fisher matrix) of $\Omega_{\rm m}$ and $\sigma_{\rm 8}$, $\sqrt{det(\mathbf{F})}$ . Other parameters are fixed throughout the MCMC and a flat prior that limits parameters in the ranges 0.05<$\Omega_{\rm m}$<0.6 and 0.5<$\sigma_{\rm 8}$<1.1 is assumed. Figure 2 demonstrates the relation between the information content retained and the number of principal components used. It shows that if we apply data compression such that we have 80% fewer data points, we lose only 23% of the information content.

3.5 Multivariate likelihood models

After the PCA transformation, the data points are linearly uncorrelated. We next continue to use parametric and non-parametric functions to describe the likelihood distributions in the PCA space.

In this paper we consider the following likelihood models for the PCA coordinates:

• •

Gaussian function

• •

non-Gaussian Edgeworth function

• •

$k$-nearest neighbors

• •

spectral series

Here we describe the two parametric and the two non-parametric models listed above. The performance and results of the models are covered in more details in Sects. 4 and 5.

Under the assumption that the underlying likelihood function measured in simulations is close to a multivariate Gaussian function, we approximate the multivariate likelihood function in the PCA coordinates as a product of parametric “marginal distribution functions”. The independence of PCA components is a strong assumption, but it is the assumption that the standard multivariate Gaussian likelihood makes. The product of one-dimensional Gaussian distributions of PCA components is identical to the multivariate Gaussian likelihood.

The Edgeworth function is a Gaussian function multiplied by correction terms constructed by its cumulants. It serves as an improvement upon the Gaussian likelihood. In this paper, we adopt Petrov’s formula of the Edgeworth expansion (Blinnikov & Moessner, 1998; Petrov, 1962) and the coefficients in the expansion are fixed by the cumulants of the simulation data. In the case where the standard deviation $\sigma$=1, the first four terms in the expansion are

[TABLE]

where $He_{n}(x)$ are Hermite polynomials given by

[TABLE]

$\Delta x=x-\mu$ and $\kappa_{n}$ are cumulants. Moments and cumulants are two different sets of quantities that can summarize a distribution. Cumulants arise naturally from Fourier transformation. In the Fourier transformation, the probability density function $f(x)$ is transformed into

[TABLE]

The cumulants are then defined as the coefficients of the power series expansion

[TABLE]

Since the Edgeworth function is not guaranteed to be positive and could have oscillatory behavior, one should be careful with anomalies (negative probabilities) and avoid its use in strongly non-Gaussian cases. In our case, we do not find negative possibilities when modeling the marginal distributions of PCs with the Edgeworth expansion. We show in Fig. 3 an example of the marginal one-dimensional distributions of $\xi_{+}$ and the two parametric models.

In addition to these parametric models, we also model the likelihoods more generally by estimating the high-dimensional density ratio $\beta(x)=f(x)/g(x)$ non-parametrically. Unlike the parametric methods that approximate the multivariate distributions as products of marginal distributions, our non-parametric methods do not assume independence. In our case, we take $g(x)$ to be the Gaussian model. Once $\beta(x)$ is fitted, we can sample from the estimated density $f(x)$ by importance sampling with $g(x)$ as the proposal distribution. In this paper, the density ratio is estimated by non-parametric methods based on the k-nearest-neighbors kernel density estimator (Lincheng & Zhijun, 1985) and the Spectral Series estimators.

The k-nearest-neighbors estimator (knn) approximates the density at a point by a kernel smoothing applied to the k nearest neighbors of that point, and the Spectral Series estimator (Izbicki et al., 2014) combines orthogonal series expansion and adaptively chosen bases to construct non-parametric likelihood functions. Let $\{\psi_{i}\}$ be an orthonormal basis with respect to the data distribution; then the density estimator for spectral analysis has the form

[TABLE]

In low-dimensional non-parametric curve estimation, the basis is fixed to the usual choices, such as a Fourier basis. In the Spectral Series method, the basis is driven by data so as to capture the intrinsic dimensionality of the data (Izbicki et al., 2014). With the density ratio, we are improving the Gaussian likelihood model using non-parametric methods. This is different from measuring the high-dimensional density with non-parametric models directly. Since the non-parametric models have more degrees of freedom compared to parametric models and make no assumption on the likelihood distributions, including the non-parametric methods make our list of models more complete.

3.6 Skewness and Kurtosis

We examine particular departures from Gaussianity of the shear two-point correlation function by calculating higher moments of the distributions: skewness and kurtosis. The skewness can be quantified as the normalized expectation value of the third central moment

[TABLE]

and it measures the asymmetry of the distribution. Increasing the number of samples $N_{s}$ does not reduce the skewness, but rather reduces the uncertainty of the skewness estimate, so repeated measurements would not remove the non-Gaussianity. Gaussian functions are symmetric; their skewness is zero. The distribution of shear correlation functions, however, is not perfectly symmetric. In Appendix A.1 we derive for Gaussian fields the general expressions for the third moment of the likelihood of the shear correlation functions, from which the skewness can be predicted. Following the same derivation in Appendix A.1, the $n^{th}$ moment scales as $1/f_{sky}^{n-1}$ in general. Besides cosmic variance, the effect of shape noise can also be included in this framework. We show that for the scales much smaller than the survey size, the third moment decreases with the survey size as ${f_{\text{sky}}^{-2}}$ and hence the skewness as defined in Eq. (20) decreases as ${f_{\text{sky}}}^{-1/2}$ ($\sigma\propto{f_{\text{sky}}}^{-1/2}$). As the scale $\theta$ approaches the survey window size, the third moment rises faster than $\sigma^{3}$ and thus the skewness will increase. This trend is consistent with expectations from the Central Limit Theorem, which explains a decreasing skewness from an increase in survey area through the fact that the number of modes that are averaged over within given bin increases.

Besides the skewness, the asymmetry in the distribution is also captured by the mean-mode difference. In Appendix B, we limit the possible range of the mean-mode difference by assuming a unimodal distribution. For larger surveys, the mean-mode difference in terms of $\sigma$ also follows the same scaling relation as the skewness, $(\tilde{\xi}-\overline{\xi})/\sigma\propto{f_{\text{sky}}}^{-1/2}$.

Additionally, we measure the kurtosis of the likelihood function as a metric for the level of non-Gaussianity. The kurtosis is defined as:

[TABLE]

measures the symmetric outliers of the distribution. Since the fourth moment of the standard normal distribution equals 3, the kurtosis (or more precisely, the excess kurtosis) is defined as the normalized fourth moment minus 3.

4 Likelihood model assessment

The goal of this section is to introduce the statistical tools that we use to construct and assess the one-dimensional and multidimensional likelihood models in PCA coordinates. To properly compare different likelihood models and avoid overfitting, we use 10-fold cross-validation. That is, the 932 realizations are partitioned into 10 non-overlapping subsets. For each experiment, 9 of the subsets are used to compute the mean and covariance of the Gaussian model or used to train the non-parametric models, while the remaining subset is used for testing. After the models are built, we draw samples from the models (90 realizations for each sample) and compare these samples to the test samples using different two-sample tests and distance metrics. This process is then repeated 10 times with each of the 10 data sets used once as the testing data.

In order to quantify the performance of the one-dimensional models for the likelihood of the data in PCA space, we perform cross-validation as described above with univariate Kolmogorov-Smirnov (KS) tests. The two-sample KS test statistic is the maximum distance between two empirical cumulative distribution functions. The p-value of the two-sample KS test statistic is defined under the null hypothesis that the two samples are drawn from the same distribution. We take p<0.05 of this KS statistic as an indication that the null hypothesis is rejected and that the two sets of samples are not likely to have been drawn from the same distribution.

Besides tests on univariate distributions, we assess how close samples from the multivariate models are to test samples, using two (non-parametric) multivariate test statistics: the maximum mean discrepancy (MMD; Gretton et al. 2012) and the energy distance (ED; Székely & Rizzo 2004; Baringhaus & Franz 2004). The energy distance is a statistical distance between two probability distributions. It is defined as the square root of

[TABLE]

where $\mathbb{E}$ is the expectation value, $X$, $Y$, $X^{\prime}$, and $Y^{\prime}$ are independent random vectors, the distribution of $X$ and $X^{\prime}$ is $p$, and the distribution of $Y$ and $Y^{\prime}$ is $q$. Here we use the Euclidean metric (and a sample estimate of the above expression). The maximum mean discrepancy can then be seen as a generalization of the energy distance to reproducing kernel Hilbert spaces. More specifically, we use a Gaussian kernel $K_{h}(x,y)$ with bandwidth $h$ to measure the similarity between two vectors $x$ and $y$, and we define our MMD test statistic as the MMD sample estimate:

[TABLE]

where $X_{1},\ldots,X_{m}\sim p$ and $Y_{1},\ldots,Y_{n}\sim q$.

The two test statistics above provide quantitative distance measures between samples drawn from two sampling distributions of weak lensing correlations functions $p$ and $q$. In the case of weak lensing analysis, $X$ and $Y$ in the above notations are multivariate shear correlation function data vectors drawn from the hold-out simulation data or the trained models. These test statistics or distance metrics111111We use the publicly available R packages “kernlab” and “energy”. We implement MMD with the Gaussian kernel and determine the bandwidth by the median heuristic method, i.e. $h^{2}=\text{median}(\lVert X_{i}-X_{j}\rVert^{2})$. are invariant to orthogonal transformations (such as a PCA rotation) and often used for comparing higher-dimensional data.

5 Results

5.1 Biases due to non-Gaussian distributions of $\xi_{\pm}$

In the previous sections, we introduced the approach of building and assessing the likelihood distributions. Here we apply the tools to examine the biases due to non-Gaussianity in distributions of the shear correlation functions $\xi_{\pm}$.

If the likelihood function of $\xi_{\pm}$ is a multivariate Gaussian function, then its marginal distributions are Gaussian by construction. However, we detect non-Gaussianity in the marginal distributions of $\xi_{\pm}$ for many values of $\theta$ and tomographic bins. The finding that the likelihood of $\xi_{\pm}$ is skewed was previously discussed in Sellentin & Heavens (2018). They use CFHTLenS Clone simulations for the tomographic analysis of the CFHTLenS data that provide 1656 semi-independent simulations for the 210 data points of CFHTLenS. For a specific data point at $\theta=35$ arcmin, they found that the most likely lensing amplitude is about 5% below the mean, so the distribution is ‘left-skewed’.

In Fig. 4 we show the non-zero skewness and kurtosis of the 1-D distributions of $\xi_{\pm}$ in selected tomographic bins, and of PCA coordinates, for data with and without shape noise. For the shape noise-free data, the skewness and kurtosis both decrease as redshift increases. The magnitude of non-Gaussianity shown in the skewness and kurtosis is statistically significant and peaks roughly at $\theta=20$ arcmin. It is difficult to gain insight into the $\theta-$ dependence of the skewness of the shear 2PCF, since the latter is an integral over the $C(l)$ values with highly oscillating filter functions $J_{0/4}$.

Most of the scales that we consider are dominated by the shape noise, which strongly suppresses the skewness and kurtosis. At scales around 100 arcmin, the skewness and the kurtosis start to reach a comparable level as in the shape noise-free case, since the shape noise is relatively less important at large scales. Note that these results for skewness and kurtosis of the marginal 1D likelihoods of $\xi_{\pm}$ values do not fully represent the level of non-Gaussianity in the multivariate observable space, since the $\xi_{\pm}$ values are highly correlated across $\theta$ and redshift bins.

Compared to $\xi_{\pm}$, the sampling distribution of the data in PCA space has far less significant skewness and kurtosis. The PCA components are linear combinations of $\xi_{\pm}$ values. During the transformation, a large number of the $\xi_{\pm}$ data points are added with positive and negative weights. As a result the linear combination then has reduced skewness and is closer to Gaussian. This is shown in the last two rows of Fig. 4, which compares the level of skewness and kurtosis in the marginal 1D distribution of $\xi$ and in PCA components. We again note that the multivariate sampling distribution is not affected by the PCA transformation, the marginal distributions can change due to rotations and the resulting marginal likelihoods depend on the specific forms of rotations. After the PCA rotation, we find that the 1-D likelihoods empirically become more Gaussian, as illustrated in the last two rows of Fig. 4.

To estimate the biases in cosmological parameter due to the non-Gaussianity in distributions of $\xi_{\pm}$, we apply the maximum likelihood method. Maximum likelihood estimation provides an intuitive way of estimating the biases in cosmological parameters due to a failure to model the non-Gaussian distributions of weak lensing two point functions. In this method, we compute the maximum likelihood estimate of cosmological parameters of each realization in the mocks with the Gaussian covariance. Given our 932 simulations, this yields 932 MLEs in cosmological parameter space. Limited by the number of realizations, we perform the maximum likelihood estimation for individual parameters, either $\Omega_{\rm m}$ or $\sigma_{\rm 8}$ (separately), with the other parameters fixed to their true values in the simulations. To avoid a bias due to a small mismatch between the mocks and the theoretical prediction (see Fig. 1), the elements of data vectors for each simulation realization are rescaled by ratios ${\xi^{ij}_{\text{theory-fid}}(\theta)}/{\left<\xi^{ij}_{\text{mock}}(\theta)\right>}$.

In order to estimate the impact of incorrect likelihood models, we also perform the maximum likelihood estimation on 10,000 Gaussian-distributed samples of two-point correlation functions. We first model the likelihood distributions of the simulation realizations with a multivariate Gaussian, and then draw samples from them. The samples share the same mean value and covariance matrix as the mocks but follow a multivariate Gaussian distribution. Hence the difference between the 1-D parameter estimates of the two sets of samples (mocks and Gaussian samples) comes solely from the incorrect assumption of likelihood functions.

We show in Fig. 5 the 1-D parameter estimates of $\Omega_{m}$ and $\sigma_{8}$ for the mock data and the 10,000 Gaussian samples with and without shape noise. In finding the maximum likelihood estimates, only one parameter ($\Omega_{\rm m}$ or $\sigma_{\rm 8}$) is explored at a time. To compare the 1-D distributions, we consider two statistical tests: the independent T-test and the KS-test. The T-test determines whether there is a significant difference between the mean values of two samples. Judging from the high p-values of the T-test for all the four cases in Fig. 5, we do not see significant shifts between the average values of the two sets of samples. This indicates that the mean values of 1-D parameter posteriors are not affected by the non-Gaussian marginal distributions of shear correlation functions that we observe in the mock data. In addition, we find that the mean values in the four panels are located around the fiducial values of parameters ($\Omega_{\rm m}=0.2905$ and $\sigma_{\rm 8}=0.826$). The deviations from the true value with respect to the uncertainties are listed in Fig. 5.

Besides the T-test, we also report the p-values of the two-sample KS-test. The p-values of the KS statistics are all above the 5% threshold, meaning that we do not detect significant difference between the 1-D distributions of parameters of the two samples. With the T-test we learn that the average cosmological parameters are not shifted by the non-Gaussian features of likelihoods in the mocks. The KS-test results suggest that the 1-D distributions of parameters are neither skewed nor distorted in other ways significantly by failure to model the left-skewed marginal distributions of shear correlation functions.

5.2 Modeling the distributions of principal components

In this subsection, we show results of attempts to model the likelihood function of PCA components.

When we perform KS tests of the univariate Gaussian and Edgeworth models on the training data set, all of the PCA components have $p$-values that are approximately uniformly distributed between 0 and 1 as shown in the bottom row of Fig. 6. This means that both of the parametric models have good univariate fitting performance on the 90% training data for data with and without shape noise. When we perform KS tests for these models against the testing data (top row of Fig. 6), however, the $p$-values concentrate around 0 for higher principal components for all the parametric models we consider, including Gaussian, Edgeworth function to the second order and Edgeworth function to the fourth order. This indicates that neither one of the models generalizes well to hold-out data for the computed low-variance principal components as a result of the PCA decomposition. Including the components that are mostly dominated by noise would overfit the data. In addition, Fig. 6 shows that the Gaussian model fits decently to the leading principal components of the current data. We could not noticeably improve the fitting performance by using more complicated models such as the 1-D Edgeworth function, but we can avoid overfitting by adopting the PCA framework and discarding the high-order principal components with additional benefit of data compression. Note that the KS test statistic is the largest difference between the two empirical CDFs. The KS-statistics take on discrete values, since the empirical CDFs are discrete due to the finite number of samples.

An analysis of variances also shows that the first 40 components capture over 99.9% of the variance of the training data without intrinsic shape noise, whereas the corresponding number of components for data with shape noise is 400.

Figure 7 shows the multivariate performance of different likelihood models including parametric models (Gaussian, Edgeworth to the 2nd order, Edgeworth to the 4th order) and non-parametric models (k-nearest-neighbor and Spectral Series) for data with shape noise. The lower the statistics (that is, the smaller the distance between samples from the model and samples from the test sets), the higher the performance. The last column of Fig. 7 shows the null distribution of the test statistics, and tells us what distance values to expect if two samples or sets (with 93 realizations each) were drawn from the true distribution of the data. We estimated this distribution by repeatedly resampling the unmodeled data without replacement, and then in each fold computing distances between these samples and the hold-out sample.

For the data without shape noise, models based on 40 components give good performance on hold-out data, whereas we need around 150 components in the presence of shape noise as seen in Fig. 7 for similar performance. We do not see significant differences in performance between the multivariate Gaussian model and the other parametric or non-parametric methods. The estimated null distribution in the last column of Fig.  7 indicates that the MMD statistic is powerful enough to discriminate between models for our data set. There is a significant difference betweeen the MMD statistics of models and the real data, implying that there is still some room for improvement. But improving the models would require an order-of-magnitude increase in the number of realizations, which is currently beyond our reach. Hence, our conclusion is that given the current number of realizations and inherent data noise, more complex multi-dimensional likelihood models do not seem to improve upon a simple Gaussian likelihood model typical of an LSST-like survey.

The number of components that retains $99.9\%$ of the variance is just an estimate. In order to find the optimal number of components to use, we apply our multivariate testing framework to a $k$-dimensional Gaussian likelihood model for different candidate values of $k$. Figure 8 indicates that we can further reduce the dimension of the model to 150 for data with shape noise without oversmoothing or “underfitting” the data. For data without shape noise, we can reduce the dimensionality to 40, the same as for our previous 99.9% variance estimate.

5.3 Biases due to non-Gaussian distributions of principal components

In Sect. 5.1, the maximum likelihood method was used to directly estimate cosmological parameter biases due to the non-Gaussian likelihood function of weak lensing shear correlation functions. This method of estimating the level of parameter biases when using the Gaussian likelihood approximation suggested that those biases are subdominant to statistical uncertainties in the mocks. In this subsection, we consider the multivariate PCA models that we have built and assessed in Section 5.2 and estimate parameter biases by comparing the MCMC chains of different models with emcee 121212https://github.com/dfm/emcee (Foreman-Mackey et al., 2013) MCMC sampler. The MCMC chains only explore $\Omega_{m}$ and $\sigma_{8}$ with the other parameters fixed, and we use flat priors for $\Omega_{m}$ and $\sigma_{8}$ with range 0.05<$\Omega_{\rm m}$<0.6 and 0.5<$\sigma_{\rm 8}$<1.1. In the MCMC sampling, we assume a constant covariance estimated from the simulations at the fiducial cosmology, regardless of the changes in $\Omega_{m}$ and $\sigma_{8}$.

Figures 9 and 10 show contour plots of posteriors for the cosmological parameter constraints derived from the mock data with shape noise. The posteriors are re-centered to the true value as a result of the rescaling by the ratio in Eq. (7). Fig. 9 compares the performance of models with different numbers of principal components and shows the effect of data compression. The 20% difference in the contour areas is due to information loss in data compression. Note that the model with all 770 PCA components and 1-D Gaussian distributions is strictly identical to the standard multivariate Gaussian likelihood of $\xi_{\pm}$. Hence the contour labeled as G770 in Figure 9 is also what we expect from the standard Gaussian likelihood analysis. Figure 10, on the other hand, compares Gaussian and non-Gaussian Edgeworth models in the PCA coordinates to demonstrate any biases due to non-Gaussian likelihood functions. Consistent with the results of the maximum likelihood estimation, no significant difference is found between the contours of the Gaussian model and the Edgeworth model for the mock data with shape noise.

6 Conclusions and Discussions

It is well known that the multivariate likelihood function of weak lensing shear correlation functions is not a Gaussian. Approximating the likelihood with a Gaussian should introduce associated biases in recovered cosmological parameters; however, whether this would cause significant biases for weak lensing shear correlation functions for upcoming lensing surveys and invalidate the mainstream multivariate Gaussian likelihood assumption has not yet been established. In this paper, we make further advances towards answering this question by modeling and testing non-Gaussian likelihood functions with simulated weak lensing data, and estimating the resulting biases on cosmological parameters.

We have measured the level of non-Gaussianity of shear correlation function likelihoods in weak lensing simulations and suggest avenues to improve corresponding cosmic shear likelihood analyses. The simulations are based on 100 deg2 lines-of-sight with the same source redshift distribution and number density as expected for LSST. A systematic approach to constructing the likelihood models was used, and biases in the parameter space ($\Omega_{m},\sigma_{8}$) due to the assumption of a multivariate Gaussian likelihood model were assessed.

We explored the non-Gaussianity in univariate distributions of our data vectors and find detectable non-Gaussianity in the marginal 1D likelihoods of shear correlation functions in $\theta$ and redshift bins. Even though the marginal distributions of $\xi_{\pm}$ provide hints on the non-Gaussianity, the linear transformation of PCA reduces the level of skewness and kurtosis in the principal components significantly.

Besides the tests on marginal 1D likelihoods, we use maximum likelihood estimation and 2-D MCMC likelihood analyses with multivariate likelihood models to estimate the parameter biases due to the skewed likelihood distributions of shear correlation functions. In both analyses, we do not observe significant parameter biases in terms of $\Omega_{\rm m}$ and $\sigma_{\rm 8}$ due to the Gaussian likelihood assumption in our data. Note that the results of the Edgeworth likelihood model provide a lower limit on the bias due to non-Gaussian likelihood functions for the 100 deg2 simulations since this model is built from marginal Edgeworth corrections of simple Gaussian functions in the PCA space based on the assumption of independence of PCs and the marginal distributions of PCA components exhibit small skewness and kurtosis. The multivariate non-parametric models that are not sensitive to coordinate rotations provide more general constraints on the bias, but their performance in terms of predicting hold-out data does not exceed that of the Edgeworth model given our data. We are not asserting that this is the final conclusion for the parameter biases due to the uncertainty of likelihoods. Measuring the high-dimensional likelihoods is extremely difficult. It is likely that our likelihood models, especially the non-parametric models, would benefit from orders of magnitude of increase in the number of simulation realizations. However, given the currently available number of simulation realizations, our results show no significant biases using the standard multivariate Gaussian likelihood.

In Appendices A.1 and B we show that the skewness and the $(\rm{mode-mean})/\sigma$ of shear correlation functions decrease with the survey area as $f_{\text{sky}}^{-1/2}$. The non-Gaussianity of the likelihood and the resulting cosmological parameter biases would therefore be even smaller for large surveys such as LSST. The scaling relation is derived by assuming that the shear fields are Gaussian, corresponding to larger scales. Future work should explore the relative contributions of the non-linear growth (non-Gaussian field) to the non-Gaussianity in likelihoods.

Therefore, given the small bias measured from the current simulations, we do not expect noticeable biases due to the non-Gaussian distribution of the weak lensing shear correlation function for future generation of large-scale structure observations. A multivariate Gaussian likelihood will continue to be a valid approximation in cosmic shear analyses.

In addition to our findings on the Gaussian approximation to the likelihood function, PCA poses a straightforward avenue to solving some of the practical problems related to covariance matrices. Several practical challenges for debiasing the covariance matrix of the weak lensing observables from simulations have been featured in the literature (Schneider et al., 2011; Dodelson & Schneider, 2013; Taylor & Joachimi, 2014; Blot et al., 2016), all of which are connected to the large number of data points expected in future cosmic shear surveys.

One way to address the problem is to identify a more efficient observable and hence compress the cosmological information into fewer data points. We have shown that the assumption of a multivariate Gaussian in PCA coordinates is valid for the mock data and that PCA ranks the data points efficiently as a function of signal-to-noise. This allows us to reduce the dimensionality of the data vector by throwing out modes with low signal-to-noise, which in turn alleviates some of the pressing problems related to covariance inversion.

Data compression approaches have been suggested in the past, such as MOPED (Heavens et al., 2000; Heavens et al., 2017) and COSEBIs (Schneider et al., 2010); it remains to be seen if the likelihood function of these summary statistics is closer to a Gaussian than, e.g. for the shear two-point correlation function. Using PCA, we have been able to reduce the number of dimensions from 770 to 40 for data without shape noise and from 770 to 150 for data with shape noise with acceptable loss in cosmological information.

In brief, our conclusions are itemized as follows:

We find significant non-Gaussianity in the marginal distributions of the correlation functions, indicating that the multivariate likelihood distributions are non-Gaussian. 2. 2.

We estimate the bias in cosmological parameters induced by ignoring the non-Gaussianity of the likelihood with maximum likelihood estimation assuming a Gaussian likelihood, and do not detect strong biases in $\Omega_{\rm m}$ and $\sigma_{\rm 8}$ in our simulated data. 3. 3.

Since the $(\rm{mode-mean})/\sigma$ of shear correlation functions scales with the survey area as $f_{\text{sky}}^{-1/2}$ under the assumption that the fields are Gaussian and that the angular power spectra follow a gamma distribution, we expect the bias due to assuming a Gaussian likelihood to be smaller for LSST than for these small-area mock catalogs. 4. 4.

Fully reconstructing high-dimensional distributions directly from simulations is very difficult, requiring simulation volumes well beyond those presently available for this study. Our results based on the described set of simulated data suggest that neglecting the non-Gaussianity of the likelihood for shear-shear correlations is not a significant source of bias for ongoing surveys or even future ones such as LSST.

While the simulated weak lensing data implies that the impact of non-Gaussianity is negligible in current and future cosmic shear surveys, it is an open question if this result is stable for different redshift distributions, different survey parameters, and different cosmologies. In Eifler et al. (2009) they show that the covariance and hence the likelihood is cosmology-dependent and that the strength of this effect depends on the specific properties of the survey. It is also unclear if some of the systematics that affect shear observables can introduce non-Gaussianity (e.g., as foregrounds do in case of the CMB), which could warrant further mitigation strategies. Finally, it will be interesting to explore the non-Gaussianity of the likelihood function for the multi-probe analysis case, e.g., when including galaxy clustering, galaxy-galaxy lensing, and perhaps even higher-order statistics of auto and cross observables of clustering and shear. A multi-probe simulation with even more realizations is required for the joint data vector.

Acknowledgments

This paper has undergone internal review in the LSST Dark Energy Science Collaboration. The authors would like to give thanks to the DESC Publication Board and the internal reviewers A. Heavens, R. Rosenfeld and J. Zuntz for the feedback and comments that improved our manuscript. CHL, JHD, TE and RM developed the concept and designed the analysis framework; CHL carried out the analysis and prepared the manuscript with the help of JHD, TE, RM, ABL, TP and SS; JHD conducted weak-lensing simulations; TP and ABL designed and performed the statistical model assessment, and fitted the nonparametric likelihood models; SS derived expressions of the third moment of shear and mean-mode difference. JHD acknowledges support from the European Commission under a Marie-Skłodowska-Curie European Fellowship (EU project 656869) and from the European Research Council under grant number 647112. Computations for the $N$-body simulations were performed in part on the Orcinus supercomputer at the WestGrid HPC consortium (www.westgrid.ca), in part on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto. This work was partially supported by NSF DMS-1520786, DOE grant DESC0011114, and Department of Energy Cosmic Frontier program. Part of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration and is supported by NASA ROSES ATP 16-ATP16-0084 grant and by NASA ROSES 16-ADAP16-0116.

The DESC acknowledges ongoing support from the Institut National de Physique Nucléaire et de Physique des Particules in France; the Science & Technology Facilities Council in the United Kingdom; and the Department of Energy, the National Science Foundation, and the LSST Corporation in the United States. DESC uses resources of the IN2P3 Computing Center (CC-IN2P3–Lyon/Villeurbanne - France) funded by the Centre National de la Recherche Scientifique; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231; STFC DiRAC HPC Facilities, funded by UK BIS National E-infrastructure capital grants; and the UK particle physics grid, supported by the GridPP Collaboration. This work was performed in part under DOE Contract DE-AC02-76SF00515.

Data availability

The data generated and analyzed in this article are available from the corresponding author upon reasonable request.

Appendix A Third Central Moment

A.1 Convergence

Here we derive the expression for the third central moment of the (tomographic) cross-correlation function of two convergence maps denoted by $X$ and $Y$. Similar calculations have also been done by Keitel & Schneider (2011). We use a different notation that is easier to generalize for cross-correlations and restrict ourselves to computing the third central moment, though calculations can be extended straightforwardly but tediously to higher moments. We show here how skewness is affected by $f_{sky}$, a relation we used in Sect. 3.6 in the main text.

For simplicity, we will work in the flat-sky approximation. We will also assume that the fields are Gaussian and ignore the non-Gaussian terms, though we will indicate where the non-Gaussian terms should enter. Using simulations, we have already shown that at the small scales, where the non-Gaussian terms are most important, the impact of non-Gaussian likelihood is small.

The cross-correlation function of two fields can be written as

[TABLE]

where $\kappa_{1}$ belongs to field $X$ and $\kappa_{2}$ to field $Y$. $W(\mathbf{\theta})$ is the survey window function. We have assumed that the noise and the $\kappa_{i}$ have zero mean and are also uncorrelated with each other on all scales. The normalization factor is the integral over window functions

[TABLE]

The third central moment ($S_{3}$) of the correlation function is given as

[TABLE]

where $\kappa_{1},\kappa_{3},\kappa_{5}$ belong to field $X$ and $\kappa_{2},\kappa_{4},\kappa_{6}$ belong to $Y$. $\text{S}_{3}$ is a function of three $\theta$ variables and is therefore a third-order tensor. The ‘perms’ denotes permutations over combinations of $(12,34,56)$, where one combination gives $\xi$ and other two give covariance. So there are 6 permutations in the last term of Eq. (26) (or 3 terms if the symmetry of covariance is considered). We use $\left\langle\kappa_{1}\kappa_{2}\kappa_{3}\kappa_{4}\kappa_{5}\kappa_{6}\right\rangle_{ijl}$ as short-hand for the first term on the right hand side in Eq. (26), which we would like to simplify:

[TABLE]

Writing the $\kappa_{i}$ in terms of its Fourier space counterpart $\tilde{\kappa}_{i}$, we get

[TABLE]

We have integrated over $d^{2}\mathbf{\theta}$, $d^{2}\mathbf{\theta}^{\prime}$, $d^{2}\mathbf{\theta}^{\prime\prime}$ and then over $d^{2}\mathbf{\ell}_{2}$, $d^{2}\mathbf{\ell}_{4}$, $d^{2}\mathbf{\ell}_{6}$ to obtain the last expression.

We now expand the six-point function into two separable parts: the connected or non-Gaussian component $\left\langle\tilde{\kappa}_{1}\tilde{\kappa}_{2}\tilde{\kappa}_{3}\tilde{\kappa}_{4}\kappa_{5}\kappa_{6}\right\rangle^{\prime}$ and the Gaussian component, which using Wick’s theorem can be expanded as the sum of the product of two-point functions:

[TABLE]

where the ‘all perms’ now denote all possible combinations over $(1,2,3,4,5,6)$. Note that a similar expansion is also required for four-point functions in the covariance terms in Eq. (26) (see for example Singh et al. 2017 for a covariance expansion in similar notation). Simplifying, the auto-correlation terms from the same theta bin, e.g. $\left\langle\kappa_{1}\kappa_{2}\right\rangle$, cancel out and we are left with the permutations that contain only cross-correlation terms involving at most one $\kappa_{i}$ from each $\mathbf{\theta}_{j}$ bins (there are ${}^{4}C_{1}\times{}^{2}C_{1}=8$ such permutations). Thus we have

[TABLE]

To further simplify the expressions, we will assume that the scales of interest are smaller than the survey window size and ignore the coupling between the window function and the power spectra. Within this assumption, the window function integrals can be carried out without the power spectra to finally give

[TABLE]

where in the second equation we have used the fact that odd-odd combinations such as ‘13’ give the auto-correlation of field $X$, $C_{XX}$, even-even combinations give the auto-correlation of $Y$, $C_{YY}$ and the even-odd combinations give the cross-correlation $C_{XY}$. Since $A_{W}$ scales with the fraction of sky covered by a survey, $f_{\text{sky}}$, Eq. (33) suggests that the third central moment scales as $f_{\text{sky}}^{-2}$. Thus for any given scale that is well within the survey window size, the third central moment of the convergence correlation function decreases faster than the covariance as the survey area increases. However, as the scale $\theta$ approaches the size of survey, the area factor in the normalization approaches zero, i.e. $A_{W}(\theta)\rightarrow 0$ for large $\theta$ and $S_{3}$ will increase. $S_{3}$ rises faster than covariance, $S_{2}$ and thus the skewness increases at large scales. This is consistent with the expectation from the central limit theorem as large scales have fewer modes within the survey and we expect them to be more skewed.

We also note that we have ignored the connected terms, which can be important if the fields are non-Gaussian. The connected terms have three factors of area, $A_{W}$ in the denominator and a factor of window functions in the numerators. Under the assumption that the window and the connected terms can be decoupled, the window term in the numerator will cancel one factor of $A_{W}$ and the connected term will also scale similarly to the Gaussian term. The decoupling of the window is not well justified in the case of coupling between large-scale (super-sample) and small-scale modes, but such coupling terms have also been shown to scale with similar factors of area in the context of super-sample covariance. We have not tested such scaling for the six-point function in this work.

A.2 Shear

For shear we begin by noting that $\xi_{+}\propto\langle\gamma_{{}_{X}}\gamma_{{}_{Y}}^{*}\rangle$ and $\xi_{-}\propto\langle\gamma_{{}_{X}}\gamma_{{}_{Y}}\rangle$. Using these relations and the expressions from the previous subsection, we get

[TABLE]

where ‘cyc’ denotes the interchanging of order of the Bessel functions for different $\theta$. We conclude from the equation above that the same scaling with $f_{sky}$ applies to shear as well.

Appendix B Difference between mean and mode

Assuming a unimodal distribution for a random variable $X$, the difference between the most likely value, mode $\tilde{X}$, and the mean $\overline{X}$ in terms of the standard deviation $\sigma$ has an upper bound (Johnson & Rogers, 1951)

[TABLE]

In our main analysis using $100$ deg2 simulations, we have already shown that the difference between the mean and mode is much smaller than this bound.

In order to address the question of how this difference scales with the survey area, we begin by noting that for angular power spectra $C_{\ell}$, which follow a gamma distribution with $\nu\approx(2\ell+1)f_{\text{sky}}$ degrees of the freedom (see e.g.  Percival & Brown, 2006, for a detailed discussion of the $C_{\ell}$ likelihood for CMB), the difference between the mean and the mode is

[TABLE]

With the variance $\sigma_{C_{\ell}}^{2}\approx 2(C_{\ell}^{2}/\nu)$, the difference in terms of $\sigma$ is

[TABLE]

Here we assume the cosmic variance-dominated regime, but the scaling in the shot noise-dominated regime is similar since the shot noise also scales as $1/f_{\text{sky}}$ for a fixed number density of tracers.

Since the correlation functions, $\xi_{\pm}$, are the Hankel transform of $C_{\ell}$ (in the flat sky approximation), the difference between the mean and the mode for the correlation functions is

[TABLE]

where $n=0$ for $\xi_{+}$ and $n=4$ for $\xi_{-}$.

Since the covariance, $\mathbf{C}$, of $\xi$ scales as $1/f_{\text{sky}}$, the difference between the mean and the mode in terms of signal-to-noise ratio scales as

[TABLE]

Eq. (40) also gives the scaling of the bias in the log likelihood (Gaussian) and the parameter covariance with $f_{\text{sky}}$.