Type Ia supernovae (SN Ia) are faraway exploding stars used as ``standardisable candles'' to determine cosmological distances, measure the accelerating expansion of the Universe, and constrain the properties of dark energy. Inferring peak luminosities of supernovae from distance-independent observables, such as the shapes and colours of their multi-wavelength time series, underpins the evidence for cosmic acceleration. Current and future optical SN surveys aim to determine the physical nature of the mysterious dark energy driving the acceleration. However, these efforts are now limited by systematic, rather than statistical, errors. Hierarchical Bayes offers a principled approach to coherently modelling the population and individual SN, and the probabilistic generative processes underlying the observable data, including multiple latent physical effects and uncertainties, such as measurement error, host galaxy dust, and intrinsic SN variations correlated across time and wavelength. I will describe applications to the statistical modelling of SN Ia optical-to-infrared multi-wavelength time series data to optimise cosmological inferences.