We have analyzed in this work the updated sample of neutron star masses derived from the study of a variety of 96 binary systems containing at least one neutron star using Bayesian methods. After updating the multimodality of the distributions found in previous works, we determined the maximum mass implied by the sample using a robust truncation technique, with the result $m_{max} sim 2.5-2.6 , M_{odot}$. We have checked that this mass is actually consistent by generating synthetic data and employing a Posterior Predictive Check. A comparison with seven published $m_{max}$ values inferred from the remnant of the NS-NS merger GW170817 was performed and the tension between the latter and the obtained $m_{max}$ value quantified. Finally, we performed a Local Outlier Factor test and verified that the result for $m_{max}$ encompasses the highest individual mass determinations with the possible exception of PSR J1748-2021B. The conclusion is that the whole distribution already points toward a high value of $m_{max}$, while several lower values derived from the NS-NS merger event are disfavored and incompatible with the higher binary system masses. A large $m_{max}$ naturally accommodates the lower mass component of the event GW190814 as a neutron star.