Intervals are used to represent imprecise numerical values. Modelling uncertain values with precise bounds without considering their probability distribution is infeasible in many applications. As a solution, this paper proposes the use of probability density functions instead of intervals; we consider evaluation of an arithmetical function of random variables. Since the result density cannot in general be solved algebraically, an interval method for determining its guaranteed bounds is developed. This possibility challenges traditional Monte Carlo methods in which only stochastic characterizations for the result distribution, such as confidence bounds for fractiles, can be determined.