Power spectrum recovery¶

Gaussian and lognormal fields¶

We begin by testing the consistency of our Gaussian and lognormal maps generation pipeline.

We show example realisations of the two fields in Fig. 9, Gaussian on left and lognormal on the right. There’s a visible difference between the two, as it can be seen clearly from the distribution plot. The main check to perform is for testing whether the generated fields recover the theoretical power spectrum. Fig. 10 shows that this is the case for the Gaussian fields. They recover the fiducial $C(l)$ within a few percent error, with larger deviations $\sim 5\%$ at the low and high ends of the $l$-range. Instead, lognormal fields present deviations $\gtrsim 10\%$. As the lognormal transformations we use have been reported by different sources Zhou et al. (2023)Boruah et al. (2022), the issue must lie with our JAX implementation of the Hankel transformation. Resolving such issues could be achieved by future iterations of this work. In this work, we restrict ourselves to the use of Gaussian fields, as it is enough to prove our thesis and show that Gaussian processes can be applied to cosmological fields.

Gaussian process priors¶

First, we test the ability of the kernels we have built in Sec. ?? to recover the power spectrum of our cosmology.

We test five models in Fig. 11: integration, FFTlog, full-range FFT, half-range FFT and sinc FFTlog. The first four methods are described in the Gaussian process kernels Sec. ??, whereas the sinc FFTlog referes to a FFTlog model on which we applied smoothing, by multiplying the power spectrum by a factor of $sinc^4(l\frac{L}{2\pi N})$. The recovered power spectra are plotted against the fiducial power spectrum, or smoothed power spectrum for the sinc FFTlog. Mean and standard deviation associated to the plots are calculated from 500 samples. As expected the integration, FFTlog and full-range FFT perform similarly, as they all contain the same ammount of information. As these models deviate so strongly from the fiducial power spectrum we tried applying smoothing, which helps to recover half of the $l$-range at large scales. The only method that seems to be consistently recovering the fiducial power spectrum is the half-range FFT. One could argue that due to the inherent discreteness and boundedness of the fields we are working with, using FFTs is the most natural choice; also, half-range FFT uses the only grid that recovers a correlation function of the same shape as the field without having to perform binning.

We have also tested the efficacy of the half-range FFT model for different cosmologies of values $\{\sigma_8, S_8\}$ equal to $\{0.4,0.2\}$, $\{1.2,1.5\}$ and, our fiducial cosmology, $\{0.8,0.8\}$. As Fig. 12 shows, the model is independent of the choice of cosmology. From here on the results will be presented assuming a kernel built with the half-range FFT model.

Gaussian process map reconstruction¶

Armed with a reliable kernel, let’s embark upon the journey of reconstructing a heavily masked cosmological field. What we will do is: create a noiseless GRF in the fiducial cosmology Tab. %s, True map; apply a mask to obtain the Data map; condition a Gaussian process which assumes the fiducial cosmology. Fig. 13 lists the result of this operation, showing the resulting mean $\mu$ and standard deviation $\sigma$ of the conditioned GP. We also plot the ratio between residuals $\Delta=\mu-$True and standard deviation squared, to test the goodness of fit of our model, the values of the map sum up to $\chi^2 \sim 2495$. With the mask covering $\nu=2353$ pixels, we obtain $\chi^2 / \nu = 1.06$. Of course, this is just a noiseless application, which is unreasonable for a real application.

Inference of cosmological parameters¶

To test the ability of Gaussian processes to recover cosmological parameters without any prior knowledge except a noisy and masked map, we perform a MCMC simulation to infer the posterior distributions of $\sigma_8$ and $S_8$. We use the convention

