  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. 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.
