Type Ia supernovae (SNe Ia), standardizable candles that allow tracing the expansion history of the Universe, are instrumental in constraining cosmological parameters, particularly dark energy. State-of-the-art likelihood-based analyses scale poorly to future large data sets, are limited to simplified probabilistic descriptions, and must explicitly sample a high-dimensional latent posterior to infer the few parameters of interest, which makes them inefficient. Marginal likelihood-free inference, on the other hand, is based on forward simulations of data, and thus can fully account for complicated redshift uncertainties, contamination from non-SN Ia sources, selection effects, and a realistic instrumental model. All latent parameters, including instrumental and survey-related ones, per object and population-level properties, are implicitly marginalized, while the cosmological parameters of interest are inferred directly. As a proof of concept, we apply truncated marginal neural ratio estimation (TMNRE), a form of marginal likelihood-free inference, to BAHAMAS, a Bayesian hierarchical model for SALT parameters. We verify that TMNRE produces unbiased and precise posteriors for cosmological parameters from up to 100 000 SNe Ia. With minimal additional effort, we train a network to infer simultaneously the ∼100 000 latent parameters of the supernovae (e.g. absolute brightnesses). In addition, we describe and apply a procedure that utilizes local amortization of the inference to convert the approximate Bayesian posteriors into frequentist confidence regions with exact coverage. Finally, we discuss the planned improvements to the model that are enabled by using a likelihood-free inference framework, like selection effects and non-Ia contamination.
