The dependence of the nuclear level density on intrinsic deformation is an important input to dynamical nuclear processes such as fission. Auxiliary-field Monte Carlo (AFMC) method is a powerful method for computing nuclear level densities. However, the statistical distribution of intrinsic shapes is not readily accessible due to the formulation of AFMC in a spherical configurationinteraction shell-model approach. Instead, theory of deformation up to now has largely relied on a mean-field approximation which breaks rotational symmetry. We show here how the distributions of the intrinsic quadrupole deformation parameters can be calculated within the AFMC method, and present results for a chain of even-mass samarium nuclei ( 148 Sm, 150 Sm, 152 Sm, 154 Sm) which includes spherical, transitional, and strongly deformed isotopes. The method relies on a Landau-like expansion of the Helmholtz free energy in invariant polynomials of the quadrupole tensor. We find that an expansion to fourth order provides an excellent description of the AFMC results.