In Bayesian inversion, the solution is characterized by its posterior probability density (PPD). A fast Gibbs sampler (FGS) has been developed to estimate the multi-dimensional integrals of the PPD, which requires solving the forward models many times and leads to intensive computation for multi-frequency or range-dependent inversion cases. This paper presents an alternative approach based on a neighborhood 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 estimating the PPD accurately with limited model samples. According to the preliminary PPD estimation from the NAB, this paper developed a multi-step inversion scheme, which adjusts the parameter search intervals flexibly, in order to improve the approximation accuracy of the NAB and obtain more complete parameter uncertainties. The prominent feature of the NAB is to approximate the PPD by incorporating all models for which the forward problem has been solved into the appraisal stage. Comparison of the FGS and NAB for noisy synthetic benchmark test cases and Mediterranean real data indicates that the NAB provides reasonable estimates of the PPD moments while requiring significantly less computation time.