In inversion of geodetic data for spatial distribution of fault slip under smoothing and positivity constraints, subjective or incorrect techniques are often employed to select a smoothing parameter that determines the relative weight placed on fitting the data versus smoothing the slip distribution. A popular objective method based on Akaike’s Bayesian information criterion (ABIC) is incorrect if positivity constraints are employed. We introduce a fully probabilistic inversion method to simultaneously estimate the slip distribution and smoothing parameter objectively and correctly in a Bayesian framework. The complete solution to the inverse problem, the joint posterior probability density function of slip and smoothing parameter, is formulated using Bayes’ theorem and sampled with a Markov chain Monte Carlo method. We apply the new method to coseismic displacement data from the 1999 Chi-Chi, Taiwan, earthquake and compare the results of our method with conventional methods. The slip distribution estimated with our method is significantly smoother than that estimated by the ABIC-based method, demonstrating that the ABIC-based method may in general generate significantly undersmoothed or oversmoothed slip distributions under positivity constraints. On the other hand, the estimated slip distribution is similar to that estimated using another conventional method based on cross validation. Multiple trade-off curve plots between slip roughness and data misfit show that trade-off curves provide no basis for selecting a single smoothing parameter. The new method is easily generalized to slip inversions of multiple data sets with unknown relative weighting and to slip inversions with unknown fault geometry parameters.