pnasresearcharticle Hahn The three-dimensional spatial distribution of
galaxies encodes key cosmological information on the nature of dark energy and
the contents of the Universe. Current analyses of the statistical clustering
of galaxies successfully extract information on large scales that are well
described by analytic models. They, however, struggle on smaller, non-linear,
scales. Here we present SimBIG, a new approach to galaxy clustering analyses
that can extract information on non-linear regimes by exploiting high-fidelity
simulations and inference based on machine learning. To demonstrate its
advantages, we apply SimBIG to 109,636 galaxies of the BOSS survey and analyze
a standard summary statistic of the galaxy distribution. Our constraints are
consistent with previous works and substantially improve their precision on
select cosmological parameters. C.H., M.E., S.H., P.L., E.M., A.M.D., B.R.B.,
and M.A. designed research; C.H. performed research; C.H., M.E., S.H., J.H.,
P.L., E.M., A.M.D., and B.R.B. contributed new reagents/analytic tools; C.H.,
M.E., S.H., J.H., P.L., E.M., C.M., A.M.D., and B.R.B. analyzed data; and
C.H., M.E., S.H., J.H., P.L., E.M., C.M., A.M.D., and B.R.B. wrote the paper.
The authors declare no conflict of interest<EMAIL_ADDRESS>
#  SimBIG: A Forward Modeling Approach To Analyzing Galaxy Clustering

ChangHoon Hahn  Department of Astrophysical Sciences, Princeton University,
Princeton NJ 08544, USA  Michael Eickenberg  Center for Computational
Mathematics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA  Shirley
Ho  Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave 5th
floor, New York, NY 10010, USA  Jiamin Hou  Department of Astronomy,
University of Florida, 211 Bryant Space Science Center, Gainesville, FL 32611,
USA  Max-Planck-Institut für Extraterrestrische Physik, Postfach 1312,
Giessenbachstrasse 1, 85748 Garching bei München, Germany  Pablo Lemos
Department of Physics, Université de Montréal, Montréal, 1375 Avenue Thérèse-
Lavoie-Roux, QC H2V 0B3, Canada  Mila - Quebec Artificial Intelligence
Institute, Montréal, 6666 Rue Saint-Urbain, QC H2S 3H1, Canada  Center for
Computational Astrophysics, Flatiron Institute, 162 5th Ave 5th floor, New
York, NY 10010, USA  Elena Massara  Waterloo Centre for Astrophysics,
University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada
Department of Physics and Astronomy, University of Waterloo, 200 University
Ave W, Waterloo, ON N2L 3G1, Canada  Chirag Modi  Center for Computational
Mathematics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA  Center
for Computational Astrophysics, Flatiron Institute, 162 5th Ave 5th floor, New
York, NY 10010, USA  Azadeh Moradinezhad Dizgah  Département de Physique
Théorique, Université de Genève, 24 quai Ernest Ansermet, 1211 Genève 4,
Switzerland  Bruno Régaldo-Saint Blancard  Center for Computational
Mathematics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA
Muntazir M. Abidi  Département de Physique Théorique, Université de Genève, 24
quai Ernest Ansermet, 1211 Genève 4, Switzerland

###### Abstract

We present the first-ever cosmological constraints from a simulation-based
inference (SBI) analysis of galaxy clustering from the new SimBIG forward
modeling framework. SimBIG leverages the predictive power of high-fidelity
simulations and provides an inference framework that can extract cosmological
information on small non-linear scales, inaccessible with standard analyses.
In this work, we apply SimBIG to the BOSS CMASS galaxy sample and analyze the
power spectrum, $P_{\ell}(k)$, to $k_{\rm max}=0.5\,h/{\rm Mpc}$. We construct
20,000 simulated galaxy samples using our forward model, which is based on
high-resolution Quijote $N$-body simulations and includes detailed survey
realism for a more complete treatment of observational systematics. We then
conduct SBI by training normalizing flows using the simulated samples and
infer the posterior distribution of $\Lambda$CDM cosmological parameters:
$\Omega_{m},\Omega_{b},h,n_{s},\sigma_{8}$. We derive significant constraints
on $\Omega_{m}$ and $\sigma_{8}$, which are consistent with previous works.
Our constraints on $\sigma_{8}$ are 27% more precise than standard analyses.
This improvement is equivalent to the statistical gain expected from analyzing
a galaxy sample that is $\sim$60% larger than CMASS with standard methods. It
results from additional cosmological information on non-linear scales beyond
the limit of current analytic models, $k>0.25\,h/{\rm Mpc}$. While we focus on
$P_{\ell}$ in this work for validation and comparison to the literature,
SimBIG provides a framework for analyzing galaxy clustering using any summary
statistic. We expect further improvements on cosmological constraints from
subsequent SimBIG analyses of summary statistics beyond $P_{\ell}$.

###### keywords:

cosmology $|$ machine learning $|$ galaxies $|$ simulation

###### doi:

www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX

This manuscript was compiled on

The three-dimensional spatial distribution of galaxies provides key
cosmological information that can be used to constrain the nature of dark
matter and dark energy and measure the contents of the Universe. The next
generation spectroscopic galaxy surveys, conducted using the Dark Energy
Spectroscopic Instrument (DESI; 1, 2, 3), Subaru Prime Focus Spectrograph
(PFS; 4, 5), the ESA Euclid satellite mission (6), and the Nancy Gracy Roman
Space Telescope (Roman; 7, 8), will probe galaxies over unprecedented cosmic
volumes out to $z\sim 3$, over 10 Gyrs of cosmic history. Combined with other
cosmological probes, they will provide the most stringent tests of the
standard $\Lambda$CDM cosmological model and potentially lead to discoveries
of new physics.

In current analyses, the power spectrum is used as the primary measurement of
galaxy clustering (_e.g._ 9, 10, 11). Furthermore, the analyses are limited to
large, linear scales where the impact of non-linear structure formation is
small. These restrictions result from the fact that standard analyses use
analytic models based on perturbation theory (PT) of large-scale structure
(see 12, 13, for a review). PT struggles to accurately model scales beyond
quasi-linear scales, especially for higher-order clustering statistics (_e.g._
bispectrum). In the (14) PT analysis, for instance, the authors restrict the
power spectrum to $k<0.2\,h/{\rm Mpc}$ and the bispectrum to $k<0.08\,h/{\rm
Mpc}$. In a recent PT analysis, (15) analyzes the bispectrum to
$k<0.23\,h/{\rm Mpc}$; however, they require 33 extra parameters for the
theoretical consistency of their model. PT also cannot be used to model
various newly proposed summary statistics (_e.g._ 16) or to exploit the full
galaxy distribution at the field level.

