We present a new catalogue of Damped Lyman-$alpha$ absorbers from SDSS DR16Q, as well as new estimates of their statistical properties. Our estimates are computed with the Gaussian process models presented in Garnett et al. (2017); Ho et al. (2020) with an improved model for marginalising uncertainty in the mean optical depth of each quasar. We compute the column density distribution function (CDDF) at $2 < z < 5$, the line density ($textrm{d} N/ textrm{d} X$), and the neutral hydrogen density ($Omega_{textrm{DLA}}$). Our Gaussian process model provides a posterior probability distribution of the number of DLAs per spectrum, thus allowing unbiased probabilistic predictions of the statistics of DLA populations even with the noisiest data. We measure a non-zero column density distribution function for $N_{textrm{HI}} < 3 times 10^{22} ,textrm{cm}^{-2}$ with $95%$ confidence limits, and $N_{textrm{HI}} lesssim 10^{22} ,textrm{cm}^{-2}$ for spectra with signal-to-noise ratios $> 4$. Our results for DLA line density and total hydrogen density are consistent with previous measurements. Despite a small bias due to the poorly measured blue edges of the spectra, we demonstrate that our new model can measure the DLA population statistics when the DLA is in the Lyman-$beta$ forest region. We verify our results are not sensitive to the signal-to-noise ratios and redshifts of the background quasars although a residual correlation remains for detections from $z_{textrm{QSO}} < 2.5$, indicating some residual systematics when applying our models on very short spectra, where the SDSS spectral observing window only covers part of the Lyman-$alpha$ forest.
Download