Intensity mapping surveys will provide access to a coarse view of the cosmic large-scale structure in unprecedented large volumes at high redshifts. Given the large fractions of the sky that can be efficiently scanned using emission from cosmic neutral hydrogen (HI), intensity mapping is ideally suited to probe a wide range of density environments and hence to constrain cosmology and fundamental physics. To efficiently extract information from 21cm intensities beyond average, one needs non-Gaussian statistics that capture large deviations from mean HI density. Counts-in-cells statistics are ideally suited for this purpose, as the statistics of matter densities in spheres can be predicted accurately on scales where their variance is below unity. We use a large state-of-the-art magneto-hydrodynamic simulation from the IllustrisTNG project to determine the relation between neutral hydrogen and matter densities in cells. We demonstrate how our theoretical knowledge about the matter PDF for a given cosmology can be used to extract a parametrisation-independent HI bias function from a measured neutral hydrogen PDF. When combining the predicted matter PDFs with a simple bias fit to the simulation, we obtain a prediction for neutral hydrogen PDFs at a few percent accuracy at scale R=5 Mpc/h from redshift z=5 to z=1. Furthermore, we find a density-dependent HI clustering signal that is consistent with theoretical expectations and could allow for joint constraints of HI bias and the amplitude of matter fluctuations or the growth of structure.