The distribution of density as a function of position within the Earth is much less well constrained than the seismic velocities. The primary information comes from the mass and moment of inertia of the Earth and this information alone requires that there be a concentration of mass towards the centre of the globe. Additional information is to be found in the frequencies of the graver normal modes of the Earth which are sensitive to density through self-gravitation effects induced in deformation. The present generation of density models has been constructed using linearized inversion techniques from earlier models, which ultimately relate back to models developed by Bullen and based in large part on physical arguments. A number of experiments in non-linear inversion have been conducted using the PREM reference model, with fixed velocity and attenuation, but with the density model constrained to lie within fixed bounds on both density and density gradient. A set of models is constructed from a uniform probability density within the bound and slope constraints. Each of the resultant density models is tested against the mass and moment of inertia of the Earth, and for successful models a comparison is made with observed normal mode frequencies. From the misfit properties of the ensemble of models the robustness of the density profile in different portions of the Earth can be assessed, which can help with the design of parametrization for future reference models. In both the lower mantle and the outer core it would be desirable to allow a more flexible representation than the single cubic polynomial employed in PREM.