The structure of the magnitude distribution of earthquakes from three different seismotectonically homogeneous areas of Greece has been investigated by means of statistical inference methods. Unlike in previous studies, a nonparametric approach, namely, the smoothed bootstrap test for multimodality applied in this work makes it possible to test the complexity of the distribution without specifying any particular probabilistic model. Two null hypotheses, the number of modes in magnitude density equals 1 and the number of bumps in magnitude density equals 1, have been considered. Their alternatives mean that the magnitude population has a multicomponent structure. In two out of three studied cases the significance of the null hypotheses is less than 10%, which indicates that the magnitude distribution follows neither the log linear nor any smoothly non--log linear law but is more complex. Consequently, when the log linear model is applied to represent the magnitude distribution of each of the cases studied, estimates of mean return periods dramatically disagree with earthquake recurrence observations. The most significant differences are for large magnitude range; that is, such return period estimates are the most erroneous for large earthquakes. It has also been shown that the hazard estimation can be improved considerably by using the model-free approach with the kernel estimator of magnitude density. This approach ensures a satisfactory agreement between the mean return period estimates and actual observations, and in most of the cases, the agreement is very good.