Quantifying diffeomorphic shape variability is an important tool for relating brain shape to disease processes and changes in cognitive and behavioral measures. The common template-based approach to statistical shape analysis is to use deformable image registration to map input images to a template, and then compute statistics of resulting transformations. In this talk, I will first present a Bayesian framework of atlas building that estimates the parameters that control the diffeomorphic transformation regularity. A Monte Carlo Expectation Maximization algorithm (MCEM) is developed for inference where the expectation step is approximated via sampling on the manifold of diffeomorphisms. Based on this setting, I will then introduce how to encode shape variability directly in a Bayesian model through diffeomorphic deformations. To overcome the challenge of high dimensionality of the deformations, combined with a relatively small number of image samples, we extract intrinsic low-dimensional, second-order statistics of anatomical shape variability. Bayesian inference with Markov Chain Monte Carlo (MCMC) is intractable due to the large computational cost of diffeomorphic image registration. Therefore, I will propose a fast geodesic shooting algorithm, which breaks through the prohibitive time and memory requirements of the Bayesian inference. This is achieved by introducing a novel finite-dimensional Lie algebra structure on the space of bandlimited velocity fields of diffeomorphisms.