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 binning of LFs. One challenge in applying KDE to LFs is how to treat the boundary bias problem, since astronomical surveys usually obtain truncated samples predominantly due to the flux-density limits of surveys. We use two solutions, the transformation KDE method (φ t ), and the transformation-reflection KDE method (φ tr ) 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 bandwidth and transformation parameters forφ t andφ tr are derived within a Markov chain Monte Carlo (MCMC) sampling procedure. The simulation result shows thatφ t andφ tr perform better than the traditional binned 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 (φ tra ). Monte Carlo simulations suggest that it has a good stability and reliability in performance, and is around an order of magnitude more accurate than using the binned 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 non-parametric methods.