An appropriate rock physics model is chosen from well-log data and/or experience. The rock physics model is then used to generate a plurality of attribute values, including compression and shear velocities, at different porosities and gas hydrate saturations. Gas hydrates are classified into different ranges of porosity and hydrate saturation based on the population of multiple attributes, and probability density functions for each individual gas hydrate class are created. Probability density distribution functions for individual attributes and a joint conditional probability density functions are created using a Bayesian function. The conditional probability of the occurrences of gas hydrate classes given a set of values of the chosen attributes derived from seismic inversion and or well measurements is inverted. Finally, a maximum a-posteriori (MAP) rule is employed to obtain the optimal porosity and hydrate saturation estimation. This information may be used to make decisions regarding management of the well site and the subsurface hydrate resources.