  We present a Bayesian Stacking technique to directly measure the HI mass
function (HIMF) and its evolution with redshift using galaxies formally below
the nominal detection threshold. We generate galaxy samples over several sky
areas given an assumed HIMF described by a Schechter function and simulate the
HI emission lines with different levels of background noise to test the
technique. We use Multinest to constrain the parameters of the HIMF in a broad
redshift bin, demonstrating that the HIMF can be accurately reconstructed,
using the simulated spectral cube far below the HI mass limit determined by the
$5\sigma$ flux-density limit, i.e. down to $M_{\rm HI} = 10^{7.5}$ M$_{\odot}$
over the redshift range $0 < z < 0.55$ for this particular simulation, with a
noise level similar to that expected for the MIGHTEE survey. We also find that
the constraints on the parameters of the Schechter function, $\phi_{\star}$,
$M_\star$ and $\alpha$ can be reliably fit, becoming tighter as the background
noise decreases as expected, although the constraints on the redshift evolution
are not significantly affected. All the parameters become better constrained as
the survey area increases. In summary, we provide an optimal method for
estimating the HI mass at cosmological distances that allows us to constrain
the HI mass function below the detection threshold in forthcoming HI surveys.
This study is a first step towards the measurement of the HIMF at high
($z>0.1$) redshifts.