Another major challenge for current analyses is accurately accounting for
observational systematics. Observations suffer from imperfections in _e.g._
targeting, imaging, completeness that can significantly impact the analysis
(17, 18). Current analyses account for these effects by applying correction
weights to the galaxies. Fiber collisions, for example, prevent galaxy surveys
that use fiber-fed spectrographs (_e.g._ DESI, PFS) from successfully
measuring redshifts from galaxies within some angular scale of one another on
the focal plane. This significantly biases the power spectrum by more than the
amplitude of cosmic variance on scales smaller than $k>0.1\,h/{\rm Mpc}$ (19,
20, 21). To correct for this effect, the weights of the “collided” galaxies
missed by survey are assigned to their nearest angular neighbors (22, 23).
Even for current analyses, these correction weights do not sufficiently
correct the measured power spectrum (20). Furthermore, they are only designed
and demonstrated for the power spectrum.

Meanwhile, additional cosmological information is available on non-linear
scales and in higher-order statistics. Recent studies have accurately
quantified the information content in these regimes using large suites of
simulations. (24) and (25) used the Quijote suite of simulations to
demonstrate that constraints on cosmological parameters,
$\Omega_{m},\Omega_{b},h,n_{s},\sigma_{8}$, improve by a factor of $\sim$2 by
including non-linear scales ($0.2<k<0.5\,h/{\rm Mpc}$) in power spectrum
analyses. Even more improvement comes from including higher-order clustering
information in the bispectrum. Similar forecasts for other summary statistics,
_e.g._ marked power spectrum (26, 27), reconstructed power spectrum (28), skew
spectra (29), wavelet statistics (30), find consistent improvements from
including non-linear scales and higher-order clustering. Despite the growing
evidence of the significant constraining power available in non-linear scales
and higher-order statistics, it cannot be exploited by standard methods with
PT.

Robustly exploiting non-linear and non-Gaussian cosmological information
requires a framework that can both accurately model non-linear structure
formation and account for detailed observational systematics. In this work, we
present SIMulation-Based Inference of Galaxies (SimBIG), a framework for
analyzing galaxy clustering that achieves these requirements by using a
forward modeling approach. Instead of analyzing galaxy clustering using
analytic models, a forward model approach uses simulations that model the full
details of the observations.

In SimBIG, our forward model is based on cosmological $N$-body simulations
that accurately models non-linear structure formation. We also use the halo
occupation framework, which provides a compact and flexible prescription for
connecting the galaxy distribution to the dark matter distribution. Our
forward model also takes advantage of the fact that many observational
systematics can be more easily included in simulations (_e.g._ 31, 32) than
corrected in the observations. With this forward model, we can rigorously
analyze galaxy clustering on non-linear scales and with higher-order
statistics.

To infer the cosmological parameters, our approach does not require sampling
the posterior using an assumed analytic likelihood. We instead use simulation-
based inference (SBI; see 33, for a review). SBI, also known as “likelihood-
free inference”, enables accurate Bayesian inference using forward models
(_e.g._ 34, 35, 36, 37). Moreover, they leverage neural density estimation
from machine learning (_e.g._ 38, 34) to more efficiently infer the posterior
without sampling or making strong assumptions on the functional form of the
likelihood.

In this work, we apply SimBIG to the CMASS galaxy sample observed by the Sloan
Digital Sky Survey SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS; 39,
40). With the main goal of demonstrating the accuracy and potential of SimBIG,
we use the power spectrum as our summary statistic. We present the
cosmological constraints inferred from our analysis and compare them to
previous constraints in the literature. In an accompanying paper (41,
hereafter H22b), we present our forward model in further detail and the mock
challenge that we conduct to rigorously validate the accuracy of SimBIG
cosmological constraints.

Figure 1:  The SimBIG forward model produces simulated galaxy samples with the
same survey geometry and observational systematics as the observed BOSS CMASS
SGC galaxy sample. We present the 3D distribution of the galaxies from three
different viewing angles. The colormap represents the redshift of the
galaxies. In the top set of panels, we present the distribution of galaxies in
the CMASS sample. In the bottom, we present the distribution of a simulated
galaxy sample, generated from our forward model. The SimBIG galaxy samples are
constructed from Quijote $N$-body dark matter simulations using an HOD model
that populates dark matter halos identified using the Rockstar algorithm. The
3D distributions illustrate that our forward model is able to generate galaxy
distributions that are difficult to statistically distinguish from
observations. For more comparisons of the 3D distributions, we refer readers
to .

## Simulation-Based Inference of Galaxies SimBIG

Modern cosmological analyses use Bayesian inference to constrain the posterior
distribution $p(\boldsymbol{\theta}\,|\,\textbf{{x}})$ of cosmological
parameters, $\boldsymbol{\theta}$, given observation x. In standard galaxy
clustering analyses, the posterior is evaluated using Bayes’ rule. The
likelihood is assumed to have a Gaussian functional form and evaluated using
an analytic PT model.

SBI offers an alternative that requires no assumptions on the form of the
likelihood. SBI only requires a forward model, _i.e._ a simulation to generate
mock observations $\textbf{{x}}^{\prime}$ given parameters
$\boldsymbol{\theta}^{\prime}$. It uses a training dataset of simulated pairs
$(\boldsymbol{\theta}^{\prime},\textbf{{x}}^{\prime})$ to estimate the
posterior. SBI has already been successfully applied to a wide range of
inference problems in astronomy and cosmology (42, 43, 44, 35, 45, 36, 46,
47).

In this work, we utilize SBI based on neural density estimation, where a
neural network $q$ with parameters $\boldsymbol{\phi}$ is trained to estimate
$p(\boldsymbol{\theta}\,|\,\textbf{{x}})\approx
q_{\boldsymbol{\phi}}(\boldsymbol{\theta}\,|\,\textbf{{x}})$. In particular,
we use “normalizing flow” models that are capable of accurately estimating
complex distributions (48, 49). Below, we briefly describe our forward model
and SBI framework.

###  0.1 Forward Model

SBI requires a forward model that is capable of generating mock observations
which are statistically indistinguishable from the observations. We start with
high-resolution $N$-body simulations from the Quijote suite (50). These
simulations follow the evolution of $1024^{3}$ cold dark matter (CDM)
particles in a volume of $(1\,h^{-1}{\rm Gpc})^{3}$ from $z=127$ to $z=0.5$
using the TreePM Gadget-III code. They accurately model the clustering of
matter down to non-linear scales beyond $k=0.5\,h/{\rm Mpc}$ (50).

To model the galaxy distribution, we identify gravitationally bound dark
matter halos and populate them with galaxies using a flexible halo occupation
framework. We identify halos using the Rockstar phase-space-based halo finder
(51), which accurately determines the location of halos and resolves their
substructure (52). We then populate the halos using Halo Occupation
Distribution (HOD) models that provide a statistical prescription for
populating halos with galaxies based on halo properties such as their mass and
concentration. In this work, we use a state-of-the-art HOD model that
supplements the standard (53) model with assembly, concentration, and velocity
biases. The extra features of our HOD model add additional flexibility that
recent works suggest may be necessary to describe galaxy clustering (_e.g._
54, 55, 56).

