Abstract: Type Ia supernovae (SNe Ia) are thermonuclear exploding stars that can be used to put constrains on the nature of our Universe. One challenge with population analyses of SNe Ia is Malmquist bias, where we preferentially observe the brighter SNe due to limitations of our telescopes. If untreated, this bias can propagate through to our posteriors on cosmological parameters. In this paper we develop a novel technique of using a normalising flow to learn the non-analytical likelihood of observing a SN Ia for a given survey from simulations, that is independent of any cosmological model. The learnt likelihood is then used in a hierarchical Bayesian model with Hamiltonian Monte Carlo sampling to put constraints on different sets of cosmological parameters conditioned on our observed data. We verify this technique on toy model simulations finding excellent agreement with analytically-derived posteriors to within $1 \sigma$.