Estimation of the covariance matrix is a pivotal step in landmark based statistical shape analysis. For high dimensional representation of the shapes, often the number of available shape examples is far too small for reliable estimation of the covariance matrix by the traditionally used Maximum Likelihood (ML) method. The ML covariance matrix is rank deficient and the eigenvectors corresponding to the small eigenvalues are arbitrary. The effect of this biasing phenomenon is the exaggeration of the importance associated with low variance subspace spanned by the eigenvectors corresponding to the smallest eigenvalues. We take a Bayesian approach to the problem and show how the prior information can be used to estimate the covariance matrix from a small number of samples in a high dimensional shape space. The performance of the proposed method is evaluated in the context of reconstructions of high resolution vertebral boundary from an incomplete and lower dimensional representation. The algorithm performs better than the ML method, especially for small number of samples in the training set. The superiority of the proposed Bayesian approach was also observed when noisy incomplete lower dimensional representation of the vertebral boundary was used in the reconstruction algorithm. Moreover, unlike other commonly used approaches, e.g., regularization, the presented method does not depend heavily on the choice of the parameters values.