Sex estimation of the individuals in a sample is fundamental for any bioarchaeological study to define a particular demographic assemblage or to classify isolated remains. Long bones are an excellent alternative for sex estimation when the most dimorphic anatomical parts are not preserved or are highly altered. Here we propose a set of discriminant functions and classification models to estimate the sex of prehistoric individuals using linear discriminant analysis and machine learning approaches. Different osteometric variables were taken from the humeri, ulnae, radii, femurs and tibias of a sample of 109 articulated skeletons buried in the collective tomb of Camino del Molino (Region of Murcia, SE-Spain), dated to the 3rd millennium BC. Sex was estimated based on standard anthropological methods and ancient DNA analysis of a control sample. Fifty-two discriminant functions with prediction thresholds higher than 0.8 on the ROC curve were obtained using independent (22) and combined variables (30). The best LDA models for sex prediction were those based on proximal epiphyseal widths or their combination with other variables, reaching values close to 0.98 on the ROC curve. The random forest-based model obtained an accuracy of 0.94 and confirmed the importance of epiphyseal widths in sex classification. This analysis is more comprehensive than univariate LDA, as it allows for ranking the importance of bones in sex discrimination and considers correlations between long bones rather than treating them as independent observations. In contrast, applying LDA to each bone makes it easier to predict the sex of other coeval collections that do not have such a complete sample. This work aims to overcome the scarcity of methods that can be applied to sex estimation of the large volume of isolated remains from Camino del Molino and for other Mediterranean skeletal series from the Late Prehistory with high biological affinity and that share similar environmental conditions.