Rock physics models are physical equations that map petrophysical properties into geophysical variables, such as elastic properties and density. These equations are generally used in quantitative log and seismic interpretation to estimate the properties of interest from measured well logs and seismic data. Such models are generally calibrated using core samples and well log data and result in accurate predictions of the unknown properties. Because the input data are often affected by measurement errors, the model predictions are often uncertain. Instead of applying rock physics models to deterministic measurements, I propose to apply the models to the probability density function of the measurements. This approach has been previously adopted in literature using Gaussian distributions, but for petrophysical properties of porous rocks, such as volumetric fractions of solid and fluid components, the standard probabilistic formulation based on Gaussian assumptions is not applicable due to the bounded nature of the properties, the multimodality, and the non-symmetric behavior. The proposed approach is based on the Kumaraswamy probability density function for continuous random variables, which allows modeling double bounded non-symmetric distributions and is analytically tractable, unlike the Beta or Dirichtlet distributions. I present a probabilistic rock physics model applied to double bounded continuous random variables distributed according to a Kumaraswamy distribution and derive the analytical solution of the posterior distribution of the rock physics model predictions. The method is illustrated for three rock physics models: Raymer’s equation, Dvorkin’s stiff sand model, and Kuster-Toksoz inclusion model.