Neutral hydrogen (HI) traces the dark matter distribution of the Universe. Upcoming surveys such as the Square Kilometre Array Observatory (SKAO) will trace neutral hydrogen up to z < 6 using several detection techniques including Intensity Mapping, which offers a unique window to explore the post-reionization Universe. Beyond two-point statistics promise to extract additional non-Gaussian information but require an accurate modelling of observational systematics such as foregrounds and the telescope beam. This work develops a theoretical model for the HI one-point probability density function (PDF) in spherical cells based on large-deviation statistics and spherical collapse for dark matter along with a nonlinear tracer bias and stochasticity parameterisation. It incorporates foreground removal and telescope beam effects that are validated against high-resolution simulations. We show that, despite these observational systematics, the HI PDF is able to capture additional non-Gaussian information from HI intensity maps compared to the power spectrum and can thus tighten constraints on cosmological parameters, breaking the degeneracy between the linear bias and the clustering amplitude.