We derive a Bayesian statistical procedure for inversion of geophysical data to rock properties. The procedure is for simplicity presented in the seismic AVO setting where rock properties influence the data through elastic parameters. The framework may however easily be extended. The procedure combines sampling based techniques and a compound Gaussian approximation to assess local approximations to marginal posterior distributions of rock properties, which the inversion is based on. The framework offers a range of approximations where inversion speed and accuracy may be balanced. The approach is also well suited for parallelisation, making it attractive for large inversion problems. We apply the procedure to a 4D CO2 monitoring case with focus on predicting saturation content. Promising results are obtained for both synthetic and real data. Finally we compare our method with regular linear Gaussian inversion for density prediction, where our method gives an improved fit.