We present a general methodology that aims to learn multi-variate statistics of high dimensional images, in order to capture the inter-individual variability of imaging data from a limited number of training images. The statistical learning procedure is used for identifying abnormalities as deviations from the normal variation. In most practical applications, learning an accurate statistical model of the observed data is a very challenging task due to the very high dimensionality of the images, and the limited number of available training samples. We attempt to overcome this problem by capturing the statistics of a large number of lower dimensional subspaces, which can be estimated more reliably. The subspaces are derived in a multi-scale fashion, and capture image characteristics ranging from fine and localized to coarser and relatively more global. The main premise is that an imaging pattern that is consistent with the statistics of a large number of subspaces, each reflecting a marginal probability density function (pdf), is likely to be consistent with the overall pdf, which hasn't been explicitly estimated. Abnormalities in a new image are identified as significant deviations from the normal variation captured by the learned subspace models, and are determined via iterative projections on these subspaces.