  We study and compare fitting methods for the Lyman-$\alpha$ (Ly$\alpha$)
forest 3D correlation function. We use the nested sampler PolyChord and the
community code picca to perform a Bayesian analysis which we compare with
previous frequentist analyses. By studying synthetic correlation functions, we
find that the frequentist profile likelihood produces results in good agreement
with a full Bayesian analysis. On the other hand, Maximum Likelihood Estimation
with the Gaussian approximation for the uncertainties is inadequate for current
data sets. We compute for the first time the full posterior distribution from
the Ly$\alpha$ forest correlation functions measured by the extended Baryon
Oscillation Spectroscopic Survey (eBOSS). We highlight the benefits of sampling
the full posterior distribution by expanding the baseline analysis to better
understand the contamination by Damped Ly$\alpha$ systems (DLAs). We make our
improvements and results publicly available as part of the picca package.
