We describe a Bayesian approach to estimating quasar black hole mass functions (BHMF) when using the broad emission lines to estimate black hole mass. We show how using the broad line mass estimates in combination with statistical techniques developed for luminosity function estimation leads to statistically biased results. We derive the likelihood function for the BHMF based on the broad line mass estimates, and derive the posterior distribution for the BHMF, given the observed data. We develop our statistical approach for a flexible model where the BHMF is modelled as a mixture of Gaussian functions. Statistical inference is performed using markov chain monte carlo (MCMC) methods. Our method has the advantage that it is able to constrain the BHMF even beyond the survey detection limits at the adopted confidence level, accounts for measurement errors and the intrinsic uncertainty in broad line mass estimates, and provides a natural way of estimating the probability distribution of any quantities derived from the BHMF. We conclude by using our method to estimate the local active BHMF using the z 10^8. Our analysis implies that z