Frequent significant deviations of the observed magnitude distribution of anthropogenic seismicity from the Gutenberg-Richter relation require alternative estimation methods for probabilistic seismic hazard assessments. We evaluate five nonparametric kernel density estimation (KDE) methods on simulated samples drawn from four magnitude distribution models: the exponential, concave and convex bi-exponential, and exponential-Gaussian distributions. The latter three represent deviations from the Gutenberg-Richter relation due to the finite thickness of the seismogenic crust and the effect of characteristic earthquakes. The assumed deviations from exponentiality are never more than those met in practice. The studied KDE methods include Silverman's and Scott's rules with Abramson's bandwidth adaptation, two diffusion-based methods (ISJ and diffKDE), and adaptiveKDE, which formulates the bandwidth estimation as an optimization problem. We assess their performance for magnitudes from 2 to 6 with sample sizes of 400 to 5000, using the mean integrated square error (MISE) over 100,000 simulations. Their suitability in hazard assessments is illustrated by the mean of the mean return period (MRP) for a sample size of 1000. Among the tested methods, diffKDE provides the most accurate cumulative distribution function estimates for larger magnitudes. Even when the data is drawn from an exponential distribution, diffKDE performs comparably to maximum likelihood estimation when the sample size is at least 1000. Given that anthropogenic seismicity often deviates from the exponential model, we recommend using diffKDE for probabilistic seismic hazard assessments whenever a sufficient sample size is available.