The focus of the present paper is the determination of the lower and upper bounds of the probability of failure under uncertain inputs by means of random set theory. Under this general framework, it is possible to model uncertainty in the form of probability boxes, fuzzy sets, cumulative distribution functions, Dempster-Shafer structures or intervals. In addition, the dependence between the input variables can be expressed using copulas. In order to speed up the calculation, a very efficient probability-based reliability method - known as "subset simulation" - will be used. This method is specially suited for finding small probabilities of failure in both low- and high-dimensional spaces, disjoint failure regions and nonlinear limit state functions. The proposed method represents a drastic reduction of the computational labor implied by plain Monte Carlo simulation for problems defined with a mixture of representations for the input variables, while delivering similar results. A numerical example shows the usefulness of the proposed approach.