We present a comprehensive method for determining stellar mass functions, and apply it to samples in the local Universe. We combine the classical 1/Vmax approach with STY, a parametric maximum likelihood method and SWML, a non-parametric maximum likelihood technique. In the parametric approach, we are assuming that the stellar mass function can be modelled by either a single or a double Schechter function and we use a likelihood ratio test to determine which model provides a better fit to the data. We discuss how the stellar mass completeness as a function of z biases the three estimators and how it can affect, especially the low mass end of the stellar mass function. We apply our method to SDSS DR7 data in the redshift range from 0.02 to 0.06. We find that the entire galaxy sample is best described by a double Schechter function with the following parameters: $log (M^{*}/M_odot) = 10.79 pm 0.01$, $log (Phi^{*}_1/mathrm{h^3 Mpc^{-3}}) = -3.31 pm 0.20$, $alpha_1 = -1.69 pm 0.10$, $log (Phi^{*}_2/mathrm{h^3 Mpc^{-3}}) = -2.01 pm 0.28$ and $alpha_2 = -0.79 pm 0.04$. We also use morphological classifications from Galaxy Zoo and halo mass, overdensity, central/satellite, colour and sSFR measurements to split the galaxy sample into over 130 subsamples. We determine and present the stellar mass functions and the best fit Schechter function parameters for each of these subsamples.