We have developed a seismic inversion method for the joint estimation of facies and elastic properties from prestack seismic data based on a geostatistical approach. The objectives of our inversion methodology are to sample from the posterior distribution of seismic properties and to simultaneously classify the lithology conditioned by seismic data. The inversion algorithm is a sequential Gaussian mixture inversion based on Bayesian linearized amplitude variation with offset inverse theory and sequential geostatistical simulations. The stochastic approach to the inversion allows generating multiple elastic models that match the seismic data. To mathematically represent the multimodal behavior of elastic properties due to their variations within different lithologies, we adopt a Gaussian mixture distribution for the prior model of the elastic properties and we use the prior probability of the facies as weights of the Gaussian components of the mixture. The solution of the inverse problem is achieved by deriving the explicit analytical expression of the posterior distribution of the elastic properties and facies and by sampling from this distribution according to a spatial correlation model. The inversion methodology has been validated using well logs and synthetic seismic data with different noise levels, and it is then applied to a real 3D seismic data set in North Sea.