This paper presents a two-step Bayesian framework for the estimation of the intrinsic single crystal elastic stiffness parameters from the measurements of spherical indentation stress-strain responses in multiple individual grains of a polycrystalline sample, whose crystal lattice orientations have been measured using electron back-scattered diffraction technique. The first step requires the establishment of the functional dependence of the indentation elastic modulus given the lattice orientation and the intrinsic single crystal elastic stiffness parameters. Previous efforts for this step required a large database of computationally expensive finite element (FE) simulations in order to establish this function with adequate accuracy. In this paper, it is shown that the introduction of a Bayesian framework can greatly reduce the number of simulations necessary to establish this function, while introducing practically useful measures of uncertainty which can guide the selection of specific additional simulations that are expected to best improve the predictive accuracy of the function. The second step involves a Markov Chain Monte-Carlo (MCMC) sampling of the distribution of possible values for the single crystal elastic stiffness parameters based on a given set of experimentally measured elastic indentation moduli in individual grains of different lattice orientations. This second step is accomplished by calibrating the available experimental data to the function established in the first step. This novel framework is presented and demonstrated in this paper for an as-cast cubic polycrystalline Fe-3% Si sample and a hexagonal polycrystalline commercially pure (CP-Ti) titanium sample.