Once we have our galaxy distribution in the simulation box, we apply survey
realism. We remap the box to a cuboid (57) and then cut out the detailed
survey geometry of the BOSS CMASS SGC sample (see Materials and Methods). This
includes masking for bright stars, centerpost, bad field, and collision
priority (17, 40, 18). We apply fiber collisions by first identifying all
pairs of galaxies within an angular scale of $62^{\prime\prime}$ then, for 60%
of the pairs, removing one of the galaxies from the sample. Lastly, we trim
the forward modeled galaxy catalog to match the $0.45<z<0.6$ redshift range
and angular range of the observations.

In total, our forward model has 14 parameters. 5 $\Lambda$CDM cosmological
parameters, $\Omega_{m},\Omega_{b},h,n_{s},\sigma_{8}$, that determine the
matter distribution and 9 HOD parameters that determine the connection between
galaxies and halos. For further details on our forward model, we refer readers
to H22b. In the bottom panels of Fig. 1, we present the three-dimensional
spatial distribution of galaxies in our forward model. We present the angular
distribution of galaxies in our forward model in Fig. 2. The forward model
accurately reproduces the survey geometry and angular footprint of the
observed BOSS sample. For additional comparisons of the 3D distributions of
galaxies in CMASS and our forward model, we refer readers to
111https://youtube.com/playlist?list=PLQk8Faa2x0twK3fgs55ednnHD2vbIzo4z.

Figure 2:  Angular distribution of galaxies from the CMASS sample (top) and a
galaxy sample generated using the SimBIG forward model (bottom). Comparison of
the angular distributions highlights the detailed CMASS angular selection that
we include in our forward model to account for observational systematics.
Figure 3:  Power spectrum, $P_{\ell}(k)$, measured from the simulated galaxy
catalogs constructed using the SimBIG forward model. We present $P_{\ell}(k)$
of 100 out of the total 20,000 catalogs for clarity. In each of the panels, we
plot the monopole, quadrupole, and hexadecapole of the power spectrum
($\ell=0,2,4$). For reference, we include $P_{\ell}(k)$ measured from the BOSS
CMASS SGC galaxy sample (black) with uncertainties estimated from H22b
simulations. $P_{\ell}$ is the most commonly used summary statistic of galaxy
distribution that measures the two-point clustering. We use $P_{\ell}$ in this
work to showcase and validate the SimBIG framework and make detailed
comparisons to previous works in the literature. The $P_{\ell}$ of the SimBIG
catalogs encompass the BOSS $P_{\ell}$ and, thus, provide a sufficiently broad
dataset to conduct SBI.

###  0.2 Training Dataset for Simulation-Based Inference

Using our forward model, we construct 20,000 simulated galaxy catalogs. They
are constructed from 2,000 Quijote $N$-body simulations with 10 different sets
of HOD parameters, sampled from a broad prior. The $N$-body simulations are
arranged in a Latin Hypercube configuration, which therefore imposes uniform
priors on the cosmological parameters that conservatively encompasses the
Planck cosmological constraints (58).

In principle, SimBIG can be directly applied to the full galaxy catalog if the
forward model is capable of accurately modeling observations at all scales.
Even with $N$-body simulations, however, this is not the case due to
limitations on mass and time resolution and inadequacies of halo occupation
models. Instead, we use summary statistics of the galaxy sample, where we can
impose cuts, _e.g.,_ based on physical scales, to which our forward model is
accurate. Since the primary goal of this work is to present and demonstrate
the SimBIG framework, we use the most commonly used summary statistic: the
galaxy power spectrum multipole, $P_{\ell}(k)$. We also include the average
galaxy number density of the sample, $\bar{n}_{g}$.

In this work, we use the redshift-space galaxy power spectrum monopole,
quadrupole, and hexadecapole ($\ell=0,2$, and $4$). We measure $P_{0}$,
$P_{2}$, and $P_{4}$ for each of the simulated galaxy catalogs using the (59)
algorithm. In this work, we impose a conservative $k<k_{\rm max}=0.5\,h/{\rm
Mpc}$ limit on the $P_{\ell}$, based on the convergence of matter clustering
of the Quijote simulations (see H22b for further details). We also measure the
power spectrum for the BOSS CMASS-SGC galaxy sample with the same algorithm.
For the observed $\hat{P}_{\ell}(k)$ we include systematics weights for
redshift failures, stellar density, and seeing conditions, which are effects
not included in our forward model but shown to be successfully accounted for
using the weights (17, 23).

By using $P_{\ell}$, we can compare the constraints inferred using SimBIG with
previous constraints (_e.g._ 10, 11) as further validation of SimBIG. To be
further consistent with previous analyses, we include a nuisance parameter,
$A_{\rm shot}$, that is typically included to account for residual shot noise
contribution (_e.g._ 9, 10, 11). In Fig. 3, we present $P_{\ell}(k)$ of our
forward modeled galaxy catalogs. We randomly select 100 out of the total
20,000 power spectra for clarity. The left, center, and right panels present
the monopole, quadrupole, and hexadecapole. These $P_{\ell}$ measurements and
$\bar{n}_{g}$ serve as the training dataset
$\\{(\boldsymbol{\theta}^{\prime},\textbf{{x}}^{\prime})\\}$ for our SBI
posterior estimation using normalizing flows.

###  0.3 Simulation-Based Inference with Normalizing Flows

Normalizing flow models use an invertible bijective transformation,
$f:\textbf{{z}}\mapsto\boldsymbol{\theta}$, to map a complex target
distribution to a simple base distribution, $\pi(\textbf{{z}})$, that is fast
to evaluate. For SBI, the target distribution is the posterior,
$p(\boldsymbol{\theta}|\textbf{{x}})$ while $\pi(\textbf{{z}})$ is typically a
multivariate Gaussian. The transformation $f$ must be invertible and have a
tractable Jacobian so that we can evaluate the target distribution from
$\pi(\textbf{{z}})$ by change of variables. Since $\pi(\textbf{{z}})$ is easy
to sample and evaluate, we can also easily sample and evaluate the target
distribution. A neural network with parameters $\boldsymbol{\phi}$ is trained
to obtain $f$.

Among the various normalizing flow-based neural density estimators now
available in the literature, we use a Masked Autoregressive Flow (MAF; 34).
MAF combines normalizing flows with an autoregressive design (60), which is
well-suited for estimating conditional probability distributions such as a
posterior. A MAF model is built by stacking multiple Masked Autoencoder for
Distribution Estimation (MADE; 38) models so that it has the autoregressive
structure of MADE models but with additional flexibility to describe complex
probability distributions. We use the MAF implementation of the sbi Python
package (61, 62).

In training, our goal is to determine the parameters, $\boldsymbol{\phi}$, of
our normalizing flow,
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}\,|\,\textbf{{x}})$, so that it
accurately estimates the posterior, $p(\boldsymbol{\theta}\,|\,\textbf{{x}})$.
We can formulate this into an optimization problem of minimizing the Kullback-
Leibler (KL) divergence between
$p(\boldsymbol{\theta},\textbf{{x}})=p(\boldsymbol{\theta}\,|\,\textbf{{x}})p(\textbf{{x}})$
and
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}\,|\,\textbf{{x}})p(\textbf{{x}})$,
which measures the difference between the two distributions.

