Usually, when inverting geodetic data to estimate the slip distributions on a fault, the area is made large enough to more than cover the rupture zone, with regularization producing regions of large slip with very small slip over the rest of the surface. We have developed a new inverse method which assumes that nonzero slip is confined to a rectangular region, and which jointly estimates, using Bayesian methods, the boundaries of this region as well as the slip distribution within it, using a smoothing parameter also determined as part of the inversion. Synthetic tests show that our method can successfully image deeper slip regions not resolved by previous methods, and does not produce spurious regions of nonzero slip. We apply our method to coseismic displacements measured by GPS for the 2009 L’Aquila earthquake, first determining the orientation of the fault assuming a simplified model with uniform slip, and then determining probability density functions for the location, length, and width of the rupture area and for the slip distribution. The standard deviation of slip is about 10 cm and describes a normal-faulting earthquake with a maximum slip of 88 ± 11 cm and seismic moment of 3.32 ± 0.30 × 1018 N m.