Abstract: In Bayesian inference theory, the solution of an inverse problem is characterized by its Posterior Probability Density (PPD), which combines prior information about the model with information from an observed data set. An efficient and Fast Gibbs Sampler (FGS) has been developed to estimate the multi-dimensional integrals of the PPD, which requires solving the forward model problems many times and leads to intensive computation for multi-frequency or range-dependent inversion cases. This paper presents an alternative approach in order to speed this estimation process based on a Neighbourhood Approximation Bayes (NAB) algorithm. For lower dimension geoacoustic inversion, the NAB can approximate the PPD very well. For higher dimensional problems and sensitive parameters, however, the NAB algorithm has difficulty to estimate the PPD accurately with limited model samples. According to the preliminary estimation of the PPD by NAB, this paper developed a multi-step inversion scheme, which adjusts the parameter search intervals flexibly, in order to improve the approximation accuracy of NAB and obtain more complete information about parameter uncertainties. The prominent feature of NAB is to approximate the PPD by incorporating all models for which the forward problem has been solved into the appraisal stage. Comparison of FGS and NAB for noisy synthetic benchmark test cases and Mediterranean real data indicates that NAB provides reasonable estimates of the PPD moments while requiring less computation time.