$\displaystyle\min_{\boldsymbol{\phi}}D_{\rm
KL}\bigl{(}p(\boldsymbol{\theta},\textbf{{x}})\,$
$\displaystyle\parallel\,q_{\boldsymbol{\phi}}(\boldsymbol{\theta}\,|\,\textbf{{x}})p(\textbf{{x}})\bigl{)}$
$\displaystyle=\min_{\boldsymbol{\phi}}\int
p(\boldsymbol{\theta},\textbf{{x}})\log\frac{p(\boldsymbol{\theta}\,|\,\textbf{{x}})}{q_{\boldsymbol{\phi}}(\boldsymbol{\theta}\,|\,\textbf{{x}})}\,{\rm
d}\boldsymbol{\theta}{\rm d}\textbf{{x}}$ (1)
$\displaystyle\approx\min_{\boldsymbol{\phi}}\sum\limits_{i}\log
p(\boldsymbol{\theta}_{i}\,|\,\textbf{{x}}_{i})-\log
q_{\phi}(\boldsymbol{\theta}_{i}\,|\,\textbf{{x}}_{i})$ (2)
$\displaystyle\approx\min_{\boldsymbol{\phi}}\sum\limits_{i}-\log
q_{\phi}(\boldsymbol{\theta}_{i}\,|\,\textbf{{x}}_{i})$ (3)
$\displaystyle\approx\max_{\boldsymbol{\phi}}\sum\limits_{i}\log
q_{\phi}(\boldsymbol{\theta}_{i}\,|\,\textbf{{x}}_{i}).$ (4)

Eq. 2 follows from the fact that the training dataset
$\\{(\boldsymbol{\theta}^{\prime},\textbf{{x}}^{\prime})\\}$ is constructed by
sampling from $p(\boldsymbol{\theta},\textbf{{x}})$ with our forward model.

We split the training data into a training and validation set with a 90/10
split, then maximize Eq. 4 over the training set. We use the Adam optimizer
(63) with a learning rate of $5\times 10^{-4}$. We prevent overfitting by
stopping the training when the log-likelihood (Eq. 4) evaluated on the
validation set fails to increase after 20 epochs. We determine the
architecture of our normalizing flow through experimentation. Our final
trained model has 6 MADE blocks, each with 9 hidden layers and 186 hidden
units. For further details on the training procedure, we refer readers to
H22b. Once trained, we estimate the posterior of our 5 cosmological, 9 HOD
parameters, and 1 nuisance parameter. for the BOSS CMASS SGC $\hat{P_{\ell}}$
and $\hat{\bar{n}}_{g}$, by sampling our normalizing flow
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}\,|\,\hat{\textbf{{x}}})$.

Figure 4:  Left: Posterior of cosmological parameters inferred from $P_{\ell}$
using SimBIG. In the diagonal panels we present the marginalized 1D posterior
of each parameter. The other panels present the 2D posteriors that illustrate
the degeneracies between two parameters. The contours mark the 68 and 95
percentiles. We accurately analyze $P_{\ell}$ down to non-linear regimes,
$k_{\rm max}=0.5\,h/{\rm Mpc}$, by using a simulation-based forward model that
includes observational systematics. Right: We focus on the posteriors of
($\Omega_{m}$, $\sigma_{8}$), the parameters that can be most significantly
constrained by galaxy clustering alone. We derive
$\Omega_{m}=0.292^{+0.055}_{-0.040}$ and $\sigma_{8}=0.812^{+0.067}_{-0.068}$.
Our $\sigma_{8}$ constraints are 27% tighter than the (10) $k_{\rm
max}=0.25\,h/{\rm Mpc}$ PT constraint (orange).

##  1 Results

We present the posterior distribution of the $\Lambda$CDM cosmological
parameters, $\Omega_{m},\Omega_{b},h,n_{s},\sigma_{8}$, inferred from
$P_{\ell}(k)$ using SimBIG in Fig. 4. The posterior is inferred from the BOSS
$P_{\ell}(k)$ down to $k_{\rm max}=0.5\,h/{\rm Mpc}$. The diagonal panels
present the marginalized one-dimensional posteriors for each parameter. The
other panels present marginalized two-dimensional posteriors of different
parameter pairs that highlight parameter degeneracies. We mark the 68 and 95
percentiles of the posteriors with the contours. We infer the posterior of HOD
and nuisance parameters; however, we do not include them in the figure for
clarity. Among the cosmological parameters, the SimBIG posterior significantly
constrains $\Omega_{m}$ and $\sigma_{8}$. This is consistent with previous
works that relied on priors from Big Bang nucleosynthesis or cosmic microwave
background (CMB) experiments for the other parameters, $\Omega_{b}$ and
$n_{s}$ (_e.g._ 10, 11). We infer $\Omega_{m}=0.292^{+0.055}_{-0.040}$ and
$\sigma_{8}=0.812^{+0.067}_{-0.068}$.

In the accompanying paper H22b, we present the validation of the SimBIG
posterior using a suite of 1,500 test simulations. We construct the test suite
using different forward models than the one used for our training data. They
are constructed using different $N$-body simulations, halo finders, and HOD
models. This is to ensure that the cosmological constraints we derive are
independent of the choices and assumptions made in our forward model.

For validation, we conduct a mock challenge where we infer posteriors of the
cosmological parameters for each of the test simulations. Since we know the
true cosmological parameter values of the test simulations, we can access the
accuracy and precision of the inferred posteriors. H22b reveals that SimBIG
produces unbiased posteriors. On the other hand, the posteriors are
conservative, _i.e._ they are broader than the true posterior. This is due to
the limited number of $N$-body simulations used to construct our training
dataset, which makes the estimate of the KL divergence (Eq. 2) noisy. Our
constraint on $\Omega_{m}$ is particularly conservative. Additional $N$-body
simulations would significantly improve the precision of our posteriors.

Despite the fact that they are conservative, the $\sigma_{8}$ posterior from
SimBIG is significantly more precise than constraints from previous works.
(10) analyzed the $P_{\ell}$ of the BOSS CMASS galaxy sample using the PT
approach with an analytic model based on effective field theory. For the CMASS
SGC sample, with uniform priors on the cosmological parameters, and with
$k_{\rm max}=0.25\,h/{\rm Mpc}$, (10) inferred
$\sigma_{8}=0.719^{+0.100}_{-0.085}$. With SimBIG, we improve $\sigma_{8}$
constraints by 27% over the standard galaxy clustering analysis. We emphasize
that this improvement is roughly equivalent to analyzing a galaxy survey
$\sim$60% larger than the original survey using PT.

