We present an analytical description of the probability distribution function (PDF) of the smoothed three-dimensional matter density field for modified gravity and dark energy. Our approach, based on the principles of Large Deviations Theory, is applicable to general extensions of the standard $Lambda$CDM cosmology. We show that late-time changes to the law of gravity and background expansion can be included through Einstein-de Sitter spherical collapse dynamics combined with linear theory calculations and a calibration measurement of the non-linear variance of the smoothed density field from a simple numerical simulation. In a comparison to $N$-body simulations for $f(R)$, DGP and evolving dark energy theories, we find percent level accuracy around the peak of the distribution for predictions in the mildly non-linear regime. A Fisher forecast of an idealised experiment with a Euclid-like survey volume demonstrates the power of combining measurements of the 3D matter PDF with the 3D matter power spectrum. This combination is shown to halve the uncertainty on parameters for an evolving dark energy model, relative to a power spectrum analysis on its own. The PDF is also found to substantially increase the detection significance for small departures from General Relativity, with improvements of up to six times compared to the power spectrum alone. This analysis is therefore very promising for future studies including non-Gaussian statistics, as it has the potential to alleviate the reliance of these analyses on expensive high resolution simulations and emulators.