Detecting primordial non-Gaussianity on mildly non-linear scales requires precise modelling of late-time structure formation. Accurately predicting the impact of non-linear gravitational collapse, non-linear tracer bias and baryonic physics on the variance and higher order moments of the cosmic density field is challenging, as they strongly depend on the tails of the probability distribution function (PDF) of density fluctuations. A way around this problem is to directly analyse the bulk of the PDF instead. For this purpose we devise a new method to predict the impact of general non-Gaussian initial conditions on the late-time density PDF. With this formalism we show that - even when marginalizing over potential ignorance of the amplitude and slope of the non-linear power spectrum - an analysis of the PDF at mildly non-linear densities can measure the amplitude of different primordial bispectrum shapes to an accuracy of $\Delta f_{\mathrm{NL}}^{\mathrm{loc}} = \pm 3.1\ ,\ \Delta f_{\mathrm{NL}}^{\mathrm{equi}} = \pm 10.0\ ,\ \Delta f_{\mathrm{NL}}^{\mathrm{ortho}} = \pm 17.0\ $. This assumes a joint analysis of the PDF on smoothing scales of $15$Mpc/$h$ and $30$Mpc/$h$ in a survey volume of $V=100(\mathrm{Gpc}/h)3$ at $z=1$, analysing only densities of $\delta (15\mathrm{Mpc}/h)\in [-0.4,0.5]$ ($\approx 87\%$ of probability) and $\delta (30\mathrm{Mpc}/h)\in [-0.3,0.4]$ ($\approx 95\%$ of probability). Note that a formalism closely related to ours was already successfully applied to observational data \citep{Gruen2018, Friedrich2018}, demonstrating that the methodology developed here can indeed be carried over to real data analysis.