  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.
