We propose a new method of inverting gravity data for the lateral density variation of the terrain above sea level. The density variation is modeled by a piecewise step function. We impose an a priori condition that the Bouguer anomaly distribution expressed by a third-order spline function is smooth relative to the topographic distribution. No a priori condition is imposed on the density distribution. The Akaike's Bayesian Information Criterion is used to obtain optimum trade-off parameters controlling the smoothness of the Bouguer anomaly surface against its fitness to the observed anomaly values. We apply the method to about 45,000 gravity observations in southwest Japan. The calculated density distribution, in general, correlates well with major geologic features: the density of sedimentary, volcanic, or granitic units increases as their age increases in a way consistent with laboratory measurements of rock density and age. In some regions, however, where the geologic features described in the literature are likely to be surficial and to be underlain by other geologic units, the estimated densities deviate from this general trend. Some active faults sharply offset topography and geology, resulting generally in a step in the Bouguer anomaly. In such cases the estimated densities tend to be unreasonably large or small. Synthetic tests indicate that a larger number of data points and spline knots are required near the offset to climinate this problem. ¿ 1997 American Geophysical Union