We propose a flexible method for estimating luminosity functions (LFs) based on kernel density estimation (KDE), the most popular nonparametric density estimation approach developed in modern statistics, to overcome issues surrounding the binning of LFs. One challenge in applying KDE to LFs is how to treat the boundary bias problem, as astronomical surveys usually obtain truncated samples predominantly due to the flux-density limits of surveys. We use two solutions, the transformation KDE method) and the transformation-reflection KDE method to reduce the boundary bias. We develop a new likelihood cross-validation criterion for selecting optimal bandwidths, based on which the posterior probability distribution of the bandwidth and transformation parameters for are derived within a Markov Chain Monte Carlo sampling procedure. The simulation result shows that perform better than the traditional binning method, especially in the sparse data regime around the flux limit of a survey or at the bright end of the LF. To further improve the performance of our KDE methods, we develop the transformation-reflection adaptive KDE approach ). Monte Carlo simulations suggest that it has good stability and reliability in performance, and is around an order of magnitude more accurate than using the binning method. By applying our adaptive KDE method to a quasar sample, we find that it achieves estimates comparable to the rigorous determination in a previous work, while making far fewer assumptions about the LF. The KDE method we develop has the advantages of both parametric and nonparametric methods.