We present a probabilistic approach for inferring the parameters of the present day power-law stellar mass function ( MF ) of a resolved young star cluster . This technique ( a ) fully exploits the information content of a given dataset ; ( b ) can account for observational uncertainties in a straightforward way ; ( c ) assigns meaningful uncertainties to the inferred parameters ; ( d ) avoids the pitfalls associated with binning data ; and ( e ) can be applied to virtually any resolved young cluster , laying the groundwork for a systematic study of the high mass stellar MF ( M \gtrsim 1 M _ { \odot } ) . Using simulated clusters and Markov chain Monte Carlo sampling of the probability distribution functions , we show that estimates of the MF slope , \alpha , are unbiased and that the uncertainty , \Delta \alpha , depends primarily on the number of observed stars and on the range of stellar masses they span , assuming that the uncertainties on individual masses and the completeness are both well-characterized . Using idealized mock data , we compute the theoretical precision , i.e. , lower limits , on \alpha , and provide an analytic approximation for \Delta \alpha as a function of the observed number of stars and mass range . Comparison with literature studies shows that \sim 3/4 of quoted uncertainties are smaller than the theoretical lower limit . By correcting these uncertainties to the theoretical lower limits , we find the literature studies yield \langle \alpha \rangle = 2.46 , with a 1- \sigma dispersion of 0.35 dex . We verify that it is impossible for a power-law MF to obtain meaningful constraints on the upper mass limit of the IMF , beyond the lower bound of the most massive star actually observed . We show that avoiding substantial biases in the MF slope requires : ( 1 ) including the MF as a prior when deriving individual stellar mass estimates ; ( 2 ) modeling the uncertainties in the individual stellar masses ; and ( 3 ) fully characterizing and then explicitly modeling the completeness for stars of a given mass . The precision on MF slope recovery in this paper are lower limits , as we do not explicitly consider all possible sources of uncertainty , including dynamical effects ( e.g. , mass segregation ) , unresolved binaries , and non-coeval populations . We briefly discuss how each of these effects can be incorporated into extensions of the present framework . Finally , we emphasize that the technique and lessons learned are applicable to more general problems involving power-law fitting .