Application of seismic attenuation estimation using sonic waveform data is limited because the estimation methods have not yet been fully developed. Although the median frequency shift method is considered to be effective and robust compared to the conventional spectral ratio method, we demonstrated that the median frequency shift methods strongly depend on reference data under lower signal-to-noise ratios. We modified an existing median frequency shift method not to depend on arbitrarily choosing a reference value and to quantify the uncertainties in attenuation estimation. Furthermore, we implied the optimum selection of receiver pairs used for more stable attenuation analysis. Our numerical experiments supported the advantage of the proposed method. Although our main findings by applying the proposed methods in methane hydrate-bearing sediments are almost consistent with past field sonic logging measurements, we find some differences in the magnitude of attenuation values compared to existing sonic attenuation measurements and discuss various possible factors. We believe that more stable and reliable attenuation results can lead to clarifying various factors affecting attenuation estimation, such as the effect of scattering, near-field effects, and source-coupling effects. Furthermore, we emphasized the importance of scattering effect caused by the heterogeneity of the formation and demonstrated the limitation of characterizing the 1D heterogeneity using the sonic logging data spatially sampled at 0.15 m to adequately estimate the effect of scattering attenuation.