  Type Ia supernovae (SNae Ia), standardisable 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 datasets, 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
marginalised, 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 SNae Ia.
With minimal additional effort, we train a network to infer simultaneously the
O(100 000) latent parameters of the supernovae (e.g. absolute brightnesses). In
addition, we describe and apply a procedure that utilises local amortisation 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.
