We use a sampling-based Markov chain Monte Carlo method to invert seismic data directly to porosity and quantify its uncertainty distribution in a hard-rock carbonate reservoir in Southwest Iran. Due to the processing of the seismic data, the remainder noise is correlated with the bandwidth in the range of the seismic wavelet. Hence, we assume the estimated seismic wavelet as a suitable proxy for capturing the coupling of noise samples and we propose a simple and pragmatic approach to account for the correlated (colored) noise in the probabilistic inversion of real seismic data. We also calibrate an empirical and a semi-empirical inclusion-based rock-physics model to characterize the rock-frame elastic moduli via lithology constrained fitting parameters of these models, i.e. the critical porosity and the pore aspect ratio. These calibrated rock-physics models are embedded in the inversion procedure to link petrophysical and elastic properties. In addition, we obtain the pointwise critical porosity and pore aspect ratio, which can potentially facilitate the interpretation of the reservoir for further studies. The inversion results are evaluated by comparing with porosity logs and an existing geological model (porosity model) constructed through a geostatistical simulation approach. We assess the consistency of the geological model through a geomodel-to-seismic modeling approach. The results confirm the performance of the probabilistic inversion in resolving some thin layers and reconstructing the observed seismic data. We also present the applicability of the proposed sampling-based approach to invert 3D seismic data for estimating the porosity distribution and its associated uncertainty for four subzones of the reservoir. The porosity time maps and the facies probabilities obtained via porosity cut-offs indicate the relative quality of the reservoir’s subzones over each other.