  The non-Gaussisan spatial distribution of galaxies traces the large-scale
structure of the Universe and therefore constitutes a prime observable to
constrain cosmological parameters. We conduct Bayesian inference of the
$\Lambda$CDM parameters $\Omega_m$, $\Omega_b$, $h$, $n_s$, and $\sigma_8$ from
the BOSS CMASS galaxy sample by combining the wavelet scattering transform
(WST) with a simulation-based inference approach enabled by the ${\rm
S{\scriptsize IM}BIG}$ forward model. We design a set of reduced WST statistics
that leverage symmetries of redshift-space data. Posterior distributions are
estimated with a conditional normalizing flow trained on 20,000 simulated ${\rm
S{\scriptsize IM}BIG}$ galaxy catalogs with survey realism. We assess the
accuracy of the posterior estimates using simulation-based calibration and
quantify generalization and robustness to the change of forward model using a
suite of 2,000 test simulations. When probing scales down to $k_{\rm
max}=0.5~h/\text{Mpc}$, we are able to derive accurate posterior estimates that
are robust to the change of forward model for all parameters, except
$\sigma_8$. We mitigate the robustness issues with $\sigma_8$ by removing the
WST coefficients that probe scales smaller than $k \sim 0.3~h/\text{Mpc}$.
Applied to the BOSS CMASS sample, our WST analysis yields seemingly improved
constraints obtained from a standard PT-based power spectrum analysis with
$k_{\rm max}=0.25~h/\text{Mpc}$ for all parameters except $h$. However, we
still raise concerns on these results. The observational predictions
significantly vary across different normalizing flow architectures, which we
interpret as a form of model misspecification. This highlights a key challenge
for forward modeling approaches when using summary statistics that are
sensitive to detailed model-specific or observational imprints on galaxy
clustering.
