  We apply statistical inference on the Pierre Auger Open Data to discern for
the first time the full mass composition of cosmic rays at different energies.
Working with longitudinal electromagnetic profiles of cosmic ray showers, in
particular their peaking depths $X_{\rm max}$, we employ central moments of the
$X_{\rm max}$ distributions as features to discriminate between different
shower compositions. We find that already the first few moments entail the most
relevant information to infer the primary cosmic ray mass spectrum. Our
approach, based on an unbinned likelihood, allows us to consistently account
for sources of statistical uncertainties due to finite datasets, both measured
and simulated, as well as systematic effects. Finally, we provide a
quantitative comparison of different high energy hadronic interaction models
available in the atmospheric shower simulation codes.
