We use a sampling-based Markov chain Monte Carlo method to invert seismic data directly for porosity and to quantify its uncertainty distribution in a hard-rock carbonate reservoir in southwest Iran. The noise that remains in seismic data after the processing flow is correlated with the bandwidth in the range of the seismic wavelet. Hence, to account for the inherent correlated nature of the band-limited seismic noise in the probabilistic inversion of real seismic data, we assume the estimated seismic wavelet as a suitable proxy for capturing the coupling of noise samples. In contrast to the common approach of inserting a delta function on the main diagonal of the covariance matrix, we insert the seismic wavelet on its main diagonal. We also calibrate an empirical and a semiempirical inclusion-based rock-physics model (RPM) to characterize the rock-frame elastic moduli via a lithology-constrained fitting of the parameters of these models, i.e., the critical porosity and the pore aspect ratio. These calibrated RPMs 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 geologic model (porosity model) constructed through a geostatistical simulation approach. We assess the consistency of the geologic 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 have evaluated the applicability of our 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 cutoffs indicate the relative quality of the reservoir’s subzones.