This article presents a unifying framework to uncertainty quantification for systems subject to several design requirements that depend polynomially on both aleatory and epistemic uncertainties. This methodology, which is based on the Bernstein expansions of polynomials, enables calculating bounding intervals for the range of means, variances and failure probabilities of response metrics corresponding to all possible epistemic realizations. Moreover, it enables finding sets that contain the critical combination of epistemic uncertainties leading to the best-case and worst-case results, e.g., cases where the failure probability attains its smallest and largest value. These bounding intervals and sets, whose analytical structure renders them free of approximation error while eliminating the possibility of convergence to non-global optima, can be made arbitrarily tight with additional computational effort. This framework enables the consideration of arbitrary and possibly dependent aleatory variables as well as the efficient accommodation for changes in the models used to describe the uncertainty.