This paper proposes an approach to effectively representing the statistics of high-dimensional deformations, when relatively few training samples are available, and conventional methods, like PCA, fail due to insufficient training. Based on previous work on scale-space decomposition of deformation fields, herein we represent the space of "valid deformations" as the intersection of three subspaces: one that satisfies constraints on deformations themselves, one that satisfies constraints on Jacobian determinants of deformations, and one that represents smooth deformations via a Markov Random Field (MRF). The first two are extensions of PCA-based statistical shape models. They are based on a wavelet packet basis decomposition that allows for more accurate estimation of the covariance structure of deformation or Jacobian fields, and they are used jointly due to their complementary strengths and limitations. The third is a nested MRF regularization aiming at eliminating potential discontinuities introduced by assumptions in the statistical models. A randomly sampled deformation field is projected onto the space of valid deformations via iterative projections on each of these subspaces until convergence, i.e. all three constraints are met. A deformation field simulator uses this process to generate random samples of deformation fields that are not only realistic but also representative of the full range of anatomical variability. These simulated deformations can be used for validation of deformable registration methods. Other potential uses of this approach include representation of shape priors in statistical shape models as well as various estimation and hypothesis testing paradigms in the general fields of computational anatomy and pattern recognition.