Recently, (11) also analyzed the $P_{\ell}$ of BOSS CMASS sample but using a
theoretical model based on a halo power spectrum emulator. Instead of using a
galaxy bias scheme used by PT to connect the galaxy and matter distributions,
(11) used halo power spectra predicted by an emulator and a halo occupation
framework, similar to the HOD model in our forward model. We note that while
the halo power spectrum emulator is trained using simulations, the approach in
(11) does not forward model observational systematics. They also make the same
assumptions on the form of the likelihood as PT analyses for their inference.
For the CMASS SGC sample, with uniform priors on all cosmological parameters,
and with $k_{\rm max}=0.25\,h/{\rm Mpc}$, (11) inferred
$\sigma_{8}=0.790^{+0.083}_{-0.072}$. The (11) constraints are tighter than
the (10) PT constraints because the halo occupation model provides a more
compact framework for modeling galaxies. Nevertheless, with SimBIG, we improve
on their $\sigma_{8}$ constraints by 13%.

SimBIG produces significantly tighter constraints on $\sigma_{8}$ because we
are able to accurately extract cosmological information available on small,
non-linear, scales. With our forward modeling approach, we can accurately
model non-linear clustering and robustly account for observational systematics
down to $k_{\rm max}=0.5\,h/{\rm Mpc}$. In both (10) and (11), they restrict
their analysis to $k_{\rm max}<0.25\,h/{\rm Mpc}$ due to the limitations of
their analyses on smaller scales.

To further verify that our improvement comes from constraining power at
$k>0.25\,h/{\rm Mpc}$, we analyze $P_{\ell}$ to $k_{\rm max}=0.25\,h/{\rm
Mpc}$ using SimBIG. In Fig. 5, we present the SimBIG $k_{\rm max}=0.25\,h/{\rm
Mpc}$ posterior (blue) along with the posteriors from (10, orange) and (11,
green). We focus our comparison on $\Omega_{m}$ and $\sigma_{8}$, the
cosmological parameters that can be most competitively constrained with galaxy
clustering alone. The contours again represent the 68 and 95 percentiles. We
find overall good agreement among the posteriors. All of the posteriors are
statistically consistent with each other. For $\sigma_{8}$, our $k_{\rm
max}=0.25\,h/{\rm Mpc}$ places a $\sigma_{8}=0.861^{+0.070}_{-0.091}$
constraint. This is significantly broader than our $k_{\rm max}=0.5\,h/{\rm
Mpc}$ constraint, which demonstrates that the constraining power is in fact
from non-linear scales. Furthermore, the precision of the $k_{\rm
max}=0.25\,h/{\rm Mpc}$ SimBIG constraint is in excellent agreement with the
(11) constraint. This serves as further validation of SimBIG, since (11) uses
a similar halo occupation framework to model the power spectrum.

For $\Omega_{m}$, we infer broader posteriors than (10) and (11). As we
discuss in H22b, this is due to the fact that the SimBIG normalizing flow is
trained using a limited number of simulations. We use 20,000 forward modeled
simulations; however, they are constructed from 2000 $N$-body simulations with
different values of cosmological parameters. To estimate the expected
improvement in the $\Omega_{m}$ constraints, we use the posterior “re-
calibration” procedure from (64). The re-calibration uses the posteriors
inferred for the test simulations and their true parameter values. We
calculate the local probability integral transform (65), a diagnostic of the
inferred posteriors, and use this quantity to derive a weighting scheme that
corrects the posteriors so that it matches the true posterior of the test
simulations.

The re-calibration uses test simulations, so we do not use it for inference.
However, it provides a bound for the SimBIG constraints if we were to have
sufficient training simulations. The re-calibrated posterior constrains
$\Omega_{m}=0.284^{+0.021}_{-0.017}$. For reference, we mark the re-calibrated
$\Omega_{m}$ constraint (black dotted) in Fig. 5. The re-calibrated
$\Omega_{m}$ is in good agreement with both the (10) and (11) constraints. It
is significantly tighter than the original SimBIG constraint and illustrates
that additional training simulations would significantly improve the precision
of the SimBIG $\Omega_{m}$ constraints.

Based on our $k_{\rm max}=0.5\,h/{\rm Mpc}$ posterior, we infer
$S_{8}=\sigma_{8}\sqrt{\Omega_{m}/0.3}=0.802^{+0.102}_{-0.092}$ (and
$0.797^{+0.078}_{-0.076}$ for our re-calibrated posterior). Multiple recent
large-scale structure studies have reported a “$S_{8}$ tension” with
constraints from the (58) CMB analysis. They find significantly lower values
of $S_{8}$ than Planck (66, 67, 68, 69, 70, 71, 72, 73). PT analyses of BOSS
also infer relatively low values of $S_{8}$. (10), for instance, infers
$S_{8}=0.737^{+0.110}_{-0.092}$. This $S_{8}$ tension between constraints from
large-scale structure and CMB analyses has motivated a number of works to
explore modifications of the standard $\Lambda$CDM cosmological model (_e.g._
74, 75, 76, 77). We do not find a significant $S_{8}$ tension with the Planck
constraints (58). However, given the statistical precision of our $S_{8}$
constraint, we refrain from more detailed comparison and discussion.

Figure 5:  Comparison of the ($\Omega_{m}$, $\sigma_{8}$) posteriors inferred
from the $P_{\ell}$ CMASS-SGC analysis for $k_{\rm max}=0.25\,h/{\rm Mpc}$
from SimBIG (blue), the (10) PT approach (orange), and the (11) emulator
approach (green). The contours represent the 68 and 95 percentiles. We find
overall good agreement among the posteriors. For $\Omega_{m}$, SimBIG infers
consistent but broader posterior due to the fact that we use a limited number
of simulations. We estimate the expected $\Omega_{m}$ constraints with more
simulations using posterior “recalibration” (black dotted). For $\sigma_{8}$,
both SimBIG and (11) derive tighter constraints than (10) due to the fact that
we use halo occupation instead of a galaxy bias prescription. Meanwhile, the
precision of our $\sigma_{8}$ constraint is in excellent agreement with (11),
which uses a similar halo occupation model. The consistency of the $k_{\rm
max}=0.25\,h/{\rm Mpc}$ posteriors demonstrate that the improvements in the
$k_{\rm max}=0.5\,h/{\rm Mpc}$ constraints come from additional cosmological
information in the non-linear regime that SimBIG can robustly extract.

##  2 Conclusions

We present SimBIG, a forward modeling framework for analyzing galaxy
clustering using SBI. As a demonstration of the framework, we apply it to the
BOSS CMASS SGC, a galaxy sample at $z\sim 0.5$. We analyze the galaxy power
spectrum multipoles ($P_{\ell}$), the most commonly used summary statistic of
the galaxy spatial distribution, to showcase and validate the SimBIG
framework.

