ABSTRACT Observations of the cosmic 21-cm power spectrum (PS) are starting to enable precision Bayesian inference of galaxy properties and physical cosmology, during the first billion years of our Universe. Here we investigate the impact of common approximations about the likelihood used in such inferences, including: (i) assuming a Gaussian functional form; (ii) estimating the mean from a single realization; and (iii) estimating the (co)variance at a single point in parameter space. We compare ‘classical’ inference that uses an explicit likelihood with simulation-based inference (SBI) that estimates the likelihood from a training set. Our forward models include: (i) realizations of the cosmic 21-cm signal computed with 21cmFAST by varying ultraviolet (UV) and X-ray galaxy parameters together with the initial conditions; (ii) realizations of the telescope noise corresponding to a $1000 \, \mathrm{h}$ integration with the low-frequency component of the Square Kilometre Array (SKA1-Low); and (iii) the excision of Fourier modes corresponding to a foreground-dominated horizon ‘wedge’. We find that the 1D PS likelihood is well described by a Gaussian accounting for covariances between wave modes and redshift bins (higher order correlations are small). However, common approaches of estimating the forward-modelled mean and (co)variance from a random realization or at a single point in parameter space result in biased and overconstrained posteriors. Our best results come from using SBI to fit a non-Gaussian likelihood with a Gaussian mixture neural density estimator. Such SBI can be performed with up to an order of magnitude fewer simulations than classical, explicit likelihood inference. Thus SBI provides accurate posteriors at a comparably low computational cost.