We present a comprehensive method for determining stellar mass functions , and apply it to samples in the local Universe . We combine the classical 1/ V _ { \mathrm { max } } 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 .