  We present the first fully simulation-based hierarchical analysis of the
light curves of a population of low-redshift type Ia supernovae (SNae Ia). Our
hardware-accelerated forward model, released in the Python package slicsim,
includes stochastic variations of each SN's spectral flux distribution (based
on the pre-trained BayeSN model), extinction from dust in the host and in the
Milky Way, redshift, and realistic instrumental noise. By utilising truncated
marginal neural ratio estimation (TMNRE), a neural network-enabled
simulation-based inference technique, we implicitly marginalise over 4000
latent variables (for a set of $\approx 100$ SNae Ia) to efficiently infer SN
Ia absolute magnitudes and host-galaxy dust properties at the population level
while also constraining the parameters of individual objects. Amortisation of
the inference procedure allows us to obtain coverage guarantees for our results
through Bayesian validation and frequentist calibration. Furthermore, we show a
detailed comparison to full likelihood-based inference, implemented through
Hamiltonian Monte Carlo, on simulated data and then apply TMNRE to the light
curves of 86 SNae Ia from the Carnegie Supernova Project, deriving marginal
posteriors in excellent agreement with previous work. Given its ability to
accommodate arbitrarily complex extensions to the forward model -- e.g.
different populations based on host properties, redshift evolution, complicated
photometric redshift estimates, selection effects, and non-Ia contamination --
without significant modifications to the inference procedure, TMNRE has the
potential to become the tool of choice for cosmological parameter inference
from future, large SN Ia samples.
