Interpretation of hydrocarbon-bearing shale is subject to great uncertainty because of pervasive heterogeneity, thin beds, and incomplete and uncertain knowledge of saturation-porosity-resistivity models. We developed a stochastic joint-inversion method specifically developed to address the quantitative petrophysical interpretation of hydrocarbon-bearing shale. The method was based on the rapid and interactive numerical simulation of resistivity and nuclear logs. Instead of property values themselves, the estimation method delivered the a posteriori probability of each property. The Markov-chain Monte Carlo algorithm was used to sample the model space to quantify the a posteriori distribution of formation properties. Additionally, the new interpretation method allows the use of fit-for-purpose statistical correlations between water saturation, salt concentration, porosity, and electrical resistivity to implement uncertain, non-Archie resistivity models derived from core data, including those affected by total organic carbon (TOC). In the case of underdetermined estimation problems, i.e., when the number of measurements was lower than the number of unknowns, the use of a priori information enabled plausible results within prespecified petrophysical and compositional bounds. The developed stochastic interpretation technique was successfully verified with data acquired in the Barnett and Haynesville Shales. Core data (including X-ray diffraction data) were combined into a priori information for interpretation of nuclear and resistivity logs. Results consisted of mineral concentrations, TOC, and porosity together with their uncertainty. Eighty percent of the core data was located within the 95% credible interval of estimated mineral/fluid concentrations.