The mass composition of high energy cosmic rays above $10^{17}$ eV is a crucial issue to solve some open questions in astrophysics such as the acceleration and propagation mechanisms. Unfortunately, the standard procedures to identify the primary particle of a cosmic ray shower have low efficiency mainly due to large fluctuations and limited experimental observables. We present a statistical method for composition studies based on several measurable features of the longitudinal development of the CR shower such as $N_{max}$, $X_{max}$, asymmetry, skewness and kurtosis. Principal component analysis (PCA) was used to evaluate the relevance of each parameter in the representation of the overall shower features and a linear discriminant analysis (LDA) was used to combine the different parameters to maximize the discrimination between different particle showers. The new parameter from LDA provides a separation between primary gammas, proton and iron nuclei better than the procedures based on $X_{max}$ only. The method proposed here was successfully tested in the energy range from $10^{17}$ to $10^{20}$ eV even when limitations of shower track length were included in order to simulate the field of view of fluorescence telescopes.