to infer deterministically a posterior for $\Omega_m$. Such a reparametrisation is needed due to the strong degeneracy between $\sigma_8$ and $\Omega_m$. Eq. (30) breaks this degeneracy, changing the geometry of the sampling space and making the sampling more consistent. The model assumes uninformed flat priors for the cosmological parameters, as shown in Tab. %s, such prior bounds are also in accordance with the jaxcosmo release Campagne et al. (2023). The likelihood of the model is given by a Gaussian process distribution conditioned on Data, with a standard deviation equal to the noise applied to the map. The analysis is coded with numpyro Phan et al. (2019) Bingham et al. (2019), using a the No-U-Turn Sampler (NUTS) method with max_tree_depth=16, target_accept_prob=0.8. We simulate 8 chains for the $32\times32$ grid and 4 chains for the $64\times64$. Each chain performs 1000 warmup steps and 3000 samples.

One parameter¶

As a first step and for a consistency check, we run the inference model for one cosmological parameter, keeping all others fixed. Using a $64\times64$ grid with $n_g=10 \text{ galaxies}/\text{arcmin}^2$. In Fig. 15 we show the inferred distribution for both $\sigma_8$ and $\Omega_m$. We find that we are able to recover the true value for both parameters within two sigmas, $\sigma_8 = 0.776\pm0.015$ and $\Omega_m = 0.284\pm0.010$. We notice a slight tendency of the inferred distribution to be biased low; a tendency we also observe next for both sampled parameters, $S_8$ and $\sigma_8$.

Two parameters¶

Effect of noise¶

We perform some tests on low resolution $32\times32$ grids to see the effect that the noise level has on the recovered parameters, see Tab. 1. Here we report the inferred cosmological parameters for one data realisation and different noise levels, corresponding respectively to $n_g = 4$, 10, 30 and $100 \text{ galaxies}/\text{arcmin}^2$, see Eq. (22). The inferred value of $S_8$ can vary as much as a full $\sigma$ between high and low noise runs. Keeping in mind that $\sigma_8$ and $\Omega_m$ are extremely unreliable due to relative uncertainties of $\sim 25-30\%$ caused by the degeneracy: as a general trend we notice $\Omega_m$ gets bigger when $\sigma_8$ gets smaller with less noise.

Inferred cosmological parameters¶

Running the model for a larger $64\times64$ grid with $n_g=10 \text{ galaxies}/\text{arcmin}^2$, gives much better constraints on the cosmological parameters. We present the values recovered by the posterior distributions, listed as follows in Tab. 2.

Fig. 16 shows the inferred posterior distributions and contours for the three cosmological parameters $\sigma_8$, $\Omega_m$ and $S_8$. Looking at the contours, we obtain the well known banana-shaped degeneracy between $\sigma_8$ and $\Omega_m$. The $S_8$ and $\Omega_m$ contour presents sharp cuts for high and low $\Omega_m$, indicating an issue with the bounds of the uniform priors imposed. Unfortunately the jaxcosmo package does not allow for the choice of priors to be wider than what shown in Tab. %s, as the model then starts to have divergent samples.

Posterior checks¶

Following the two parameter inference model, we perform some posterior checks at the map level Porqueres et al. (2021). Fig. 18 sums up the ability of the model to recover the true map, noiseless and unmasked. Here we present the run with noise level $\sigma_\text{noise}\sim 0.0088$ and a $64\times64$ grid. We show the mean and standard deviation for the sample with highest likelihood out of the 12000. The mean field $\mu$ is visibly different to the true field in the masked regions and it seems to be of overall lower amplitude. The sample map is comparable to the noisy data; which is to be expected, as the internal noise given to the Gaussian process is the same as the noise level of the data. The standard deviation map $\sigma$ presents an overall amplitude comparable to the noise level $\sim 0.010$, with higher values for the masked regions. Summing up the map values of the residuals divided by standard deviation squared, we obtain a $\chi^2\sim 1297.4$. Compared to the number of free parameters $\nu$ in our inference model, which for a $10\%$ mask and a $64\times64$ grid, is $\nu=3689$. The value of $\chi^2$ therefore seems to be low, indicating that the noise level assumed by the GP is overestimated. This is supported by the fact that the sample map looks just as noisy as the data, according to Fig. 17, its distribution is in fact just as wide as the noise.