SimBIG utilizes a full forward model of the CMASS sample, unlike standard
approaches that use analytic models of the summary statistic. The forward
model is based on high-resolution Quijote $N$-body simulations that can
accurately model the matter distribution on small scales. It uses halo
modeling and a state-of-the-art HOD model with assembly, concentration, and
velocity biases that provide a flexible mapping between the matter and galaxy
distributions. The forward model also includes realistic observational
systematics such as survey geometry and fiber collisions. With this forward
modeling approach, we can leverage the predictive power of simulations to
analyze small, non-linear, scales as well as higher-order clustering. It also
provides a framework to account for systematics for any summary statistic.

Using the forward model, we construct 20,000 simulated CMASS-like samples that
span a wide range of cosmological and HOD parameters. We measure $P_{\ell}$
and $\bar{n}_{g}$ for each of these samples and use the measurements as the
training dataset for SBI. To estimate the posterior, we use neural density
estimation based on normalizing flows. Using the training dataset, we train
our normalizing flows by minimizing the KL divergence between its posterior
estimate and the true posterior. Once trained, we apply our normalizing flow
to the observed summary statistics to infer the posterior of 5 cosmological, 9
HOD, and 1 nuisance parameter.

Focusing on the cosmological parameters, we derive significant constraints on:
$\Omega_{m}=0.316^{+0.040}_{-0.036}$ and
$\sigma_{8}=0.668^{+0.0324}_{-0.0284}$. Our $\sigma_{8}$ constraints are 27%
tighter than the (10) constraints using a standard PT approach on the same
galaxy sample. This improvement is roughly equivalent to increasing the volume
of the galaxy survey by $\sim$60% for a standard PT analysis. The SimBIG
constraints are inferred from $P_{\ell}$ out to $k_{\rm max}=0.5\,h/{\rm Mpc}$
while the PT constraints are limited to $k_{\rm max}=0.25\,h/{\rm Mpc}$ The
improvement is driven by the additional cosmological information on non-linear
scales that SimBIG can robustly extract.

We also infer the posterior using SimBIG from $P_{\ell}$ with $k_{\rm
max}=0.25\,h/{\rm Mpc}$ and compare it to posteriors in the literature. In
particular, the SimBIG $\sigma_{8}$ constraint for $k_{\rm max}=0.25\,h/{\rm
Mpc}$ are in excellent agreement with the constraint from the recent halo
model based emulator analysis of (11). Since they use a similar halo
occupation framework as SimBIG, this comparison firmly verifies the robustness
of our constraints. The comparison also confirms that the improvement in
SimBIG $k_{\rm max}=0.5\,h/{\rm Mpc}$ constraints come from the non-linear
regime. In the accompanying H22b, we present additional tests of SimBIG
through a mock challenge. The tests use a suite of 1,500 test simulations
constructed with different forward models to demonstrate that SimBIG produces
unbiased cosmological constraints. H22b also presents further details on our
forward model and discusses posterior constraints on HOD parameters.

SimBIG can also extract higher-order cosmological information. Standard galaxy
clustering analyses primarily focus on two-point clustering statistics.
Analyses of higher-order statistics have been limited to, _e.g._ the
bispectrum (78, 14, 15), and even these analyses extract only limited
cosmological information beyond linear scales. In subsequent work, we will use
SimBIG to analyse the BOSS CMASS galaxies using higher-order statistics (the
bispectrum) and non-standard observables that contain additional cosmological
information: _e.g._ marked power spectrum, skew spectra, void probability
functions, and wavelet- scattering-like statistics. We will also apply SimBIG
to analyze field-level summary statistics that capture all the information in
the galaxy field using convolutional and graph neural networks.

SimBIG can also be extended to upcoming spectroscopic galaxy surveys observed
using DESI, PFS, Euclid, and Roman will probe unprecedented cosmic volumes
over the next decade. They will produce the largest and most detailed three-
dimensional maps of galaxies in the Universe. These surveys are already
expected to provide the most precise constraints on cosmological parameters
using standard analyses. SimBIG can further exploit the statistical power of
these surveys to place even tighter constraints on cosmological parameters and
produce the most stringent tests of the standard $\Lambda$CDM cosmological
model and beyond.

### Observations: SDSS-III BOSS

In this work, we analyze observations from the Sloan Digital Sky Survey SDSS-
III (39, 40) Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 12.
In particular, we use the CMASS galaxy sample, which selects high stellar mass
Luminous Red Galaxies (LRGs) over the redshift $0.43<z<0.7$ (79). We restrict
our analysis to CMASS galaxies in the Southern Galactic Cap (SGC) and impose a
redshift cut of $0.45<z<0.6$ and the following angular cuts: ${\rm Dec}>-6$
and $-25<{\rm RA}<28$. In upper panels in Fig. 1, we present the three-
dimensional distributions of our CMASS SGC galaxy sample at three different
viewing angles. We also present the angular distribution of the sample in Fig.
2. In total, our CMASS SGC galaxy sample contains 109,636 galaxies.

It is a pleasure to thank Peter Melchior, Uros̆ Seljak, David Spergel, Licia
Verde, and Benjamin D. Wandelt for valuable discussions. We also thank Mikhail
M. Ivanov and Yosuke Kobayashi for providing us with the posteriors used in
the comparison. This work was supported by the AI Accelerator program of the
Schmidt Futures Foundation. This work was also supported by NASA ROSES grant
12-EUCLID12-0004 NASA grant 15-WFIRST15-0008. JH has received funding from the
European Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie grant agreement No 101025187. AMD acknowledges funding from
Tomalla Foundation for Research in Gravity and Boninchi Foundation.

