  The distribution of depth in which a cosmic ray air shower reaches its
maximum number of particles ($X_{max}$) is studied and parametrized. Three
functions are studied for proton, carbon, silicon, and iron primary particles
with energies ranging from $10^{17}$ eV to $10^{20}$ eV for three hadronic
interaction models: EPOS-LHC, QGSJetII.04, and Sibyll2.3c. The function which
best describes the $X_{max}$ distribution of a mixed composition is also
studied. A very large number of simulated showers and a detailed analysis
procedure were used to guarantee negligible effects of undersampling and of
fitting in the final results. For the first time, a comparison of several
functions is presented under the same assumption with the intention of
selecting the best functional form to describe the $X_{max}$ distribution. The
Generalized Gumbel distribution is shown to be the best option for a general
description of all cases. The Log-normal distribution is also a good choice for
some cases while the Exponentially Modified Gaussian distribution has shown to
be the worst choice in almost all cases studied. All three functions are
parametrized as a function of energy and primary mass.
