Estimating Cosmological Constraints from Galaxy Cluster Abundance using Simulation-Based Inference

Simulation-based inference (SBI) (Cranmer et al., 2020) methods can incorporate complex physical processes and observational effects in forward simulations (Tam et al., 2022). SBI is a machine-learning approach to learn the sampling distribution of data as a function of the model parameters. The main goal of simulation-based inference is to identify parameter sets which are both compatible with prior knowledge and match empirical observations. The outputs of SBI are not point estimates; rather all probabilistic values of parameter space consistent with the inputs are identified (Lueckmann et al., 2021). We use SBI (Tejero-Cantero et al., 2020), a neural network-based PyTorch package (Ketkar & Moolayil, 2021). The process is explained by the block diagram shown in Figure 1. The prior is first sampled to obtain an initial set of parameters, which are used to create synthetic data using forward simulations. Data generated by forward simulations are then passed into a gaussian mixture density network. The output of this network is the probability distribution of the cosmological variables as a function of the observables (number counts and mean masses of dark matter halos). The posterior (probability distribution function) is then sampled to generate the parameter space of the cosmological variables for specific observables. One of the major advantages of SBI over traditional likelihood-based methods is that it can be implemented in stages – simulations, training, inference. It provides us with an opportunity to more efficiently perform complicated analysis or analysis of large volume of data.

The Quijote simulation suite contains subsets of simulations which are run for different sets of cosmological parameters (Villaescusa-Navarro et al., 2020). One such example is the latin-hypercube simulation set which contains 2000 simulations and their cosmological parameters are varied from the standard fiducial cosmology model with $\Omega_{m}$ in the range of $[0.1-0.5]$, $\Omega_{b}=[0.03-0.07]$, $h=[0.5-0.9]$, $n_{s}=[0.8-1.2]$ and $\sigma_{8}=[0.6-1.0]$. We use the same summary statistics (described in Section 2) for the galaxy cluster obervables from these simulations to train the machine learning model.

The galaxy cluster observables can also be calculated using analytical formulas as:

In these equations, $N(\Delta_{M}|z,\Theta)$ and $M(\Delta_{M},z|\Theta)$ are the mock observables calculated from analytical models. $V(z)$ is the cosmic volume of the Quijote simulations that the problem is based upon. $n(M,z|\Theta)$ is the theoretical halo mass function that describes the number density of dark matter halos at redshift $z$, which has an analytical form depending on cosmology, indicated by $\Theta$. In this analysis $\Theta$ refers to the five cosmological parameters, $\Omega_{m}$, $\Omega_{b}$, $h$, $n_{s}$ and $\sigma_{8}$. We make use of the Bhattacharya et al. (2011) halo mass function for the FOF halo mass definitions implemented in Colossus, with a linear correction (Costanzi et al., 2019) to account for differences with the Quijote simulations. Furthermore, $\delta_{N}$ and $\delta_{M}$ represent noises added to the analytical models, caused by cosmic variance. We use the Quijote fixed cosmology simulations to estimate the cosmic variances for $N(\Delta_{M}|z,\Theta)$ and $M(\Delta_{M}|z,\Theta)$, and then randomly draw a gaussian uncertainty according to those variances, $\delta_{N}$ and $\delta_{M}$, for each set of mock observables. In the end, we generate those analytical model simulations for over 10,000 sets of cosmological parameters. Specifically, we use a total of 11378 simulations to train the model. It is to be noted that the fast analytical simulations give valid insights and generating these simulations is feasible for this particular study.

Later in this analysis, the cosmic variances used in this fast analytical methods are used as the covariance matrices in a multi-dimensional Gaussian likelihood implemented in a Markov Chain Monte Carlo method as a comparison analysis.

We derive the posterior distributions of the five cosmological parameters with SBI method using the Quijote simulations (blue dashed lines), and the analytical models (green solid lines) as the training samples. The black dashed lines represent the truth values of those parameters. For comparison, we show the posterior results from Markov Chain Monte Carlo (MCMC) process, a state-of-the-art default method for sampling the posterior parameter distributions in such an analysis. In this comparison, the MCMC results are shown by red dashdot lines.

As shown in Figure 2 and listed in Table 1, the truth values of the parameters are reasonably recovered by the MCMC method, and the SBI method applied to both the Quijote simulations and the analytical-model based simulations. For $\Omega_{m}$, the truth value is recovered at 0.1$\sigma$, 1.6$\sigma$ and 2.2$\sigma$ level by the Quijote simulations, analytical models and MCMC method respectively. For the other four parameters, the truth values are recovered within the 1$\sigma$ range by all three methods. For $h$ and $n_{s}$, the SBI bias (both Quijote and analytical simulations) is much smaller than the corresponding MCMC bias.

In cosmological analysis, it is also important to accurately estimate the posterior uncertainties of the cosmological parameters, so as to correctly evaluate consistency and tensions between different cosmological models. While the MCMC method (state-of-the-art model) and the analytical-model based SBI method yield similar levels of uncertainties, the Quijote simulations based SBI method results in larger uncertainties than the other two methods for all parameters except $\Omega_{b}$. The histograms in Figure 3 show the bias and uncertainty distributions for 100 test samples for the analytical simulations (red) and the Quijote simulations (blue) for $\Omega_{m}$ and $h$. Both methods result in symmetric bias distributions. For the analytical simulations, the bias tends to be distributed over a smaller range of values, and uncertainties are smaller (same result as Figure 2) than the Quijote simulations. We suspect that this might be due to the training sample size being too small in the latin-hypercube (Quijote) based results, and explore further in the next subsection.

We further evaluate how the bias and uncertainty level of the posterior constraints vary with different training sample sizes used in analytical-model based SBI method. We define bias and uncertainty ($\sigma$) as:

$bias=\bar{\theta}_{Posterior}-\theta_{truth}$

$\sigma=\sqrt{\frac{\sum_{n=1}^{N}(\theta_{Posterior,i}-\bar{\theta}_{Posterior})^{2}}{N}}$

In Figure 4, we plot the uncertainty (black circles) and bias (black squares) for different sizes of training sets using data generated from the analytical models for $\Omega_{m}$ and $h$ . We start with 1000 samples and increase the training sample size in steps of 1000 until we reach 10000 samples. Uncertainty and bias corresponding to the MCMC (blue dashed lines and blue circles) and latin-hypercube simulations (red stars) are also plotted on the same axes. The uncertainty plots (upper panel) show that apart from some random fluctuations, uncertainty reduces with an increase in the number of training samples for $\Omega_{m}$ and stablizes when the sample size reaches $\sim$5000, for which the MCMC bias is about 33% lower than the analytical-simulations SBI bias. A similar trend is also observed for $\sigma_{8}$ (not shown). There is no observed correlation between the training sample size and the uncertainty of $h$, $\Omega_{b}$ (not shown), and $n_{s}$(not shown), indicating that the contraints are not affected by the training sample size we tested here. The convergent value of analytical SBI uncertainty for these parameters is comparable to that of the MCMC method.

According to the results shown in the lower panel of Figure 4, on average, the SBI bias based on analytical simulations diminishes as the size of training sample increases for $\Omega_{m}$. For $h$, the bias is nearly constant and independent of training sample size. Like the uncertainty, the bias for $\sigma_{8}$ shows a trend similar to $\Omega_{m}$, and the bias for $\Omega_{b}$ and $n_{s}$ follow the trend similar to h. MCMC bias is comparable to the SBI analytical bias for $\Omega_{m}$, and is higher than the analytical bias for $h$.