## References

  * (1) D Collaboration, et al., The DESI Experiment Part I: Science,Targeting, and Survey Design.  arXiv:1611.00036 [astro-ph] (2016).
  * (2) D Collaboration, et al., The DESI Experiment Part II: Instrument Design.  arXiv:1611.00037 [astro-ph] (2016).
  * (3) B Abareshi, et al., Overview of the Instrumentation for the Dark Energy Spectroscopic Instrument (2022).
  * (4) M Takada, et al., Extragalactic science, cosmology, and Galactic archaeology with the Subaru Prime Focus Spectrograph.  Publications of the Astronomical Society of Japan 66, R1 (2014).
  * (5) N Tamura, et al., Prime Focus Spectrograph (PFS) for the Subaru telescope: Overview, recent progress, and future perspectives in Ground-Based and Airborne Instrumentation for Astronomy VI.  (eprint: arXiv:1608.01075), Vol. 9908, p. 99081M (2016).
  * (6) R Laureijs, et al., Euclid Definition Study Report.  arXiv e-prints p. arXiv:1110.3193 (2011).
  * (7) D Spergel, et al., Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report (2015).
  * (8) Y Wang, et al., The High Latitude Spectroscopic Survey on the Nancy Grace Roman Space Telescope.  The Astrophysical Journal 928, 1 (2022).
  * (9) F Beutler, et al., The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: Anisotropic galaxy clustering in Fourier space.  Monthly Notices of the Royal Astronomical Society 466, 2242–2260 (2017).
  * (10) MM Ivanov, M Simonović, M Zaldarriaga, Cosmological parameters from the BOSS galaxy power spectrum.  Journal of Cosmology and Astroparticle Physics 2020, 042 (2020).
  * (11) Y Kobayashi, T Nishimichi, M Takada, H Miyatake, Full-shape cosmology analysis of SDSS-III BOSS galaxy power spectrum using emulator-based halo model: A $5\%$ determination of $\sigma_8$.  arXiv:2110.06969 [astro-ph] (2021).
  * (12) F Bernardeau, S Colombi, E Gaztanaga, R Scoccimarro, Large-Scale Structure of the Universe and Cosmological Perturbation Theory.  Physics Reports 367, 1–248 (2002).
  * (13) V Desjacques, D Jeong, F Schmidt, Large-Scale Galaxy Bias.  arXiv:1611.09787 [astro-ph, physics:gr-qc, physics:hep-ph] (2016).
  * (14) OHE Philcox, MM Ivanov, The BOSS DR12 Full-Shape Cosmology: $\Lambda$CDM Constraints from the Large-Scale Galaxy Power Spectrum and Bispectrum Monopole.  arXiv:2112.04515 [astro-ph, physics:hep-ex] (2021).
  * (15) G D’Amico, Y Donath, M Lewandowski, L Senatore, P Zhang, The BOSS bispectrum analysis at one loop from the Effective Field Theory of Large-Scale Structure (2022).
  * (16) K Naidoo, E Massara, O Lahav, Cosmology and neutrino mass with the Minimum Spanning Tree.  Monthly Notices of the Royal Astronomical Society (2022).
  * (17) AJ Ross, et al., The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: Analysis of potential systematics.  Monthly Notices of the Royal Astronomical Society 424, 564–590 (2012).
  * (18) AJ Ross, et al., The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: Observational systematics and baryon acoustic oscillations in the correlation function.  Monthly Notices of the Royal Astronomical Society 464, 1168–1191 (2017).
  * (19) H Guo, I Zehavi, Z Zheng, A New Method to Correct for Fiber Collisions in Galaxy Two-point Statistics.  The Astrophysical Journal 756, 127 (2012).
  * (20) C Hahn, R Scoccimarro, MR Blanton, JL Tinker, SA Rodríguez-Torres, The Effect of Fiber Collisions on the Galaxy Power Spectrum Multipoles.  Monthly Notices of the Royal Astronomical Society 467, 1940–1956 (2017).
  * (21) D Bianchi, et al., Unbiased clustering estimates with the DESI fibre assignment.  Monthly Notices of the Royal Astronomical Society 481, 2338–2348 (2018).
  * (22) I Zehavi, et al., Galaxy Clustering in Early Sloan Digital Sky Survey Redshift Data.  The Astrophysical Journal 571, 172–190 (2002).
  * (23) L Anderson, et al., The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: Baryon acoustic oscillations in the Data Releases 10 and 11 Galaxy samples.  Monthly Notices of the Royal Astronomical Society 441, 24–62 (2014).
  * (24) C Hahn, F Villaescusa-Navarro, E Castorina, R Scoccimarro, Constraining M$\nu$ with the bispectrum. Part I. Breaking parameter degeneracies.  Journal of Cosmology and Astroparticle Physics 03, 040 (2020).
  * (25) C Hahn, F Villaescusa-Navarro, Constraining M$\nu$ with the bispectrum. Part II. The information content of the galaxy bispectrum monopole.  Journal of Cosmology and Astroparticle Physics 2021, 029 (2021).
  * (26) E Massara, F Villaescusa-Navarro, S Ho, N Dalal, DN Spergel, Using the Marked Power Spectrum to Detect the Signature of Neutrinos in Large-Scale Structure.  arXiv:2001.11024 [astro-ph] (2020).
  * (27) E Massara, et al., Cosmological Information in the Marked Power Spectrum of the Galaxy Field (2022).
  * (28) Y Wang, et al., Extracting high-order cosmological information in galaxy surveys with power spectra (2022).
  * (29) J Hou, A Moradinezhad Dizgah, C Hahn, E Massara, Cosmological Information in Skew Spectra of Biased Tracers in Redshift Space.  arXiv e-prints p. arXiv:2210.12743 (2022).
  * (30) M Eickenberg, et al., Wavelet Moments for Cosmological Parameter Estimation (2022).
  * (31) SA Rodríguez-Torres, et al., The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: Modelling the clustering and halo occupation distribution of BOSS CMASS galaxies in the Final Data Release.  Monthly Notices of the Royal Astronomical Society 460, 1173–1187 (2016).
  * (32) G Rossi, et al., The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: N-body Mock Challenge for Galaxy Clustering Measurements.  Monthly Notices of the Royal Astronomical Society 505, 377–407 (2021).
  * (33) K Cranmer, J Brehmer, G Louppe, The frontier of simulation-based inference.  Proceedings of the National Academy of Sciences 117, 30055–30062 (2020).
  * (34) G Papamakarios, T Pavlakou, I Murray, Masked Autoregressive Flow for Density Estimation.  arXiv e-prints 1705, arXiv:1705.07057 (2017).
  * (35) J Alsing, T Charnock, S Feeney, B Wandelt, Fast likelihood-free cosmology with neural density estimators and active learning.  Monthly Notices of the Royal Astronomical Society 488, 4440–4458 (2019).
  * (36) N Jeffrey, J Alsing, F Lanusse, Likelihood-free inference with neural compression of DES SV weak lensing map statistics.  Monthly Notices of the Royal Astronomical Society 501, 954–969 (2021).
  * (37) L Tortorelli, et al., The PAU Survey: Measurement of Narrow-band galaxy properties with Approximate Bayesian Computation.  arXiv:2106.02651 [astro-ph] (2021).
  * (38) M Germain, K Gregor, I Murray, H Larochelle, MADE: Masked Autoencoder for Distribution Estimation.  Proceedings of the 32nd International Conference on Machine Learning 37, 881–889 (2015).
  * (39) DJ Eisenstein, et al., SDSS-III: Massive Spectroscopic Surveys of the Distant Universe, the Milky Way, and Extra-Solar Planetary Systems.  The Astronomical Journal 142, 72 (2011).
  * (40) KS Dawson, et al., The Baryon Oscillation Spectroscopic Survey of SDSS-III.  The Astronomical Journal 145, 10 (2013).
  * (41) C Hahn, et al., SimBIG: Mock Challenge for a Forward Modeling Approach to Galaxy Clustering.  arXiv (2022).
  * (42) E Cameron, AN Pettitt, Approximate Bayesian Computation for astronomical model analysis: A case study in galaxy demographics and morphological transformation at high redshift.  Monthly Notices of the Royal Astronomical Society 425, 44–65 (2012).
  * (43) A Weyant, C Schafer, WM Wood-Vasey, Likelihood-free Cosmological Inference with Type Ia Supernovae: Approximate Bayesian Computation for a Complete Treatment of Uncertainty.  The Astrophysical Journal 764, 116 (2013).
  * (44) C Hahn, et al., Approximate Bayesian Computation in Large Scale Structure: Constraining the galaxy-halo connection.  Monthly Notices of the Royal Astronomical Society 469, 2791–2805 (2017).
  * (45) D Huppenkothen, M Bachetti, Accurate X-ray Timing in the Presence of Systematic Biases With Simulation-Based Inference (2021).
  * (46) K Zhang, et al., Real-time Likelihood-free Inference of Roman Binary Microlensing Events with Amortized Neural Posterior Estimation.  Astronomical Journal 161, 262 (2021).
  * (47) C Hahn, P Melchior, Accelerated Bayesian SED Modeling using Amortized Neural Posterior Estimation (2022).
  * (48) EG Tabak, E Vanden-Eijnden, Density estimation by dual ascent of the log-likelihood.  Communications in Mathematical Sciences 8, 217–233 (2010).
  * (49) EG Tabak, CV Turner, A Family of Nonparametric Density Estimation Algorithms.  Communications on Pure and Applied Mathematics 66, 145–164 (2013).
  * (50) F Villaescusa-Navarro, et al., The Quijote Simulations.  The Astrophysical Journal Supplement Series 250, 2 (2020).
  * (51) PS Behroozi, RH Wechsler, HY Wu, The ROCKSTAR Phase-space Temporal Halo Finder and the Velocity Offsets of Cluster Cores.  The Astrophysical Journal 762, 109 (2013).
  * (52) A Knebe, et al., Haloes gone MAD: The Halo-Finder Comparison Project.  Monthly Notices of the Royal Astronomical Society 415, 2293–2318 (2011).
  * (53) Z Zheng, AL Coil, I Zehavi, Galaxy Evolution from Halo Occupation Distribution Modeling of DEEP2 and SDSS Galaxy Clustering.  The Astrophysical Journal 667, 760–779 (2007).
  * (54) AR Zentner, A Hearin, FC van den Bosch, JU Lange, A Villarreal, Constraints on Assembly Bias from Galaxy Clustering.  arXiv:1606.07817 [astro-ph] (2016).
  * (55) M Vakili, C Hahn, How Are Galaxies Assigned to Halos? Searching for Assembly Bias in the SDSS Galaxy Clustering.  The Astrophysical Journal 872, 115 (2019).
  * (56) B Hadzhiyska, et al., Galaxy assembly bias and large-scale distribution: A comparison between IllustrisTNG and a semi-analytic model.  Monthly Notices of the Royal Astronomical Society 508, 698–718 (2021).
  * (57) J Carlson, M White, Embedding realistic surveys in simulations through volume remapping.  The Astrophysical Journal Supplement Series 190, 311–314 (2010).
  * (58) Planck Collaboration, et al., Planck 2018 results. VI. Cosmological parameters.  arXiv:1807.06209 [astro-ph] (2018).
  * (59) N Hand, Y Li, Z Slepian, U Seljak, An optimal FFT-based anisotropic power spectrum estimator.  Journal of Cosmology and Astro-Particle Physics 07, 002 (2017).
  * (60) B Uria, MA Côté, K Gregor, I Murray, H Larochelle, Neural Autoregressive Distribution Estimation.  arXiv:1605.02226 [cs] (2016).
  * (61) DS Greenberg, M Nonnenmacher, JH Macke, Automatic Posterior Transformation for Likelihood-Free Inference (2019).
  * (62) A Tejero-Cantero, et al., Sbi: A toolkit for simulation-based inference.  Journal of Open Source Software 5, 2505 (2020).
  * (63) DP Kingma, J Ba, Adam: A Method for Stochastic Optimization.  arXiv:1412.6980 [cs] (2017).
  * (64) B Dey, et al., Calibrated Predictive Distributions via Diagnostics for Conditional Coverage (2022).
  * (65) D Zhao, N Dalmasso, R Izbicki, AB Lee, Diagnostics for conditional density models and Bayesian inference algorithms in Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence.  (PMLR), pp. 1830–1840 (2021).
  * (66) M Cacciato, FC van den Bosch, S More, H Mo, X Yang, Cosmological constraints from a combination of galaxy clustering and lensing - III. Application to SDSS data.  Monthly Notices of the Royal Astronomical Society 430, 767–786 (2013).
  * (67) R Mandelbaum, et al., Cosmological parameter constraints from galaxy-galaxy lensing and galaxy clustering with the SDSS DR7.  Monthly Notices of the Royal Astronomical Society 432, 1544–1575 (2013).
  * (68) A Leauthaud, et al., Lensing is low: Cosmology, galaxy formation or new physics?  Monthly Notices of the Royal Astronomical Society 467, 3024–3047 (2017).
  * (69) C Hikage, et al., Cosmology from cosmic shear power spectra with Subaru Hyper Suprime-Cam first-year data.  Publications of the Astronomical Society of Japan 71, 43 (2019).
  * (70) M Asgari, et al., KiDS-1000 cosmology: Cosmic shear constraints and comparison between two point statistics.  Astronomy & Astrophysics, Volume 645, id.A104, $<$NUMPAGES$>$31$<$/NUMPAGES$>$ pp. 645, A104 (2021).
  * (71) A Krolewski, S Ferraro, M White, Cosmological constraints from unWISE and Planck CMB lensing tomography.  Journal of Cosmology and Astroparticle Physics 2021, 028 (2021).
  * (72) A Amon, et al., Dark Energy Survey Year 3 results: Cosmology from cosmic shear and robustness to data calibration.  Physical Review D 105, 023514 (2022).
  * (73) JU Lange, et al., Five per cent measurements of the growth rate from simulation-based modelling of redshift-space clustering in BOSS LOWZ.  Monthly Notices of the Royal Astronomical Society 509, 1779–1804 (2022).
  * (74) PD Meerburg, Alleviating the tension at low $l$ through axion monodromy.  Physical Review D 90, 063529 (2014).
  * (75) A Chudaykin, D Gorbunov, I Tkachev, Dark matter component decaying after recombination: Sensitivity to baryon acoustic oscillation and redshift space distortion probes.  Physical Review D 97, 083508 (2018).
  * (76) E Di Valentino, A Melchiorri, O Mena, S Vagnozzi, Nonminimal dark sector physics and cosmological tensions.  Physical Review D 101, 063502 (2020).
  * (77) GF Abellán, R Murgia, V Poulin, J Lavalle, Implications of the S8 tension for decaying dark matter with warm decay products.  Physical Review D 105, 063525 (2022).
  * (78) H Gil-Marín, et al., The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the power spectrum and bispectrum of the DR12 BOSS galaxies.  Monthly Notices of the Royal Astronomical Society 465, 1757–1788 (2017).
  * (79) B Reid, et al., SDSS-III Baryon Oscillation Spectroscopic Survey Data Release 12: Galaxy target selection and large-scale structure catalogues.  Monthly Notices of the Royal Astronomical Society 455, 1553–1573 (2016).